An `A <-> B` reaction¶

with 1st-order kinetics in both directions, taken to equilibrium, using a simple, coarse fixed-timestep simulation.

Afterwards, perform some analysis of the results

See also the experiment "1D/reactions/reaction_1" for a multi-compartment version.

This experiment gets continued in "react_2_b" , with a more sophisticated approach, involving adaptive variable time steps.

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

import numpy as np
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
from src.modules.visualization.graphic_log import GraphicLog

In [3]:

# Initialize the HTML logging (for the graphics)
log_file = get_notebook_basename() + ".log.htm"    # Use the notebook base filename for the log file

# Set up the use of some specified graphic (Vue) components
GraphicLog.config(filename=log_file,
                  components=["vue_cytoscape_2"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

-> Output will be LOGGED into the file 'react_2_a.log.htm'

In [ ]:

Initialize the System¶

In [4]:

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

# Reaction A <-> B , with 1st-order kinetics in both directions
dynamics.add_reaction(reactants=["A"], products=["B"], 
                       forward_rate=3., reverse_rate=2.)

print("Number of reactions: ", dynamics.number_of_reactions())

Number of reactions:  1

In [5]:

dynamics.describe_reactions()

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

In [6]:

# Send a plot of the network of reactions to the HTML log file
dynamics.plot_reaction_network("vue_cytoscape_2")

[GRAPHIC ELEMENT SENT TO LOG FILE `react_2_a.log.htm`]

In [7]:

# Initial concentrations of all the chemicals, in index order
dynamics.set_conc([10., 50.])

In [8]:

dynamics.describe_state()

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

In [9]:

dynamics.get_history()

Out[9]:

	SYSTEM TIME	A	B	caption
0	0.0	10.0	50.0	Initialized state

Test your intuition:¶

given that this reaction operates mostly in the forward direction (kF = 3 , kR = 2 , K = 1.5),¶

do you think that A will be consumed and B will be produced??¶

We can take a sneak preview at the final equilibrium concentrations without actually running the simulation:

In [10]:

dynamics.find_equilibrium_conc(rxn_index=0)    # This is an EXACT solution

Out[10]:

{'A': 24.0, 'B': 36.0}

The reaction will actually proceed IN REVERSE, because of the large initial concentration of B (which we had set to 50), relative to the small initial concentration of A (10)¶

Now, let's see the reaction in action!

In [ ]:

Run the reaction¶

In [11]:

# First step of reaction
dynamics.single_compartment_react(initial_step=0.1, n_steps=1, variable_steps=False, 
                                  snapshots={"initial_caption": "first reaction step"})

1 total step(s) taken

In [12]:

dynamics.get_history()

Out[12]:

	SYSTEM TIME	A	B	caption
0	0.0	10.0	50.0	Initialized state
1	0.1	17.0	43.0	first reaction step

We can already see the reaction proceeding in reverse...

In [13]:

# Numerous more fixed steps
dynamics.single_compartment_react(initial_step=0.1, n_steps=10, variable_steps=False, 
                                  snapshots={"initial_caption": "2nd reaction step",
                                             "final_caption": "last reaction step"})

10 total step(s) taken

In [14]:

dynamics.get_history()

Out[14]:

	SYSTEM TIME	A	B	caption
0	0.0	10.000000	50.000000	Initialized state
1	0.1	17.000000	43.000000	first reaction step
2	0.2	20.500000	39.500000	2nd reaction step
3	0.3	22.250000	37.750000
4	0.4	23.125000	36.875000
5	0.5	23.562500	36.437500
6	0.6	23.781250	36.218750
7	0.7	23.890625	36.109375
8	0.8	23.945312	36.054688
9	0.9	23.972656	36.027344
10	1.0	23.986328	36.013672
11	1.1	23.993164	36.006836	last reaction step

NOTE: for demonstration purposes, we're using FIXED time steps...¶

Typically, one would use the option for adaptive variable time steps (see experiment `react_2_b`)¶

Check the final equilibrium¶

In [15]:

dynamics.get_system_conc()   # The current concentrations, in the order the chemicals were added 

Out[15]:

array([23.99316406, 36.00683594])

NOTE: Consistent with the 3/2 ratio of forward/reverse rates (and the 1st order of the reactions), the systems settles in the following equilibrium:

[A] = 23.99316406

[B] = 36.00683594

In [16]:

# Verify that the reaction has reached equilibrium
dynamics.is_in_equilibrium()

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

Out[16]:

True

As noted earlier, because of the high initial concentration of B relative to A, the overall reaction has proceeded IN REVERSE¶

Plots of changes of concentration with time¶

In [17]:

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

Note the raggedness of the left-side (early times) of the curves.¶

In experiment `react_2_b` this simulation gets repeated with an adaptive variable time resolution that takes smaller steps at the beginning, when the reaction is proceeding faster¶

By contrast, here we used FIXED time steps (shown below), which generally gives poor results, unless taking a very large number of very small steps!¶

In [18]:

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

In [19]:

df = dynamics.get_history()         # Revisited from earlier
df

Out[19]:

	SYSTEM TIME	A	B	caption
0	0.0	10.000000	50.000000	Initialized state
1	0.1	17.000000	43.000000	first reaction step
2	0.2	20.500000	39.500000	2nd reaction step
3	0.3	22.250000	37.750000
4	0.4	23.125000	36.875000
5	0.5	23.562500	36.437500
6	0.6	23.781250	36.218750
7	0.7	23.890625	36.109375
8	0.8	23.945312	36.054688
9	0.9	23.972656	36.027344
10	1.0	23.986328	36.013672
11	1.1	23.993164	36.006836	last reaction step

Now investigate A_dot, i.e. d[A]/dt¶

In [20]:

A = list(df.A)
A

Out[20]:

[10.0,
 17.0,
 20.5,
 22.25,
 23.125,
 23.5625,
 23.78125,
 23.890625,
 23.9453125,
 23.97265625,
 23.986328125,
 23.9931640625]

In [21]:

len(A)

Out[21]:

In [22]:

A_dot = np.gradient(A, 0.1)      # 0.1 is the constant step size

In [23]:

A_dot

Out[23]:

array([7.00000000e+01, 5.25000000e+01, 2.62500000e+01, 1.31250000e+01,
       6.56250000e+00, 3.28125000e+00, 1.64062500e+00, 8.20312500e-01,
       4.10156250e-01, 2.05078125e-01, 1.02539062e-01, 6.83593750e-02])

In [24]:

df['A_dot'] = A_dot     # Add a column to the Pandas dataframe

In [25]:

df

Out[25]:

	SYSTEM TIME	A	B	caption	A_dot
0	0.0	10.000000	50.000000	Initialized state	70.000000
1	0.1	17.000000	43.000000	first reaction step	52.500000
2	0.2	20.500000	39.500000	2nd reaction step	26.250000
3	0.3	22.250000	37.750000		13.125000
4	0.4	23.125000	36.875000		6.562500
5	0.5	23.562500	36.437500		3.281250
6	0.6	23.781250	36.218750		1.640625
7	0.7	23.890625	36.109375		0.820312
8	0.8	23.945312	36.054688		0.410156
9	0.9	23.972656	36.027344		0.205078
10	1.0	23.986328	36.013672		0.102539
11	1.1	23.993164	36.006836	last reaction step	0.068359

In [26]:

dynamics.plot_history(chemicals=["A", "A_dot"], colors=['navy', 'brown'], 
                      ylabel="concentration (darkturquoise) /<br> concentration change per unit time (brown)",
                      title="Concentration of A with time (darkturquoise), and its rate of change (brown)")

At t=0 :¶

[A]=10 and [A] has a high rate of change (70)

As the system approaches equilibrium :¶

[A] approaches a value of 24, and its rate of change decays to zero.

NOTE: The curves are jagged because of limitations of numerically estimating derivatives, as well as the large time steps taken (especially in the early times, when there's a lot of change.)¶

In experiment "react_2_b", we revisit the same reaction using a better simulator that employs adaptive variable time steps¶

In [ ]:

An A <-> B reaction¶

Initialize the System¶

Test your intuition:¶

given that this reaction operates mostly in the forward direction (kF = 3 , kR = 2 , K = 1.5),¶

do you think that A will be consumed and B will be produced??¶

The reaction will actually proceed IN REVERSE, because of the large initial concentration of B (which we had set to 50), relative to the small initial concentration of A (10)¶

Run the reaction¶

NOTE: for demonstration purposes, we're using FIXED time steps...¶

Typically, one would use the option for adaptive variable time steps (see experiment react_2_b)¶

Check the final equilibrium¶

As noted earlier, because of the high initial concentration of B relative to A, the overall reaction has proceeded IN REVERSE¶

Plots of changes of concentration with time¶

Note the raggedness of the left-side (early times) of the curves.¶

In experiment react_2_b this simulation gets repeated with an adaptive variable time resolution that takes smaller steps at the beginning, when the reaction is proceeding faster¶

By contrast, here we used FIXED time steps (shown below), which generally gives poor results, unless taking a very large number of very small steps!¶

Now investigate A_dot, i.e. d[A]/dt¶

At t=0 :¶

As the system approaches equilibrium :¶

NOTE: The curves are jagged because of limitations of numerically estimating derivatives, as well as the large time steps taken (especially in the early times, when there's a lot of change.)¶

In experiment "react_2_b", we revisit the same reaction using a better simulator that employs adaptive variable time steps¶

An `A <-> B` reaction¶

Typically, one would use the option for adaptive variable time steps (see experiment `react_2_b`)¶

In experiment `react_2_b` this simulation gets repeated with an adaptive variable time resolution that takes smaller steps at the beginning, when the reaction is proceeding faster¶