#!/usr/bin/env python # coding: utf-8 # ## 2 COUPLED reactions of different speeds, forming a "cascade": # ### `A <-> B` (fast) and `B <-> C` (slow) # All 1st order. Taken to equilibrium. Both reactions are mostly forward. # The concentration of the intermediate product B manifests 1 oscillation (transient "overshoot") # # (Adaptive variable time resolution is used, with extensive diagnostics.) # # LAST REVISED: May 23, 2023 # ## 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 red 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.reactions.reaction_data import ReactionData as chem from src.modules.reactions.reaction_dynamics import ReactionDynamics import numpy as np import plotly.express as px 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_1"], 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]: # Specify the chemicals chem_data = chem(names=["A", "B", "C"]) # Reaction A <-> B (fast) chem_data.add_reaction(reactants=["A"], products=["B"], forward_rate=64., reverse_rate=8.) # Reaction B <-> C (slow) chem_data.add_reaction(reactants=["B"], products=["C"], forward_rate=12., reverse_rate=2.) print("Number of reactions: ", chem_data.number_of_reactions()) # In[5]: chem_data.describe_reactions() # In[6]: # Send a plot of the network of reactions to the HTML log file graph_data = chem_data.prepare_graph_network() GraphicLog.export_plot(graph_data, "vue_cytoscape_1") # # Start the simulation # In[7]: dynamics = ReactionDynamics(reaction_data=chem_data) # ### Set the initial concentrations of all the chemicals, in their index order # In[8]: dynamics.set_conc([50., 0, 0.], snapshot=True) # In[9]: dynamics.describe_state() # In[10]: dynamics.get_history() # ## Run the reaction # In[11]: dynamics.set_diagnostics() # To save diagnostic information about the call to single_compartment_react() # These settings can be tweaked to make the time resolution finer or coarser dynamics.set_thresholds(norm="norm_A", low=0.5, high=1.0, abort=1.44) dynamics.set_thresholds(norm="norm_B", low=0.2, high=0.5, abort=1.5) dynamics.set_step_factors(upshift=1.4, downshift=0.5, abort=0.5) dynamics.set_error_step_factor(0.333) 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 # In[12]: dynamics.plot_curves(title="Coupled reactions A <-> B and B <-> C", colors=['blue', 'red', 'green']) # In[13]: dynamics.curve_intersection("A", "B", t_start=0, t_end=0.05) # In[14]: dynamics.curve_intersection("A", "C", t_start=0, t_end=0.05) # In[15]: dynamics.curve_intersection("B", "C", t_start=0.05, t_end=0.1) # In[16]: dynamics.get_history() # In[17]: dynamics.explain_time_advance() # ### Check the final equilibrium # In[18]: # Verify that all the reactions have reached equilibrium dynamics.is_in_equilibrium() # In[ ]: # # EVERYTHING BELOW IS FOR DIAGNOSTIC INSIGHT # # ### Perform some verification # ### Take a peek at the diagnostic data saved during the earlier reaction simulation # In[19]: # Concentration increments due to reaction 0 (A <-> B) # Note that [C] is not affected dynamics.get_diagnostic_rxn_data(rxn_index=0) # In[20]: # 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[ ]: # ### Provide a detailed explanation of all the steps/substeps of the reactions, from the saved diagnostic data # In[ ]: # # Re-run with very small constanst steps # 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[21]: dynamics2 = ReactionDynamics(reaction_data=chem_data) # In[22]: dynamics2.set_conc([50., 0, 0.], snapshot=True) # In[23]: dynamics2.single_compartment_react(initial_step=0.0005, reaction_duration=0.4, snapshots={"initial_caption": "1st reaction step", "final_caption": "last reaction step"}, ) # In[24]: fig = px.line(data_frame=dynamics2.get_history(), x="SYSTEM TIME", y=["A", "B", "C"], title="Coupled reactions A <-> B and B <-> C (run with constant steps)", color_discrete_sequence = ['blue', 'red', 'green'], labels={"value":"concentration", "variable":"Chemical"}) fig.show() # In[25]: dynamics2.curve_intersection(t_start=0, t_end=0.05, chem1="A", chem2="B") # In[26]: dynamics2.curve_intersection(t_start=0, t_end=0.05, chem1="A", chem2="C") # In[27]: dynamics2.curve_intersection(t_start=0.05, t_end=0.1, chem1="B", chem2="C") # In[28]: df2 = dynamics2.get_history() df2 # ## Notice that we now did 801 steps - vs. the 56 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[29]: # Earlier run (using variable time steps) fig1 = px.line(data_frame=dynamics.get_history(), x="SYSTEM TIME", y=["B"], color_discrete_sequence = ['red'], labels={"value":"concentration", "variable":"Variable-step run"}) fig1.show() # In[30]: # Latest run (high-precision result from fine fixed-resolution run) fig2 = px.line(data_frame=dynamics2.get_history(), x="SYSTEM TIME", y=["B"], color_discrete_sequence = ['orange'], labels={"value":"concentration", "variable":"Fine fixed-step run"}) fig2.show() # In[31]: import plotly.graph_objects as go all_fig = go.Figure(data=fig1.data + fig2.data) # Note that the + is concatenating lists all_fig.update_layout(title="The 2 runs contrasted") all_fig['data'][0]['name']="B (variable steps)" all_fig['data'][1]['name']="B (fixed high precision)" all_fig.show() # In[ ]: