One-bin `2A <-> B` reaction,¶

COMPARING 1st-order and 2nd-order kinetics in forward direction; reverse direction 1-st order¶

Diffusion not applicable (just 1 bin)

See also the experiment "reactions_single_compartment/react_4"

LAST REVISED: June 23, 2024 (using v. 1.0 beta34.1)

In [1]:

import set_path      # Importing this module will add the project's home directory to sys.path

Added 'D:\Docs\- MY CODE\BioSimulations\life123-Win7' to sys.path

In [2]:

from experiments.get_notebook_info import get_notebook_basename

from life123 import ChemData as chem
from life123 import BioSim1D

import plotly.express as px
from life123 import HtmlLog as log
from life123 import GraphicLog

In [3]:

# Initialize the HTML logging
log_file = get_notebook_basename() + ".log.htm"    # Use the notebook base filename for the log file

# Set up the use of some specified graphic (Vue) components
GraphicLog.config(filename=log_file,
                  components=["vue_cytoscape_2"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

-> Output will be LOGGED into the file 'reaction_7.log.htm'

INITIALLY, with 1st-order kinetics in both directions¶

In [4]:

# Initialize the system
chem_data = chem(names=["A", "B"])     # NOTE: Diffusion not applicable (just 1 bin)

# Reaction  2A <-> B , FOR NOW with 1st-order kinetics in both directions
chem_data.add_reaction(reactants=[(2, "A", 1)], products=["B"], forward_rate=5., reverse_rate=2.)

bio = BioSim1D(n_bins=1, chem_data=chem_data)

bio.set_all_uniform_concentrations( [3., 5.] )

bio.describe_state()

SYSTEM STATE at Time t = 0:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [3.]
  Species 1 (B). Diff rate: None. Conc: [5.]

In [5]:

# Save the state of the concentrations of all species at bin 0
bio.add_snapshot(bio.bin_snapshot(bin_address = 0), caption="Initial state")
bio.get_history()

Out[5]:

	SYSTEM TIME	A	B	caption
0	0	3.0	5.0	Initial state

In [6]:

chem_data.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: 2 A <-> B  (kF = 5 / kR = 2 / delta_G = -2,271.4 / K = 2.5) | 1st order in all reactants & products
Set of chemicals involved in the above reactions: {'B', 'A'}

In [7]:

# Send a header and a plot to the HTML log file
log.write("Reaction 2A <-> B is 1st order in all species:",
          style=log.h2)

chem_data.plot_reaction_network("vue_cytoscape_2")

Reaction 2A  B is 1st order in all species:
[GRAPHIC ELEMENT SENT TO LOG FILE `reaction_7.log.htm`]

In [8]:

# First step
bio.react(time_step=0.02, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()

SYSTEM STATE at Time t = 0.02:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [2.8]
  Species 1 (B). Diff rate: None. Conc: [5.1]

Small conc. changes so far: [A] = 2.8 , [B] = 5.1

In [9]:

bio.get_history()

Out[9]:

	SYSTEM TIME	A	B	caption
0	0.00	3.0	5.0	Initial state
1	0.02	2.8	5.1

In [10]:

# Numerous more steps, to equilibrium
bio.react(time_step=0.02, n_steps=20, snapshots={"sample_bin": 0})

bio.describe_state()

SYSTEM STATE at Time t = 0.42:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [2.16928427]
  Species 1 (B). Diff rate: None. Conc: [5.41535786]

Consistent with the 5/2 ratio of forward/reverse rates (and the 1st order reactions), the systems settles in the following equilibrium:
[A] = 2.16928427 , [B] = 5.41535786

In [11]:

# Verify that the reaction has reached equilibrium
bio.reaction_dynamics.is_in_equilibrium(rxn_index=0, conc=bio.bin_snapshot(bin_address = 0))

2 A <-> B
Final concentrations: [A] = 2.169 ; [B] = 5.415
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2.49638
    Formula used:  [B] / [A]
2. Ratio of forward/reverse reaction rates: 2.5
Discrepancy between the two values: 0.1448 %
Reaction IS in equilibrium (within 1% tolerance)

Out[11]:

True

In [12]:

df = bio.get_history()
df

Out[12]:

	SYSTEM TIME	A	B	caption
0	0.00	3.000000	5.000000	Initial state
1	0.02	2.800000	5.100000
2	0.04	2.648000	5.176000
3	0.06	2.532480	5.233760
4	0.08	2.444685	5.277658
5	0.10	2.377960	5.311020
6	0.12	2.327250	5.336375
7	0.14	2.288710	5.355645
8	0.16	2.259420	5.370290
9	0.18	2.237159	5.381421
10	0.20	2.220241	5.389880
11	0.22	2.207383	5.396309
12	0.24	2.197611	5.401194
13	0.26	2.190184	5.404908
14	0.28	2.184540	5.407730
15	0.30	2.180251	5.409875
16	0.32	2.176990	5.411505
17	0.34	2.174513	5.412744
18	0.36	2.172630	5.413685
19	0.38	2.171199	5.414401
20	0.40	2.170111	5.414945
21	0.42	2.169284	5.415358

In [13]:

fig = px.line(data_frame=bio.get_history(), x="SYSTEM TIME", y=["A", "B"], 
              title="2A <-> B : changes in concentrations with time",
              color_discrete_sequence = ['navy', 'orange'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()

A gets depleted, while B gets produced.

Let's verify that the stoichiometry is being respected¶

In [14]:

# We'll check the first two arrays of concentrations, from the run's history
arr0 = bio.reaction_dynamics.get_historical_concentrations(row=0, df=df)
arr1 = bio.reaction_dynamics.get_historical_concentrations(row=1, df=df)
arr0, arr1

Out[14]:

(array([3., 5.], dtype=float32), array([2.8, 5.1], dtype=float32))

In [15]:

bio.reaction_dynamics.stoichiometry_checker(rxn_index=0, 
                               conc_arr_before = arr0, 
                               conc_arr_after = arr1)

Out[15]:

True

STARTING OVER, this time with 2nd-order kinetics in the forward reaction¶

In [16]:

bio.reaction_dynamics.clear_reactions()

In [17]:

# Reaction  2A <-> B , NOW WITH 2nd-order kinetics in the forward direction
chem_data.add_reaction(reactants=[(2, "A", 2)], products=["B"], forward_rate=5., reverse_rate=2.)

Out[17]:

In [18]:

# RESET the concentrations to their original values
bio.set_all_uniform_concentrations( [3., 5.] )

bio.describe_state()

SYSTEM STATE at Time t = 0.42:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [3.]
  Species 1 (B). Diff rate: None. Conc: [5.]

In [19]:

# Save the state of the concentrations of all species at bin 0
bio.add_snapshot(bio.bin_snapshot(bin_address = 0),
                 caption = "RESET all concentrations to initial values")
bio.get_history()

Out[19]:

	SYSTEM TIME	A	B	caption
0	0.00	3.000000	5.000000	Initial state
1	0.02	2.800000	5.100000
2	0.04	2.648000	5.176000
3	0.06	2.532480	5.233760
4	0.08	2.444685	5.277658
5	0.10	2.377960	5.311020
6	0.12	2.327250	5.336375
7	0.14	2.288710	5.355645
8	0.16	2.259420	5.370290
9	0.18	2.237159	5.381421
10	0.20	2.220241	5.389880
11	0.22	2.207383	5.396309
12	0.24	2.197611	5.401194
13	0.26	2.190184	5.404908
14	0.28	2.184540	5.407730
15	0.30	2.180251	5.409875
16	0.32	2.176990	5.411505
17	0.34	2.174513	5.412744
18	0.36	2.172630	5.413685
19	0.38	2.171199	5.414401
20	0.40	2.170111	5.414945
21	0.42	2.169284	5.415358
22	0.42	3.000000	5.000000	RESET all concentrations to initial values

In [20]:

chem_data.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: 2 A <-> B  (kF = 5 / kR = 2 / delta_G = -2,271.4 / K = 2.5) | 2-th order in reactant A
Set of chemicals involved in the above reactions: {'B', 'A'}

In [21]:

# Send a header and a plot to the HTML log file
log.write("Reaction 2A <-> B is 2nd order in A, and 1st order in B:",
          style=log.h2)
chem_data.plot_reaction_network("vue_cytoscape_2")

Reaction 2A  B is 2nd order in A, and 1st order in B:
[GRAPHIC ELEMENT SENT TO LOG FILE `reaction_7.log.htm`]

In [22]:

# First step
bio.react(time_step=0.02, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()

SYSTEM STATE at Time t = 0.44:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [1.6]
  Species 1 (B). Diff rate: None. Conc: [5.7]

[A] = 1.6 , [B] = 5.7 (Contrast with the counterpart in the 1st order kinetics: [A] = 2.8 , [B] = 5.1)

In [23]:

bio.get_history()

Out[23]:

	SYSTEM TIME	A	B	caption
0	0.00	3.000000	5.000000	Initial state
1	0.02	2.800000	5.100000
2	0.04	2.648000	5.176000
3	0.06	2.532480	5.233760
4	0.08	2.444685	5.277658
5	0.10	2.377960	5.311020
6	0.12	2.327250	5.336375
7	0.14	2.288710	5.355645
8	0.16	2.259420	5.370290
9	0.18	2.237159	5.381421
10	0.20	2.220241	5.389880
11	0.22	2.207383	5.396309
12	0.24	2.197611	5.401194
13	0.26	2.190184	5.404908
14	0.28	2.184540	5.407730
15	0.30	2.180251	5.409875
16	0.32	2.176990	5.411505
17	0.34	2.174513	5.412744
18	0.36	2.172630	5.413685
19	0.38	2.171199	5.414401
20	0.40	2.170111	5.414945
21	0.42	2.169284	5.415358
22	0.42	3.000000	5.000000	RESET all concentrations to initial values
23	0.44	1.600000	5.700000

In [24]:

# Numerous more steps
bio.react(time_step=0.02, n_steps=20, snapshots={"sample_bin": 0})

bio.describe_state()

SYSTEM STATE at Time t = 0.84:
1 bins and 2 species:
  Species 0 (A). Diff rate: None. Conc: [1.51554944]
  Species 1 (B). Diff rate: None. Conc: [5.74222528]

The systems settles in the following equilibrium: [A] = 1.51554944 , [B] = 5.74222528

In [25]:

# Verify that the reaction has reached equilibrium
bio.reaction_dynamics.is_in_equilibrium(rxn_index=0, conc=bio.bin_snapshot(bin_address = 0))

2 A <-> B
Final concentrations: [A] = 1.516 ; [B] = 5.742
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2.5
    Formula used:  [B] / [A]^2 
2. Ratio of forward/reverse reaction rates: 2.5
Discrepancy between the two values: 1.041e-08 %
Reaction IS in equilibrium (within 1% tolerance)

Out[25]:

True

In [26]:

df2 = bio.get_history()
df2

Out[26]:

	SYSTEM TIME	A	B	caption
0	0.00	3.000000	5.000000	Initial state
1	0.02	2.800000	5.100000
2	0.04	2.648000	5.176000
3	0.06	2.532480	5.233760
4	0.08	2.444685	5.277658
5	0.10	2.377960	5.311020
6	0.12	2.327250	5.336375
7	0.14	2.288710	5.355645
8	0.16	2.259420	5.370290
9	0.18	2.237159	5.381421
10	0.20	2.220241	5.389880
11	0.22	2.207383	5.396309
12	0.24	2.197611	5.401194
13	0.26	2.190184	5.404908
14	0.28	2.184540	5.407730
15	0.30	2.180251	5.409875
16	0.32	2.176990	5.411505
17	0.34	2.174513	5.412744
18	0.36	2.172630	5.413685
19	0.38	2.171199	5.414401
20	0.40	2.170111	5.414945
21	0.42	2.169284	5.415358
22	0.42	3.000000	5.000000	RESET all concentrations to initial values
23	0.44	1.600000	5.700000
24	0.46	1.544000	5.728000
25	0.48	1.525453	5.737274
26	0.50	1.519033	5.740483
27	0.52	1.516780	5.741610
28	0.54	1.515984	5.742008
29	0.56	1.515703	5.742148
30	0.58	1.515604	5.742198
31	0.60	1.515569	5.742216
32	0.62	1.515556	5.742222
33	0.64	1.515552	5.742224
34	0.66	1.515550	5.742225
35	0.68	1.515550	5.742225
36	0.70	1.515550	5.742225
37	0.72	1.515549	5.742225
38	0.74	1.515549	5.742225
39	0.76	1.515549	5.742225
40	0.78	1.515549	5.742225
41	0.80	1.515549	5.742225
42	0.82	1.515549	5.742225
43	0.84	1.515549	5.742225

In [27]:

fig = px.line(data_frame=bio.get_history(), x="SYSTEM TIME", y=["A", "B"], 
              title="2A <-> B : changes in concentrations (the jump at 0.42 is the concentration reset)",
              color_discrete_sequence = ['navy', 'orange'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()

Compared to first-order kinetics in A, the (2nd order in A) reaction now takes place much more quickly, and proceeds to almost complete depletion of A

Let's verify that the stoichiometry is still being respected¶

In [28]:

# We'll check the first two arrays of concentrations, from the run's history
arr0 = bio.reaction_dynamics.get_historical_concentrations(row=22, df=df2)  # Row 22 is the conc. reset
arr1 = bio.reaction_dynamics.get_historical_concentrations(row=23, df=df2)
arr0, arr1

Out[28]:

(array([3., 5.], dtype=float32), array([1.6, 5.7], dtype=float32))

In [29]:

bio.reaction_dynamics.stoichiometry_checker(rxn_index=0, 
                               conc_arr_before = arr0, 
                               conc_arr_after = arr1)

Out[29]:

True

In [ ]:

One-bin 2A <-> B reaction,¶