Reaction-Diffusion in 1-D: `A + B <-> C` in 1-D, taken to equilibrium¶

Mostly forward reaction, with 1st-order kinetics for each species¶

Initial concentrations of A and B are spatially separated to the opposite ends of the system; as a result, no C is being generated.

But, as soon as A and B, from their respective distant originating points at the edges, diffuse into the middle - and into each other - the reaction starts, consuming both A and B, until an equilibrium is reached in both diffusion and reactions.

TAGS : "reactions 1D", "diffusion 1D", "quick-start"¶

In [1]:

LAST_REVISED = "Jan. 24, 2024"
LIFE123_VERSION = "1.0.0rc2"        # 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 check_version, BioSim1D, ChemData, UniformCompartment

In [4]:

check_version(LIFE123_VERSION)

OK

In [ ]:

In [5]:

# Initialize the system
chem_data = ChemData(names=["A", "B", "C"], diffusion_rates=[50., 50., 1.],   # `A` and `B` diffuse fast; `C` diffuses slowly
                     plot_colors=["red", "blue", "purple"])                   # Color choice is a reminder that red + blue = purple

In [6]:

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

In [7]:

reactions = bio.get_reactions()

# Reaction A + B <-> C , with 1st-order kinetics for each species; note that it's mostly in the forward direction
reactions.add_reaction(reactants=["A", "B"], products="C", forward_rate=20., reverse_rate=2.)
reactions.describe_reactions()

Number of reactions: 1 (at temp. 25 C)
0: A + B <-> C  (kF = 20 / kR = 2 / delta_G = -5,708 / K = 10) | 1st order in all reactants & products
Set of chemicals involved in the above reactions: {"B" (blue), "A" (red), "C" (purple)}

In [ ]:

TIME 0 : Inject initial concentrations of `A` and `B` at opposite ends of the system¶

In [8]:

bio.set_bin_conc(bin_address=0, species_name="A", conc=20.)
bio.set_bin_conc(bin_address=6, species_name="B", conc=20.)

bio.describe_state()

SYSTEM STATE at Time t = 0:
7 bins and 3 chemical species:

Out[8]:

	Species	Diff rate	Bin 0	Bin 6
0	A	50.0	20.0	0.0
1	B	50.0	0.0	20.0
2	C	1.0	0.0	0.0

In [9]:

bio.show_system_snapshot()     # A more streamlined alternate view

SYSTEM SNAPSHOT at time 0:

Out[9]:

	A	B
0	20.0	0.0
1	0.0	0.0
2	0.0	0.0
3	0.0	0.0
4	0.0	0.0
5	0.0	0.0
6	0.0	20.0

In [10]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [11]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Enable History¶

In [12]:

# Let's take a peek at the current concentrations of all chemicals in the bins with the initial concentration injections,
# as well as at the bin in the very center
bio.selected_concentrations(bins=[0, 6, 3])

Out[12]:

{0: {'A': 20.0, 'B': 0.0, 'C': 0.0},
 6: {'A': 0.0, 'B': 20.0, 'C': 0.0},
 3: {'A': 0.0, 'B': 0.0, 'C': 0.0}}

In [13]:

# Let's enable history for those same 3 bins
bio.enable_history(bins=[0, 6, 3], frequency=2, take_snapshot=True)     # Taking a snapshot to include the current initial state in the history

History enabled for bins [0, 6, 3] and chemicals None (None means 'all')

In [ ]:

First step : advance to time t=0.002 ¶

In [14]:

delta_t = 0.002   # This will be our time "quantum" for this experiment

In [15]:

# First step
bio.react_diffuse(time_step=delta_t, n_steps=1)
bio.describe_state()

SYSTEM STATE at Time t = 0.002:
7 bins and 3 chemical species:

Out[15]:

	Species	Diff rate	Bin 0	Bin 1	Bin 5	Bin 6
0	A	50.0	18.0	2.0	0.0	0.0
1	B	50.0	0.0	0.0	2.0	18.0
2	C	1.0	0.0	0.0	0.0	0.0

In [16]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0.002:

Out[16]:

	A	B
0	18.0	0.0
1	2.0	0.0
2	0.0	0.0
3	0.0	0.0
4	0.0	0.0
5	0.0	2.0
6	0.0	18.0

In [17]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [18]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Several more steps : advance to time t=0.016 ¶

In [19]:

# Continue with several delta_t steps
for _ in range(7):
    bio.react_diffuse(time_step=delta_t, n_steps=1)
    bio.describe_state(concise=True)

SYSTEM STATE at Time t = 0.004:
[[16.4  3.4  0.2  0.   0.   0.   0. ]
 [ 0.   0.   0.   0.   0.2  3.4 16.4]
 [ 0.   0.   0.   0.   0.   0.   0. ]]
SYSTEM STATE at Time t = 0.006:
[[15.1   4.38  0.5   0.02  0.    0.    0.  ]
 [ 0.    0.    0.    0.02  0.5   4.38 15.1 ]
 [ 0.    0.    0.    0.    0.    0.    0.  ]]
SYSTEM STATE at Time t = 0.008:
[[1.4028e+01 5.0640e+00 8.4000e-01 6.5984e-02 2.0000e-03 0.0000e+00
  0.0000e+00]
 [0.0000e+00 0.0000e+00 2.0000e-03 6.5984e-02 8.4000e-01 5.0640e+00
  1.4028e+01]
 [0.0000e+00 0.0000e+00 0.0000e+00 1.6000e-05 0.0000e+00 0.0000e+00
  0.0000e+00]]
SYSTEM STATE at Time t = 0.01:
[[1.31316000e+01 5.53800000e+00 1.18493120e+00 1.36813108e-01
  8.13120000e-03 2.00000000e-04 0.00000000e+00]
 [0.00000000e+00 2.00000000e-04 8.13120000e-03 1.36813108e-01
  1.18493120e+00 5.53800000e+00 1.31316000e+01]
 [0.00000000e+00 0.00000000e+00 6.72320000e-05 1.90027530e-04
  6.72320000e-05 0.00000000e+00 0.00000000e+00]]
SYSTEM STATE at Time t = 0.012:
[[1.23722400e+01 5.86200882e+00 1.51504114e+00 2.28008774e-01
  1.98211433e-02 9.28816000e-04 2.00000000e-05]
 [2.00000000e-05 9.28816000e-04 1.98211433e-02 2.28008774e-01
  1.51504114e+00 5.86200882e+00 1.23722400e+01]
 [0.00000000e+00 4.44384640e-05 4.52470702e-04 9.37489304e-04
  4.52470702e-04 4.44384640e-05 0.00000000e+00]]
SYSTEM STATE at Time t = 0.014:
[[1.17212070e+01 6.07811756e+00 1.81983529e+00 3.33817478e-01
  3.75512896e-02 2.50955578e-03 1.00983808e-04]
 [1.00983808e-04 2.50955578e-03 3.75512896e-02 3.33817478e-01
  1.81983529e+00 6.07811756e+00 1.17212070e+01]
 [9.98666893e-06 2.62777001e-04 1.65200869e-03 3.01131931e-03
  1.65200869e-03 2.62777001e-04 9.98666893e-06]]
SYSTEM STATE at Time t = 0.016:
[[1.11568507e+01 6.21598919e+00 2.09433486e+00 4.48347321e-01
  6.09468566e-02 5.16378807e-03 2.94534867e-04]
 [2.94534867e-04 5.16378807e-03 6.09468566e-02 4.48347321e-01
  2.09433486e+00 6.21598919e+00 1.11568507e+01]
 [5.77983875e-05 8.74133777e-04 4.37882730e-03 7.45120113e-03
  4.37882730e-03 8.74133777e-04 5.77983875e-05]]

In [20]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0.016:

Out[20]:

	A	B	C
0	11.156851	0.000295	0.000058
1	6.215989	0.005164	0.000874
2	2.094335	0.060947	0.004379
3	0.448347	0.448347	0.007451
4	0.060947	2.094335	0.004379
5	0.005164	6.215989	0.000874
6	0.000295	11.156851	0.000058

In [21]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

A is continuing to diffuse from the left.
B is continuing to diffuse from the right.
They're finally beginning to overlap in the middle bins!

In [22]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Several groups of longer runs : advance to time t=0.096 ¶

In [23]:

# Now, do several group of longer runs
for _ in range(4):
    print("\n\n+ 10 steps later:")
    bio.react_diffuse(time_step=delta_t, n_steps=10)
    bio.describe_state(concise=True)


+ 10 steps later:
SYSTEM STATE at Time t = 0.036:
[[7.92576526 5.95824446 3.31917619 1.31581241 0.35918764 0.07038175
  0.01251229]
 [0.01251229 0.07038175 0.35918764 1.31581241 3.31917619 5.95824446
  7.92576526]
 [0.01555131 0.07914514 0.24409632 0.36133444 0.24409632 0.07914514
  0.01555131]]


+ 10 steps later:
SYSTEM STATE at Time t = 0.056:
[[6.29600683 5.06118233 3.15970431 1.44600532 0.47532616 0.12019317
  0.03045075]
 [0.03045075 0.12019317 0.47532616 1.44600532 3.15970431 5.06118233
  6.29600683]
 [0.07498099 0.28735559 0.78250579 1.12144639 0.78250579 0.28735559
  0.07498099]]


+ 10 steps later:
SYSTEM STATE at Time t = 0.076:
[[5.13609493 4.21433746 2.7454954  1.35328847 0.4992179  0.14748071
  0.04559633]
 [0.04559633 0.14748071 0.4992179  1.35328847 2.7454954  4.21433746
  5.13609493]
 [0.1615343  0.52706799 1.31954666 1.8421909  1.31954666 0.52706799
  0.1615343 ]]


+ 10 steps later:
SYSTEM STATE at Time t = 0.096:
[[4.2209458  3.50170661 2.34873485 1.23042685 0.50075376 0.16821567
  0.0606064 ]
 [0.0606064  0.16821567 0.50075376 1.23042685 2.34873485 3.50170661
  4.2209458 ]
 [0.25982593 0.75464804 1.76534285 2.40897641 1.76534285 0.75464804
  0.25982593]]

In [24]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0.096:

Out[24]:

	A	B	C
0	4.220946	0.060606	0.259826
1	3.501707	0.168216	0.754648
2	2.348735	0.500754	1.765343
3	1.230427	1.230427	2.408976
4	0.500754	2.348735	1.765343
5	0.168216	3.501707	0.754648
6	0.060606	4.220946	0.259826

In [25]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

A is continuing to diffuse from the left.
B is continuing to diffuse from the right.
The overlap of A and B, especially in the central bins, is by now extensive, and the reaction is proceeding in earnest.
Notice the continue symmetry about the center of the system

In [26]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Advance to time t=0.336 ¶

In [27]:

# Continue the simulation
for _ in range(4):
    print("\n\n+++ 30 steps later:")
    bio.react_diffuse(time_step=delta_t, n_steps=30)
    bio.describe_state(concise=True)


+++ 30 steps later:
SYSTEM STATE at Time t = 0.156:
[[2.39736516 2.05453684 1.50075745 0.92781069 0.48529917 0.22472198
  0.11560109]
 [0.11560109 0.22472198 0.48529917 0.92781069 1.50075745 2.05453684
  2.39736516]
 [0.57202182 1.29674075 2.59761288 3.36115673 2.59761288 1.29674075
  0.57202182]]


+++ 30 steps later:
SYSTEM STATE at Time t = 0.216:
[[1.43652347 1.28435794 1.02979899 0.73906751 0.47257834 0.28028088
  0.18445971]
 [0.18445971 0.28028088 0.47257834 0.73906751 1.02979899 1.28435794
  1.43652347]
 [0.8597085  1.64384498 2.94653088 3.67276444 2.94653088 1.64384498
  0.8597085 ]]


+++ 30 steps later:
SYSTEM STATE at Time t = 0.276:
[[0.94369275 0.88396666 0.77535511 0.63019458 0.46959173 0.33300009
  0.25664867]
 [0.25664867 0.33300009 0.46959173 0.63019458 0.77535511 0.88396666
  0.94369275]
 [1.09382006 1.85282552 3.05530325 3.70365274 3.05530325 1.85282552
  1.09382006]]


+++ 30 steps later:
SYSTEM STATE at Time t = 0.336:
[[0.69798039 0.68111864 0.64213556 0.57196145 0.47435422 0.37878712
  0.32097915]
 [0.32097915 0.37878712 0.47435422 0.57196145 0.64213556 0.68111864
  0.69798039]
 [1.27482053 1.97696102 3.05404907 3.62102222 3.05404907 1.97696102
  1.27482053]]

In [28]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0.336:

Out[28]:

	A	B	C
0	0.697980	0.320979	1.274821
1	0.681119	0.378787	1.976961
2	0.642136	0.474354	3.054049
3	0.571961	0.571961	3.621022
4	0.474354	0.642136	3.054049
5	0.378787	0.681119	1.976961
6	0.320979	0.697980	1.274821

In [29]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [30]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Advance to time t=0.736 ¶

In [31]:

# Continue the simulation
for _ in range(4):
    print("\n+++++ 50 steps later:")
    bio.react_diffuse(time_step=delta_t, n_steps=50)
    bio.describe_state(concise=True)

+++++ 50 steps later:
SYSTEM STATE at Time t = 0.436:
[[0.53513401 0.5437576  0.54867719 0.53174198 0.48696931 0.43256384
  0.39647767]
 [0.39647767 0.43256384 0.48696931 0.53174198 0.54867719 0.5437576
  0.53513401]
 [1.49508439 2.09028492 2.96849022 3.41695934 2.96849022 2.09028492
  1.49508439]]

+++++ 50 steps later:
SYSTEM STATE at Time t = 0.536:
[[0.48751897 0.50165856 0.51768297 0.51798994 0.49498758 0.46076288
  0.4365229 ]
 [0.4365229  0.46076288 0.49498758 0.51798994 0.51768297 0.50165856
  0.48751897]
 [1.65797359 2.15520413 2.86719131 3.22213812 2.86719131 2.15520413
  1.65797359]]

+++++ 50 steps later:
SYSTEM STATE at Time t = 0.636:
[[0.47470241 0.48869324 0.50591479 0.51120986 0.49756771 0.47365233
  0.45594669]
 [0.45594669 0.47365233 0.49756771 0.51120986 0.50591479 0.48869324
  0.47470241]
 [1.78774781 2.20070508 2.77778668 3.05983384 2.77778668 2.20070508
  1.78774781]]

+++++ 50 steps later:
SYSTEM STATE at Time t = 0.736:
[[0.47257504 0.48488523 0.50040062 0.50652836 0.49733067 0.47935337
  0.46567693]
 [0.46567693 0.47935337 0.49733067 0.50652836 0.50040062 0.48488523
  0.47257504]
 [1.89367665 2.23536347 2.70338342 2.9284027  2.70338342 2.23536347
  1.89367665]]

In [32]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0.736:

Out[32]:

	A	B	C
0	0.472575	0.465677	1.893677
1	0.484885	0.479353	2.235363
2	0.500401	0.497331	2.703383
3	0.506528	0.506528	2.928403
4	0.497331	0.500401	2.703383
5	0.479353	0.484885	2.235363
6	0.465677	0.472575	1.893677

In [33]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [34]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Advance to time t=1.936 ¶

In [35]:

# Continue the simulation
for _ in range(4):
    print("\n+++++++++++++++ 150 steps later:")
    bio.react_diffuse(time_step=delta_t, n_steps=150)
    bio.describe_state(concise=True)

+++++++++++++++ 150 steps later:
SYSTEM STATE at Time t = 1.036:
[[0.47744061 0.4845086  0.49346999 0.49747162 0.49331726 0.48423339
  0.47709743]
 [0.47709743 0.48423339 0.49331726 0.49747162 0.49346999 0.4845086
  0.47744061]
 [2.11014111 2.3007984  2.55138084 2.66782039 2.55138084 2.3007984
  2.11014111]]

+++++++++++++++ 150 steps later:
SYSTEM STATE at Time t = 1.336:
[[0.48168804 0.48554921 0.49041386 0.49259511 0.49040626 0.48553552
  0.48167096]
 [0.48167096 0.48553552 0.49040626 0.49259511 0.49041386 0.48554921
  0.48168804]
 [2.2288825  2.33377315 2.46847456 2.52988061 2.46847456 2.33377315
  2.2288825 ]]

+++++++++++++++ 150 steps later:
SYSTEM STATE at Time t = 1.636:
[[0.48405962 0.48615505 0.48878337 0.48995858 0.48878299 0.48615437
  0.48405877]
 [0.48405877 0.48615437 0.48878299 0.48995858 0.48878337 0.48615505
  0.48405962]
 [2.29360822 2.350868   2.42345244 2.45618991 2.42345244 2.350868
  2.29360822]]

+++++++++++++++ 150 steps later:
SYSTEM STATE at Time t = 1.936:
[[0.48534444 0.4864795  0.48789956 0.48853323 0.48789955 0.48647947
  0.4853444 ]
 [0.4853444  0.48647947 0.48789955 0.48853323 0.48789956 0.4864795
  0.48534444]
 [2.32876222 2.35988613 2.39905402 2.4166151  2.39905402 2.35988613
  2.32876222]]

In [36]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 1.936:

Out[36]:

	A	B	C
0	0.485344	0.485344	2.328762
1	0.486480	0.486479	2.359886
2	0.487900	0.487900	2.399054
3	0.488533	0.488533	2.416615
4	0.487900	0.487900	2.399054
5	0.486479	0.486480	2.359886
6	0.485344	0.485344	2.328762

In [37]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [38]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Advance to time t=5.936 ¶

In [39]:

# Continue the simulation
for _ in range(2):
    print("\n++++++++++ ... ++++++++++ 1,000 steps later:")
    bio.react_diffuse(time_step=delta_t, n_steps=1000)
    bio.describe_state(concise=True)

++++++++++ ... ++++++++++ 1,000 steps later:
SYSTEM STATE at Time t = 3.936:
[[0.48683089 0.48684974 0.48687325 0.48688372 0.48687325 0.48684974
  0.48683089]
 [0.48683089 0.48684974 0.48687325 0.48688372 0.48687325 0.48684974
  0.48683089]
 [2.36959744 2.37011659 2.37076408 2.37105229 2.37076408 2.37011659
  2.36959744]]

++++++++++ ... ++++++++++ 1,000 steps later:
SYSTEM STATE at Time t = 5.936:
[[0.48685551 0.48685582 0.48685621 0.48685639 0.48685621 0.48685582
  0.48685551]
 [0.48685551 0.48685582 0.48685621 0.48685639 0.48685621 0.48685582
  0.48685551]
 [2.37027551 2.37028411 2.37029484 2.37029961 2.37029484 2.37028411
  2.37027551]]

In [40]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 5.936:

Out[40]:

	A	B	C
0	0.486856	0.486856	2.370276
1	0.486856	0.486856	2.370284
2	0.486856	0.486856	2.370295
3	0.486856	0.486856	2.370300
4	0.486856	0.486856	2.370295
5	0.486856	0.486856	2.370284
6	0.486856	0.486856	2.370276

In [41]:

bio.visualize_system(title_prefix="Reaction-Diffusion A + B <-> C")   # Line curve view

In [42]:

bio.system_heatmap(title_prefix="Reaction-Diffusion A + B <-> C")

In [ ]:

Equilibrium in both diffusion and reaction¶

All bins now have essentially uniform concentration, for each of the chemicals

In each bin at equilibrium, [A] = 0.49, [B] = 0.49, [C] = 2.37
Let's verify that it's consistent with our (mostly-forward) equation A + B <-> C:

In [43]:

bio.reaction_in_equilibrium(bin_address=0, rxn_index=0, explain=True)  # Choice of bin is immaterial now, because they have all equilibrated

A + B <-> C
Current concentrations: [A] = 0.4869 ; [B] = 0.4869 ; [C] = 2.37
1. Ratio of reactant/product concentrations, adjusted for reaction orders: 9.99997
    Formula used:  [C] / ([A][B])
2. Ratio of forward/reverse reaction rates: 10
Discrepancy between the two values: 0.0003116 %
Reaction IS in equilibrium (within 1% tolerance)

Out[43]:

True

In [ ]:

Mass conservation: The initial "40 units of concentration" (20 from A and 20 from B) can be located by adding up the final concentrations in the 7 bins.
Notice that the concentration of C needs to be multiplied by 2, because the production of each molecule of C consumed 2 of the original molecules (one A and one B)

In [44]:

(2.37 * 2 + 0.49 + 0.49) * 7

Out[44]:

40.040000000000006

In [ ]:

Plots of changes of concentration with time¶

In [45]:

bio.plot_history_single_bin(bin_address=0)

In [46]:

bio.plot_history_single_bin(bin_address=6)

In [47]:

bio.plot_history_single_bin(bin_address=3)

A and B overlap on the plot, due to the symmetry of the system.
Initially, in the middle bin, neither A nor B are present; over time they diffuse there... but then they react and get consumed (producing C), to an equilibrium value.
Meanwhile, C gradually diffuses to uniformity.

In [ ]:

Reaction-Diffusion in 1-D: A + B <-> C in 1-D, taken to equilibrium¶

Mostly forward reaction, with 1st-order kinetics for each species¶

TAGS : "reactions 1D", "diffusion 1D", "quick-start"¶

TIME 0 : Inject initial concentrations of A and B at opposite ends of the system¶

Enable History¶

First step : advance to time t=0.002¶

Several more steps : advance to time t=0.016¶

Several groups of longer runs : advance to time t=0.096¶

Advance to time t=0.336¶

Advance to time t=0.736¶

Advance to time t=1.936¶

Advance to time t=5.936¶

Equilibrium in both diffusion and reaction¶

Plots of changes of concentration with time¶

Reaction-Diffusion in 1-D: `A + B <-> C` in 1-D, taken to equilibrium¶

TIME 0 : Inject initial concentrations of `A` and `B` at opposite ends of the system¶

First step : advance to time t=0.002 ¶

Several more steps : advance to time t=0.016 ¶

Several groups of longer runs : advance to time t=0.096 ¶

Advance to time t=0.336 ¶

Advance to time t=0.736 ¶

Advance to time t=1.936 ¶

Advance to time t=5.936 ¶