Merge pull request #69 from BrainAnnex/dev

Dev beta 38

Author: BrainAnnex

experiments/reactions_single_compartment/cascade_1.ipynb

@@ -21,10 +21,9 @@
 # Adaptive variable time resolution is used, with extensive diagnostics, 
 # and finally compared to a new run using fixed time intervals, with the same initial data.
 #
-# In part2, some diagnotic insight is explored.   
-# In part3, two identical runs ("adaptive variable steps" and "fixed small steps") are compared. 
+# In part2, some diagnostic insight is explored.  
 #
-# LAST REVISED: June 23, 2024 (using v. 1.0 beta34.1)
+# In part3, two identical runs ("adaptive variable steps" vs. "fixed small steps") are compared. 
 
 # %% [markdown]
 # ## Bathtub analogy:
@@ -41,20 +40,28 @@
 # ![2 Coupled Reactions](../../docs/2_coupled_reactions.png)
 
 # %%
-import set_path      # Importing this module will add the project's home directory to sys.path
+LAST_REVISED = "July 26, 2024"
+LIFE123_VERSION = "1.0.0.beta.38"      # Version this experiment is based on
+
+# %%
+#import set_path            # Using MyBinder?  Uncomment this before running the next cell!
+                            # Importing this local file will add the project's home directory to sys.path
 
 # %% tags=[]
-from experiments.get_notebook_info import get_notebook_basename
+#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 ipynbname
 
-from life123 import ChemData
-from life123 import UniformCompartment
-from life123.visualization.plotly_helper import PlotlyHelper
+from life123 import check_version, ChemData, UniformCompartment, PlotlyHelper, GraphicLog
 
-from life123 import GraphicLog
+# %%
+check_version(LIFE123_VERSION)
 
 # %% tags=[]
 # Initialize the HTML logging (for the graphics)
-log_file = get_notebook_basename() + ".log.htm"    # Use the notebook base filename for the log file
+log_file = ipynbname.name() + ".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,
@@ -66,15 +73,16 @@
 # Specify the chemicals and the reactions
 
 # %% tags=[]
-# Instantiate the simulator and specify the chemicals
-chem = ChemData(names=["A", "B", "C"])
+# Specify the chemicals and the reactions; this data structure will get re-used in 
+# the various simulations below
+chem = ChemData()
 
 # Reaction A <-> B (fast)
-chem.add_reaction(reactants=["A"], products=["B"],
+chem.add_reaction(reactants="A", products="B",
                        forward_rate=64., reverse_rate=8.) 
 
 # Reaction B <-> C (slow)
-chem.add_reaction(reactants=["B"], products=["C"],
+chem.add_reaction(reactants="B", products="C",
                        forward_rate=12., reverse_rate=2.) 
 
 print("Number of reactions: ", chem.number_of_reactions())
@@ -86,6 +94,8 @@
 # Send a plot of the network of reactions to the HTML log file
 chem.plot_reaction_network("vue_cytoscape_2")
 
+# %%
+
 # %% [markdown]
 # ## Run the simulation
 
@@ -102,7 +112,7 @@
 # ## Run the reaction
 
 # %%
-dynamics.enable_diagnostics()         # To save diagnostic information about the call to single_compartment_react()
+dynamics.enable_diagnostics()    # To save diagnostic information about the call to single_compartment_react()
 
 dynamics.single_compartment_react(initial_step=0.02, duration=0.4,
                                   snapshots={"initial_caption": "1st reaction step",
@@ -115,7 +125,7 @@
 
 # %%
 dynamics.plot_history(title="Coupled reactions A <-> B and B <-> C",
-                      colors=['blue', 'orange', 'green'], show_intervals=True)
+                      colors=['darkturquoise', 'orange', 'green'], show_intervals=True)
 
 # %%
 dynamics.curve_intersect("A", "B", t_start=0, t_end=0.05)
@@ -130,7 +140,7 @@
 dynamics.get_history()
 
 # %%
-dynamics.explain_time_advance()
+dynamics.diagnostics.explain_time_advance()
 
 # %% [markdown]
 # ### Check the final equilibrium
@@ -174,14 +184,13 @@
 # %%
 # Concentration increments due to reaction 0 (A <-> B)
 # Note that [C] is not affected
-dynamics.get_diagnostic_rxn_data(rxn_index=0)
+dynamics.diagnostics.get_diagnostic_rxn_data(rxn_index=0)
 
 # %%
-# Concentration increments due to reaction 1 (B <-> C)
-# Note that [A] is not affected.  
-# Also notice that the 0-th row from the A <-> B reaction isn't seen here, because that step was aborted
-# early on, before even getting to THIS reaction
-dynamics.get_diagnostic_rxn_data(rxn_index=1)
+# Concentration increments due to reaction 1 (B <-> C) 
+# Also notice that the 0-th row from the A <-> B reaction isn't seen here (start time 0 and step 0.02), 
+# because that step was aborted early on, BEFORE even getting to THIS reaction
+dynamics.diagnostics.get_diagnostic_rxn_data(rxn_index=1)
 
 # %%
 
@@ -209,10 +218,10 @@
 
 # %%
 dynamics2.plot_history(title="Coupled reactions A <-> B and B <-> C , re-run with CONSTANT STEPS",
-                       colors=['blue', 'orange', 'green'], show_intervals=True)
+                       colors=['darkturquoise', 'orange', 'green'], show_intervals=True)
 
 # %% [markdown]
-# _(Notice that the vertical steps are now equally spaced - and that there are so many of them that we're only showing some)_
+# _(Notice that the vertical steps are now equally spaced - and that there are so many of them that we're only showing 1 every 6)_
 
 # %%
 dynamics2.curve_intersect(t_start=0, t_end=0.05, chem1="A", chem2="B")

experiments/reactions_single_compartment/cascade_1.py

@@ -13,7 +13,7 @@
 # ---
 
 # %% [markdown]
-# ## A complex (multistep) reaction `A <-> C` derived from 2 coupled elementary reactions: 
+# ## A complex (composite/multistep) reaction `A <-> C` derived from 2 coupled elementary reactions: 
 # ## `A <-> B` and `B <-> C`  
 # We are given the time evolution of the complex reaction,  
 # and want to determine whether it can be modeled as an elementary reaction.  
@@ -22,17 +22,24 @@
 # In PART 2, the time functions generated in Part 1 are taken as a _starting point,_ to explore how to model the composite reaction `A <-> C`  
 #
 # **Background**: please see experiments `cascade_1` and `mystery_reaction_1`
-#
-# LAST REVISED: June 23, 2024 (using v. 1.0 beta36)
 
 # %%
-import set_path      # Importing this module will add the project's home directory to sys.path
+LAST_REVISED = "July 26, 2024"
+LIFE123_VERSION = "1.0.0.beta.38"    # Version this experiment is based on
+
+# %%
+#import set_path            # Using MyBinder?  Uncomment this before running the next cell!
+                            # Importing this module will add the project's home directory to sys.path
 
 # %% tags=[]
-from experiments.get_notebook_info import get_notebook_basename
+#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 UniformCompartment
-from life123.visualization.plotly_helper import PlotlyHelper
+from life123 import check_version, UniformCompartment, PlotlyHelper
+
+# %%
+check_version(LIFE123_VERSION)
 
 # %%
 
@@ -44,7 +51,7 @@
 
 # %% tags=[]
 # Instantiate the simulator and specify the chemicals
-dynamics = UniformCompartment(names=["A", "B", "C"], preset="mid")
+dynamics = UniformCompartment(preset="mid")
 
 # Reaction A <-> B (slower, and with a smaller K)
 dynamics.add_reaction(reactants="A", products="B",
@@ -73,7 +80,7 @@
 dynamics.plot_history(colors=['darkturquoise', 'orange', 'green'], show_intervals=True)
 
 # %%
-dynamics.is_in_equilibrium(tolerance=15)
+dynamics.is_in_equilibrium(tolerance=12)
 
 # %%
 
@@ -87,7 +94,8 @@
 # Let's start by taking stock of the actual data (saved during the simulation of part 1):
 
 # %%
-df = dynamics.get_history(columns=["SYSTEM TIME", "A", "C", "caption"])   # We're NOT given the intermediary B
+# For this analysis, we're NOT given the intermediary B
+df = dynamics.get_history(columns=["SYSTEM TIME", "A", "C", "caption"])
 df
 
 # %% [markdown]
@@ -111,11 +119,11 @@
 # Let's see what happens if we try to do such a linear fit!
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr, reactant_conc=A_conc, product_conc=C_conc, reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr, A_conc=A_conc, B_conc=C_conc, reactant_name="A", product_name="C")
 
 # %% [markdown]
 # ### The least-square fit is awful : the complex reaction `A <-> C` doesn't seem to be amenable to being modeled as a simple reaction with some suitable rate constants
-# Probably not too surprising given our "secret" knowledge from Part 1 that the complex reaction originates from 2 elementary reactions where one doesn't dominate the other one in terms of reaction kinetics
+# Probably not too surprising given our "secret" knowledge from Part 1 that the complex reaction originates from 2 elementary reactions where the intermediate product builds up at one point
 
 # %% [markdown]
 # ### A glance at the above diagram reveals much-better linear fits, if split into 2 portions, one where A(t) ranges from 0 to about 24, and one from about 24 to 50   
@@ -135,8 +143,8 @@
 
 # %% [markdown]
 # ### Let's split the `A_conc` and `C_conc` arrays we extracted earlier (with the entire time evolution of, respectively, [A] and [C]) into two parts:  
-# 1) points numbered 0-47   
-# 2) points 48-end
+# 1) points numbered 0 thru 47   
+# 2) points 48 - end
 
 # %%
 A_conc_early = A_conc[:48]
@@ -161,12 +169,12 @@
 # ### I. Let's start with the EARLY region, when t < 0.1
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_early, reactant_conc=A_conc_early, product_conc=C_conc_early, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_early, A_conc=A_conc_early, B_conc=C_conc_early,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # Trying to fit an elementary reaction to that region leads to a **negative** reverse rate constant!  
-# It's not surprise that an elementary reaction is a good fit, if one observes what happens to the time evolution of the concentrations.  Repeating the earlier plot, but only showing `A` and `C` (i.e. hiding the intermediary `B`):
+# It's no surprise that an elementary reaction is a good fit, if one observes what happens to the time evolution of the concentrations.  Repeating the earlier plot, but only showing `A` and `C` (i.e. hiding the intermediary `B`):
 
 # %%
 dynamics.plot_history(colors=['darkturquoise', 'green'], xrange=[0, 0.4], vertical_lines=[0.1], 
@@ -177,17 +185,19 @@
 # when the reactant `A` is plentiful, the rate of change (gradient) of the product `C` is low - and vice versa.  
 # Does that look like an elementary reaction in its kinetics?  Nope!
 
+# %%
+
 # %% [markdown]
 # ### II. And now let's consider the LATE region, when t > 0.1
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_late, reactant_conc=A_conc_late, product_conc=C_conc_late, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_late, A_conc=A_conc_late, B_conc=C_conc_late,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # This time we have an adequate linear fit AND meaningful rate constants : kF of about 8 and kR of about 0.  Do those numbers sound familiar?  A definite resemblance to the kF=8, kR=2 of the SLOWER elementary reaction `A <-> B`!  
 #
-# #### The slower `A <-> B` reaction dominates the kinetics from about t=0.1  
+# #### The slower `A <-> B` reaction dominates the kinetics, from about t=0.1 on  
 #
 # Let's see the graph again:
 
@@ -204,11 +214,16 @@
 # The slow link (`A <-> B`) largely determines the kinetics of the supply line.
 
 # %% [markdown]
-# ### While it's a well-known Chemistry notion that the slower reaction is the rate-determining step in a chain, we saw in this experiment  that the complex reaction could be roughly modeled with the rate constants of the slower reaction only after some time.  
+# ### While it's a well-known Chemistry notion that the slower reaction is the rate-determining step in a chain, we saw in this experiment  that **the complex reaction could be roughly modeled with the rate constants of the slower reaction ONLY AFTER SOME TIME**.  
 
 # %% [markdown]
 # If we were interested in early transients (for example, if diffusion quickly intervened), we couldn't use that model.
 
+# %% [markdown]
+# #### Is that surprising?  At early times, compare the inflection of the final product, C, of the composite reaction vs. the inflection of the product of a simple reaction (such as B in experiment `react1`, both appearing in green.)
+
+# %%
+
 # %% [markdown]
 # #### In the continuation experiment, `cascade_2_b`, we explore the scenario where the 2 elementary reactions are much more different from each other
 

experiments/reactions_single_compartment/cascade_2_a.ipynb

@@ -111,7 +111,7 @@
 # Let's see what happens if we try to do such a linear fit!
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr, reactant_conc=A_conc, product_conc=C_conc, reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr, A_conc=A_conc, B_conc=C_conc, reactant_name="A", product_name="C")
 
 # %% [markdown]
 # ### The least-square fit is awful : the complex reaction `A <-> C` doesn't seem to be amenable to being modeled as a simple reaction with some suitable rate constants
@@ -153,8 +153,8 @@
 # ### I. Let's start with the EARLY region, when t < 0.028
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_early, reactant_conc=A_conc_early, product_conc=C_conc_early, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_early, A_conc=A_conc_early, B_conc=C_conc_early,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # Just as we saw in experiment `cascade_2_a`, trying to fit an elementary reaction to that region leads to a **negative** reverse rate constant!   
@@ -164,8 +164,8 @@
 # ### II. And now let's consider the LATE region, when t > 0.028
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_late, reactant_conc=A_conc_late, product_conc=C_conc_late, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_late, A_conc=A_conc_late, B_conc=C_conc_late,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # This time we have an adequate linear fit AND meaningful rate constants : kF of about 8 and kR of about 0.  Do those numbers sound familiar?  A definite resemblance to the kF=8, kR=2 of the SLOWER elementary reaction `A <-> B`!  

experiments/reactions_single_compartment/cascade_2_a.py

@@ -111,7 +111,7 @@
 # Let's see what happens if we try to do such a linear fit!
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr, reactant_conc=A_conc, product_conc=C_conc, reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr, A_conc=A_conc, B_conc=C_conc, reactant_name="A", product_name="C")
 
 # %% [markdown]
 # ### A terrible fit!  But it looks like we'll do much better if split it into 3 portions, one where A(t) ranges from 0 to about 0.04, one from about 0.04 to 5, and one from about 5 to 50   
@@ -159,8 +159,8 @@
 # ### I. Let's start with the EARLY region, when t < 0.028
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_early, reactant_conc=A_conc_early, product_conc=C_conc_early, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_early, A_conc=A_conc_early, B_conc=C_conc_early,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # Just as we saw in experiment `cascade_2_b`, trying to fit an elementary reaction to that region leads to a **negative** reverse rate constant!  
@@ -170,8 +170,8 @@
 # ### II. Let's now consider the MID region (not seen in earlier experiments in this series), when 0.028 < t < 0.1
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_mid, reactant_conc=A_conc_mid, product_conc=C_conc_mid, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_mid, A_conc=A_conc_mid, B_conc=C_conc_mid,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # For this region, too, trying to fit an elementary reaction to that region leads to a **negative** reverse rate!  
@@ -181,8 +181,8 @@
 # ### III. And now let's consider the LATE region, when t > 0.1
 
 # %%
-dynamics.estimate_rate_constants(t=t_arr_late, reactant_conc=A_conc_late, product_conc=C_conc_late, 
-                                 reactant_name="A", product_name="C")
+dynamics.estimate_rate_constants_simple(t=t_arr_late, A_conc=A_conc_late, B_conc=C_conc_late,
+                                        reactant_name="A", product_name="C")
 
 # %% [markdown]
 # This time we have an adequate linear fit AND positive rate constants, but a **huge value of kF** relative to that of any of the elementary reactions!