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

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

A repeat of experiment cascade_2_a, with more DISSIMILAR elementary reactions.

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_a

In [1]:

LAST_REVISED = "Nov. 12, 2024"
LIFE123_VERSION = "1.0.0.rc.0"      # 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 

from life123 import UniformCompartment, ReactionKinetics, 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 [4]:

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

# Reaction A <-> B (much slower, and with a much smaller K)
dynamics.add_reaction(reactants="A", products="B",
                      forward_rate=8., reverse_rate=2.) 

# Reaction B <-> C (much faster, and with a much larger K)
dynamics.add_reaction(reactants="B", products="C",
                      forward_rate=80., reverse_rate=0.1)   # <===== THIS IS THE KEY DIFFERENCE FROM THE EARLIER EXPERIMENT `cascade_2_a`
                                   
dynamics.describe_reactions()

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

Run the simulation¶

In [5]:

dynamics.set_conc({"A": 50.}, snapshot=True)  # Set the initial concentrations of the chemicals
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: {'A', 'B', 'C'}

In [6]:

dynamics.single_compartment_react(initial_step=0.01, duration=0.8,
                                  variable_steps=True)

102 total step(s) taken in 0.200 sec
Number of step re-do's because of elective soft aborts: 1
Norm usage: {'norm_A': 15, 'norm_B': 16, 'norm_C': 14, 'norm_D': 14}
System Time is now: 0.81029

In [7]:

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

In [8]:

dynamics.is_in_equilibrium()

0: A <-> B
Final concentrations: [A] = 0.08481 ; [B] = 0.06976
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 0.822525
    Formula used:  [B] / [A]
2. Ratio of forward/reverse reaction rates: 4
Discrepancy between the two values: 79.44 %
Reaction is NOT in equilibrium (not within 1% tolerance)

1: B <-> C
Final concentrations: [B] = 0.06976 ; [C] = 49.85
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 714.525
    Formula used:  [C] / [B]
2. Ratio of forward/reverse reaction rates: 800
Discrepancy between the two values: 10.68 %
Reaction is NOT in equilibrium (not within 1% tolerance)

Out[8]:

{False: [0, 1]}

We didn't quite advance to equilibrium this time. A separate run (not shown) displayed far-better values if we had progressed the simulation just a little more in time.

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

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

Out[9]:

	SYSTEM TIME	A	C	caption
0	0.000000	50.000000	0.000000	Set concentration
1	0.004000	48.400000	0.000000	1st reaction step
2	0.008000	46.864000	0.512000
3	0.010000	46.124672	0.931738
4	0.011000	45.761562	1.167132
...	...	...	...	...
98	0.758938	0.121082	49.805302
99	0.771777	0.110536	49.816969
100	0.784616	0.101044	49.827470
101	0.797455	0.092501	49.836921
102	0.810294	0.084812	49.845428	last reaction step

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

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

In [11]:

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

In [12]:

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

ReactionKinetics.estimate_rate_constants_simple(t=t_arr, A_conc=A_conc, B_conc=C_conc, reactant_name="A", product_name="C")

Reaction A <-> C
Total REACTANT + PRODUCT has a median of 49.19, 
    with standard deviation 1.519 (ideally should be zero)
The sum of the time derivatives of the reactant and the product 
    has a median of 1.073 (ideally should be zero)
Least square fit to model as elementary reaction: C'(t) = kF * A(t) - kR * C(t)

-> ESTIMATED RATE CONSTANTS: kF = 6.549 , kR = -0.2912

The least-square fit is awful : the complex reaction `A <-> C` doesn't seem to be amenable to being modeled as a simple reaction with some suitable rate constants¶

But it looks like we'll do much better if splitting into 2 portions, one where A(t) ranges from 0 to about 40, and one from about 40 to 50¶

Indeed, revisiting the early portion of the time plot from Part 1, one can see an inflection in the [C] green curve roughly around time t=0.028, which is when [A] is around 40 (turquoise). That's around the peak of the mystery intermediate B (orange).

We'll pick time t=0.028 as the divider between the 2 domains of the A <-> C time evolution that we want to model.

Note that this is a much smaller time than we saw in experiment cascade_2_a

In [15]:

dynamics.plot_history(colors=['darkturquoise', 'orange', 'green'], range_x=[0, 0.4],
                      vertical_lines_to_add=[0.028])

Let's locate where the t = 0.028 point occurs in the data¶

In [16]:

dynamics.get_history(t=0.028)

Out[16]:

	search_value	SYSTEM TIME	A	B	C	caption
19	0.028	0.02728	40.260772	3.8859	5.853329

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

points numbered 0-19
points 19-end

In [17]:

A_conc_early = A_conc[:20]
A_conc_late = A_conc[20:]

C_conc_early = C_conc[:20]
C_conc_late = C_conc[20:]

t_arr_early = t_arr[:20]
t_arr_late = t_arr[20:]

In [ ]:

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

In [18]:

ReactionKinetics.estimate_rate_constants_simple(t=t_arr_early, A_conc=A_conc_early, B_conc=C_conc_early,
                                        reactant_name="A", product_name="C")

Reaction A <-> C
Total REACTANT + PRODUCT has a median of 46.44, 
    with standard deviation 0.9145 (ideally should be zero)
The sum of the time derivatives of the reactant and the product 
    has a median of -62.86 (ideally should be zero)
Least square fit to model as elementary reaction: C'(t) = kF * A(t) - kR * C(t)

-> ESTIMATED RATE CONSTANTS: kF = 2.946 , kR = -43.72

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

In [ ]:

II. And now let's consider the LATE region, when t > 0.028¶

In [19]:

ReactionKinetics.estimate_rate_constants_simple(t=t_arr_late, A_conc=A_conc_late, B_conc=C_conc_late,
                                        reactant_name="A", product_name="C")

Reaction A <-> C
Total REACTANT + PRODUCT has a median of 49.59, 
    with standard deviation 1.353 (ideally should be zero)
The sum of the time derivatives of the reactant and the product 
    has a median of 2.347 (ideally should be zero)
Least square fit to model as elementary reaction: C'(t) = kF * A(t) - kR * C(t)

-> ESTIMATED RATE CONSTANTS: kF = 8.427 , kR = -0.04512

This time we have an adequate linear fit AND meaningful rate constants : kF of about 8 and kR of about 0. Do those numbers sound familiar? A definite resemblance to the kF=8, kR=2 of the SLOWER elementary reaction A <-> B!

The slower `A <-> B` reaction dominates the kinetics from about t=0.028 on.¶

Let's see the concentrations over time again:

In [20]:

dynamics.plot_history(colors=['darkturquoise', 'orange', 'green'], range_x=[0, 0.4],
                      vertical_lines_to_add=[0.028])

Just as we concluded in experiment cascade_2_a, the earlier part of the complex (compound) reaction A <-> C cannot be modeled by an elementary reaction, while the later part can indeed be modeled by a 1st order elementary reaction, with kinetics similar to the slower A <-> B reaction. This time, with a greater disparity between the two elementary reaction, the transition happens much sooner.

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, and in the previous one, that the complex reaction could be roughly modeled with the rate constants of the slower reaction only after some time - especially if the 2 elementary reactions are relatively similar (as in the previous experiment `cascade_2_a`).¶

If we were interested in early transients (for example, if diffusion quickly intervened), we couldn't use that model.

In [ ]:

NEXT: In the continuation experiment, `cascade_2_c`, we explore the scenario where the 2 elementary reactions are reversed in their relative speeds¶

In [ ]:

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

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

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¶

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

The least-square fit is awful : the complex reaction A <-> C doesn't seem to be amenable to being modeled as a simple reaction with some suitable rate constants¶

But it looks like we'll do much better if splitting into 2 portions, one where A(t) ranges from 0 to about 40, and one from about 40 to 50¶

Let's locate where the t = 0.028 point occurs in the data¶

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

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

II. And now let's consider the LATE region, when t > 0.028¶

The slower A <-> B reaction dominates the kinetics from about t=0.028 on.¶

NEXT: In the continuation experiment, cascade_2_c, we explore the scenario where the 2 elementary reactions are reversed in their relative speeds¶

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

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

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

The least-square fit is awful : the complex reaction `A <-> C` doesn't seem to be amenable to being modeled as a simple reaction with some suitable rate constants¶

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

The slower `A <-> B` reaction dominates the kinetics from about t=0.028 on.¶

NEXT: In the continuation experiment, `cascade_2_c`, we explore the scenario where the 2 elementary reactions are reversed in their relative speeds¶