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 22, 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 ReactionData 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_3.log.htm'

Initialize the System¶

Specify the chemicals and the reactions

In [4]:

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

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

chem_data.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: A + B <-> C  (kF = 5 / kR = 2 / Delta_G = -2,271.45 / K = 2.5) | 1st order in all reactants & products

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_3.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([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

In [9]:

dynamics.get_history()

Out[9]:

	SYSTEM TIME	A	B	C	caption
0	0.0	10.0	50.0	20.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=0.5, high=0.8, abort=1.44)
dynamics.set_thresholds(norm="norm_B", low=0.08, high=0.5, abort=1.5)
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.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.5 (set to 0.002) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.002) 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.001) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.001) 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.0005) [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
39 total step(s) taken

In [11]:

dynamics.get_history()

Out[11]:

	SYSTEM TIME	A	B	C	caption
0	0.000000	10.000000	50.000000	20.000000	Initial state
1	0.000500	8.770000	48.770000	21.230000	1st reaction step
2	0.001000	7.721948	47.721948	22.278052
3	0.001500	6.822960	46.822960	23.177040
4	0.002000	6.047459	46.047459	23.952541
5	0.002500	5.375236	45.375236	24.624764
6	0.003000	4.790104	44.790104	25.209896
7	0.003500	4.278941	44.278941	25.721059
8	0.004000	3.830995	43.830995	26.169005
9	0.004500	3.437373	43.437373	26.562627
10	0.005000	3.090659	43.090659	26.909341
11	0.005500	2.784622	42.784622	27.215378
12	0.006000	2.513990	42.513990	27.486010
13	0.006500	2.274277	42.274277	27.725723
14	0.007000	2.061644	42.061644	27.938356
15	0.007500	1.872792	41.872792	28.127208
16	0.008000	1.704872	41.704872	28.295128
17	0.008500	1.555413	41.555413	28.444587
18	0.009000	1.422268	41.422268	28.577732
19	0.009500	1.303562	41.303562	28.696438
20	0.010000	1.197654	41.197654	28.802346
21	0.010500	1.103105	41.103105	28.896895
22	0.011250	0.976421	40.976421	29.023579
23	0.012000	0.869918	40.869918	29.130082
24	0.012750	0.780288	40.780288	29.219712
25	0.013500	0.704791	40.704791	29.295209
26	0.014250	0.641153	40.641153	29.358847
27	0.015000	0.587476	40.587476	29.412524
28	0.015750	0.542179	40.542179	29.457821
29	0.016875	0.484816	40.484816	29.515184
30	0.018000	0.440819	40.440819	29.559181
31	0.019125	0.407050	40.407050	29.592950
32	0.020813	0.368149	40.368149	29.631851
33	0.022500	0.342763	40.342763	29.657237
34	0.025031	0.317892	40.317892	29.682108
35	0.028828	0.299973	40.299973	29.700027
36	0.034523	0.294024	40.294024	29.705976
37	0.043066	0.295518	40.295518	29.704482
38	0.055881	0.293836	40.293836	29.706164
39	0.075103	0.297939	40.297939	29.702061	last reaction step

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

In [12]:

dynamics.explain_time_advance()

From time 0 to 0.0105, in 21 steps of 0.0005
From time 0.0105 to 0.01575, in 7 steps of 0.00075
From time 0.01575 to 0.01913, in 3 steps of 0.00113
From time 0.01913 to 0.0225, in 2 steps of 0.00169
From time 0.0225 to 0.02503, in 1 step of 0.00253
From time 0.02503 to 0.02883, in 1 step of 0.0038
From time 0.02883 to 0.03452, in 1 step of 0.0057
From time 0.03452 to 0.04307, in 1 step of 0.00854
From time 0.04307 to 0.05588, in 1 step of 0.0128
From time 0.05588 to 0.0751, in 1 step of 0.0192
(39 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 [13]:

dynamics.plot_step_sizes(show_intervals=True)

Plots of changes of concentration with time¶

In [14]:

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

Check the final equilibrium¶

In [15]:

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

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

Out[15]:

True

In [ ]:

Everthing below is just for diagnostic insight¶

into the adaptive variable time steps¶

In [16]:

dynamics.get_diagnostic_decisions_data()

Out[16]:

	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.5	0.004000	excessive norm value(s)
1	0.000000	-4.920000	-4.920000	4.920000	8.068800e+00	NaN	ABORT	0.5	0.002000	excessive norm value(s)
2	0.000000	-2.460000	-2.460000	2.460000	2.017200e+00	NaN	ABORT	0.5	0.001000	excessive norm value(s)
3	0.000000	-1.230000	-1.230000	1.230000	5.043000e-01	0.123000	OK (stay)	1.0	0.000500
4	0.000500	-1.048052	-1.048052	1.048052	3.661378e-01	0.119504	OK (stay)	1.0	0.000500
5	0.001000	-0.898988	-0.898988	0.898988	2.693931e-01	0.116420	OK (stay)	1.0	0.000500
6	0.001500	-0.775501	-0.775501	0.775501	2.004672e-01	0.113660	OK (stay)	1.0	0.000500
7	0.002000	-0.672223	-0.672223	0.672223	1.506278e-01	0.111158	OK (stay)	1.0	0.000500
8	0.002500	-0.585132	-0.585132	0.585132	1.141264e-01	0.108857	OK (stay)	1.0	0.000500
9	0.003000	-0.511163	-0.511163	0.511163	8.709597e-02	0.106712	OK (stay)	1.0	0.000500
10	0.003500	-0.447946	-0.447946	0.447946	6.688532e-02	0.104686	OK (stay)	1.0	0.000500
11	0.004000	-0.393622	-0.393622	0.393622	5.164603e-02	0.102747	OK (stay)	1.0	0.000500
12	0.004500	-0.346714	-0.346714	0.346714	4.007008e-02	0.100866	OK (stay)	1.0	0.000500
13	0.005000	-0.306037	-0.306037	0.306037	3.121956e-02	0.099020	OK (stay)	1.0	0.000500
14	0.005500	-0.270632	-0.270632	0.270632	2.441392e-02	0.097188	OK (stay)	1.0	0.000500
15	0.006000	-0.239713	-0.239713	0.239713	1.915417e-02	0.095352	OK (stay)	1.0	0.000500
16	0.006500	-0.212633	-0.212633	0.212633	1.507090e-02	0.093495	OK (stay)	1.0	0.000500
17	0.007000	-0.188852	-0.188852	0.188852	1.188836e-02	0.091603	OK (stay)	1.0	0.000500
18	0.007500	-0.167920	-0.167920	0.167920	9.399084e-03	0.089663	OK (stay)	1.0	0.000500
19	0.008000	-0.149459	-0.149459	0.149459	7.445949e-03	0.087666	OK (stay)	1.0	0.000500
20	0.008500	-0.133145	-0.133145	0.133145	5.909198e-03	0.085601	OK (stay)	1.0	0.000500
21	0.009000	-0.118706	-0.118706	0.118706	4.697054e-03	0.083463	OK (stay)	1.0	0.000500
22	0.009500	-0.105908	-0.105908	0.105908	3.738831e-03	0.081245	OK (stay)	1.0	0.000500
23	0.010000	-0.094549	-0.094549	0.094549	2.979837e-03	0.078945	OK (low)	1.5	0.000500
24	0.010500	-0.126684	-0.126684	0.126684	5.349575e-03	0.114843	OK (stay)	1.0	0.000750
25	0.011250	-0.106503	-0.106503	0.106503	3.780970e-03	0.109075	OK (stay)	1.0	0.000750
26	0.012000	-0.089630	-0.089630	0.089630	2.677874e-03	0.103033	OK (stay)	1.0	0.000750
27	0.012750	-0.075497	-0.075497	0.075497	1.899922e-03	0.096755	OK (stay)	1.0	0.000750
28	0.013500	-0.063639	-0.063639	0.063639	1.349957e-03	0.090294	OK (stay)	1.0	0.000750
29	0.014250	-0.053676	-0.053676	0.053676	9.603761e-04	0.083718	OK (stay)	1.0	0.000750
30	0.015000	-0.045297	-0.045297	0.045297	6.839364e-04	0.077104	OK (low)	1.5	0.000750
31	0.015750	-0.057364	-0.057364	0.057364	1.096869e-03	0.105802	OK (stay)	1.0	0.001125
32	0.016875	-0.043996	-0.043996	0.043996	6.452305e-04	0.090749	OK (stay)	1.0	0.001125
33	0.018000	-0.033769	-0.033769	0.033769	3.801200e-04	0.076606	OK (low)	1.5	0.001125
34	0.019125	-0.038901	-0.038901	0.038901	5.044334e-04	0.095569	OK (stay)	1.0	0.001688
35	0.020813	-0.025386	-0.025386	0.025386	2.148211e-04	0.068957	OK (low)	1.5	0.001688
36	0.022500	-0.024871	-0.024871	0.024871	2.061837e-04	0.072560	OK (low)	1.5	0.002531
37	0.025031	-0.017919	-0.017919	0.017919	1.070309e-04	0.056368	OK (low)	1.5	0.003797
38	0.028828	-0.005948	-0.005948	0.005948	1.179398e-05	0.019829	OK (low)	1.5	0.005695
39	0.034523	0.001493	0.001493	-0.001493	7.433416e-07	0.005079	OK (low)	1.5	0.008543
40	0.043066	-0.001682	-0.001682	0.001682	9.429553e-07	0.005691	OK (low)	1.5	0.012814
41	0.055881	0.004103	0.004103	-0.004103	5.611315e-06	0.013963	OK (low)	1.5	0.019222

In [17]:

dynamics.get_diagnostic_rxn_data(rxn_index=0)

Reaction:  A + B <-> C

Out[17]:

	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	-4.920000	-4.920000	4.920000	0.002000	aborted: excessive norm value(s)
2	0.000000	-2.460000	-2.460000	2.460000	0.001000	aborted: excessive norm value(s)
3	0.000000	-1.230000	-1.230000	1.230000	0.000500
4	0.000500	-1.048052	-1.048052	1.048052	0.000500
5	0.001000	-0.898988	-0.898988	0.898988	0.000500
6	0.001500	-0.775501	-0.775501	0.775501	0.000500
7	0.002000	-0.672223	-0.672223	0.672223	0.000500
8	0.002500	-0.585132	-0.585132	0.585132	0.000500
9	0.003000	-0.511163	-0.511163	0.511163	0.000500
10	0.003500	-0.447946	-0.447946	0.447946	0.000500
11	0.004000	-0.393622	-0.393622	0.393622	0.000500
12	0.004500	-0.346714	-0.346714	0.346714	0.000500
13	0.005000	-0.306037	-0.306037	0.306037	0.000500
14	0.005500	-0.270632	-0.270632	0.270632	0.000500
15	0.006000	-0.239713	-0.239713	0.239713	0.000500
16	0.006500	-0.212633	-0.212633	0.212633	0.000500
17	0.007000	-0.188852	-0.188852	0.188852	0.000500
18	0.007500	-0.167920	-0.167920	0.167920	0.000500
19	0.008000	-0.149459	-0.149459	0.149459	0.000500
20	0.008500	-0.133145	-0.133145	0.133145	0.000500
21	0.009000	-0.118706	-0.118706	0.118706	0.000500
22	0.009500	-0.105908	-0.105908	0.105908	0.000500
23	0.010000	-0.094549	-0.094549	0.094549	0.000500
24	0.010500	-0.126684	-0.126684	0.126684	0.000750
25	0.011250	-0.106503	-0.106503	0.106503	0.000750
26	0.012000	-0.089630	-0.089630	0.089630	0.000750
27	0.012750	-0.075497	-0.075497	0.075497	0.000750
28	0.013500	-0.063639	-0.063639	0.063639	0.000750
29	0.014250	-0.053676	-0.053676	0.053676	0.000750
30	0.015000	-0.045297	-0.045297	0.045297	0.000750
31	0.015750	-0.057364	-0.057364	0.057364	0.001125
32	0.016875	-0.043996	-0.043996	0.043996	0.001125
33	0.018000	-0.033769	-0.033769	0.033769	0.001125
34	0.019125	-0.038901	-0.038901	0.038901	0.001688
35	0.020813	-0.025386	-0.025386	0.025386	0.001688
36	0.022500	-0.024871	-0.024871	0.024871	0.002531
37	0.025031	-0.017919	-0.017919	0.017919	0.003797
38	0.028828	-0.005948	-0.005948	0.005948	0.005695
39	0.034523	0.001493	0.001493	-0.001493	0.008543
40	0.043066	-0.001682	-0.001682	0.001682	0.012814
41	0.055881	0.004103	0.004103	-0.004103	0.019222

In [18]:

dynamics.get_diagnostic_conc_data()

Out[18]:

	TIME	A	B	C
0	0.000000	10.000000	50.000000	20.000000
1	0.000500	8.770000	48.770000	21.230000
2	0.001000	7.721948	47.721948	22.278052
3	0.001500	6.822960	46.822960	23.177040
4	0.002000	6.047459	46.047459	23.952541
5	0.002500	5.375236	45.375236	24.624764
6	0.003000	4.790104	44.790104	25.209896
7	0.003500	4.278941	44.278941	25.721059
8	0.004000	3.830995	43.830995	26.169005
9	0.004500	3.437373	43.437373	26.562627
10	0.005000	3.090659	43.090659	26.909341
11	0.005500	2.784622	42.784622	27.215378
12	0.006000	2.513990	42.513990	27.486010
13	0.006500	2.274277	42.274277	27.725723
14	0.007000	2.061644	42.061644	27.938356
15	0.007500	1.872792	41.872792	28.127208
16	0.008000	1.704872	41.704872	28.295128
17	0.008500	1.555413	41.555413	28.444587
18	0.009000	1.422268	41.422268	28.577732
19	0.009500	1.303562	41.303562	28.696438
20	0.010000	1.197654	41.197654	28.802346
21	0.010500	1.103105	41.103105	28.896895
22	0.011250	0.976421	40.976421	29.023579
23	0.012000	0.869918	40.869918	29.130082
24	0.012750	0.780288	40.780288	29.219712
25	0.013500	0.704791	40.704791	29.295209
26	0.014250	0.641153	40.641153	29.358847
27	0.015000	0.587476	40.587476	29.412524
28	0.015750	0.542179	40.542179	29.457821
29	0.016875	0.484816	40.484816	29.515184
30	0.018000	0.440819	40.440819	29.559181
31	0.019125	0.407050	40.407050	29.592950
32	0.020813	0.368149	40.368149	29.631851
33	0.022500	0.342763	40.342763	29.657237
34	0.025031	0.317892	40.317892	29.682108
35	0.028828	0.299973	40.299973	29.700027
36	0.034523	0.294024	40.294024	29.705976
37	0.043066	0.295518	40.295518	29.704482
38	0.055881	0.293836	40.293836	29.706164
39	0.075103	0.297939	40.297939	29.702061