2 COUPLED reactions: `A + B <-> C` and `C + D <-> E` ,¶

with 1st-order kinetics for each species, taken to equilibrium¶

Both reactions are stronger in their respective forward rates. For the most part, "C" is produced by the 1st reaction, and consumed by the 2nd one

Diffusion not applicable (just 1 bin)

LAST REVISED: June 4, 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_data import ChemData as chem
from src.modules.reactions.reaction_dynamics import ReactionDynamics
from src.life_1D.bio_sim_1d import BioSim1D

import plotly.express as px
from src.modules.html_log.html_log import HtmlLog as log
from src.modules.visualization.graphic_log import GraphicLog

In [3]:

# Initialize the HTML logging
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_1"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")

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

In [4]:

# Initialize the system
chem_data = chem(names=["A", "B", "C", "D", "E"])     # NOTE: Diffusion not applicable (just 1 bin)

# Specify the reactions


# Reactions A + B <-> C  and  C + D <-> E , with 1st-order kinetics for each species
chem_data.add_reaction(reactants=["A", "B"], products=["C"], forward_rate=5., reverse_rate=2.)
chem_data.add_reaction(reactants=["C", "D"], products=["E"], forward_rate=8., reverse_rate=4.)

bio = BioSim1D(n_bins=1, chem_data=chem_data)

bio.set_all_uniform_concentrations( [3., 5., 1., 0.4, 0.1] )

bio.describe_state()

SYSTEM STATE at Time t = 0:
1 bins and 5 species:
  Species 0 (A). Diff rate: None. Conc: [3.]
  Species 1 (B). Diff rate: None. Conc: [5.]
  Species 2 (C). Diff rate: None. Conc: [1.]
  Species 3 (D). Diff rate: None. Conc: [0.4]
  Species 4 (E). Diff rate: None. Conc: [0.1]

In [5]:

# Save the state of the concentrations of all species at bin 0
bio.add_snapshot(bio.bin_snapshot(bin_address = 0), caption="Initial state")
bio.get_history()

Out[5]:

	SYSTEM TIME	A	B	C	D	E	caption
0	0	3.0	5.0	1.0	0.4	0.1	Initial state

In [6]:

chem_data.describe_reactions()

Number of reactions: 2 (at temp. 25 C)
0: A + B <-> C  (kF = 5 / kR = 2 / Delta_G = -2,271.45 / K = 2.5) | 1st order in all reactants & products
1: C + D <-> E  (kF = 8 / kR = 4 / Delta_G = -1,718.28 / K = 2) | 1st order in all reactants & products

In [7]:

# Send a header and a plot to the HTML log file
log.write("2 COUPLED reactions:  A + B <-> C  and  C + D <-> E",
          style=log.h2)
# Send the plot to the HTML log file
graph_data = chem_data.prepare_graph_network()
GraphicLog.export_plot(graph_data, "vue_cytoscape_1")

2 COUPLED reactions:  A + B  C  and  C + D  E
[GRAPHIC ELEMENT SENT TO LOG FILE `reaction_8.log.htm`]

First step ¶

In [8]:

# First step
bio.react(time_step=0.01, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()

SYSTEM STATE at Time t = 0.01:
1 bins and 5 species:
  Species 0 (A). Diff rate: None. Conc: [2.27]
  Species 1 (B). Diff rate: None. Conc: [4.27]
  Species 2 (C). Diff rate: None. Conc: [1.702]
  Species 3 (D). Diff rate: None. Conc: [0.372]
  Species 4 (E). Diff rate: None. Conc: [0.128]

1 bins and 5 species: [[2.27 ] [4.27 ] [1.702] [0.372] [0.128]]

In [9]:

bio.get_history()

Out[9]:

	SYSTEM TIME	A	B	C	D	E	caption
0	0.00	3.00	5.00	1.000	0.400	0.100	Initial state
1	0.01	2.27	4.27	1.702	0.372	0.128

In [10]:

# Identical 2nd step
bio.react(time_step=0.01, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()

SYSTEM STATE at Time t = 0.02:
1 bins and 5 species:
  Species 0 (A). Diff rate: None. Conc: [1.819395]
  Species 1 (B). Diff rate: None. Conc: [3.819395]
  Species 2 (C). Diff rate: None. Conc: [2.10707348]
  Species 3 (D). Diff rate: None. Conc: [0.32646848]
  Species 4 (E). Diff rate: None. Conc: [0.17353152]

1 bins and 5 species: [[1.819395 ] [3.819395 ] [2.10707348] [0.32646848] [0.17353152]]

In [11]:

bio.get_history()

Out[11]:

	SYSTEM TIME	A	B	C	D	E	caption
0	0.00	3.000000	5.000000	1.000000	0.400000	0.100000	Initial state
1	0.01	2.270000	4.270000	1.702000	0.372000	0.128000
2	0.02	1.819395	3.819395	2.107073	0.326468	0.173532

Numerous more steps ¶

In [12]:

# Numerous more identical steps, to equilibrium
bio.react(time_step=0.01, n_steps=200, snapshots={"sample_bin": 0, "frequency": 10})
bio.describe_state()

SYSTEM STATE at Time t = 2.02:
1 bins and 5 species:
  Species 0 (A). Diff rate: None. Conc: [0.50508029]
  Species 1 (B). Diff rate: None. Conc: [2.50508029]
  Species 2 (C). Diff rate: None. Conc: [3.16316668]
  Species 3 (D). Diff rate: None. Conc: [0.06824696]
  Species 4 (E). Diff rate: None. Conc: [0.43175304]

1 bins and 5 species: [[0.50508029] [2.50508029] [3.16316668] [0.06824696] [0.43175304]]

Equilibrium ¶

In [13]:

# Verify that each reaction has reached equilibrium
bio.reaction_dynamics.is_in_equilibrium(rxn_index=0, conc=bio.bin_snapshot(bin_address = 0))

A + B <-> C
Final concentrations:  [C] = 3.163 ; [A] = 0.5051 ; [B] = 2.505
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2.5
    Formula used:  [C] / ([A][B])
2. Ratio of forward/reverse reaction rates: 2.5
Discrepancy between the two values: 8.882e-14 %
Reaction IS in equilibrium (within 1% tolerance)

Out[13]:

True

In [14]:

bio.reaction_dynamics.is_in_equilibrium(rxn_index=1, conc=bio.bin_snapshot(bin_address = 0))

C + D <-> E
Final concentrations:  [E] = 0.4318 ; [C] = 3.163 ; [D] = 0.06825
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2
    Formula used:  [E] / ([C][D])
2. Ratio of forward/reverse reaction rates: 2.0
Discrepancy between the two values: 4.441e-14 %
Reaction IS in equilibrium (within 1% tolerance)

Out[14]:

True

In [15]:

# Do a consistent check with the equilibrium concentrations:

A_eq = bio.bin_concentration(0, 0)
B_eq = bio.bin_concentration(0, 1)
C_eq = bio.bin_concentration(0, 2)
D_eq = bio.bin_concentration(0, 3)
E_eq = bio.bin_concentration(0, 4)

Rf0 = chem_data.get_forward_rate(0)
Rb0 = chem_data.get_reverse_rate(0)

Rf1 = chem_data.get_forward_rate(1)
Rb1 = chem_data.get_reverse_rate(1)

equil = -(Rf0 * A_eq * B_eq - Rf1 * C_eq * D_eq) + (Rb0 * C_eq - Rb1 * E_eq)

print("\nAt equilibrium: ", equil, " (this should be close to 0 at equilibrium)")

At equilibrium:  -4.440892098500626e-15  (this should be close to 0 at equilibrium)

In [16]:

bio.get_history()

Out[16]:

	SYSTEM TIME	A	B	C	D	E	caption
0	0.00	3.000000	5.000000	1.000000	0.400000	0.100000	Initial state
1	0.01	2.270000	4.270000	1.702000	0.372000	0.128000
2	0.02	1.819395	3.819395	2.107073	0.326468	0.173532
3	0.12	0.654654	2.654654	3.032120	0.086774	0.413226
4	0.22	0.527627	2.527627	3.141994	0.069620	0.430380
5	0.32	0.508578	2.508578	3.159830	0.068408	0.431592
6	0.42	0.505625	2.505625	3.162646	0.068270	0.431730
7	0.52	0.505165	2.505165	3.163085	0.068251	0.431749
8	0.62	0.505094	2.505094	3.163154	0.068248	0.431752
9	0.72	0.505082	2.505082	3.163165	0.068247	0.431753
10	0.82	0.505081	2.505081	3.163166	0.068247	0.431753
11	0.92	0.505080	2.505080	3.163167	0.068247	0.431753
12	1.02	0.505080	2.505080	3.163167	0.068247	0.431753
13	1.12	0.505080	2.505080	3.163167	0.068247	0.431753
14	1.22	0.505080	2.505080	3.163167	0.068247	0.431753
15	1.32	0.505080	2.505080	3.163167	0.068247	0.431753
16	1.42	0.505080	2.505080	3.163167	0.068247	0.431753
17	1.52	0.505080	2.505080	3.163167	0.068247	0.431753
18	1.62	0.505080	2.505080	3.163167	0.068247	0.431753
19	1.72	0.505080	2.505080	3.163167	0.068247	0.431753
20	1.82	0.505080	2.505080	3.163167	0.068247	0.431753
21	1.92	0.505080	2.505080	3.163167	0.068247	0.431753
22	2.02	0.505080	2.505080	3.163167	0.068247	0.431753

Plots of changes of concentration with time¶

In [17]:

fig = px.line(data_frame=bio.get_history(), x="SYSTEM TIME", y=["A", "B", "C", "D", "E"], 
              title="2 COUPLED reactions:  A + B <-> C  and  C + D <-> E . Changes in concentrations",
              color_discrete_sequence = ['navy', 'cyan', 'red', 'orange', 'green'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()

A and B get consumed.
C gets produced by the 1st reaction more quickly than consumed by the 2nd one. D gets consumed, while E gets produced.

In [ ]:

2 COUPLED reactions: A + B <-> C and C + D <-> E ,¶

with 1st-order kinetics for each species, taken to equilibrium¶

First step¶

Numerous more steps¶

Equilibrium¶

Plots of changes of concentration with time¶

2 COUPLED reactions: `A + B <-> C` and `C + D <-> E` ,¶

First step ¶

Numerous more steps ¶

Equilibrium ¶