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 and the use of diagnostic tools for insight into the details of the simulation.

Background: please see experiment react_2_a

TAGS : "basic", "under-the-hood", "uniform compartment"¶

In [1]:

LAST_REVISED = "Feb. 16, 2026"
LIFE123_VERSION = "1.0.0rc7"        # Library version this experiment is based on

In [2]:

#import set_path                    # Using MyBinder?  Uncomment this before running the next cell!

In [3]:

#import sys
#sys.path.append("C:/some_path/my_env_or_install")   # CHANGE to the folder containing your venv or libraries installation!
# NOTE: If any of the imports below can't find a module, uncomment the lines above, or try:  import set_path   

import ipynbname

from life123 import check_version, UniformCompartment, PlotlyHelper

In [4]:

check_version(LIFE123_VERSION)

OK

In [5]:

# Initialize the HTML logging (for the graphics)
log_file = ipynbname.name() + ".log.htm"    # Use the notebook base filename for the log file
                                            # IN CASE OF PROBLEMS, set manually to any desired name

In [ ]:

PART 1 - RUN THE SIMULATION¶

Initialize the System¶

In [6]:

# Instantiate the simulator and specify the chemicals
# The diagnostics will be give insight into the inner workings of the simulation
uc = UniformCompartment(names=["A", "B"], preset="mid", enable_diagnostics=True)

# Reaction A <-> B , with 1st-order kinetics in both directions
uc.add_reaction(reactants="A", products="B", kF=3., kR=2.)

uc.describe_reactions()

add_reaction(): detected reaction type `ReactionUnimolecular`
Number of reactions: 1
0: A <-> B  Elementary Unimolecular reaction  (kF = 3 / kR = 2 / delta_G = -1,005.1 / K = 1.5 / Temp = 25 C)
Chemicals involved in the above reactions: ['A', 'B']

In [7]:

# Send a plot of the network of reactions to the HTML log file
uc.plot_reaction_network(log_file=log_file)

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

In [8]:

# Set the initial concentrations of all the chemicals
uc.set_conc({"A": 10., "B": 50.})
uc.describe_state()

SYSTEM STATE at Time t = 0:
2 species:
  Species 0 (A). Conc: 10.0
  Species 1 (B). Conc: 50.0
Chemicals involved in reactions: ['B', 'A']

In [9]:

uc.get_history()

Out[9]:

	SYSTEM TIME	A	B	step	caption
0	0.0	10.0	50.0		Set concentration

In [10]:

# For experiment repeatability, we specified, when instantiating the "UniformCompartment" class, 
# a particular PRESET applicable to the adaptive time steps; 
# let's take a look at the value assigned by that preset
uc.adaptive_steps.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, 'error': 0.25}

In [ ]:

Run the reaction¶

Passing True to variable_steps automatically adjusts up or down the time steps¶

In [11]:

uc.single_compartment_react(initial_step=0.1, target_end_time=1.2,
                            variable_steps=True)

19 total variable step(s) taken in 0.119 sec
Number of step re-do's because of elective soft aborts: 2
Norm usage: {'norm_A': 17, 'norm_B': 15, 'norm_C': 15, 'norm_D': 15}
System Time is now: 1.2268

In [12]:

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

Out[12]:

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

In [13]:

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

0: A <-> B
Current 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[13]:

True

In [ ]:

PART 2 - Visualize the Results¶

Plots of changes of concentration with time¶

In [14]:

uc.plot_history(colors=['darkturquoise', 'green'], 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 [15]:

uc.plot_step_sizes(show_intervals=True)

In [ ]:

PART 2 - Scrutinizing the inner workings of the step-size changes¶

NOTE: this part is NOT meant for typically end users. It's for debugging, and for anyone interested in taking an "under the hood" look

In [16]:

diagnostics = uc.get_diagnostics()      # Available because we enabled the diagnostics

The "Diagnostics" object contains a treasure trove of data and methods yo get insights into the inner workings of the simulation.
Diagnostic data was saved because of the call to enable_diagnostics() prior to running the reaction

In [17]:

type(diagnostics)    # this will show an object

Out[17]:

life123.diagnostics.Diagnostics

In [18]:

# Let's revisit the variables steps taken
diagnostics.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)

In [ ]:

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

In [19]:

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

Reaction:  A <-> B

Out[19]:

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

Note that diagnostic data with the DELTA Concentrations - in the above listing - also records the values that were considered (but not actually used) during ABORTED steps. For example, in steps 0-2, above, the START_TIME remains the same, as the time_step gets progressively reduced until the resulting changes are deemed acceptable.¶

Line 2, above, shows that the concentration of the product [B], which was set by us to an initial value of 50, gets affected by a rate of change of -70, sustained over a delta_time of 0.016 .
The reaction rate is negative because the product is decreasing.
The new value for [B] is:

In [20]:

50. - 70. * 0.016

Out[20]:

48.88

That's indeed the value we saw in the history of the product [B] at time t=0.016 :

In [21]:

uc.get_history(t=0.016)

Out[21]:

	search_value	SYSTEM TIME	A	B	step	caption
1	0.016	0.016	11.12	48.88	1	1st reaction step

In [ ]:

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

history[1:4]

Out[22]:

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

In [23]:

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

Out[23]:

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

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

In [24]:

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

Out[24]:

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

In [25]:

# Computes some measures of how large delta_concentrations is, and propose a course of action
uc.adaptive_steps.adjust_timestep(delta_conc=delta_concentrations, baseline_conc=baseline_conc,
                                  n_chems=2, indexes_of_active_chemicals=uc.get_reactions().indexes_of_active_chemicals())  

Out[25]:

{'action': 'stay',
 'step_factor': 1,
 'norms': {'norm_A': 0.5308620929718018, 'norm_B': 0.09266187},
 'applicable_norms': 'ALL'}

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

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

Out[26]:

0.016000000000000004

In [27]:

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

Out[27]:

0.016000000000000007

In [28]:

next_step / original_step

Out[28]:

1.0000000000000002

In [ ]:

Example 2: very late in the run¶

In [29]:

history[17:20]

Out[29]:

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

In [30]:

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

Out[30]:

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

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

In [31]:

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

Out[31]:

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

In [32]:

# Computes a measure of how large delta_concentrations is, and propose a course of action
uc.adaptive_steps.adjust_timestep(delta_conc=delta_concentrations, baseline_conc=baseline_conc,
                                  n_chems=2, indexes_of_active_chemicals=uc.get_reactions().indexes_of_active_chemicals())  

Out[32]:

{'action': 'low',
 'step_factor': 1.2,
 'norms': {'norm_A': 0.0008898575906641781, 'norm_B': 0.0017613951},
 'applicable_norms': 'ALL'}

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

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

Out[33]:

0.17118912860651525

In [34]:

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

Out[34]:

0.2054269543278182

In [35]:

next_step / original_step

Out[35]:

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

uc.adaptive_steps.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, 'error': 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 instantiated the class UniformCompartment with the option to enable the diagnostics)

In [37]:

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

Out[37]:

	TIME	A	B	caption
0	0.000000	10.000000	50.000000	Set concentration
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 [38]:

diagnostics.get_decisions_data()

Out[38]:

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

In [ ]:

PART 3 - Investigate A_dot, i.e. d[A]/dt¶

In [39]:

rates_df = uc.get_rate_history()
rates_df

Out[39]:

	SYSTEM TIME	rxn0_rate	step
0	0.000000	-70.000000	0
1	0.016000	-64.400000	1
2	0.032000	-59.248000	2
3	0.048000	-54.508160	3
4	0.067200	-49.275377	4
5	0.086400	-44.544940	5
6	0.109440	-39.413363	6
7	0.132480	-34.872944	7
8	0.160128	-30.052108	8
9	0.193306	-25.066824	9
10	0.233119	-20.076882	10
11	0.280894	-15.280942	11
12	0.338225	-10.900592	12
13	0.407022	-7.150948	13
14	0.489579	-4.199162	14
15	0.588647	-2.119154	15
16	0.707528	-0.859515	16
17	0.850186	-0.246433	17
18	1.021375	-0.035500	18

Note that reaction rates are defined for the reaction products; since A is a reactant (in reaction 0, our only reaction), we must flip its sign; since the stoichiometry of A is simply 1, no further adjustment needed.

In [40]:

rates_df['A_dot'] = - rates_df['rxn0_rate']    # Add a column to the Pandas dataframe
rates_df

Out[40]:

	SYSTEM TIME	rxn0_rate	step	A_dot
0	0.000000	-70.000000	0	70.000000
1	0.016000	-64.400000	1	64.400000
2	0.032000	-59.248000	2	59.248000
3	0.048000	-54.508160	3	54.508160
4	0.067200	-49.275377	4	49.275377
5	0.086400	-44.544940	5	44.544940
6	0.109440	-39.413363	6	39.413363
7	0.132480	-34.872944	7	34.872944
8	0.160128	-30.052108	8	30.052108
9	0.193306	-25.066824	9	25.066824
10	0.233119	-20.076882	10	20.076882
11	0.280894	-15.280942	11	15.280942
12	0.338225	-10.900592	12	10.900592
13	0.407022	-7.150948	13	7.150948
14	0.489579	-4.199162	14	4.199162
15	0.588647	-2.119154	15	2.119154
16	0.707528	-0.859515	16	0.859515
17	0.850186	-0.246433	17	0.246433
18	1.021375	-0.035500	18	0.035500

In [41]:

p2 = PlotlyHelper.plot_pandas(df=rates_df, x_var="SYSTEM TIME", fields="A_dot", colors='brown',
                              y_label="dA/dt",
                              title="Rate of change of A with time")
p2

Let's create a combined plot like we had in experiment react_2_a:

In [42]:

p1 = uc.plot_history(chemicals="A", colors='darkturquoise')   # The plot of [A] vs. time that we saw earlier

In [43]:

PlotlyHelper.combine_plots([p1, p2],
                           title="Concentration of A with time, and its rate of change (A_dot)",
                           y_label="[A] (turquoise) /<br> A_dot (brown)",
                           legend_title="Plot",
                           curve_labels=["A", "A_dot"]
                           )

Notice how much smoother the lines are, compared to what we had in experiment `react_2_a` !¶

In [ ]:

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

TAGS : "basic", "under-the-hood", "uniform compartment"¶

PART 1 - RUN THE SIMULATION¶

Initialize the System¶

Run the reaction¶

Passing True to variable_steps automatically adjusts up or down the time steps¶

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 - Visualize the Results¶

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¶

PART 2 - Scrutinizing the inner workings of the step-size changes¶

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

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¶

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 - Investigate A_dot, i.e. d[A]/dt¶

Notice how much smoother the lines are, compared to what we had in 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`¶

Notice how much smoother the lines are, compared to what we had in experiment `react_2_a` !¶