An `A <-> B` reaction¶

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

In Part 2, below, perform some analysis of the results: in particular, examine the reaction rates.

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

This experiment gets repeated in "react_2_b" , with a more sophisticated approach,¶

involving adaptive variable time steps¶

In [1]:

LAST_REVISED = "Aug. 29, 2024"
LIFE123_VERSION = "1.0.0rc5"      # 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   

import numpy as np
import ipynbname

from life123 import check_version, UniformCompartment, PlotlyHelper, GraphicLog

In [4]:

check_version(LIFE123_VERSION)

OK

In [5]:

# Initialize the HTML logging (for the graphics)
log_file = ipynbname.name() + ".log.htm"    # Use the notebook base filename for the log file
                                            # IN CASE OF PROBLEMS, set manually to any desired name

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

PART 1 - simulation of the reaction¶

In [6]:

# Instantiate the simulator and specify the chemicals
uc = UniformCompartment()

# Reaction A <-> B , with 1st-order kinetics in both directions
uc.add_reaction(reactants="A", products="B", kF=3., kR=2.)

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

Number of reactions:  1

In [7]:

uc.describe_reactions()

Number of reactions: 1
0: A <-> B  (Elementary Unimolecular reaction)  (kF = 3 / kR = 2 / delta_G = -1,005.1 / K = 1.5 / Temp = 25 C)
Chemicals involved in the above reactions: ['A', 'B']

In [8]:

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

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

In [ ]:

In [9]:

# Initial concentrations of all the chemicals
uc.set_conc({"A": 10., "B": 50.})

In [10]:

uc.describe_state()

SYSTEM STATE at Time t = 0:
2 species:
  Species 0 (A). Conc: 10.0
  Species 1 (B). Conc: 50.0
Chemicals involved in reactions: ['B', 'A']

In [11]:

uc.get_history()

Out[11]:

	SYSTEM TIME	A	B	step	caption
0	0.0	10.0	50.0		Set concentration

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

uc.find_equilibrium_conc(rxn_index=0)    # This is an EXACT equilibrium solution, 
                                         # for 1 reaction (the only reaction)

Out[12]:

{'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 [13]:

# First step of reaction
uc.single_compartment_react(initial_step=0.1, n_steps=1, variable_steps=False) # NOT using variable steps!

1 total fixed step(s) taken in 0.00 min

In [14]:

uc.get_history()

Out[14]:

	SYSTEM TIME	A	B	step	caption
0	0.0	10.0	50.0		Set concentration
1	0.1	17.0	43.0	1	last reaction step

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

In [15]:

# Numerous more fixed steps
uc.single_compartment_react(initial_step=0.1, n_steps=10, variable_steps=False)  # Again, fixed steps used

10 total fixed step(s) taken in 0.00 min

In [16]:

uc.get_history()

Out[16]:

	SYSTEM TIME	A	B	step	caption
0	0.0	10.000000	50.000000		Set concentration
1	0.1	17.000000	43.000000	1	last reaction step
2	0.2	20.500000	39.500000	1	1st reaction step
3	0.3	22.250000	37.750000	2
4	0.4	23.125000	36.875000	3
5	0.5	23.562500	36.437500	4
6	0.6	23.781250	36.218750	5
7	0.7	23.890625	36.109375	6
8	0.8	23.945312	36.054688	7
9	0.9	23.972656	36.027344	8
10	1.0	23.986328	36.013672	9
11	1.1	23.993164	36.006836	10	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 [17]:

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

Out[17]:

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

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

0: A <-> B
Current 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[18]:

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

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

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 in dashed lines), which generally gives poor results, unless taking a very large number of very small steps! In particular, the early steps we took were too large, and the later steps were unnecessarily small¶

In [ ]:

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

There's no need to compute this; it is automatically saved during the reaction simulations

In [20]:

rates_df = uc.get_rate_history()
rates_df

Out[20]:

	SYSTEM TIME	rxn0_rate	step
0	0.0	-70.000000	0
1	0.1	-35.000000	0
2	0.2	-17.500000	1
3	0.3	-8.750000	2
4	0.4	-4.375000	3
5	0.5	-2.187500	4
6	0.6	-1.093750	5
7	0.7	-0.546875	6
8	0.8	-0.273438	7
9	0.9	-0.136719	8
10	1.0	-0.068359	9

Note that reaction rates are defined for the reaction products; since A is a reactant (in reaction 0, our only reaction), we must flip its sign; since the stoichiometry of A is simply 1, no further adjustment needed.

In [21]:

rates_df['A_dot'] = - rates_df['rxn0_rate']    # Add a column to the Pandas dataframe
rates_df

Out[21]:

	SYSTEM TIME	rxn0_rate	step	A_dot
0	0.0	-70.000000	0	70.000000
1	0.1	-35.000000	0	35.000000
2	0.2	-17.500000	1	17.500000
3	0.3	-8.750000	2	8.750000
4	0.4	-4.375000	3	4.375000
5	0.5	-2.187500	4	2.187500
6	0.6	-1.093750	5	1.093750
7	0.7	-0.546875	6	0.546875
8	0.8	-0.273438	7	0.273438
9	0.9	-0.136719	8	0.136719
10	1.0	-0.068359	9	0.068359

In [22]:

uc.get_history()         # Revisited from earlier

Out[22]:

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

Notice that we lack a rate for the last time value, in the above table, because no reaction simulation starting at that time has been performed

In [23]:

p1 = uc.plot_history(chemicals="A", colors="darkturquoise")   # The plot of [A] from the system history

In [24]:

p2 = PlotlyHelper.plot_pandas(df=rates_df, x_var="SYSTEM TIME", fields="A_dot", colors="brown")   # The plot of A_dot, from rates_df

In [25]:

PlotlyHelper.combine_plots([p1, p2], 
                           title="Concentration of A with time, and its rate of change (A_dot)",
                           y_label="[A] (turquoise) /<br> A_dot (brown)",
                           legend_title="Plot",
                           curve_labels=["A", "A_dot"])

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 above curves are jagged because of 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 approach that employs adaptive variable time steps.¶

In [ ]:

An A <-> B reaction¶

This experiment gets repeated in "react_2_b" , with a more sophisticated approach,¶

involving adaptive variable time steps¶

PART 1 - simulation of the reaction¶

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 in dashed lines), which generally gives poor results, unless taking a very large number of very small steps! In particular, the early steps we took were too large, and the later steps were unnecessarily small¶

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

At t=0 :¶

As the system approaches equilibrium :¶

NOTE: The above curves are jagged because of 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 approach 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¶

In experiment `react_2_b`, we revisit the same reaction using a better approach that employs adaptive variable time steps.¶