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 substeps are used)

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

For the 0-th reaction (the only reaction in our case)¶

In [1]:

# Extend the sys.path variable, to contain the project's root directory
import set_path
set_path.add_ancestor_dir_to_syspath(2)  # The number of levels to go up 
                                         # to reach the project's home, from the folder containing this notebook

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.)

print("Number of reactions: ", chem_data.number_of_reactions())

Number of reactions:  1

In [5]:

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 [6]:

# 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 [7]:

dynamics = ReactionDynamics(reaction_data=chem_data)

In [8]:

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

In [9]:

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 [10]:

dynamics.get_history()

Out[10]:

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

Run the reaction¶

In [11]:

dynamics.set_diagnostics()       # To save diagnostic information about the call to single_compartment_react()
#dynamics.verbose_list = [1]      # Uncomment for detailed run information (meant for debugging the adaptive variable time step)

dynamics.single_compartment_react(time_step=0.004, reaction_duration=0.06,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  dynamic_substeps=2, rel_fast_threshold=1)

single_compartment_react(): setting abs_fast_threshold to 2.5
15 total step(s) taken

Note: the argument dynamic_step=2 splits the time steps in 2 whenever the reaction is "fast" (as determined using the given value of fast_threshold )¶

In [12]:

df = dynamics.get_history()
df

Out[12]:

	SYSTEM TIME	A	B	C	caption
0	0.000	10.000000	50.000000	20.000000	Initial state
1	0.002	5.080000	45.080000	24.920000	Interm. step, due to the fast rxns: [0]
2	0.004	2.889616	42.889616	27.110384	1st reaction step
3	0.006	1.758712	41.758712	28.241288	Interm. step, due to the fast rxns: [0]
4	0.008	1.137262	41.137262	28.862738
5	0.010	0.784874	40.784874	29.215126	Interm. step, due to the fast rxns: [0]
6	0.012	0.581625	40.581625	29.418375
7	0.014	0.463266	40.463266	29.536734	Interm. step, due to the fast rxns: [0]
8	0.016	0.393960	40.393960	29.606040
9	0.018	0.353248	40.353248	29.646752	Interm. step, due to the fast rxns: [0]
10	0.020	0.329288	40.329288	29.670712
11	0.022	0.315171	40.315171	29.684829	Interm. step, due to the fast rxns: [0]
12	0.024	0.306849	40.306849	29.693151
13	0.026	0.301940	40.301940	29.698060	Interm. step, due to the fast rxns: [0]
14	0.028	0.299045	40.299045	29.700955
15	0.032	0.295628	40.295628	29.704372
16	0.036	0.295013	40.295013	29.704987
17	0.040	0.294902	40.294902	29.705098
18	0.044	0.294882	40.294882	29.705118
19	0.048	0.294878	40.294878	29.705122
20	0.052	0.294878	40.294878	29.705122
21	0.056	0.294877	40.294877	29.705123
22	0.060	0.294877	40.294877	29.705123	last reaction step

In [13]:

dynamics.explain_time_advance()

From time 0 to 0.028, in 14 substeps of 0.002 (each 1/2 of full step)
From time 0.028 to 0.06, in 8 FULL steps of 0.004

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

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

Check the final equilibrium¶

In [14]:

dynamics.get_system_conc()

Out[14]:

array([ 0.29487741, 40.29487741, 29.70512259])

NOTE: Consistent with the 3/2 ratio of forward/reverse rates (and the 1st order reactions),
the systems settles in the following equilibrium:

[A] = 0.29487741 , [B] = 40.29487741 , [C] = 29.70512259

In [15]:

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

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

Out[15]:

True

Plots of changes of concentration with time¶

In [16]:

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

In [ ]:

Everthing below is just for diagnostic insight¶

into the adaptive variable time steps¶

WARNING: The explanations below are based on an older system of using RELATIVE concentration changes, and no longer applicable to the newer approach; it'll be corrected in future versions (COMMENTED OUT FOR NOW)¶

In [17]:

# This approach, from the run data, is only usable with single-reaction runs
#dynamics.examine_run(df=df, time_step=0.004, rel_fast_threshold=1)
# The arguments step MUST match the value used in call to single_compartment_react()

Take a peek at internal diagnostic data from the reactions¶

In [18]:

# This approach, from internal diagnostic data, 
# is more generally applicable also to runs with multiple reactions
#dynamics.diagnose_variable_time_steps()

The above diagnostics are based on the following diagnostic data:¶

In [19]:

dynamics.get_diagnostic_data(rxn_index=0)

Reaction:  A + B <-> C

Out[19]:

	TIME	Delta A	Delta B	Delta C	substep	time_subdivision	delta_time
0	0.000	-4.920000e+00	-4.920000e+00	4.920000e+00	0	2	0.002
1	0.002	-2.190384e+00	-2.190384e+00	2.190384e+00	1	2	0.002
2	0.004	-1.130904e+00	-1.130904e+00	1.130904e+00	0	2	0.002
3	0.006	-6.214505e-01	-6.214505e-01	6.214505e-01	1	2	0.002
4	0.008	-3.523874e-01	-3.523874e-01	3.523874e-01	0	2	0.002
5	0.010	-2.032495e-01	-2.032495e-01	2.032495e-01	1	2	0.002
6	0.012	-1.183593e-01	-1.183593e-01	1.183593e-01	0	2	0.002
7	0.014	-6.930543e-02	-6.930543e-02	6.930543e-02	1	2	0.002
8	0.016	-4.071193e-02	-4.071193e-02	4.071193e-02	0	2	0.002
9	0.018	-2.396011e-02	-2.396011e-02	2.396011e-02	1	2	0.002
10	0.020	-1.411669e-02	-1.411669e-02	1.411669e-02	0	2	0.002
11	0.022	-8.322570e-03	-8.322570e-03	8.322570e-03	1	2	0.002
12	0.024	-4.908484e-03	-4.908484e-03	4.908484e-03	0	2	0.002
13	0.026	-2.895574e-03	-2.895574e-03	2.895574e-03	1	2	0.002
14	0.028	-3.416720e-03	-3.416720e-03	3.416720e-03	0	1	0.004
15	0.032	-6.153737e-04	-6.153737e-04	6.153737e-04	0	1	0.004
16	0.036	-1.108825e-04	-1.108825e-04	1.108825e-04	0	1	0.004
17	0.040	-1.998121e-05	-1.998121e-05	1.998121e-05	0	1	0.004
18	0.044	-3.600700e-06	-3.600700e-06	3.600700e-06	0	1	0.004
19	0.048	-6.488634e-07	-6.488634e-07	6.488634e-07	0	1	0.004
20	0.052	-1.169284e-07	-1.169284e-07	1.169284e-07	0	1	0.004
21	0.056	-2.107106e-08	-2.107106e-08	2.107106e-08	0	1	0.004