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 functions 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

LAST REVISED: June 14, 2024 (using v. 1.0 beta33)

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.uniform_compartment import UniformCompartment
from src.modules.visualization.plotly_helper import PlotlyHelper

import numpy as np

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

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

# Reaction A <-> B
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: {'B', 'A'}

Run the simulation¶

In [4]:

dynamics.set_conc([40., 10.], snapshot=True)  # Set the initial concentrations of all the chemicals, in their index order
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: {'B', 'A'}

In [5]:

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

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 negative concentrations: 0
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 ¶

Notice the variable time steps (vertical dashed lines)

In [6]:

dynamics.plot_history(title="Reaction A <-> B",
                      colors=['darkturquoise', 'green'], show_intervals=True)

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

df = dynamics.get_history() 
df

Out[7]:

	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

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

In [8]:

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

In [9]:

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

In [10]:

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

dynamics.estimate_rate_constants(t=t_arr, reactant_conc=A_conc, product_conc=B_conc, reactant_name="A", product_name="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 reactant and product 
    has a median of 0 (ideally should be zero)
Least square fit: Y = -96.25 + 14.11 X
    where X is the array [A] and Y is the time gradient of B

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

len(B_conc)

Out[12]:

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

In [13]:

t_arr  # Time points in our data set

Out[13]:

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

np.gradient(t_arr)

Out[14]:

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

t_arr    # The independent variable : Time

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

A_conc

Out[16]:

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

B_conc

Out[17]:

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

A_conc + B_conc

Out[18]:

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! We'll call that constant value, "TOT_conc"¶

In [19]:

TOT_conc = 50.

In [20]:

# 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.stoichiometry_checker_entire_run() 

Out[20]:

True

In [ ]:

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

In [21]:

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

# 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

Out[22]:

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

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

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 [B] is the difference of the forward rate (producing B) and the reverse rate (consuming it):

B'(t) = kF * A(t) - kR * B(t) (Eqn. 1)

We also know that A(t) + B(t) = TOT_conc (a CONSTANT), i.e.
B(t) = TOT_conc - A(t) (Eqn. 2)

Replacing B(t) from Eqn. 2 into Eqn. 1:

B'(t) = kF * A(t) - kR * [TOT_conc - A(t)]

Simplifying and rearranging:

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

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

B'(t) = [- kR * TOT_conc] + [kF + kR] * A(t) (Eqn. 3)

TOT_conc is a known constant; 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), as:

B'(t) = a + b * A(t) , for some constants a, b

then, comparing with Eqn. 3, we get the following system of equations:

- kR * TOT_conc = a
kF + kR = b

which can be immediately solved as:

kR = - a / TOT_conc (Eqn. 4)
kF = b - kR

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

In [24]:

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, and the rate of change in the product B is higher when the concentration of the reactant A is larger.

If we fit a linear model (least-square fit straight line), we can estimate B'(t) = a + b * A(t) , for some numbers a and b.
I.e. we want to fit: Y = a + b * X , for some numbers a and b
where Y is Deriv_B and X is A_conc, the Numpy arrays we computed earlier:

In [25]:

Y = Deriv_B    # The dependent variable (the time-gradient of the PRODUCT concentrations)
Y

Out[25]:

array([ 4.64683636e+02,  4.55316364e+02,  4.43652086e+02,  4.32261395e+02,
        4.19601598e+02,  4.04133402e+02,  3.88490530e+02,  3.73453149e+02,
        3.58997824e+02,  3.45721185e+02,  3.30434659e+02,  3.15086419e+02,
        3.00451083e+02,  2.87114799e+02,  2.71889118e+02,  2.56734462e+02,
        2.43056308e+02,  2.27599659e+02,  2.12376419e+02,  1.98794743e+02,
        1.83636872e+02,  1.68897580e+02,  1.55931743e+02,  1.41678422e+02,
        1.28621423e+02,  1.14530202e+02,  1.01857698e+02,  8.90860531e+01,
        7.56977131e+01,  6.14464933e+01,  4.96702996e+01,  3.89563341e+01,
        2.88956253e+01,  1.99661621e+01,  1.25889323e+01,  7.03327418e+00,
        3.33147049e+00,  1.24469169e+00,  3.19817889e-01, -1.71634173e-01])

In [26]:

X = A_conc    # The independent variable (the REACTANT concentrations)
X

Out[26]:

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

Let's do the least-square fit:¶

In [27]:

M = np.vstack([np.ones(len(Y)), X]).T
# M is an nx2 matrix , where n is the number of data points.  
# The 2nd column contains the values of X

M[:10, :]   # Show the first 10 rows

Out[27]:

array([[ 1.        , 40.        ],
       [ 1.        , 39.264     ],
       [ 1.        , 38.40058368],
       [ 1.        , 37.56037599],
       [ 1.        , 36.5792285 ],
       [ 1.        , 35.43982899],
       [ 1.        , 34.34453244],
       [ 1.        , 33.29163175],
       [ 1.        , 32.27948591],
       [ 1.        , 31.30651739]])

In [28]:

a, b = np.linalg.lstsq(M, Y, rcond=None)[0]  # Carry out the least-square fit  as: Y = a + b X
a, b

Out[28]:

(-96.24760283632102, 14.110812144902573)

Visually verify the least-square fit¶

In [29]:

# Plot both Y and its least-square fit, as functions of X

PlotlyHelper.plot_curves(x=A_conc, y=[Deriv_B , a + b*A_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!

Finally, from equations 4, repeated here:

kR = - a / TOT_conc
kF = b - kR

In [30]:

kR = - a / TOT_conc
kR

Out[30]:

1.9249520567264204

In [31]:

kF = b - kR
kF

Out[31]:

12.185860088176153

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¶

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

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! We'll call that constant value, "TOT_conc"¶

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:¶

Visually verify the least-square fit¶

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¶

A simple `A <-> B` reaction whose rate constants are to be estimated¶

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¶

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¶