Association/Dissociation reaction `2A <-> C`¶

with 2nd-order kinetics for `A`,¶

and 1-st order kinetics for `C`¶

Taken to equilibrium. (Adaptive variable time teps are used)

See also the experiment "1D/reactions/reaction_7"

LAST REVISED: June 4, 2023

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 src.modules.reactions.reaction_data import ChemData as chem
from src.modules.reactions.reaction_dynamics import ReactionDynamics

import numpy as np
import plotly.express as px
from src.modules.visualization.graphic_log import GraphicLog

In [3]:

# Initialize the HTML logging (for the graphics)
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_1"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

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

Initialize the System¶

Specify the chemicals and the reactions

In [4]:

# Specify the chemicals
chem_data = chem(names=["A", "C"])

# Reaction 2A <-> C , with 2nd-order kinetics for A, and 1st-order kinetics for C
chem_data.add_reaction(reactants=[(2, "A", 2)], products=["C"],
                       forward_rate=3., reverse_rate=2.)   
# Note: the first 2 in (2, "A", 2) is the stoichiometry coefficient, while the other one is the order

chem_data.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: 2 A <-> C  (kF = 3 / kR = 2 / Delta_G = -1,005.13 / K = 1.5) | 2-th order in reactant A

In [5]:

# Send a plot of the network of reactions to the HTML log file
graph_data = chem_data.prepare_graph_network()
GraphicLog.export_plot(graph_data, "vue_cytoscape_1")

[GRAPHIC ELEMENT SENT TO LOG FILE `react_4.log.htm`]

Start the simulation¶

In [6]:

dynamics = ReactionDynamics(reaction_data=chem_data)

In [7]:

# Initial concentrations of all the chemicals, in index order
dynamics.set_conc([200., 40.], snapshot=True)

In [8]:

dynamics.describe_state()

SYSTEM STATE at Time t = 0:
2 species:
  Species 0 (A). Conc: 200.0
  Species 1 (C). Conc: 40.0

In [9]:

dynamics.get_history()

Out[9]:

	SYSTEM TIME	A	C	caption
0	0.0	200.0	40.0	Initial state

Run the reaction¶

In [10]:

dynamics.set_diagnostics()       # To save diagnostic information about the call to single_compartment_react()

# All of these settings are currently close to the default values... but subject to change; set for repeatability
dynamics.set_thresholds(norm="norm_A", low=1.0, high=2.0, abort=3.5)
dynamics.set_thresholds(norm="norm_B", low=0.1, high=0.5, abort=3.0)
dynamics.set_step_factors(upshift=1.5, downshift=0.5, abort=0.5)
dynamics.set_error_step_factor(0.5)

dynamics.single_compartment_react(initial_step=0.002, reaction_duration=0.04,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  variable_steps=True, explain_variable_steps=False)

*** CAUTION: negative concentration in chemical `A` in step starting at t=0.  It will be AUTOMATICALLY CORRECTED with a reduction in time step size, as follows:
    INFO: the tentative time step (0.002) leads to a NEGATIVE concentration of `A` from reaction 2 A <-> C (rxn # 0): 
      Baseline value: 200 ; delta conc: -479.68
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 0.001) [Step started at t=0, and will rewind there]

*** CAUTION: negative concentration in chemical `A` in step starting at t=0.  It will be AUTOMATICALLY CORRECTED with a reduction in time step size, as follows:
    INFO: the tentative time step (0.001) leads to a NEGATIVE concentration of `A` from reaction 2 A <-> C (rxn # 0): 
      Baseline value: 200 ; delta conc: -239.84
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 0.0005) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.0005) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 0.00025) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.00025) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 0.000125) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.000125) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 6.25e-05) [Step started at t=0, and will rewind there]
INFO: the tentative time step (6.25e-05) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 3.125e-05) [Step started at t=0, and will rewind there]
INFO: the tentative time step (3.125e-05) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 1.5625e-05) [Step started at t=0, and will rewind there]
INFO: the tentative time step (1.5625e-05) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.5 (set to 7.8125e-06) [Step started at t=0, and will rewind there]
Some steps were backtracked and re-done, to prevent negative concentrations or excessively large concentration changes
103 total step(s) taken

Note how the (tentative) original time step that we provide, 0.002, turned out to be so large that the simulation backtracks several times, because of "hard" aborts (negative concentrations) or "soft" aborts (concentration changes surpassing the thresholds we provided)¶

In [11]:

dynamics.get_history()

Out[11]:

	SYSTEM TIME	A	C	caption
0	0.000000	200.000000	40.000000	Initial state
1	0.000008	198.126250	40.936875	1st reaction step
2	0.000016	196.287498	41.856251
3	0.000023	194.482769	42.758616
4	0.000031	192.711126	43.644437
...	...	...	...	...
99	0.019385	10.796656	134.601672
100	0.023088	10.200460	134.899770
101	0.028643	9.729999	135.135000
102	0.036975	9.500894	135.249553
103	0.049473	9.493349	135.253325	last reaction step

104 rows × 4 columns

In [12]:

dynamics.plot_curves(colors=['red', 'green'],
                     title="Reaction 2A <-> C  (2nd order in A).  Changes in concentrations with time")

Note: "A" (now largely depleted) is the limiting reagent¶

Let's take a look at time t=0.002, which in our simulation run had proposed as the first step:¶

In [13]:

# Locate the value closest to the original time step we had requested
dynamics.get_history(t=0.002)

Out[13]:

	search_value	SYSTEM TIME	A	C	caption
74	0.002	0.002044	57.328404	111.335798

Because of the very large changes happening between t=0 and 0.002, the simulation automatically slowed down and opted to actually take 74 steps in lieu of the 1 step we had proposed.¶

The number of variable steps actually taken can be modulated by changing the "norm thresholds" and "step factors" that we can optionally specify

Notice how, later in the simulation, the step sizes get BIGGER than the 0.002 we had originally proposed:¶

In [14]:

dynamics.explain_time_advance()

From time 0 to 3.125e-05, in 4 steps of 7.81e-06
From time 3.125e-05 to 4.297e-05, in 1 step of 1.17e-05
From time 4.297e-05 to 4.883e-05, in 1 step of 5.86e-06
From time 4.883e-05 to 8.398e-05, in 4 steps of 8.79e-06
From time 8.398e-05 to 9.717e-05, in 1 step of 1.32e-05
From time 9.717e-05 to 0.0001038, in 1 step of 6.59e-06
From time 0.0001038 to 0.0001433, in 4 steps of 9.89e-06
From time 0.0001433 to 0.0001581, in 1 step of 1.48e-05
From time 0.0001581 to 0.0001656, in 1 step of 7.42e-06
From time 0.0001656 to 0.0001989, in 3 steps of 1.11e-05
From time 0.0001989 to 0.0002156, in 1 step of 1.67e-05
From time 0.0002156 to 0.000224, in 1 step of 8.34e-06
From time 0.000224 to 0.0002615, in 3 steps of 1.25e-05
From time 0.0002615 to 0.0002803, in 1 step of 1.88e-05
From time 0.0002803 to 0.0002897, in 1 step of 9.39e-06
From time 0.0002897 to 0.0003319, in 3 steps of 1.41e-05
From time 0.0003319 to 0.000353, in 1 step of 2.11e-05
From time 0.000353 to 0.0003636, in 1 step of 1.06e-05
From time 0.0003636 to 0.0003952, in 2 steps of 1.58e-05
From time 0.0003952 to 0.000419, in 1 step of 2.38e-05
From time 0.000419 to 0.0004309, in 1 step of 1.19e-05
From time 0.0004309 to 0.0004665, in 2 steps of 1.78e-05
From time 0.0004665 to 0.0004932, in 1 step of 2.67e-05
From time 0.0004932 to 0.0005066, in 1 step of 1.34e-05
From time 0.0005066 to 0.0005467, in 2 steps of 2e-05
From time 0.0005467 to 0.0005768, in 1 step of 3.01e-05
From time 0.0005768 to 0.0005918, in 1 step of 1.5e-05
From time 0.0005918 to 0.0006369, in 2 steps of 2.26e-05
From time 0.0006369 to 0.0006707, in 1 step of 3.38e-05
From time 0.0006707 to 0.0006876, in 1 step of 1.69e-05
From time 0.0006876 to 0.0007384, in 2 steps of 2.54e-05
From time 0.0007384 to 0.001081, in 9 steps of 3.81e-05
From time 0.001081 to 0.001138, in 1 step of 5.71e-05
From time 0.001138 to 0.001166, in 1 step of 2.85e-05
From time 0.001166 to 0.001209, in 1 step of 4.28e-05
From time 0.001209 to 0.001659, in 7 steps of 6.42e-05
From time 0.001659 to 0.002237, in 6 steps of 9.63e-05
From time 0.002237 to 0.002959, in 5 steps of 0.000144
From time 0.002959 to 0.003826, in 4 steps of 0.000217
From time 0.003826 to 0.004801, in 3 steps of 0.000325
From time 0.004801 to 0.006264, in 3 steps of 0.000488
From time 0.006264 to 0.007727, in 2 steps of 0.000731
From time 0.007727 to 0.009922, in 2 steps of 0.0011
From time 0.009922 to 0.01321, in 2 steps of 0.00165
From time 0.01321 to 0.01568, in 1 step of 0.00247
From time 0.01568 to 0.02309, in 2 steps of 0.0037
From time 0.02309 to 0.02864, in 1 step of 0.00555
From time 0.02864 to 0.03697, in 1 step of 0.00833
From time 0.03697 to 0.04947, in 1 step of 0.0125
(103 steps total)

Notice how the reaction proceeds in far-smaller steps in the early times, when the concentrations are changing much more rapidly¶

In [15]:

# Let's look at the first two arrays of concentrations, from the run's history
arr0 = dynamics.get_historical_concentrations(0)   # The initial concentrations
arr1 = dynamics.get_historical_concentrations(1)   # After the first actual step
arr0, arr1

Out[15]:

(array([200.,  40.], dtype=float32),
 array([198.12625 ,  40.936874], dtype=float32))

In [16]:

# Let's verify that the reaction's stoichiometry is being respected
dynamics.stoichiometry_checker(rxn_index=0, 
                               conc_arr_before = arr0, 
                               conc_arr_after = arr1)

Out[16]:

True

Indeed, it can be easy checked that the drop in [A] is twice the increase in [C], as dictated by the stoichiometry¶

In [17]:

dynamics.get_diagnostic_rxn_data(rxn_index=0, head=15)    # Easily seen in the diagnostic data

Reaction:  2 A <-> C

Out[17]:

	START_TIME	Delta A	Delta C	time_step	caption
0	0.000000	NaN	NaN	0.002000	aborted: neg. conc. in `A`
1	0.000000	NaN	NaN	0.001000	aborted: neg. conc. in `A`
2	0.000000	-119.920000	59.960000	0.000500	aborted: excessive norm value(s)
3	0.000000	-59.960000	29.980000	0.000250	aborted: excessive norm value(s)
4	0.000000	-29.980000	14.990000	0.000125	aborted: excessive norm value(s)
5	0.000000	-14.990000	7.495000	0.000063	aborted: excessive norm value(s)
6	0.000000	-7.495000	3.747500	0.000031	aborted: excessive norm value(s)
7	0.000000	-3.747500	1.873750	0.000016	aborted: excessive norm value(s)
8	0.000000	-1.873750	0.936875	0.000008
9	0.000008	-1.838752	0.919376	0.000008
10	0.000016	-1.804729	0.902364	0.000008
11	0.000023	-1.771643	0.885821	0.000008
12	0.000031	-2.609190	1.304595	0.000012
13	0.000043	-1.269449	0.634725	0.000006
14	0.000049	-1.878784	0.939392	0.000009

From the diagnostic data, it can be seen that the first step had several false starts - and was automatically repeatedly shrunk - but finally happened. `Delta A` indeed equals - 2 * `Delta C`, satisfying the stoichiometry¶

Check the final equilibrium¶

In [18]:

# Verify that the reaction has reached equilibrium
dynamics.is_in_equilibrium()

2 A <-> C
Final concentrations:  [C] = 135.3 ; [A] = 9.493
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 1.50075
    Formula used:  [C] / [A]^2 
2. Ratio of forward/reverse reaction rates: 1.5
Discrepancy between the two values: 0.05016 %
Reaction IS in equilibrium (within 1% tolerance)

Out[18]:

True

In [19]:

dynamics.stoichiometry_checker_entire_run()

Out[19]:

True

Display the variable time steps¶

In [20]:

dynamics.plot_curves(colors=['red', 'green'], show_intervals=True,
                     title="Reaction 2A <-> C  (2nd order in A).  Changes in concentrations with time")

The intersection of the two lines may be found as follows:¶

In [21]:

dynamics.curve_intersection('A', 'C', t_start=0, t_end=0.01)

Min abs distance found at data row: 56

Out[21]:

(0.0009405253857178681, 93.33333333333333)

In [ ]:

For diagnostic insight, uncomment the following lines:¶

In [22]:

#dynamics.get_diagnostic_decisions_data()

#dynamics.get_diagnostic_rxn_data(rxn_index=0)

#dynamics.get_diagnostic_conc_data()