A complex (multistep) reaction `A <-> C` derived from 2 coupled elementary reactions:¶

`A <-> B` (fast) and `B <-> C` (slow)¶

A repeat of experiment cascade_2_b, but with the fast/slow elementary reactions in reverse order.

In PART 1, a time evolution of [A], [B] and [C] is obtained by simulation
In PART 2, the time functions generated in Part 1 are taken as a starting point, to explore how to model the composite reaction A <-> C

Background: please see experiment cascade_2_b

LAST REVISED: Dec. 6, 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_dynamics import ReactionDynamics
from src.modules.visualization.plotly_helper import PlotlyHelper

In [ ]:

PART 1 - We'll generate the time evolution of [A] and [C] by simulating coupled elementary reactions of KNOWN rate constants...¶

but pretend you don't see this section! (because we later assume that those time evolutions are just GIVEN to us)¶

In [3]:

# Instantiate the simulator and specify the chemicals
dynamics = ReactionDynamics(names=["A", "B", "C"])

# Reaction A <-> B (much faster, and with a much larger K)
dynamics.add_reaction(reactants="A", products="B",
                       forward_rate=80., reverse_rate=0.1)   # <===== SPEEDS of 2 reactions ARE REVERSED, relative to experiment `cascade_2_b`

# Reaction B <-> C (much slower, and with a much smaller K)
dynamics.add_reaction(reactants="B", products="C",                     
                      forward_rate=8., reverse_rate=2.)      # <===== SPEEDS of 2 reactions ARE REVERSED, relative to experiment `cascade_2_b`
                                   
dynamics.describe_reactions()

Number of reactions: 2 (at temp. 25 C)
0: A <-> B  (kF = 80 / kR = 0.1 / delta_G = -16,571 / K = 800) | 1st order in all reactants & products
1: B <-> C  (kF = 8 / kR = 2 / delta_G = -3,436.6 / K = 4) | 1st order in all reactants & products
Set of chemicals involved in the above reactions: {'B', 'A', 'C'}

Run the simulation¶

In [4]:

dynamics.set_conc([50., 0., 0.], snapshot=True)  # Set the initial concentrations of all the chemicals, in their index order
dynamics.describe_state()

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

In [5]:

# These settings can be tweaked to make the time resolution finer or coarser.  
# Here we use a "mid" heuristic: neither too fast nor too prudent
dynamics.use_adaptive_preset(preset="mid")

dynamics.single_compartment_react(initial_step=0.01, reaction_duration=0.8,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  variable_steps=True, explain_variable_steps=False)

* INFO: the tentative time step (0.01) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.4 (set to 0.004) [Step started at t=0, and will rewind there]
* INFO: the tentative time step (0.004) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.4 (set to 0.0016) [Step started at t=0, and will rewind there]
* INFO: the tentative time step (0.0016) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.4 (set to 0.00064) [Step started at t=0, and will rewind there]
* INFO: the tentative time step (0.00064) leads to a least one norm value > its ABORT threshold:
      -> will backtrack, and re-do step with a SMALLER delta time, multiplied by 0.4 (set to 0.000256) [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
152 total step(s) taken

In [6]:

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

plot_curves() WARNING: Excessive number of vertical lines (153) - only showing 1 every 2 lines

In [7]:

dynamics.is_in_equilibrium(tolerance=3)

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

1: B <-> C
Final concentrations: [B] = 9.998 ; [C] = 39.99
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 3.99992
    Formula used:  [C] / [B]
2. Ratio of forward/reverse reaction rates: 4.0
Discrepancy between the two values: 0.001889 %
Reaction IS in equilibrium (within 3% tolerance)

Out[7]:

True

A profoundly different plot than we got in experiment `cascade_2_b`¶

In [ ]:

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

We're given the data of the time evolution of `A` and `C`, and we want to try to model the complex reaction `A <-> C`¶

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

In [8]:

df = dynamics.get_history(columns=["SYSTEM TIME", "A", "C", "caption"])    # We're NOT given the intermediary B
df

Out[8]:

	SYSTEM TIME	A	C	caption
0	0.000000	50.000000	0.000000	Initial state
1	0.000256	48.976000	0.000000	1st reaction step
2	0.000563	47.772397	0.002517
3	0.000717	47.185404	0.005250
4	0.000794	46.895519	0.006975
...	...	...	...	...
148	0.450771	0.012844	39.746776
149	0.516094	0.012619	39.905478
150	0.594482	0.012517	39.971659
151	0.688547	0.012538	39.988899
152	0.801426	0.012177	39.990107	last reaction step

153 rows × 4 columns

Column B is NOT given to us. For example, `B` might be an intermediary we can't measure.¶

Only [A] and [C] are given to us, on some variably-spaced time grid¶

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

In [9]:

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

In [10]:

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

In [11]:

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

If the composite reaction `A <-> C` could be modeled as an elementary reaction, we'd expect the rate of change of [C] to be proportional to [A]¶

Let's see what happens if we try to do such a linear fit!

In [12]:

dynamics.estimate_rate_constants(t=t_arr, reactant_conc=A_conc, product_conc=C_conc, reactant_name="A", product_name="C")

Total REACTANT + PRODUCT has a median of 24.85, 
    with standard deviation 12.08 (ideally should be zero)
The sum of the time derivatives of reactant and product 
    has a median of -234.3 (ideally should be zero)
Least square fit: Y = 214.9 + -3.109 X
    where X is the array [A] and Y is the time gradient of C

-> ESTIMATED RATE CONSTANTS: kF = 5.54 , kR = -8.649

A terrible fit! But it looks like we'll do much better if split it into 3 portions, one where A(t) ranges from 0 to about 0.04, one from about 0.04 to 5, and one from about 5 to 50¶

This is a larger numer of splits than we did in experiments cascade_2_a and cascade_2_b

Let's visually locate at what times those [A] values occur:

In [23]:

dynamics.plot_history(chemicals='A', colors='blue', xrange=[0, 0.15], 
                      vertical_lines=[0.028, 0.1], title="[A] as a function of time")

In [14]:

dynamics.get_history(columns=["SYSTEM TIME", "A"], t_start=0.02, t_end=0.03) 

Out[14]:

	SYSTEM TIME	A
72	0.020462	9.428805
73	0.021645	8.540614
74	0.022828	7.736565
75	0.024012	7.008682
76	0.025195	6.349747
77	0.026378	5.753223
78	0.027561	5.213195
79	0.028745	4.724310
80	0.029928	4.281718

[A] assumes the value 5 around t=0.028

In [15]:

dynamics.get_history(columns=["SYSTEM TIME", "A"], t_start=0.09, t_end=0.11) 

Out[15]:

	SYSTEM TIME	A
126	0.090918	0.062461
127	0.093371	0.057137
128	0.095825	0.052750
129	0.098769	0.048391
130	0.101714	0.044909
131	0.105247	0.041540
132	0.109487	0.038396

[A] assumes the value 0.05 around t=0.1

Let's split the `A_conc` and `C_conc` arrays we extracted earlier (with the entire time evolution of, respectively, [A] and [C]) into 3 parts:¶

points numbered 0-78
points 79-128
points 128-end

In [16]:

A_conc_early = A_conc[:79]
A_conc_mid = A_conc[79:129]
A_conc_late = A_conc[129:]

C_conc_early = C_conc[:79]
C_conc_mid = C_conc[79:129]
C_conc_late = C_conc[129:]

t_arr_early = t_arr[:79]
t_arr_mid = t_arr[79:129]
t_arr_late = t_arr[129:]

I. Let's start with the EARLY region, when t < 0.028¶

In [17]:

dynamics.estimate_rate_constants(t=t_arr_early, reactant_conc=A_conc_early, product_conc=C_conc_early, 
                                 reactant_name="A", product_name="C")

Total REACTANT + PRODUCT has a median of 35.2, 
    with standard deviation 11.66 (ideally should be zero)
The sum of the time derivatives of reactant and product 
    has a median of -2,699 (ideally should be zero)
Least square fit: Y = 354.9 + -6.951 X
    where X is the array [A] and Y is the time gradient of C

-> ESTIMATED RATE CONSTANTS: kF = 3.132 , kR = -10.08

Just as we saw in experiment cascade_2_b, trying to fit an elementary reaction to that region leads to a negative reverse rate constant!
We won't discuss this part any further.

II. Let's now consider the MID region (not seen in earlier experiments in this series), when 0.028 < t < 0.1¶

In [18]:

dynamics.estimate_rate_constants(t=t_arr_mid, reactant_conc=A_conc_mid, product_conc=C_conc_mid, 
                                 reactant_name="A", product_name="C")

Total REACTANT + PRODUCT has a median of 14.9, 
    with standard deviation 3.612 (ideally should be zero)
The sum of the time derivatives of reactant and product 
    has a median of 199 (ideally should be zero)
Least square fit: Y = 224.4 + 24.68 X
    where X is the array [A] and Y is the time gradient of C

-> ESTIMATED RATE CONSTANTS: kF = 39.74 , kR = -15.06

For this region, too, trying to fit an elementary reaction to that region leads to a negative reverse rate!
We won't discuss this part any further.

III. And now let's consider the LATE region, when t > 0.1¶

In [19]:

dynamics.estimate_rate_constants(t=t_arr_late, reactant_conc=A_conc_late, product_conc=C_conc_late, 
                                 reactant_name="A", product_name="C")

Total REACTANT + PRODUCT has a median of 34.37, 
    with standard deviation 6.092 (ideally should be zero)
The sum of the time derivatives of reactant and product 
    has a median of 62.73 (ideally should be zero)
Least square fit: Y = -55.81 + 5,404 X
    where X is the array [A] and Y is the time gradient of C

-> ESTIMATED RATE CONSTANTS: kF = 5,402 , kR = 1.624

This time we have an adequate linear fit AND positive rate constants, but a huge value of kF relative to that of any of the elementary reactions!

Let's see the time evolution again, but just for A and C:

In [22]:

dynamics.plot_history(chemicals=['A', 'C'], colors=['blue',  'green'], xrange=[0, 0.25], 
                      vertical_lines=[0.028, 0.1], title='A <-> C compound reaction')

For t > 0.1, the product [C] appears to have a lot of response to nearly non-existing changes in [A] !

That doesn't seem like a good fit to an elementary reaction...

This time, modeling the compound reaction `A <-> C` as an elementary reaction is not a good fit at any time¶

EXPANDED conclusion :¶

While it's a well-known Chemistry notion that the slower reaction is the rate-determining step in a chain, we saw in this experiment that the complex reaction CANNOT always be modeled, not even roughly, with the rate constants of the slower reaction.¶

We saw this in a scenario where the slower reaction is the later one.

In [ ]:

A complex (multistep) reaction A <-> C derived from 2 coupled elementary reactions:¶

A <-> B (fast) and B <-> C (slow)¶

PART 1 - We'll generate the time evolution of [A] and [C] by simulating coupled elementary reactions of KNOWN rate constants...¶

but pretend you don't see this section! (because we later assume that those time evolutions are just GIVEN to us)¶

Run the simulation¶

A profoundly different plot than we got in experiment cascade_2_b¶

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

We're given the data of the time evolution of A and C, and we want to try to model the complex reaction A <-> C¶

Column B is NOT given to us. For example, B might be an intermediary we can't measure.¶

Only [A] and [C] are given to us, on some variably-spaced time grid¶

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

If the composite reaction A <-> C could be modeled as an elementary reaction, we'd expect the rate of change of [C] to be proportional to [A]¶

A terrible fit! But it looks like we'll do much better if split it into 3 portions, one where A(t) ranges from 0 to about 0.04, one from about 0.04 to 5, and one from about 5 to 50¶

Let's split the A_conc and C_conc arrays we extracted earlier (with the entire time evolution of, respectively, [A] and [C]) into 3 parts:¶

I. Let's start with the EARLY region, when t < 0.028¶

II. Let's now consider the MID region (not seen in earlier experiments in this series), when 0.028 < t < 0.1¶

III. And now let's consider the LATE region, when t > 0.1¶

This time, modeling the compound reaction A <-> C as an elementary reaction is not a good fit at any time¶

EXPANDED conclusion :¶

While it's a well-known Chemistry notion that the slower reaction is the rate-determining step in a chain, we saw in this experiment that the complex reaction CANNOT always be modeled, not even roughly, with the rate constants of the slower reaction.¶

A complex (multistep) reaction `A <-> C` derived from 2 coupled elementary reactions:¶

`A <-> B` (fast) and `B <-> C` (slow)¶

A profoundly different plot than we got in experiment `cascade_2_b`¶

We're given the data of the time evolution of `A` and `C`, and we want to try to model the complex reaction `A <-> C`¶

Column B is NOT given to us. For example, `B` might be an intermediary we can't measure.¶

If the composite reaction `A <-> C` could be modeled as an elementary reaction, we'd expect the rate of change of [C] to be proportional to [A]¶

Let's split the `A_conc` and `C_conc` arrays we extracted earlier (with the entire time evolution of, respectively, [A] and [C]) into 3 parts:¶

This time, modeling the compound reaction `A <-> C` as an elementary reaction is not a good fit at any time¶