#!/usr/bin/env python # coding: utf-8 # ## 2 COUPLED reactions of different speeds, forming a "cascade": # ### `A <-> B` (fast) and `B <-> C` (slow) # Taken to equilibrium. Both reactions are mostly forward. All 1st order. # **The concentration of the intermediate product B manifests 1 oscillation (transient "overshoot")** # # 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. # # LAST REVISED: May 27, 2024 (using v. 1.0beta32) # ## Bathtub analogy: # A is initially full, while B and C are empty. # Tubs are progressively lower (reactions are mostly forward.) # A BIG pipe connects A and B: fast kinetics. A small pipe connects B and C: slow kinetics. # # INTUITION: B, unable to quickly drain into C while at the same time being blasted by a hefty inflow from A, # will experience a transient surge, in excess of its final equilibrium level. # # * **[Compare with the final reaction plot (the orange line is B)](#cascade_1_plot)** # ![2 Coupled Reactions](../../docs/2_coupled_reactions.png) # In[1]: import set_path # Importing this module will add the project's home directory to sys.path # In[2]: from experiments.get_notebook_info import get_notebook_basename from src.modules.chemicals.chem_data import ChemData from src.modules.reactions.reaction_dynamics import ReactionDynamics from src.modules.visualization.plotly_helper import PlotlyHelper from src.modules.visualization.graphic_log import GraphicLog # In[3]: # Initialize the HTML logging (for the graphics) 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_2"], extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js") # # Initialize the System # Specify the chemicals and the reactions # In[4]: # Instantiate the simulator and specify the chemicals chem = ChemData(names=["A", "B", "C"]) # Reaction A <-> B (fast) chem.add_reaction(reactants=["A"], products=["B"], forward_rate=64., reverse_rate=8.) # Reaction B <-> C (slow) chem.add_reaction(reactants=["B"], products=["C"], forward_rate=12., reverse_rate=2.) print("Number of reactions: ", chem.number_of_reactions()) # In[5]: chem.describe_reactions() # In[6]: # Send a plot of the network of reactions to the HTML log file chem.plot_reaction_network("vue_cytoscape_2") # ## Run the simulation # In[7]: dynamics = ReactionDynamics(chem_data=chem, preset="fast") dynamics.set_conc([50., 0, 0.], snapshot=True) # Set the initial concentrations of all the chemicals, in their index order dynamics.describe_state() # In[8]: dynamics.get_history() # ## Run the reaction # In[9]: dynamics.set_diagnostics() # To save diagnostic information about the call to single_compartment_react() dynamics.single_compartment_react(initial_step=0.02, reaction_duration=0.4, snapshots={"initial_caption": "1st reaction step", "final_caption": "last reaction step"}, variable_steps=True, explain_variable_steps=False) # ### Plots of changes of concentration with time # Notice the variable time steps (vertical dashed lines) # In[10]: dynamics.plot_history(title="Coupled reactions A <-> B and B <-> C", colors=['blue', 'orange', 'green'], show_intervals=True) # In[11]: dynamics.curve_intersect("A", "B", t_start=0, t_end=0.05) # In[12]: dynamics.curve_intersect("A", "C", t_start=0, t_end=0.05) # In[13]: dynamics.curve_intersect("B", "C", t_start=0.05, t_end=0.1) # In[14]: dynamics.get_history() # In[15]: dynamics.explain_time_advance() # ### Check the final equilibrium # In[16]: # Verify that all the reactions have reached equilibrium dynamics.is_in_equilibrium() # ### Let's look at the final concentrations of `A` and `C` (i.e., the reactant and product of the composite reaction) # In[17]: A_final = dynamics.get_chem_conc("A") A_final # In[18]: C_final = dynamics.get_chem_conc("C") C_final # #### Their ratio: # In[19]: C_final / A_final # ### As expected the equilibrium constant for the overall reaction `A <-> C` (approx. 48) is indeed the product of the equilibrium constants of the two elementary reactions (K = 8 and K = 6, respectively) that we saw earlier. # In[ ]: # In[ ]: # # PART 2 - DIAGNOSTIC INSIGHT # # Perform some verification # ### Take a peek at the diagnostic data saved during the earlier reaction simulation # In[20]: # Concentration increments due to reaction 0 (A <-> B) # Note that [C] is not affected dynamics.get_diagnostic_rxn_data(rxn_index=0) # In[21]: # 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) # In[ ]: # In[ ]: # # PART 3 : Re-run with very small constant steps, and compare with original run # We'll use **constant steps of size 0.0005** , which is 1/4 of the smallest steps (the "substep" size) previously used in the variable-step run # In[22]: dynamics2 = ReactionDynamics(chem_data=chem) # Re-use the same chemicals and reactions of the previous simulation # In[23]: dynamics2.set_conc([50., 0, 0.], snapshot=True) # In[24]: # Notice that we're using FIXED steps this time dynamics2.single_compartment_react(initial_step=0.0005, reaction_duration=0.4, variable_steps=False, snapshots={"initial_caption": "1st reaction step", "final_caption": "last reaction step"}, ) # In[25]: dynamics2.plot_history(title="Coupled reactions A <-> B and B <-> C , re-run with CONSTANT STEPS", colors=['blue', 'orange', 'green'], show_intervals=True) # _(Notice that the vertical steps are now equally spaced - and that there are so many of them that we're only showing some)_ # In[26]: dynamics2.curve_intersect(t_start=0, t_end=0.05, chem1="A", chem2="B") # In[27]: dynamics2.curve_intersect(t_start=0, t_end=0.05, chem1="A", chem2="C") # In[29]: dynamics2.curve_intersect(t_start=0.05, t_end=0.1, chem1="B", chem2="C") # In[30]: df2 = dynamics2.get_history() df2 # ## Notice that we now did 800 steps - vs. the 48 of the earlier variable-resolution run! # ## Let's compare some entries with the coarser previous variable-time run # #### Let's compare the plots of [B] from the earlier (variable-step) run, and the latest (high-precision, fixed-step) one: # In[35]: # Earlier run (using variable time steps) fig1 = dynamics.plot_history(chemicals='B', colors='orange', title="Adaptive variable-step run", show=True) # In[36]: # Latest run (high-precision result from fine fixed-resolution run) fig2 = dynamics2.plot_history(chemicals='B', colors=['violet'], title="Fine fixed-step run", show=True) # In[37]: PlotlyHelper.combine_plots(fig_list=[fig1, fig2], title="The 2 runs, contrasted together", curve_labels=["B (adaptive variable steps)", "B (fixed small steps)"]) # #### They overlap fairly well! The 800 fixed-timestep points vs. the 48 adaptable variable-timestep ones. # The adaptive algorithms avoided 752 extra steps of limited benefit... # In[ ]: