Adaptive time substeps (variable time resolution) for reaction `A <-> B`,¶

with 1st-order kinetics in both directions, taken to equilibrium

Same as the experiment "react_1" , but with adaptive time substeps

LAST REVISED: Feb. 5, 2023

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_2.log.htm'

Initialize the System¶

Specify the chemicals and the reactions

In [4]:

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

# Reaction A <-> B , with 1st-order kinetics in both directions
chem_data.add_reaction(reactants=["A"], products=["B"], forward_rate=3., 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  (kF = 3 / kR = 2 / Delta_G = -1,005.13 / K = 1.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_2.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.], snapshot=True)

dynamics.describe_state()

SYSTEM STATE at Time t = 0:
2 species:
  Species 0 (A). Conc: 10.0
  Species 1 (B). Conc: 50.0

In [9]:

dynamics.history.get()

Out[9]:

	SYSTEM TIME	A	B	caption
0	0.0	10.0	50.0	Initial state

Run the reaction¶

In [10]:

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

dynamics.single_compartment_react(time_step=0.1, n_steps=11,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  dynamic_substeps=2, rel_fast_threshold=100)

single_compartment_react(): setting abs_fast_threshold to 10.0
12 total step(s) taken

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

In [11]:

df = dynamics.get_history()   # The system's history, saved during the run of single_compartment_react()
df

Out[11]:

	SYSTEM TIME	A	B	caption
0	0.00	10.000000	50.000000	Initial state
1	0.05	13.500000	46.500000	Interm. step, due to the fast rxns: [0]
2	0.10	16.125000	43.875000	1st reaction step
3	0.15	18.093750	41.906250	Interm. step, due to the fast rxns: [0]
4	0.20	19.570312	40.429688
5	0.25	20.677734	39.322266	Interm. step, due to the fast rxns: [0]
6	0.30	21.508301	38.491699
7	0.35	22.131226	37.868774	Interm. step, due to the fast rxns: [0]
8	0.40	22.598419	37.401581
9	0.50	23.299210	36.700790
10	0.60	23.649605	36.350395
11	0.70	23.824802	36.175198
12	0.80	23.912401	36.087599
13	0.90	23.956201	36.043799
14	1.00	23.978100	36.021900
15	1.10	23.989050	36.010950
16	1.20	23.994525	36.005475	last reaction step

In [12]:

dynamics.explain_time_advance()

From time 0 to 0.4, in 8 substeps of 0.05 (each 1/2 of full step)
From time 0.4 to 1.2, in 8 FULL steps of 0.1

Notice how the reaction proceeds in smaller steps (half steps) in the early times, when [A] and [B] are changing much more rapidly¶

That resulted from passing the argument dynamic_steps=2 to single_compartment_react()¶

For example, upon completing the half step to t=0.30, i.e. going from 0.25 to 0.30, the last change in [A] was (21.508301 - 20.677734) = 0.830567

The threshold we specified for a reaction to be considered fast is 100% per full step, for any of the involved chemicals. For a half step, that corresponds to 50%, i.e. 0.5
Since abs(0.830567) > 0.5 , the reaction is therefore marked "FAST" (as it has been so far), and the simulation then proceeds in a half step, to t=0.35

(Note: at t=0, in the absence of any simulation data, ALL reactions are assumed to be fast)

By contrast, upon completing the half step to t=0.40, i.e. going from 0.35 to 0.40, the following changes occur in [A] and [B]:

In [13]:

df.iloc[7]

Out[13]:

SYSTEM TIME                                       0.35
A                                            22.131226
B                                            37.868774
caption        Interm. step, due to the fast rxns: [0]
Name: 7, dtype: object

In [14]:

df.iloc[8]

Out[14]:

SYSTEM TIME          0.4
A              22.598419
B              37.401581
caption                 
Name: 8, dtype: object

In [15]:

s_0_35 = df.iloc[7][['A', 'B']].to_numpy()
s_0_35     # Concentrations of A and B at t=0.35

Out[15]:

array([22.1312255859375, 37.8687744140625], dtype=object)

In [16]:

s_0_40 = df.iloc[8][['A', 'B']].to_numpy()
s_0_40     # Concentrations of A and B at t=0.40

Out[16]:

array([22.598419189453125, 37.401580810546875], dtype=object)

In [17]:

(s_0_40 - s_0_35)

Out[17]:

array([0.467193603515625, -0.467193603515625], dtype=object)

Again, the threshold we specified for a reaction to be considered fast is 100% per full step, for any of the involved chemicals.
For a half step, that corresponds to 50%, i.e. 0.5
BOTH A's change of abs(0.46) AND B's change of abs(-0.46) are SMALLER than that.
The reaction is therefore marked "SLOW", and the simulation then proceeds in a full time step of 0.1 (i.e. a more relaxed time resolution), to t=0.50

Check the final equilibrium¶

In [18]:

dynamics.get_system_conc()

Out[18]:

array([23.99452507, 36.00547493])

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

[A] = 23.99316406

[B] = 36.00683594

The presence of equilibrium may be automatically checked withe the following handy function:

In [19]:

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

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

Out[19]:

True

Plots of changes of concentration with time¶

In [20]:

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

Note how the left-hand side of this plot is much smoother than it was in experiment "react_1", where no adaptive time substeps were used!¶

In [ ]:

Diagnostics of the run may be investigated as follows:¶

In [21]:

dynamics.diagnostic_data_baselines.get()

Out[21]:

	TIME	A	B	is_primary	primary_timestep	n_substeps	substep_number
0	0.00	10.000000	50.000000	True	0.1	2	0
1	0.05	13.500000	46.500000	False	0.1	2	1
2	0.10	16.125000	43.875000	True	0.1	2	0
3	0.15	18.093750	41.906250	False	0.1	2	1
4	0.20	19.570312	40.429688	True	0.1	2	0
5	0.25	20.677734	39.322266	False	0.1	2	1
6	0.30	21.508301	38.491699	True	0.1	2	0
7	0.35	22.131226	37.868774	False	0.1	2	1
8	0.40	22.598419	37.401581	True	0.1	2	0
9	0.50	23.299210	36.700790	True	0.1	2	0
10	0.60	23.649605	36.350395	True	0.1	2	0
11	0.70	23.824802	36.175198	True	0.1	2	0
12	0.80	23.912401	36.087599	True	0.1	2	0
13	0.90	23.956201	36.043799	True	0.1	2	0
14	1.00	23.978100	36.021900	True	0.1	2	0
15	1.10	23.989050	36.010950	True	0.1	2	0
16	1.20	23.994525	36.005475	True	0.1	2	0

In [22]:

dynamics.get_diagnostic_data(rxn_index=0)      # For the 0-th reaction (the only reaction in our case)

Reaction:  A <-> B

Out[22]:

	TIME	Delta A	Delta B	substep	time_subdivision	delta_time
0	0.00	3.500000	-3.500000	0	2	0.05
1	0.05	2.625000	-2.625000	1	2	0.05
2	0.10	1.968750	-1.968750	0	2	0.05
3	0.15	1.476562	-1.476562	1	2	0.05
4	0.20	1.107422	-1.107422	0	2	0.05
5	0.25	0.830566	-0.830566	1	2	0.05
6	0.30	0.622925	-0.622925	0	2	0.05
7	0.35	0.467194	-0.467194	1	2	0.05
8	0.40	0.700790	-0.700790	0	1	0.10
9	0.50	0.350395	-0.350395	0	1	0.10
10	0.60	0.175198	-0.175198	0	1	0.10
11	0.70	0.087599	-0.087599	0	1	0.10
12	0.80	0.043799	-0.043799	0	1	0.10
13	0.90	0.021900	-0.021900	0	1	0.10
14	1.00	0.010950	-0.010950	0	1	0.10
15	1.10	0.005475	-0.005475	0	1	0.10