Adaptive time steps (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 variable time steps

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 pandas as pd
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.)

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

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

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

# 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")  # This has the effect of turning off use of "norm_B"
dynamics.set_step_factors(upshift=2.0, downshift=0.5, abort=0.5)
dynamics.set_error_step_factor(0.5)

In [11]:

dynamics.thresholds

Out[11]:

[{'norm': 'norm_A', 'low': 0.5, 'high': 0.8, 'abort': 1.44}]

Note how we UNSET the defaults for "norm_B" (i.e. it won't be used)¶

In [12]:

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.5 (set to 0.05) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.05) 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.025) [Step started at t=0, and will rewind there]
INFO: the tentative time step (0.025) 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.0125) [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
16 total step(s) taken

The flag variable_steps automatically adjusts up or down the time step, whenever the changes of concentrations are, respectively, "slow" or "fast" (as determined using the specified thresholds )¶

In [13]:

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

Out[13]:

	SYSTEM TIME	A	B	caption
0	0.0000	10.000000	50.000000	Initial state
1	0.0125	10.875000	49.125000	1st reaction step
2	0.0375	12.515625	47.484375
3	0.0500	13.233398	46.766602
4	0.0750	14.579224	45.420776
5	0.0875	15.168022	44.831978
6	0.1125	16.272019	43.727981
7	0.1375	17.238017	42.761983
8	0.1875	18.928513	41.071487
9	0.2125	19.562449	40.437551
10	0.2625	20.671836	39.328164
11	0.3125	21.503877	38.496123
12	0.4125	22.751939	37.248061
13	0.5125	23.375969	36.624031
14	0.7125	24.000000	36.000000
15	1.1125	24.000000	36.000000
16	1.9125	24.000000	36.000000	last reaction step

In [14]:

dynamics.explain_time_advance()

From time 0 to 0.0125, in 1 step of 0.0125
From time 0.0125 to 0.0375, in 1 step of 0.025
From time 0.0375 to 0.05, in 1 step of 0.0125
From time 0.05 to 0.075, in 1 step of 0.025
From time 0.075 to 0.0875, in 1 step of 0.0125
From time 0.0875 to 0.1375, in 2 steps of 0.025
From time 0.1375 to 0.1875, in 1 step of 0.05
From time 0.1875 to 0.2125, in 1 step of 0.025
From time 0.2125 to 0.3125, in 2 steps of 0.05
From time 0.3125 to 0.5125, in 2 steps of 0.1
From time 0.5125 to 0.7125, in 1 step of 0.2
From time 0.7125 to 1.112, in 1 step of 0.4
From time 1.112 to 1.912, in 1 step of 0.8
(16 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()¶

Scrutinizing some instances of step-size changes:¶

Detailed Example 1: going from 0.1375 to 0.1875¶

In [15]:

lookup = dynamics.get_history(t_start=0.1375, t_end=0.1875)
lookup

Out[15]:

	SYSTEM TIME	A	B	caption
7	0.1375	17.238017	42.761983
8	0.1875	18.928513	41.071487

In [16]:

delta_concentrations = dynamics.extract_delta_concentrations(lookup, 7, 8, ['A', 'B'])
delta_concentrations

Out[16]:

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

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

The above values coud also be looked up from the diagnostic data, since we only have 1 reaction:

In [17]:

rxn_data = dynamics.get_diagnostic_rxn_data(rxn_index=0)

Reaction:  A <-> B

In [18]:

rxn_data[8:12]

Out[18]:

	START_TIME	Delta A	Delta B	time_step
8	0.0875	1.103997	-1.103997	0.025
9	0.1125	0.965998	-0.965998	0.025
10	0.1375	1.690496	-1.690496	0.050
11	0.1875	0.633936	-0.633936	0.025

In [19]:

delta_row = dynamics.get_diagnostic_rxn_data(rxn_index=0, t=0.1375) # Locate the row with the interval's start time
delta_row

Reaction:  A <-> B

Out[19]:

	search_value	START_TIME	Delta A	Delta B	time_step	caption
10	0.1375	0.1375	1.690496	-1.690496	0.05

In [20]:

delta_row[["Delta A", "Delta B"]].to_numpy()   # Gives same value as delta_concentrations, above

Out[20]:

array([[ 1.69049576, -1.69049576]])

In [21]:

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

Out[21]:

('high', 0.5, {'norm_A': 1.428887963294983})

The above conclusion is that the time step is on the "high" side, and should be HALVED at the next round : that's because the computed norm is > than the "high" value previously given in the argument to set_thresholds() (but doesn't exceed the "abort" threshold)¶

In [22]:

dynamics.thresholds   # Consult the previously-set threshold values

Out[22]:

[{'norm': 'norm_A', 'low': 0.5, 'high': 0.8, 'abort': 1.44}]

Detailed Example 2: going from 0.1875 to 0.2125¶

In [23]:

lookup = dynamics.get_history(t_start=0.1875, t_end=0.2125)
lookup

Out[23]:

	SYSTEM TIME	A	B	caption
8	0.1875	18.928513	41.071487
9	0.2125	19.562449	40.437551

In [24]:

delta_concentrations = dynamics.extract_delta_concentrations(lookup, 8, 9, ['A', 'B'])
delta_concentrations

Out[24]:

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

Note how substantially smaller delta_concentrations is, compared to the previous example

In [25]:

dynamics.adjust_speed(delta_concentrations)

Out[25]:

('low', 2.0, {'norm_A': 0.20093737542629242})

The above conclusion is that the time step is on the "low" side, and should be DOUBLED at the next round : that's because the computed norm is < than the "low" value previously given in the argument to set_thresholds()¶

In [26]:

dynamics.thresholds   # Consult the previously-set threshold values

Out[26]:

[{'norm': 'norm_A', 'low': 0.5, 'high': 0.8, 'abort': 1.44}]

Check the final equilibrium¶

In [27]:

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

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

Out[27]:

True

Plots of changes of concentration with time¶

In [28]:

dynamics.plot_curves(colors=['blue', 'orange'])

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

In [29]:

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

Compare the above with the fixed step sizes (always 0.1) of experiment `react_1`¶

In [30]:

dynamics.plot_step_sizes(show_intervals=True)

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

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

Out[31]:

	TIME	A	B
0	0.0000	10.000000	50.000000
1	0.0125	10.875000	49.125000
2	0.0375	12.515625	47.484375
3	0.0500	13.233398	46.766602
4	0.0750	14.579224	45.420776
5	0.0875	15.168022	44.831978
6	0.1125	16.272019	43.727981
7	0.1375	17.238017	42.761983
8	0.1875	18.928513	41.071487
9	0.2125	19.562449	40.437551
10	0.2625	20.671836	39.328164
11	0.3125	21.503877	38.496123
12	0.4125	22.751939	37.248061
13	0.5125	23.375969	36.624031
14	0.7125	24.000000	36.000000
15	1.1125	24.000000	36.000000
16	1.9125	24.000000	36.000000

In [32]:

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

Reaction:  A <-> B

Out[32]:

	START_TIME	Delta A	Delta B	time_step	caption
0	0.0000	7.000000	-7.000000	0.1000	aborted: excessive norm value(s)
1	0.0000	3.500000	-3.500000	0.0500	aborted: excessive norm value(s)
2	0.0000	1.750000	-1.750000	0.0250	aborted: excessive norm value(s)
3	0.0000	0.875000	-0.875000	0.0125
4	0.0125	1.640625	-1.640625	0.0250
5	0.0375	0.717773	-0.717773	0.0125
6	0.0500	1.345825	-1.345825	0.0250
7	0.0750	0.588799	-0.588799	0.0125
8	0.0875	1.103997	-1.103997	0.0250
9	0.1125	0.965998	-0.965998	0.0250
10	0.1375	1.690496	-1.690496	0.0500
11	0.1875	0.633936	-0.633936	0.0250
12	0.2125	1.109388	-1.109388	0.0500
13	0.2625	0.832041	-0.832041	0.0500
14	0.3125	1.248061	-1.248061	0.1000
15	0.4125	0.624031	-0.624031	0.1000
16	0.5125	0.624031	-0.624031	0.2000
17	0.7125	0.000000	0.000000	0.4000
18	1.1125	0.000000	0.000000	0.8000

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

dynamics.get_diagnostic_decisions_data()

Out[33]:

	START_TIME	Delta A	Delta B	norm_A	norm_B	action	step_factor	time_step	caption
0	0.0000	7.000000	-7.000000	24.500000	None	ABORT	0.5	0.1000	excessive norm value(s)
1	0.0000	3.500000	-3.500000	6.125000	None	ABORT	0.5	0.0500	excessive norm value(s)
2	0.0000	1.750000	-1.750000	1.531250	None	ABORT	0.5	0.0250	excessive norm value(s)
3	0.0000	0.875000	-0.875000	0.382812	None	OK (low)	2.0	0.0125
4	0.0125	1.640625	-1.640625	1.345825	None	OK (high)	0.5	0.0250
5	0.0375	0.717773	-0.717773	0.257599	None	OK (low)	2.0	0.0125
6	0.0500	1.345825	-1.345825	0.905623	None	OK (high)	0.5	0.0250
7	0.0750	0.588799	-0.588799	0.173342	None	OK (low)	2.0	0.0125
8	0.0875	1.103997	-1.103997	0.609405	None	OK (stay)	1.0	0.0250
9	0.1125	0.965998	-0.965998	0.466576	None	OK (low)	2.0	0.0250
10	0.1375	1.690496	-1.690496	1.428888	None	OK (high)	0.5	0.0500
11	0.1875	0.633936	-0.633936	0.200937	None	OK (low)	2.0	0.0250
12	0.2125	1.109388	-1.109388	0.615371	None	OK (stay)	1.0	0.0500
13	0.2625	0.832041	-0.832041	0.346146	None	OK (low)	2.0	0.0500
14	0.3125	1.248061	-1.248061	0.778829	None	OK (stay)	1.0	0.1000
15	0.4125	0.624031	-0.624031	0.194707	None	OK (low)	2.0	0.1000
16	0.5125	0.624031	-0.624031	0.194707	None	OK (low)	2.0	0.2000
17	0.7125	0.000000	0.000000	0.000000	None	OK (low)	2.0	0.4000
18	1.1125	0.000000	0.000000	0.000000	None	OK (low)	2.0	0.8000

In [ ]:

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