Association/Dissociation reaction `A + B <-> C`¶

with 1st-order kinetics for each species, taken to equilibrium.¶

Exploration of debugging and diagnostics options¶

(Adaptive variable time steps are used)

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

LAST REVISED: May 5, 2024

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_dynamics import ReactionDynamics
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_2"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

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

In [ ]:

Initialize the System¶

Specify the chemicals, the reactions, and the initial concentrations

In [4]:

# Instantiate the simulator and specify the chemicals
dynamics = ReactionDynamics(names=["A", "B", "C"])

In [5]:

# Reaction A + B <-> C , with 1st-order kinetics for each species
dynamics.add_reaction(reactants=["A" , "B"], products="C",
                      forward_rate=5., reverse_rate=2.)

dynamics.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: A + B <-> C  (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: {'A', 'C', 'B'}

In [6]:

# Send a plot of the network of reactions to the HTML log file
dynamics.plot_reaction_network("vue_cytoscape_2")

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

In [7]:

# Set the initial concentrations of all the chemicals, in their index order
dynamics.set_conc([10., 50., 20.], snapshot=True)

In [8]:

dynamics.describe_state()

SYSTEM STATE at Time t = 0:
3 species:
  Species 0 (A). Conc: 10.0
  Species 1 (B). Conc: 50.0
  Species 2 (C). Conc: 20.0
Set of chemicals involved in reactions: {'A', 'C', 'B'}

In [9]:

dynamics.get_history()

Out[9]:

	SYSTEM TIME	A	B	C	caption
0	0.0	10.0	50.0	20.0	Initial state

In [ ]:

Sneak preview of eventual equilibrum:¶

we can preview the final equilibrium concentrations without actually running the simulation

In [10]:

dynamics.find_equilibrium_conc(rxn_index=0)    # This is an EXACT solution

Out[10]:

{'A': 0.2948774087575341, 'B': 40.294877408757536, 'C': 29.705122591242464}

The reaction will proceed forward, with A and B being consumed, and C being produced

In [ ]:

Run the reaction¶

In [11]:

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

# For repeatibility, we avoid the defaults, and instead specify a particular group of preset parameters 
# applicable to the adaptive time steps.
# Here we use a "fast" heuristic: advance quickly thru time
dynamics.use_adaptive_preset(preset="fast")

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

* INFO: the tentative time step (0.004) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.6 (set to 0.0024) [Step started at t=0, and will rewind there]
* INFO: the tentative time step (0.0024) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.6 (set to 0.00144) [Step started at t=0, and will rewind there]
* INFO: the tentative time step (0.00144) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.6 (set to 0.000864) [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
23 total step(s) taken

In [12]:

dynamics.get_history()

Out[12]:

	SYSTEM TIME	A	B	C	caption
0	0.000000	10.000000	50.000000	20.000000	Initial state
1	0.000864	7.874560	47.874560	22.125440	1st reaction step
2	0.001555	6.602265	46.602265	23.397735
3	0.002246	5.571266	45.571266	24.428734
4	0.002938	4.727594	44.727594	25.272406
5	0.003629	4.031746	44.031746	25.968254
6	0.004666	3.165305	43.165305	26.834695
7	0.005702	2.512652	42.512652	27.487348
8	0.006739	2.015898	42.015898	27.984102
9	0.007776	1.634842	41.634842	28.365158
10	0.008813	1.340804	41.340804	28.659196
11	0.009850	1.112883	41.112883	28.887117
12	0.010886	0.935595	40.935595	29.064405
13	0.011923	0.797321	40.797321	29.202679
14	0.013478	0.635211	40.635211	29.364789
15	0.015034	0.525833	40.525833	29.474167
16	0.016589	0.451805	40.451805	29.548195
17	0.018922	0.376490	40.376490	29.623510
18	0.021254	0.337393	40.337393	29.662607
19	0.024754	0.306871	40.306871	29.693129
20	0.030002	0.293965	40.293965	29.706035
21	0.037876	0.295437	40.295437	29.704563
22	0.049685	0.294082	40.294082	29.705918
23	0.067400	0.296969	40.296969	29.703031	last reaction step

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

In [13]:

dynamics.explain_time_advance()

From time 0 to 0.000864, in 1 step of 0.000864
From time 0.000864 to 0.003629, in 4 steps of 0.000691
From time 0.003629 to 0.01192, in 8 steps of 0.00104
From time 0.01192 to 0.01659, in 3 steps of 0.00156
From time 0.01659 to 0.02125, in 2 steps of 0.00233
From time 0.02125 to 0.02475, in 1 step of 0.0035
From time 0.02475 to 0.03, in 1 step of 0.00525
From time 0.03 to 0.03788, in 1 step of 0.00787
From time 0.03788 to 0.04969, in 1 step of 0.0118
From time 0.04969 to 0.0674, in 1 step of 0.0177
(23 steps total)

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

The argument argument variable_steps=True dynamically adjusts the initial_step (which is initially found to be too large, leading to some backtracking)¶

In [14]:

dynamics.plot_step_sizes(show_intervals=True)

Plots changes of concentration with time¶

In [15]:

dynamics.plot_history(colors=['red', 'violet', 'green'], show_intervals=True)

Check the final equilibrium¶

In [16]:

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

0: A + B <-> C
Final concentrations: [A] = 0.297 ; [B] = 40.3 ; [C] = 29.7
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2.48209
    Formula used:  [C] / ([A][B])
2. Ratio of forward/reverse reaction rates: 2.5
Discrepancy between the two values: 0.7163 %
Reaction IS in equilibrium (within 1% tolerance)

Out[16]:

True

Compare with the values we saw earlier for the exact solution of the equilibrium values:
{'A': 0.2948774087575341, 'B': 40.294877408757536, 'C': 29.705122591242464}

It's instructive to compare the exact values with the last few points from the simulation:

In [17]:

dynamics.get_history(tail=3)

Out[17]:

	SYSTEM TIME	A	B	C	caption
21	0.037876	0.295437	40.295437	29.704563
22	0.049685	0.294082	40.294082	29.705918
23	0.067400	0.296969	40.296969	29.703031	last reaction step

The 2nd-to-last simulation point, rather than the last one, is actually closer to the exact equilibrium values.
That's because by that time the variable steps are getting so large that they introduce some error.
If we were to run the simulation longer (not shown), we'd see the variable steps continuing to grow, and then suddenly being reduced;
then continued cycles of growth and reduction ("hitting the brakes whenever getting too fast")

In [ ]:

Everthing below is just for diagnostic insight¶

into the adaptive variable time steps¶

This information is available because we made a call to dynamics.set_diagnostics() prior to running the simulation

In [18]:

dynamics.get_diagnostic_decisions_data()   # diagnostic data about concentration changes at every step - EVEN aborted ones

Out[18]:

	START_TIME	Delta A	Delta B	Delta C	norm_A	norm_B	action	step_factor	time_step	caption
0	0.000000	-9.840000	-9.840000	9.840000	3.227520e+01	NaN	ABORT	0.6	0.004000	excessive norm value(s)
1	0.000000	-5.904000	-5.904000	5.904000	1.161907e+01	NaN	ABORT	0.6	0.002400	excessive norm value(s)
2	0.000000	-3.542400	-3.542400	3.542400	4.182866e+00	NaN	ABORT	0.6	0.001440	excessive norm value(s)
3	0.000000	-2.125440	-2.125440	2.125440	1.505832e+00	NaN	OK (high)	0.8	0.000864
4	0.000864	-1.272295	-1.272295	1.272295	5.395782e-01	0.161570	OK (stay)	1.0	0.000691
5	0.001555	-1.030999	-1.030999	1.030999	3.543195e-01	0.156158	OK (stay)	1.0	0.000691
6	0.002246	-0.843672	-0.843672	0.843672	2.372610e-01	0.151433	OK (stay)	1.0	0.000691
7	0.002938	-0.695848	-0.695848	0.695848	1.614015e-01	0.147189	OK (low)	1.5	0.000691
8	0.003629	-0.866441	-0.866441	0.866441	2.502399e-01	0.214905	OK (stay)	1.0	0.001037
9	0.004666	-0.652653	-0.652653	0.652653	1.419851e-01	0.206189	OK (stay)	1.0	0.001037
10	0.005702	-0.496755	-0.496755	0.496755	8.225505e-02	0.197701	OK (stay)	1.0	0.001037
11	0.006739	-0.381056	-0.381056	0.381056	4.840115e-02	0.189025	OK (stay)	1.0	0.001037
12	0.007776	-0.294038	-0.294038	0.294038	2.881949e-02	0.179857	OK (stay)	1.0	0.001037
13	0.008813	-0.227921	-0.227921	0.227921	1.731599e-02	0.169988	OK (stay)	1.0	0.001037
14	0.009850	-0.177288	-0.177288	0.177288	1.047695e-02	0.159305	OK (stay)	1.0	0.001037
15	0.010886	-0.138275	-0.138275	0.138275	6.373314e-03	0.147793	OK (low)	1.5	0.001037
16	0.011923	-0.162110	-0.162110	0.162110	8.759877e-03	0.203318	OK (stay)	1.0	0.001555
17	0.013478	-0.109377	-0.109377	0.109377	3.987793e-03	0.172190	OK (stay)	1.0	0.001555
18	0.015034	-0.074029	-0.074029	0.074029	1.826757e-03	0.140784	OK (low)	1.5	0.001555
19	0.016589	-0.075315	-0.075315	0.075315	1.890774e-03	0.166698	OK (stay)	1.0	0.002333
20	0.018922	-0.039097	-0.039097	0.039097	5.095229e-04	0.103846	OK (low)	1.5	0.002333
21	0.021254	-0.030522	-0.030522	0.030522	3.105277e-04	0.090464	OK (low)	1.5	0.003499
22	0.024754	-0.012906	-0.012906	0.012906	5.551890e-05	0.042056	OK (low)	1.5	0.005249
23	0.030002	0.001472	0.001472	-0.001472	7.220164e-07	0.005007	OK (low)	1.5	0.007873
24	0.037876	-0.001355	-0.001355	0.001355	6.116239e-07	0.004585	OK (low)	1.5	0.011810
25	0.049685	0.002886	0.002886	-0.002886	2.776401e-06	0.009814	OK (low)	1.5	0.017715

In [19]:

dynamics.get_diagnostic_rxn_data(rxn_index=0)   # diagnostic run data of the requested SINGLE reaction

Reaction:  A + B <-> C

Out[19]:

	START_TIME	Delta A	Delta B	Delta C	time_step	caption
0	0.000000	-9.840000	-9.840000	9.840000	0.004000	aborted: excessive norm value(s)
1	0.000000	-5.904000	-5.904000	5.904000	0.002400	aborted: excessive norm value(s)
2	0.000000	-3.542400	-3.542400	3.542400	0.001440	aborted: excessive norm value(s)
3	0.000000	-2.125440	-2.125440	2.125440	0.000864
4	0.000864	-1.272295	-1.272295	1.272295	0.000691
5	0.001555	-1.030999	-1.030999	1.030999	0.000691
6	0.002246	-0.843672	-0.843672	0.843672	0.000691
7	0.002938	-0.695848	-0.695848	0.695848	0.000691
8	0.003629	-0.866441	-0.866441	0.866441	0.001037
9	0.004666	-0.652653	-0.652653	0.652653	0.001037
10	0.005702	-0.496755	-0.496755	0.496755	0.001037
11	0.006739	-0.381056	-0.381056	0.381056	0.001037
12	0.007776	-0.294038	-0.294038	0.294038	0.001037
13	0.008813	-0.227921	-0.227921	0.227921	0.001037
14	0.009850	-0.177288	-0.177288	0.177288	0.001037
15	0.010886	-0.138275	-0.138275	0.138275	0.001037
16	0.011923	-0.162110	-0.162110	0.162110	0.001555
17	0.013478	-0.109377	-0.109377	0.109377	0.001555
18	0.015034	-0.074029	-0.074029	0.074029	0.001555
19	0.016589	-0.075315	-0.075315	0.075315	0.002333
20	0.018922	-0.039097	-0.039097	0.039097	0.002333
21	0.021254	-0.030522	-0.030522	0.030522	0.003499
22	0.024754	-0.012906	-0.012906	0.012906	0.005249
23	0.030002	0.001472	0.001472	-0.001472	0.007873
24	0.037876	-0.001355	-0.001355	0.001355	0.011810
25	0.049685	0.002886	0.002886	-0.002886	0.017715

In [20]:

dynamics.get_diagnostic_conc_data()   # diagnostic concentration data saved during the run, regardless of how much history we requested to save

Out[20]:

	TIME	A	B	C
0	0.000000	10.000000	50.000000	20.000000
1	0.000864	7.874560	47.874560	22.125440
2	0.001555	6.602265	46.602265	23.397735
3	0.002246	5.571266	45.571266	24.428734
4	0.002938	4.727594	44.727594	25.272406
5	0.003629	4.031746	44.031746	25.968254
6	0.004666	3.165305	43.165305	26.834695
7	0.005702	2.512652	42.512652	27.487348
8	0.006739	2.015898	42.015898	27.984102
9	0.007776	1.634842	41.634842	28.365158
10	0.008813	1.340804	41.340804	28.659196
11	0.009850	1.112883	41.112883	28.887117
12	0.010886	0.935595	40.935595	29.064405
13	0.011923	0.797321	40.797321	29.202679
14	0.013478	0.635211	40.635211	29.364789
15	0.015034	0.525833	40.525833	29.474167
16	0.016589	0.451805	40.451805	29.548195
17	0.018922	0.376490	40.376490	29.623510
18	0.021254	0.337393	40.337393	29.662607
19	0.024754	0.306871	40.306871	29.693129
20	0.030002	0.293965	40.293965	29.706035
21	0.037876	0.295437	40.295437	29.704563
22	0.049685	0.294082	40.294082	29.705918
23	0.067400	0.296969	40.296969	29.703031

In [ ]:

Association/Dissociation reaction A + B <-> C¶

with 1st-order kinetics for each species, taken to equilibrium.¶

Exploration of debugging and diagnostics options¶

Initialize the System¶

Sneak preview of eventual equilibrum:¶

Run the reaction¶

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

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

The argument argument variable_steps=True dynamically adjusts the initial_step (which is initially found to be too large, leading to some backtracking)¶

Plots changes of concentration with time¶

Check the final equilibrium¶

Everthing below is just for diagnostic insight¶

into the adaptive variable time steps¶

Association/Dissociation reaction `A + B <-> C`¶