A simple `A <-> B` reaction whose rate constants are to be estimated from a given time evolution of [A] and [B]¶

(values given on a variable-time grid.)¶

Assume the reaction is known to be 1st order (won't verify that.)

In PART 1, a time evolution of [A] and [B] is obtained by simulation

In PART 2, the time evolutions generated in Part 1 are taken as a starting point, to estimate the rate constants of A <-> B

In PART 3, we'll repeat what we did in Part 2, but this time showing the full details of how the answer is arrived at

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

In [1]:

LAST_REVISED = "Sep. 8, 2024"
LIFE123_VERSION = "1.0.0.beta.38"      # Version this experiment is based on

In [2]:

#import set_path            # Using MyBinder?  Uncomment this before running the next cell!
                            # Importing this local file will add the project's home directory to sys.path

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
import numpy as np

from life123 import check_version, UniformCompartment, PlotlyHelper, Numerical

In [4]:

check_version(LIFE123_VERSION)

OK

In [ ]:

PART 1 - We'll generate the time evolution of [A] and [B] by simulating a reaction of KNOWN rate constants...¶

but pretend you don't see this section! (because we later want to estimate those rate constants)¶

In [5]:

# Instantiate the simulator and specify the accuracy preset
dynamics = UniformCompartment(preset="mid")

# Reaction A <-> B (mostly in the forward direction)
dynamics.add_reaction(reactants="A", products="B",
                      forward_rate=12., reverse_rate=2.) 
 
dynamics.describe_reactions()

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

Run the simulation¶

In [6]:

dynamics.set_conc({"A": 40., "B": 10.}, snapshot=True)  # Set the initial concentrations
dynamics.describe_state()

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

In [7]:

dynamics.enable_diagnostics()         # To save diagnostic information for the simulation run, below

dynamics.single_compartment_react(initial_step=0.01, duration=0.5,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  variable_steps=True)

39 total step(s) taken
Number of step re-do's because of elective soft aborts: 2
Norm usage: {'norm_A': 24, 'norm_B': 22, 'norm_C': 22, 'norm_D': 22}

Plots of changes of concentration with time ¶

In [8]:

dynamics.plot_history(colors=['darkturquoise', 'green'], show_intervals=True)

Notice the variable time steps (vertical dashed lines), more frequent when there's more change

In [ ]:

PART 2 - This is the starting point of fitting the data from part 1¶

We're given the data of the above curves - i.e. the system history, and we want to estimate the rate constants (forward and reverse) of the reaction `A <-> B`¶

Let's start by taking stock of the actual data (saved during the simulation of part 1):

In [9]:

df = dynamics.get_history() 
df

Out[9]:

	SYSTEM TIME	A	B	caption
0	0.000000	40.000000	10.000000	Initialized state
1	0.001600	39.264000	10.736000	1st reaction step
2	0.003520	38.400584	11.599416
3	0.005440	37.560376	12.439624
4	0.007744	36.579229	13.420771
5	0.010509	35.439829	14.560171
6	0.013274	34.344532	15.655468
7	0.016038	33.291632	16.708368
8	0.018803	32.279486	17.720514
9	0.021568	31.306517	18.693483
10	0.024886	30.184148	19.815852
11	0.028204	29.113912	20.886088
12	0.031521	28.093386	21.906614
13	0.034839	27.120262	22.879738
14	0.038820	26.006754	23.993246
15	0.042802	24.955312	25.044688
16	0.046783	23.962474	26.037526
17	0.051561	22.837477	27.162523
18	0.056338	21.787726	28.212274
19	0.061116	20.808189	29.191811
20	0.066849	19.711365	30.288635
21	0.072582	18.702575	31.297425
22	0.078315	17.774755	32.225245
23	0.085195	16.750734	33.249266
24	0.092074	15.825343	34.174657
25	0.100330	14.821829	35.178171
26	0.108586	13.934301	36.065699
27	0.118492	12.992362	37.007638
28	0.130381	12.018806	37.981194
29	0.144646	11.044979	38.955021
30	0.158912	10.265644	39.734356
31	0.176031	9.517222	40.482778
32	0.196574	8.834360	41.165640
33	0.221225	8.250593	41.749407
34	0.250806	7.791835	42.208165
35	0.286304	7.469313	42.530687
36	0.328902	7.274627	42.725373
37	0.380018	7.180328	42.819672
38	0.441359	7.148149	42.851851
39	0.514967	7.142696	42.857304	last reaction step

The reaction is mostly forward; the reactant A gets consumed, while the product B gets produced

Let's extract some columns, as Numpy arrays:¶

In [10]:

t_arr = df["SYSTEM TIME"].to_numpy()   # The independent variable : Time

In [11]:

A_conc = df["A"].to_numpy()

In [12]:

B_conc = df["B"].to_numpy()

Here, we take the easy way out, using a specialized Life123 function!¶

(in Part 3, we'll do a step-by-step derivation, to see how it works)

In [13]:

dynamics.estimate_rate_constants_simple(t=t_arr,
                                        A_conc=A_conc, B_conc=B_conc,
                                        reactant_name="A", product_name="B")

Reaction A <-> B
Total REACTANT + PRODUCT has a median of 50, 
    with standard deviation 5.617e-15 (ideally should be zero)
The sum of the time derivatives of the reactant and the product 
    has a median of 0 (ideally should be zero)
Least square fit to B'(t) = kF * A(t) + kR * (- B(t) )

-> ESTIMATED RATE CONSTANTS: kF = 12.19 , kR = 1.925

The least-square fit is good... and the values estimated from the data for kF and kR are in good agreement with the values we used in the simulation to get that data, respectively 12 and 2 (see PART 1, above)¶

Note that our data set is quite skimpy in the number of points:

In [14]:

len(B_conc)

Out[14]:

and that it uses a variable grid, with more points where there's more change, such as in the early times:

In [15]:

t_arr  # Time points in our data set

Out[15]:

array([0.        , 0.0016    , 0.00352   , 0.00544   , 0.007744  ,
       0.0105088 , 0.0132736 , 0.0160384 , 0.0188032 , 0.021568  ,
       0.02488576, 0.02820352, 0.03152128, 0.03483904, 0.03882035,
       0.04280166, 0.04678298, 0.05156055, 0.05633812, 0.0611157 ,
       0.06684879, 0.07258188, 0.07831497, 0.08519467, 0.09207438,
       0.10033003, 0.10858568, 0.11849246, 0.13038059, 0.14464635,
       0.15891211, 0.17603102, 0.19657372, 0.22122495, 0.25080644,
       0.28630421, 0.32890155, 0.38001835, 0.44135851, 0.5149667 ])

In [16]:

np.gradient(t_arr)   # Notice how it gets larger in later times, as bigger steps get taken

Out[16]:

array([0.0016    , 0.00176   , 0.00192   , 0.002112  , 0.0025344 ,
       0.0027648 , 0.0027648 , 0.0027648 , 0.0027648 , 0.00304128,
       0.00331776, 0.00331776, 0.00331776, 0.00364954, 0.00398131,
       0.00398131, 0.00437944, 0.00477757, 0.00477757, 0.00525533,
       0.00573309, 0.00573309, 0.0063064 , 0.00687971, 0.00756768,
       0.00825565, 0.00908121, 0.01089746, 0.01307695, 0.01426576,
       0.01569234, 0.0188308 , 0.02259696, 0.02711636, 0.03253963,
       0.03904756, 0.04685707, 0.05622848, 0.06747418, 0.07360819])

The variable time grid, and the skimpy number of data points, are best seen in the plot that was shown at the end of PART 1¶

In [ ]:

PART 3 - investigate how the `estimate_rate_constants()` function used in part 2 works¶

Again, the starting point are the time evolutions of [A] and [B] , that is the system history that was given to us¶

Let's revisit the Numpy arrays that we had set up at the beginning of Part 2

In [17]:

t_arr    # The independent variable : Time

Out[17]:

array([0.        , 0.0016    , 0.00352   , 0.00544   , 0.007744  ,
       0.0105088 , 0.0132736 , 0.0160384 , 0.0188032 , 0.021568  ,
       0.02488576, 0.02820352, 0.03152128, 0.03483904, 0.03882035,
       0.04280166, 0.04678298, 0.05156055, 0.05633812, 0.0611157 ,
       0.06684879, 0.07258188, 0.07831497, 0.08519467, 0.09207438,
       0.10033003, 0.10858568, 0.11849246, 0.13038059, 0.14464635,
       0.15891211, 0.17603102, 0.19657372, 0.22122495, 0.25080644,
       0.28630421, 0.32890155, 0.38001835, 0.44135851, 0.5149667 ])

In [18]:

A_conc

Out[18]:

array([40.        , 39.264     , 38.40058368, 37.56037599, 36.5792285 ,
       35.43982899, 34.34453244, 33.29163175, 32.27948591, 31.30651739,
       30.18414823, 29.1139116 , 28.093386  , 27.12026243, 26.00675446,
       24.95531161, 23.96247447, 22.83747684, 21.78772586, 20.80818856,
       19.71136465, 18.70257539, 17.77475483, 16.75073404, 15.82534273,
       14.82182904, 13.93430053, 12.992362  , 12.01880624, 11.04497851,
       10.2656443 ,  9.5172222 ,  8.83436019,  8.25059325,  7.7918346 ,
        7.46931299,  7.27462691,  7.18032783,  7.14814942,  7.14269565])

In [19]:

B_conc

Out[19]:

array([10.        , 10.736     , 11.59941632, 12.43962401, 13.4207715 ,
       14.56017101, 15.65546756, 16.70836825, 17.72051409, 18.69348261,
       19.81585177, 20.8860884 , 21.906614  , 22.87973757, 23.99324554,
       25.04468839, 26.03752553, 27.16252316, 28.21227414, 29.19181144,
       30.28863535, 31.29742461, 32.22524517, 33.24926596, 34.17465727,
       35.17817096, 36.06569947, 37.007638  , 37.98119376, 38.95502149,
       39.7343557 , 40.4827778 , 41.16563981, 41.74940675, 42.2081654 ,
       42.53068701, 42.72537309, 42.81967217, 42.85185058, 42.85730435])

Let's verify that the stoichiometry is satified. From the reaction `A <-> B` we can infer that any drop in [A] corresponds to an equal increase in [B]. Their sum will remain constants:¶

In [20]:

A_conc + B_conc

Out[20]:

array([50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
       50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
       50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50., 50.,
       50.])

Just as expected!¶

In [21]:

# Incidentally, there's a function to verify that the stoichiometry 
# of a single reaction holds true across the entire simulation run 
# (overkill in this case!)
dynamics.diagnostics.stoichiometry_checker_entire_run() 

Out[21]:

True

In [ ]:

Now, let's investigate the rates of change of [A] and [B]¶

In [22]:

# The rate of change of [A] with time
Deriv_A = np.gradient(A_conc, t_arr, edge_order=2)

# The rate of change of [B] with time
Deriv_B = np.gradient(B_conc, t_arr, edge_order=2)

In [23]:

# As expected from the stoichiometry, the two derivatives are opposites: 
# when [A] increases by a certain amount, [B] decreases by that same amount
Deriv_A + Deriv_B   # Will be very close to zero throughout

Out[23]:

array([-6.82121026e-12,  1.36424205e-12,  1.36424205e-12, -4.54747351e-13,
        6.82121026e-13,  4.54747351e-13,  0.00000000e+00, -4.54747351e-13,
        4.54747351e-13,  0.00000000e+00,  4.54747351e-13,  4.54747351e-13,
        4.54747351e-13, -4.54747351e-13,  4.54747351e-13,  0.00000000e+00,
        4.54747351e-13,  4.54747351e-13,  0.00000000e+00,  6.82121026e-13,
       -2.27373675e-13,  0.00000000e+00, -4.54747351e-13,  0.00000000e+00,
       -3.41060513e-13, -1.13686838e-13, -2.27373675e-13,  0.00000000e+00,
       -5.68434189e-14, -1.70530257e-13,  3.97903932e-13, -1.13686838e-13,
        2.27373675e-13, -4.26325641e-14, -4.26325641e-14, -9.94759830e-14,
       -5.68434189e-14, -8.52651283e-14,  2.84217094e-14,  1.42108547e-13])

In [24]:

PlotlyHelper.plot_curves(x=t_arr, y=[Deriv_A , Deriv_B], title="d/dt A(t) and d/dt B(t) as a function of time",
                         xlabel="t", ylabel="Time derivatives", curve_labels=["A'(t)", "B'(t)"], 
                         legend_title="Derivative", colors=['aqua', 'greenyellow'])

The rate of changes of both [A] and [B] get smaller as the reaction marches towards equilibrium

In [ ]:

Now, let's determine what kF and kR rate constants for `A <-> B` will yield the above data¶

Assuming that A <-> B is an elementary chemical reaction (i.e. occuring in a single step),
OR THAT IT CAN BE APPROXIMATED AS ONE, then the rate of change of the reaction product concentration B(t) is the difference of the forward rate (producing B) and the reverse rate (consuming it):

B'(t) = kF * A(t) - kR * B(t)

We can re-write it as:
B'(t) = kF * {A(t)} + kR * {- B(t)} (Eqn. 1)

A(t), B(t) are given to us; B'(t) is a gradient we already computed numerically; kF and kR are the rate constants that we are trying to estimate.

If we can do a satisfactory Least Square Fit to express B'(t) as a linear function of {A(t)} and {- B(t)}, as:

B'(t) = a * {A(t)} + b * {- B(t)} , for some constants a and b,

then, comparing with Eqn. 1, it immediately follows that

kF = a
kR = b

Let's carry it out! First, let's verify that B'(t) is indeed a linear function of A(t) and of B(t). We already have, from our data, B'(t) as the Numpy array Deriv_B , and we also have A(t) and B(t) as, respectively, the Numpy arrays A_conc and B_conc

In [25]:

PlotlyHelper.plot_curves(x=A_conc, y=Deriv_B, title="d/dt B(t) as a function of A(t)",
                         xlabel="A(t)", ylabel="B'(t)", colors="green")

As expected, it appears to be a straight line (green), and the rate of change in the product B is higher when the concentration of the reactant A is larger.

In [26]:

PlotlyHelper.plot_curves(x=B_conc, y=Deriv_B, title="d/dt B(t) as a function of B(t)",
                         xlabel="B(t)", ylabel="B'(t)", colors="blue")

Let's do the least-square fit we had set out to do: `B'(t) = a * {A(t)} + b * {- B(t)}` , for some a, b¶

In [27]:

a, b = Numerical.two_vector_least_square(V = A_conc, W = -B_conc, Y = Deriv_B)
a, b

Out[27]:

(12.185860088176147, 1.924952056726416)

Voila', those are, respectively, our estimated kF and kR!¶

We just obtained the same values of the estimated kF and kR as were computed by a call to `estimate_rate_constants()` in Part 2¶

Visually verify the least-square fit¶

In [28]:

# Plot B'(t) vs. its least-square approx
PlotlyHelper.plot_curves(x=Deriv_B , y=[Deriv_B, (a * A_conc - b * B_conc)],
                         title="d/dt B(t) as a function of A(t), alongside its least-square fit",
                         xlabel="A(t)", ylabel="B'(t)", 
                         curve_labels=["B'(t)", "Linear Fit"], legend_title="Curve vs Fit:", colors=['green', 'red'])

Virtually indistinguishable lines! And the same plot we saw in Part 2!

In [ ]:

A simple A <-> B reaction whose rate constants are to be estimated from a given time evolution of [A] and [B]¶

(values given on a variable-time grid.)¶

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

PART 1 - We'll generate the time evolution of [A] and [B] by simulating a reaction of KNOWN rate constants...¶

but pretend you don't see this section! (because we later want to estimate those rate constants)¶

Run the simulation¶

Plots of changes of concentration with time¶

PART 2 - This is the starting point of fitting the data from part 1¶

We're given the data of the above curves - i.e. the system history, and we want to estimate the rate constants (forward and reverse) of the reaction A <-> B¶

Let's extract some columns, as Numpy arrays:¶

Here, we take the easy way out, using a specialized Life123 function!¶

The least-square fit is good... and the values estimated from the data for kF and kR are in good agreement with the values we used in the simulation to get that data, respectively 12 and 2 (see PART 1, above)¶

The variable time grid, and the skimpy number of data points, are best seen in the plot that was shown at the end of PART 1¶

PART 3 - investigate how the estimate_rate_constants() function used in part 2 works¶

Again, the starting point are the time evolutions of [A] and [B] , that is the system history that was given to us¶

Let's verify that the stoichiometry is satified. From the reaction A <-> B we can infer that any drop in [A] corresponds to an equal increase in [B]. Their sum will remain constants:¶

Just as expected!¶

Now, let's investigate the rates of change of [A] and [B]¶

Now, let's determine what kF and kR rate constants for A <-> B will yield the above data¶

Let's do the least-square fit we had set out to do: B'(t) = a * {A(t)} + b * {- B(t)} , for some a, b¶

Voila', those are, respectively, our estimated kF and kR!¶

We just obtained the same values of the estimated kF and kR as were computed by a call to estimate_rate_constants() in Part 2¶

Visually verify the least-square fit¶

A simple `A <-> B` reaction whose rate constants are to be estimated from a given time evolution of [A] and [B]¶

Plots of changes of concentration with time ¶

We're given the data of the above curves - i.e. the system history, and we want to estimate the rate constants (forward and reverse) of the reaction `A <-> B`¶

PART 3 - investigate how the `estimate_rate_constants()` function used in part 2 works¶

Let's verify that the stoichiometry is satified. From the reaction `A <-> B` we can infer that any drop in [A] corresponds to an equal increase in [B]. Their sum will remain constants:¶

Now, let's determine what kF and kR rate constants for `A <-> B` will yield the above data¶

Let's do the least-square fit we had set out to do: `B'(t) = a * {A(t)} + b * {- B(t)}` , for some a, b¶

We just obtained the same values of the estimated kF and kR as were computed by a call to `estimate_rate_constants()` in Part 2¶