A simple `A + B <-> C` 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.)

This is the counterpart of experiment mystery_reaction_1 for a more complex reaction; we'll also proceed faster.

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", "B"], products="C",
                      forward_rate=5., reverse_rate=1.) 
 
dynamics.describe_reactions()

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

Run the simulation¶

In [6]:

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

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

In [7]:

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

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

Some steps were backtracked and re-done, to prevent negative concentrations or excessively large concentration changes
53 total step(s) taken
Number of step re-do's because of elective soft aborts: 1
Norm usage: {'norm_A': 19, 'norm_B': 17, 'norm_C': 17, 'norm_D': 17}

Plots of changes of concentration with time ¶

In [8]:

dynamics.plot_history(colors=['darkturquoise', 'violet', '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 <-> C`¶

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	C	caption
0	0.000000	40.000000	20.000000	5.000000	Initialized state
1	0.000400	38.402000	18.402000	6.598000	1st reaction step
2	0.000600	37.696646	17.696646	7.303354
3	0.000800	37.031002	17.031002	7.968998
4	0.001000	36.401921	16.401921	8.598079
5	0.001240	35.687511	15.687511	9.312489
6	0.001480	35.017928	15.017928	9.982072
7	0.001768	34.263512	14.263512	10.736488
8	0.002114	33.422717	13.422717	11.577283
9	0.002528	32.497253	12.497253	12.502747
10	0.003026	31.492903	11.492903	13.507097
11	0.003524	30.598990	10.598990	14.401010
12	0.004121	29.639181	9.639181	15.360819
13	0.004718	28.795266	8.795266	16.204734
14	0.005315	28.048707	8.048707	16.951293
15	0.005912	27.384727	7.384727	17.615273
16	0.006510	26.791395	6.791395	18.208605
17	0.007107	26.258967	6.258967	18.741033
18	0.007823	25.683487	5.683487	19.316513
19	0.008540	25.174287	5.174287	19.825713
20	0.009257	24.721753	4.721753	20.278247
21	0.009973	24.318020	4.318020	20.681980
22	0.010690	23.956587	3.956587	21.043413
23	0.011407	23.632031	3.632031	21.367969
24	0.012123	23.339792	3.339792	21.660208
25	0.012840	23.076005	3.076005	21.923995
26	0.013700	22.789650	2.789650	22.210350
27	0.014560	22.535388	2.535388	22.464612
28	0.015420	22.309033	2.309033	22.690967
29	0.016280	22.107053	2.107053	22.892947
30	0.017140	21.926451	1.926451	23.073549
31	0.018000	21.764669	1.764669	23.235331
32	0.018860	21.619505	1.619505	23.380495
33	0.019720	21.489062	1.489062	23.510938
34	0.020580	21.371693	1.371693	23.628307
35	0.021612	21.244815	1.244815	23.755185
36	0.022644	21.132875	1.132875	23.867125
37	0.023675	21.033975	1.033975	23.966025
38	0.024707	20.946489	0.946489	24.053511
39	0.025739	20.869015	0.869015	24.130985
40	0.026771	20.800342	0.800342	24.199658
41	0.028010	20.727233	0.727233	24.272767
42	0.029248	20.663960	0.663960	24.336040
43	0.030486	20.609146	0.609146	24.390854
44	0.031725	20.561619	0.561619	24.438381
45	0.033211	20.512134	0.512134	24.487866
46	0.034697	20.470471	0.470471	24.529529
47	0.036183	20.435365	0.435365	24.564635
48	0.037966	20.399844	0.399844	24.600156
49	0.039749	20.370985	0.370985	24.629015
50	0.041889	20.342829	0.342829	24.657171
51	0.044457	20.316603	0.316603	24.683397
52	0.047538	20.293560	0.293560	24.706440
53	0.051236	20.274774	0.274774	24.725226	last reaction step

The reaction is mostly forward; the reactants A and B gets consumed, while the product C 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()

In [13]:

C_conc = df["C"].to_numpy()

In [ ]:

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

dynamics.estimate_rate_constants_association(t=t_arr, 
                                             A_conc=A_conc, B_conc=B_conc, C_conc=C_conc, 
                                             reactants=["A", "B"], product="C")

Least square fit to C'(t) = kF * A(t) * B(t) + kR * (- C(t) )

-> ESTIMATED RATE CONSTANTS: kF = 5.233 , kR = 0.8841

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

len(B_conc)

Out[15]:

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

In [16]:

t_arr  # Time points in our data set

Out[16]:

array([0.        , 0.0004    , 0.0006    , 0.0008    , 0.001     ,
       0.00124   , 0.00148   , 0.001768  , 0.0021136 , 0.00252832,
       0.00302598, 0.00352365, 0.00412084, 0.00471804, 0.00531524,
       0.00591244, 0.00650963, 0.00710683, 0.00782346, 0.0085401 ,
       0.00925674, 0.00997337, 0.01069001, 0.01140665, 0.01212328,
       0.01283992, 0.01369988, 0.01455984, 0.01541981, 0.01627977,
       0.01713974, 0.0179997 , 0.01885966, 0.01971963, 0.02057959,
       0.02161154, 0.0226435 , 0.02367546, 0.02470741, 0.02573937,
       0.02677133, 0.02800967, 0.02924802, 0.03048637, 0.03172471,
       0.03321073, 0.03469675, 0.03618276, 0.03796598, 0.0397492 ,
       0.04188907, 0.04445691, 0.04753831, 0.051236  ])

In [17]:

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

Out[17]:

array([0.0004    , 0.0003    , 0.0002    , 0.0002    , 0.00022   ,
       0.00024   , 0.000264  , 0.0003168 , 0.00038016, 0.00045619,
       0.00049766, 0.00054743, 0.0005972 , 0.0005972 , 0.0005972 ,
       0.0005972 , 0.0005972 , 0.00065692, 0.00071664, 0.00071664,
       0.00071664, 0.00071664, 0.00071664, 0.00071664, 0.00071664,
       0.0007883 , 0.00085996, 0.00085996, 0.00085996, 0.00085996,
       0.00085996, 0.00085996, 0.00085996, 0.00085996, 0.00094596,
       0.00103196, 0.00103196, 0.00103196, 0.00103196, 0.00103196,
       0.00113515, 0.00123835, 0.00123835, 0.00123835, 0.00136218,
       0.00148602, 0.00148602, 0.00163462, 0.00178322, 0.00196154,
       0.00235385, 0.00282462, 0.00338954, 0.00369769])

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], [B] and [C] , 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 [18]:

t_arr    # The independent variable : Time

Out[18]:

array([0.        , 0.0004    , 0.0006    , 0.0008    , 0.001     ,
       0.00124   , 0.00148   , 0.001768  , 0.0021136 , 0.00252832,
       0.00302598, 0.00352365, 0.00412084, 0.00471804, 0.00531524,
       0.00591244, 0.00650963, 0.00710683, 0.00782346, 0.0085401 ,
       0.00925674, 0.00997337, 0.01069001, 0.01140665, 0.01212328,
       0.01283992, 0.01369988, 0.01455984, 0.01541981, 0.01627977,
       0.01713974, 0.0179997 , 0.01885966, 0.01971963, 0.02057959,
       0.02161154, 0.0226435 , 0.02367546, 0.02470741, 0.02573937,
       0.02677133, 0.02800967, 0.02924802, 0.03048637, 0.03172471,
       0.03321073, 0.03469675, 0.03618276, 0.03796598, 0.0397492 ,
       0.04188907, 0.04445691, 0.04753831, 0.051236  ])

In [19]:

A_conc

Out[19]:

array([40.        , 38.402     , 37.696646  , 37.03100247, 36.40192117,
       35.68751098, 35.01792811, 34.26351166, 33.42271749, 32.49725273,
       31.4929025 , 30.59898987, 29.6391806 , 28.79526613, 28.04870718,
       27.38472713, 26.79139514, 26.25896664, 25.68348706, 25.17428674,
       24.72175309, 24.31802047, 23.95658748, 23.63203138, 23.33979186,
       23.07600524, 22.78964984, 22.53538845, 22.30903301, 22.10705297,
       21.92645145, 21.76466854, 21.61950517, 21.48906248, 21.37169309,
       21.24481542, 21.13287483, 21.03397486, 20.94648874, 20.86901508,
       20.80034206, 20.7272334 , 20.66396015, 20.60914571, 20.56161917,
       20.51213389, 20.47047055, 20.43536453, 20.39984363, 20.37098474,
       20.34282925, 20.31660287, 20.29355989, 20.27477404])

In [20]:

B_conc

Out[20]:

array([20.        , 18.402     , 17.696646  , 17.03100247, 16.40192117,
       15.68751098, 15.01792811, 14.26351166, 13.42271749, 12.49725273,
       11.4929025 , 10.59898987,  9.6391806 ,  8.79526613,  8.04870718,
        7.38472713,  6.79139514,  6.25896664,  5.68348706,  5.17428674,
        4.72175309,  4.31802047,  3.95658748,  3.63203138,  3.33979186,
        3.07600524,  2.78964984,  2.53538845,  2.30903301,  2.10705297,
        1.92645145,  1.76466854,  1.61950517,  1.48906248,  1.37169309,
        1.24481542,  1.13287483,  1.03397486,  0.94648874,  0.86901508,
        0.80034206,  0.7272334 ,  0.66396015,  0.60914571,  0.56161917,
        0.51213389,  0.47047055,  0.43536453,  0.39984363,  0.37098474,
        0.34282925,  0.31660287,  0.29355989,  0.27477404])

In [21]:

C_conc

Out[21]:

array([ 5.        ,  6.598     ,  7.303354  ,  7.96899753,  8.59807883,
        9.31248902,  9.98207189, 10.73648834, 11.57728251, 12.50274727,
       13.5070975 , 14.40101013, 15.3608194 , 16.20473387, 16.95129282,
       17.61527287, 18.20860486, 18.74103336, 19.31651294, 19.82571326,
       20.27824691, 20.68197953, 21.04341252, 21.36796862, 21.66020814,
       21.92399476, 22.21035016, 22.46461155, 22.69096699, 22.89294703,
       23.07354855, 23.23533146, 23.38049483, 23.51093752, 23.62830691,
       23.75518458, 23.86712517, 23.96602514, 24.05351126, 24.13098492,
       24.19965794, 24.2727666 , 24.33603985, 24.39085429, 24.43838083,
       24.48786611, 24.52952945, 24.56463547, 24.60015637, 24.62901526,
       24.65717075, 24.68339713, 24.70644011, 24.72522596])

Let's verify that the stoichiometry is satified. From the reaction `A + B <-> C` we can infer that and any drop in [A] corresponds to an equal drop in [B] , and an equal increase in [C]¶

In [22]:

A_conc - B_conc    # Constant differences will show that they change in lockstep

Out[22]:

array([20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.,
       20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.,
       20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.,
       20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20., 20.,
       20., 20.])

In [23]:

A_conc + C_conc    # Constant sums will show that any drop in [A] corresponds to an equal increase in [C]

Out[23]:

array([45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45.,
       45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45.,
       45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45.,
       45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45., 45.,
       45., 45.])

In [ ]:

In [24]:

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

True

In [ ]:

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

In [25]:

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

# The rate of change of [C] with time
Deriv_C = np.gradient(C_conc, t_arr, edge_order=2)

In [26]:

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

Out[26]:

array([ 8.00355338e-11, -1.81898940e-11, -7.27595761e-12,  7.27595761e-12,
        2.91038305e-11, -7.27595761e-12, -2.54658516e-11,  9.09494702e-12,
        3.63797881e-12,  5.45696821e-12,  0.00000000e+00, -5.45696821e-12,
        0.00000000e+00,  0.00000000e+00, -1.81898940e-12, -1.81898940e-12,
        0.00000000e+00, -1.81898940e-12,  0.00000000e+00,  1.81898940e-12,
        0.00000000e+00, -1.81898940e-12, -1.81898940e-12,  1.81898940e-12,
        0.00000000e+00, -3.63797881e-12,  1.81898940e-12,  0.00000000e+00,
       -1.81898940e-12,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        1.81898940e-12,  0.00000000e+00,  0.00000000e+00,  1.81898940e-12,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        1.81898940e-12, -9.09494702e-13,  0.00000000e+00, -2.72848411e-12,
        9.09494702e-13,  9.09494702e-13,  0.00000000e+00,  1.36424205e-12,
        9.09494702e-13, -5.45696821e-12])

In [27]:

# Similarly:
Deriv_B + Deriv_C   # Will be very close to zero throughout

Out[27]:

array([ 4.36557457e-11, -1.09139364e-11,  0.00000000e+00,  0.00000000e+00,
        3.63797881e-12,  0.00000000e+00, -1.45519152e-11,  5.45696821e-12,
        3.63797881e-12, -1.81898940e-12,  0.00000000e+00, -3.63797881e-12,
        2.72848411e-12,  2.72848411e-12,  9.09494702e-13,  9.09494702e-13,
        0.00000000e+00, -4.54747351e-13, -4.54747351e-13,  9.09494702e-13,
        2.27373675e-12,  9.09494702e-13, -9.09494702e-13, -1.81898940e-12,
       -2.27373675e-12, -1.13686838e-12,  6.82121026e-13,  6.82121026e-13,
        4.54747351e-13,  0.00000000e+00, -4.54747351e-13,  1.25055521e-12,
        1.36424205e-12, -1.25055521e-12, -1.93267624e-12,  1.02318154e-12,
       -5.68434189e-13, -5.68434189e-13, -4.54747351e-13, -6.82121026e-13,
        2.84217094e-13, -6.82121026e-13, -8.52651283e-14, -4.54747351e-13,
        2.27373675e-13, -8.52651283e-14,  0.00000000e+00, -1.81898940e-12,
        1.25055521e-12,  1.35003120e-12,  1.35003120e-13,  5.40012479e-13,
        4.05009359e-13, -3.29691829e-12])

In [28]:

# As expected from the stoichiometry, [A] and [B] vary in unison in time
Deriv_A - Deriv_B   # Will be very close to zero throughout

Out[28]:

array([ 3.63797881e-11, -7.27595761e-12, -7.27595761e-12,  7.27595761e-12,
        2.54658516e-11, -7.27595761e-12, -1.09139364e-11,  3.63797881e-12,
        0.00000000e+00,  7.27595761e-12,  0.00000000e+00, -1.81898940e-12,
       -2.72848411e-12, -2.72848411e-12, -2.72848411e-12, -2.72848411e-12,
        0.00000000e+00, -1.36424205e-12,  4.54747351e-13,  9.09494702e-13,
       -2.27373675e-12, -2.72848411e-12, -9.09494702e-13,  3.63797881e-12,
        2.27373675e-12, -2.50111043e-12,  1.13686838e-12, -6.82121026e-13,
       -2.27373675e-12,  0.00000000e+00,  4.54747351e-13, -1.25055521e-12,
        4.54747351e-13,  1.25055521e-12,  1.93267624e-12,  7.95807864e-13,
        5.68434189e-13,  5.68434189e-13,  4.54747351e-13,  6.82121026e-13,
       -2.84217094e-13,  6.82121026e-13,  8.52651283e-14,  4.54747351e-13,
        1.59161573e-12, -8.24229573e-13,  0.00000000e+00, -9.09494702e-13,
       -3.41060513e-13, -4.40536496e-13, -1.35003120e-13,  8.24229573e-13,
        5.04485342e-13, -2.16004992e-12])

In [29]:

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

The rate of changes of both [A] and [C] get smaller as the reaction marches towards equilibrium.
B'(t) not shown, because essentially identical to A'(t)

In [ ]:

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

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

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

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

A(t) and 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 C'(t) as a linear function of {A(t) * B(t)} and {- C(t)}, as:

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

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

kF = a
kR = b

Let's carry out the least-square fit:¶

In [30]:

a, b = Numerical.two_vector_least_square(V = A_conc * B_conc, W = - C_conc, Y = Deriv_C)
a, b

Out[30]:

(5.2325434338059145, 0.8841249639223787)

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¶

In [ ]:

A simple A + B <-> C 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 <-> C¶

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], [B] and [C] , that is the system history that was given to us¶

Let's verify that the stoichiometry is satified. From the reaction A + B <-> C we can infer that and any drop in [A] corresponds to an equal drop in [B] , and an equal increase in [C]¶

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

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

Let's carry out the least-square fit:¶

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¶

A simple `A + B <-> C` 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 <-> C`¶

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 <-> C` we can infer that and any drop in [A] corresponds to an equal drop in [B] , and an equal increase in [C]¶

Now, let's determine what kF and kR rate constants for `A + B <-> C` 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¶