One-bin `2A + 5B <-> 4C + 3D`, with 1st-order kinetics for each species, taken to equilibrium¶

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_6.log.htm'

In [4]:

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

# Specify the reaction

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

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

bio.set_all_uniform_concentrations( [4., 7., 5., 2.] )

bio.describe_state()

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

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	caption
0	0	4.0	7.0	5.0	2.0	Initial state

In [6]:

chem_data.describe_reactions()

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

In [7]:

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

Reaction 2 A + 5 B  4 C + 3 D
[GRAPHIC ELEMENT SENT TO LOG FILE `reaction_6.log.htm`]

First step ¶

In [8]:

# First step
bio.react(time_step=0.001, n_steps=1)
bio.describe_state()

SYSTEM STATE at Time t = 0.001:
1 bins and 4 species:
  Species 0 (A). Diff rate: None. Conc: [3.76]
  Species 1 (B). Diff rate: None. Conc: [6.4]
  Species 2 (C). Diff rate: None. Conc: [5.48]
  Species 3 (D). Diff rate: None. Conc: [2.36]

Early in the reaction : [A] = 3.76 , [B] = 6.4 , [C] = 5.48 , [D] = 2.36

In [9]:

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

Out[9]:

	SYSTEM TIME	A	B	C	D	caption
0	0.000	4.00	7.0	5.00	2.00	Initial state
1	0.001	3.76	6.4	5.48	2.36

Numerous more steps ¶

In [10]:

# Numerous more steps
bio.react(time_step=0.001, n_steps=40, snapshots={"sample_bin": 0})

bio.describe_state()

SYSTEM STATE at Time t = 0.041:
1 bins and 4 species:
  Species 0 (A). Diff rate: None. Conc: [2.80284552]
  Species 1 (B). Diff rate: None. Conc: [4.00711381]
  Species 2 (C). Diff rate: None. Conc: [7.39430896]
  Species 3 (D). Diff rate: None. Conc: [3.79573172]

Consistent with the 5/2 ratio of forward/reverse rates (and the 1st order reactions), the systems settles in the following equilibrium:
[A] = 2.80284552 , [B] = 4.00711381 , [C] = 7.39430896 , [D] = 3.79573172

In [11]:

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

2 A + 5 B <-> 4 C + 3 D
Final concentrations:  [C] = 7.394 ; [D] = 3.796 ; [A] = 2.803 ; [B] = 4.007
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 2.49898
    Formula used:  ([C][D]) / ([A][B])
2. Ratio of forward/reverse reaction rates: 2.5
Discrepancy between the two values: 0.04092 %
Reaction IS in equilibrium (within 1% tolerance)

Out[11]:

True

In [12]:

df = bio.get_history()
df

Out[12]:

	SYSTEM TIME	A	B	C	D	caption
0	0.000	4.000000	7.000000	5.000000	2.000000	Initial state
1	0.001	3.760000	6.400000	5.480000	2.360000
2	0.002	3.571091	5.927728	5.857818	2.643363
3	0.003	3.421344	5.553360	6.157312	2.867984
4	0.004	3.301981	5.254952	6.396039	3.047029
5	0.005	3.206419	5.016047	6.587162	3.190372
6	0.006	3.129645	4.824113	6.740709	3.305532
7	0.007	3.067794	4.669486	6.864411	3.398309
8	0.008	3.017854	4.544634	6.964293	3.473220
9	0.009	2.977457	4.443643	7.045085	3.533814
10	0.010	2.944734	4.361834	7.110532	3.582899
11	0.011	2.918195	4.295487	7.163611	3.622708
12	0.012	2.896651	4.241627	7.206699	3.655024
13	0.013	2.879148	4.197870	7.241704	3.681278
14	0.014	2.864920	4.162300	7.270160	3.702620
15	0.015	2.853348	4.133370	7.293304	3.719978
16	0.016	2.843932	4.109831	7.312135	3.734101
17	0.017	2.836269	4.090672	7.327463	3.745597
18	0.018	2.830029	4.075073	7.339942	3.754956
19	0.019	2.824948	4.062370	7.350104	3.762578
20	0.020	2.820809	4.052024	7.358381	3.768786
21	0.021	2.817438	4.043596	7.365123	3.773843
22	0.022	2.814692	4.036729	7.370617	3.777962
23	0.023	2.812454	4.031135	7.375092	3.781319
24	0.024	2.810630	4.026576	7.378739	3.784054
25	0.025	2.809144	4.022861	7.381711	3.786283
26	0.026	2.807933	4.019834	7.384133	3.788100
27	0.027	2.806947	4.017366	7.386107	3.789580
28	0.028	2.806142	4.015355	7.387716	3.790787
29	0.029	2.805487	4.013717	7.389027	3.791770
30	0.030	2.804952	4.012381	7.390095	3.792572
31	0.031	2.804517	4.011292	7.390966	3.793225
32	0.032	2.804162	4.010405	7.391676	3.793757
33	0.033	2.803872	4.009681	7.392255	3.794191
34	0.034	2.803637	4.009092	7.392727	3.794545
35	0.035	2.803444	4.008611	7.393111	3.794833
36	0.036	2.803288	4.008219	7.393424	3.795068
37	0.037	2.803160	4.007900	7.393680	3.795260
38	0.038	2.803056	4.007640	7.393888	3.795416
39	0.039	2.802971	4.007428	7.394058	3.795543
40	0.040	2.802902	4.007255	7.394196	3.795647
41	0.041	2.802846	4.007114	7.394309	3.795732

A and B get depleted, while C and D get produced.

2A + 5B <-> 4C + 3D

Let's verify that the stoichiometry is being respected¶

In [13]:

# We'll check the first two arrays of concentrations, from the run's history
arr0 = bio.reaction_dynamics.get_historical_concentrations(row=0, df=df)
arr1 = bio.reaction_dynamics.get_historical_concentrations(row=1, df=df)
arr0, arr1

Out[13]:

(array([4., 7., 5., 2.], dtype=float32),
 array([3.76, 6.4 , 5.48, 2.36], dtype=float32))

In [14]:

bio.reaction_dynamics.stoichiometry_checker(rxn_index=0, 
                               conc_arr_before = arr0, 
                               conc_arr_after = arr1)

Out[14]:

True

Indeed, the change in [A] is -2 x 0.12, and the change in [B] is -5 X 0.12,
while the change in [C] is 4 x 0.12, and the change in [D] is 3 X 0.12

In [15]:

(arr1 - arr0) / 0.12

Out[15]:

array([-2.0000002, -4.9999995,  4.0000005,  2.9999993], dtype=float32)

Plots of changes of concentration with time¶

In [16]:

fig = px.line(data_frame=bio.get_history(), x="SYSTEM TIME", y=["A", "B", "C", "D"], 
              title="Changes in concentrations with time",
              color_discrete_sequence = ['navy', 'cyan', 'red', 'orange'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()

In [ ]:

One-bin 2A + 5B <-> 4C + 3D, with 1st-order kinetics for each species, taken to equilibrium¶

First step¶

Numerous more steps¶

Equilibrium¶

Let's verify that the stoichiometry is being respected¶

Plots of changes of concentration with time¶

One-bin `2A + 5B <-> 4C + 3D`, with 1st-order kinetics for each species, taken to equilibrium¶

First step ¶

Numerous more steps ¶

Equilibrium ¶