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

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

This is a repeat of the experiment "react_2_a" , but with adaptive variable time steps

LAST REVISED: Dec. 3, 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_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_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_b.log.htm'

In [ ]:

PART 1 - RUN THE SIMULATION¶

Initialize the System¶

In [4]:

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

# Reaction A <-> B , with 1st-order kinetics in both directions
dynamics.add_reaction(reactants=["A"], products=["B"], 
                       forward_rate=3., reverse_rate=2.)

dynamics.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: A <-> B  (kF = 3 / kR = 2 / delta_G = -1,005.1 / K = 1.5) | 1st order in all reactants & products
Set of chemicals involved in the above reactions: {'B', 'A'}

In [5]:

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

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

In [6]:

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

dynamics.describe_state()

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

In [7]:

dynamics.get_history()

Out[7]:

	SYSTEM TIME	A	B	caption
0	0.0	10.0	50.0	Initial state

In [8]:

dynamics.set_diagnostics()      # To save diagnostic information about the call to single_compartment_react()
                                # Useful for insight into the inner workings of the simulation

In [9]:

# For experiment repeatability, we specify a particular group of preset parameters applicable to the adaptive time steps
dynamics.use_adaptive_preset(preset="mid")       # A "middle-of-the road" heuristic: somewhat "conservative" but not overly so

dynamics.show_adaptive_parameters()              # Details may vary across different versions of Life123

Parameters used for the automated adaptive time step sizes -
    THRESHOLDS:  [{'norm': 'norm_A', 'low': 0.5, 'high': 0.8, 'abort': 1.44}, {'norm': 'norm_B', 'low': 0.08, 'high': 0.5, 'abort': 1.5}]
    STEP FACTORS:  {'upshift': 1.2, 'downshift': 0.5, 'abort': 0.4}
    STEP FACTOR in case of 'ERROR ABORTS':  0.25

Run the reaction¶

In [10]:

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

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

The flag variable_steps automatically adjusts up or down the time step¶

In [11]:

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

Out[11]:

	SYSTEM TIME	A	B	caption
0	0.000000	10.000000	50.000000	Initial state
1	0.016000	11.120000	48.880000	1st reaction step
2	0.032000	12.150400	47.849600
3	0.048000	13.098368	46.901632
4	0.067200	14.144925	45.855075
5	0.086400	15.091012	44.908988
6	0.109440	16.117327	43.882673
7	0.132480	17.025411	42.974589
8	0.160128	17.989578	42.010422
9	0.193306	18.986635	41.013365
10	0.233119	19.984624	40.015376
11	0.280894	20.943812	39.056188
12	0.338225	21.819882	38.180118
13	0.407022	22.569810	37.430190
14	0.489579	23.160168	36.839832
15	0.588647	23.576169	36.423831
16	0.707528	23.828097	36.171903
17	0.850186	23.950713	36.049287
18	1.021375	23.992900	36.007100
19	1.226802	24.000193	35.999807	last reaction step

In [12]:

dynamics.explain_time_advance()

From time 0 to 0.048, in 3 steps of 0.016
From time 0.048 to 0.0864, in 2 steps of 0.0192
From time 0.0864 to 0.1325, in 2 steps of 0.023
From time 0.1325 to 0.1601, in 1 step of 0.0276
From time 0.1601 to 0.1933, in 1 step of 0.0332
From time 0.1933 to 0.2331, in 1 step of 0.0398
From time 0.2331 to 0.2809, in 1 step of 0.0478
From time 0.2809 to 0.3382, in 1 step of 0.0573
From time 0.3382 to 0.407, in 1 step of 0.0688
From time 0.407 to 0.4896, in 1 step of 0.0826
From time 0.4896 to 0.5886, in 1 step of 0.0991
From time 0.5886 to 0.7075, in 1 step of 0.119
From time 0.7075 to 0.8502, in 1 step of 0.143
From time 0.8502 to 1.021, in 1 step of 0.171
From time 1.021 to 1.227, in 1 step of 0.205
(19 steps total)

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

That resulted from passing the flag variable_steps=True to single_compartment_react()¶

In [ ]:

PART 2 - Scrutinizing some instances of step-size changes¶

The Delta-concentration values for all the individual reaction time steps, as contribued by a single reaction, may all be inspected at once from the diagnostic data:

In [13]:

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

Reaction:  A <-> B

Out[13]:

	START_TIME	Delta A	Delta B	time_step	caption
0	0.000000	7.000000	-7.000000	0.100000	aborted: excessive norm value(s)
1	0.000000	2.800000	-2.800000	0.040000	aborted: excessive norm value(s)
2	0.000000	1.120000	-1.120000	0.016000
3	0.016000	1.030400	-1.030400	0.016000
4	0.032000	0.947968	-0.947968	0.016000
5	0.048000	1.046557	-1.046557	0.019200
6	0.067200	0.946087	-0.946087	0.019200
7	0.086400	1.026315	-1.026315	0.023040
8	0.109440	0.908084	-0.908084	0.023040
9	0.132480	0.964167	-0.964167	0.027648
10	0.160128	0.997057	-0.997057	0.033178
11	0.193306	0.997988	-0.997988	0.039813
12	0.233119	0.959188	-0.959188	0.047776
13	0.280894	0.876070	-0.876070	0.057331
14	0.338225	0.749929	-0.749929	0.068797
15	0.407022	0.590357	-0.590357	0.082556
16	0.489579	0.416002	-0.416002	0.099068
17	0.588647	0.251928	-0.251928	0.118881
18	0.707528	0.122616	-0.122616	0.142658
19	0.850186	0.042187	-0.042187	0.171189
20	1.021375	0.007293	-0.007293	0.205427

Note that diagnostic data with the DELTA Concentrations - above and below - also record the values that were considered (but not actually used) during ABORTED steps¶

In the examples below, we'll re-compute Delta values for individual steps, directly from the system history.

Example 1: very early in the run¶

In [14]:

history[1:4]

Out[14]:

	SYSTEM TIME	A	B	caption
1	0.016	11.120000	48.880000	1st reaction step
2	0.032	12.150400	47.849600
3	0.048	13.098368	46.901632

In [15]:

delta_concentrations = dynamics.extract_delta_concentrations(history, 1, 2, ['A', 'B'])
delta_concentrations

Out[15]:

array([ 1.0304, -1.0304], dtype=float32)

As expected by the 1:1 stoichiometry, delta_A = - delta_B

In [16]:

# Get all the concentrations at the start of the above steps
baseline_conc = dynamics.get_historical_concentrations(row=1)
#dynamics.get_historical_concentrations(t=0.016)   # Alternate way
baseline_conc

Out[16]:

array([11.12, 48.88], dtype=float32)

In [17]:

# Computes some measures of how large delta_concentrations is, and propose a course of action
dynamics.adjust_speed(delta_conc=delta_concentrations, baseline_conc=baseline_conc)  

Out[17]:

('stay', 1, {'norm_A': 0.5308620929718018, 'norm_B': 0.09266187})

The above analysis indicates that the time step is just about right, and the simulations should STAY on that course : that's based on the shown computed norms (indicating the extent of the change taking place.)¶

Indeed, the simulator maintains the same time step :

In [18]:

original_step = history["SYSTEM TIME"][2] - history["SYSTEM TIME"][1]
original_step

Out[18]:

0.016000000000000004

In [19]:

next_step = history["SYSTEM TIME"][3] - history["SYSTEM TIME"][2]
next_step

Out[19]:

0.016000000000000007

In [20]:

next_step / original_step

Out[20]:

1.0000000000000002

Example 2: very late in the run¶

In [21]:

history[17:20]

Out[21]:

	SYSTEM TIME	A	B	caption
17	0.850186	23.950713	36.049287
18	1.021375	23.992900	36.007100
19	1.226802	24.000193	35.999807	last reaction step

In [22]:

delta_concentrations = dynamics.extract_delta_concentrations(history, 17, 18, ['A', 'B'])
delta_concentrations

Out[22]:

array([ 0.04218667, -0.04218667], dtype=float32)

Notice the far less change now that the system is approaching equilibrium¶

In [23]:

# Get all the concentrations at the start of the above steps
baseline_conc = dynamics.get_historical_concentrations(row=17)
#dynamics.get_historical_concentrations(t=0.850186)   # Alternate way
baseline_conc

Out[23]:

array([23.950714, 36.049286], dtype=float32)

In [24]:

# Computes a measure of how large delta_concentrations is, and propose a course of action
dynamics.adjust_speed(delta_conc=delta_concentrations, baseline_conc=baseline_conc)  

Out[24]:

('low', 1.2, {'norm_A': 0.0008898575906641781, 'norm_B': 0.0017613951})

The above analysis indicates that the time step is on the "LOW" side, and the simulations should increase it by a factor 1.2 : again, that's based on the shown computed norms (indicating the extent of the change taking place.)¶

Indeed, the simulator increases the time step x1.2:

In [25]:

original_step = history["SYSTEM TIME"][18] - history["SYSTEM TIME"][17]
original_step

Out[25]:

0.17118912860651525

In [26]:

next_step = history["SYSTEM TIME"][19] - history["SYSTEM TIME"][18]
next_step

Out[26]:

0.2054269543278182

In [27]:

next_step / original_step

Out[27]:

1.1999999999999995

Where does that x1.2 factor come from? It's one of the parameters that we passed to the simulator; they can be seen as follows:

In [28]:

dynamics.show_adaptive_parameters()

Parameters used for the automated adaptive time step sizes -
    THRESHOLDS:  [{'norm': 'norm_A', 'low': 0.5, 'high': 0.8, 'abort': 1.44}, {'norm': 'norm_B', 'low': 0.08, 'high': 0.5, 'abort': 1.5}]
    STEP FACTORS:  {'upshift': 1.2, 'downshift': 0.5, 'abort': 0.4}
    STEP FACTOR in case of 'ERROR ABORTS':  0.25

1.2 is stored as the "step factor" (for the time steps to take) in case an 'upshift' (in step size) is the decided course of action

In [ ]:

Diagnostics of the run may be investigated as follows:¶

(note - this is possible because we make a call to set_diagnostics() prior to running the simulation)

In [29]:

dynamics.get_diagnostic_conc_data()   # This will be complete, even if we only saved part of the history during the run

Out[29]:

	TIME	A	B
0	0.000000	10.000000	50.000000
1	0.016000	11.120000	48.880000
2	0.032000	12.150400	47.849600
3	0.048000	13.098368	46.901632
4	0.067200	14.144925	45.855075
5	0.086400	15.091012	44.908988
6	0.109440	16.117327	43.882673
7	0.132480	17.025411	42.974589
8	0.160128	17.989578	42.010422
9	0.193306	18.986635	41.013365
10	0.233119	19.984624	40.015376
11	0.280894	20.943812	39.056188
12	0.338225	21.819882	38.180118
13	0.407022	22.569810	37.430190
14	0.489579	23.160168	36.839832
15	0.588647	23.576169	36.423831
16	0.707528	23.828097	36.171903
17	0.850186	23.950713	36.049287
18	1.021375	23.992900	36.007100
19	1.226802	24.000193	35.999807

In [30]:

dynamics.get_diagnostic_decisions_data()

Out[30]:

	START_TIME	Delta A	Delta B	norm_A	norm_B	action	step_factor	time_step	caption
0	0.000000	7.000000	-7.000000	24.500000	NaN	ABORT	0.4	0.100000	excessive norm value(s)
1	0.000000	2.800000	-2.800000	3.920000	NaN	ABORT	0.4	0.040000	excessive norm value(s)
2	0.000000	1.120000	-1.120000	0.627200	0.112000	OK (stay)	1.0	0.016000
3	0.016000	1.030400	-1.030400	0.530862	0.092662	OK (stay)	1.0	0.016000
4	0.032000	0.947968	-0.947968	0.449322	0.078019	OK (low)	1.2	0.016000
5	0.048000	1.046557	-1.046557	0.547640	0.079900	OK (stay)	1.0	0.019200
6	0.067200	0.946087	-0.946087	0.447541	0.066885	OK (low)	1.2	0.019200
7	0.086400	1.026315	-1.026315	0.526662	0.068008	OK (stay)	1.0	0.023040
8	0.109440	0.908084	-0.908084	0.412308	0.056342	OK (low)	1.2	0.023040
9	0.132480	0.964167	-0.964167	0.464809	0.056631	OK (low)	1.2	0.027648
10	0.160128	0.997057	-0.997057	0.497061	0.055424	OK (low)	1.2	0.033178
11	0.193306	0.997988	-0.997988	0.497990	0.052563	OK (low)	1.2	0.039813
12	0.233119	0.959188	-0.959188	0.460021	0.047996	OK (low)	1.2	0.047776
13	0.280894	0.876070	-0.876070	0.383749	0.041830	OK (low)	1.2	0.057331
14	0.338225	0.749929	-0.749929	0.281197	0.034369	OK (low)	1.2	0.068797
15	0.407022	0.590357	-0.590357	0.174261	0.026157	OK (low)	1.2	0.082556
16	0.489579	0.416002	-0.416002	0.086529	0.017962	OK (low)	1.2	0.099068
17	0.588647	0.251928	-0.251928	0.031734	0.010686	OK (low)	1.2	0.118881
18	0.707528	0.122616	-0.122616	0.007517	0.005146	OK (low)	1.2	0.142658
19	0.850186	0.042187	-0.042187	0.000890	0.001761	OK (low)	1.2	0.171189
20	1.021375	0.007293	-0.007293	0.000027	0.000304	OK (low)	1.2	0.205427

In [ ]:

PART 3 - Analyze the reaction Dynamics¶

In [31]:

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

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

Out[31]:

True

Plots of changes of concentration with time¶

In [32]:

dynamics.plot_history(colors=['blue', 'orange'], show_intervals=True)

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

Compare the above with the fixed step sizes of experiment `react_2_a`¶

To see the sizes of the steps taken:

In [33]:

dynamics.plot_step_sizes(show_intervals=True)

In [ ]:

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

PART 1 - RUN THE SIMULATION¶

Initialize the System¶

Run the reaction¶

The flag variable_steps automatically adjusts up or down the time step¶

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

That resulted from passing the flag variable_steps=True to single_compartment_react()¶

PART 2 - Scrutinizing some instances of step-size changes¶

Note that diagnostic data with the DELTA Concentrations - above and below - also record the values that were considered (but not actually used) during ABORTED steps¶

Example 1: very early in the run¶

The above analysis indicates that the time step is just about right, and the simulations should STAY on that course : that's based on the shown computed norms (indicating the extent of the change taking place.)¶

Example 2: very late in the run¶

Notice the far less change now that the system is approaching equilibrium¶

The above analysis indicates that the time step is on the "LOW" side, and the simulations should increase it by a factor 1.2 : again, that's based on the shown computed norms (indicating the extent of the change taking place.)¶

Diagnostics of the run may be investigated as follows:¶

PART 3 - Analyze the reaction Dynamics¶

Plots of changes of concentration with time¶

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

Compare the above with the fixed step sizes of experiment react_2_a¶

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

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

Compare the above with the fixed step sizes of experiment `react_2_a`¶