{ "cells": [ { "cell_type": "markdown", "id": "49bcb5b0-f19d-4b96-a5f1-e0ae30f66d8f", "metadata": {}, "source": [ "## Excessively large time steps - and remediation\n", "### Based on experiment `cycles_1`\n", "#### A cycle of reactions `A <-> B <-> C <-> A` and `C + E_High <-> A + E_Low`\n", "\n", "### [RUN 1](#large_time_steps_1_run1) - Initially, very large time steps are taken - at the edge of pushing some concentrations into negative values.\n", "The system automatic detects those problems, intercepts the problematic steps and re-runs them with 1/2 the time step. \n", "Negative concentrations are automatically avoided, but nonetheless the plots are ragged... and the solutions are eventually unstable.\n", "\n", "### [RUN 2](#large_time_steps_1_run2) - Same primary steps as for run #1, but with the option of using 1/2 substeps as needed, \n", "with thresholds that lead to those substeps being actually utilized a fair part of the time. \n", "The raggedness and instabilities are now eliminated.\n", "(Note: the 1/2 substeps are on a per-reaction basis)\n", "\n", "LAST REVISED: Dec. 3, 2023 **THIS IS AN ARCHIVED EXPERIMENT**" ] }, { "cell_type": "markdown", "id": "ee568e86-88ba-4942-ae40-e09dc9d7e8d6", "metadata": {}, "source": [ "# IMPORTANT: DO NOT ATTEMPT TO RUN THIS NOTEBOOK! \n", "## This is a **\"frozen run\"** that depends on an old version of Life123, for demonstration purposes \n", "#### (newer versions of Life123 tend to recover from instability more gracefully!!) \n", "If you bypass the execution exit in the first cell, and run the other cells, you WILL NOT REPLICATE the results below!" ] }, { "cell_type": "code", "execution_count": 1, "id": "d00543ae-6a16-4e67-bfce-284f0234706e", "metadata": {}, "outputs": [ { "ename": "StopExecution", "evalue": "", "output_type": "error", "traceback": [] } ], "source": [ "# To stop the current and subsequent cells: USED TO PREVENT ACCIDENTAL RUNS OF THIS NOTEBOOK!\n", "\n", "class StopExecution(Exception):\n", " def _render_traceback_(self):\n", " return []\n", "\n", "raise StopExecution # See: https://stackoverflow.com/a/56953105/5478830" ] }, { "cell_type": "code", "execution_count": null, "id": "292706b9-743f-4f2a-a3af-4427d68f82ed", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "d3ef9936-020b-4ab3-b762-8a6cffb963b6", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "id": "53eebce8", "metadata": { "tags": [] }, "outputs": [], "source": [ "from experiments.get_notebook_info import get_notebook_basename\n", "\n", "from src.modules.chemicals.chem_data import ChemData as chem\n", "from src.modules.reactions.reaction_dynamics import ReactionDynamics\n", "\n", "from src.modules.visualization.graphic_log import GraphicLog" ] }, { "cell_type": "code", "execution_count": 3, "id": "cc53849f-351d-49e0-bfa8-22f8d8e22f8e", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-> Output will be LOGGED into the file 'large_time_steps_1.log.htm'\n" ] } ], "source": [ "# Initialize the HTML logging\n", "log_file = get_notebook_basename() + \".log.htm\" # Use the notebook base filename for the log file\n", "\n", "# Set up the use of some specified graphic (Vue) components\n", "GraphicLog.config(filename=log_file,\n", " components=[\"vue_cytoscape_1\"],\n", " extra_js=\"https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js\")" ] }, { "cell_type": "markdown", "id": "d6d3ca49-589d-49b7-8424-37c7b01bcacf", "metadata": {}, "source": [ "### Initialize the system" ] }, { "cell_type": "code", "execution_count": 4, "id": "23c15e66-52e4-495b-aa3d-ecddd8d16942", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Number of reactions: 3 (at temp. 25 C)\n", "0: A <-> B (kF = 9 / kR = 3 / Delta_G = -2,723.41 / K = 3) | 1st order in all reactants & products\n", "1: B <-> C (kF = 8 / kR = 4 / Delta_G = -1,718.28 / K = 2) | 1st order in all reactants & products\n", "2: C + E_high <-> A + E_low (kF = 1 / kR = 0.2 / Delta_G = -3,989.73 / K = 5) | 1st order in all reactants & products\n", "[GRAPHIC ELEMENT SENT TO LOG FILE `large_time_steps_1.log.htm`]\n" ] } ], "source": [ "# Initialize the system\n", "chem_data = chem(names=[\"A\", \"B\", \"C\", \"E_high\", \"E_low\"])\n", "\n", "# Reaction A <-> B, mostly in forward direction (favored energetically)\n", "# Note: all reactions in this experiment have 1st-order kinetics for all species\n", "chem_data.add_reaction(reactants=\"A\", products=\"B\",\n", " forward_rate=9., reverse_rate=3.)\n", "\n", "# Reaction B <-> C, also favored energetically\n", "chem_data.add_reaction(reactants=\"B\", products=\"C\",\n", " forward_rate=8., reverse_rate=4.)\n", "\n", "# Reaction C + E_High <-> A + E_Low, also favored energetically, but kinetically slow\n", "# Note that, thanks to the energy donation from E, we can go \"upstream\" from C, to the higher-energy level of \"A\"\n", "chem_data.add_reaction(reactants=[\"C\" , \"E_high\"], products=[\"A\", \"E_low\"],\n", " forward_rate=1., reverse_rate=0.2)\n", "\n", "chem_data.describe_reactions()\n", "\n", "# Send the plot of the reaction network to the HTML log file\n", "graph_data = chem_data.prepare_graph_network()\n", "GraphicLog.export_plot(graph_data, \"vue_cytoscape_1\")" ] }, { "cell_type": "markdown", "id": "d1d0eabb-b5b1-4e15-846d-5e483a5a24a7", "metadata": {}, "source": [ "### Set the initial concentrations of all the chemicals" ] }, { "cell_type": "code", "execution_count": 5, "id": "e4ff6a84-f5d5-4645-9c56-d9e981c108df", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'A': 100.0, 'B': 0.0, 'C': 0.0, 'E_high': 1000.0, 'E_low': 0.0}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "initial_conc = {\"A\": 100., \"B\": 0., \"C\": 0., \"E_high\": 1000., \"E_low\": 0.} # Note the abundant energy source \"E_high\"\n", "initial_conc" ] }, { "cell_type": "markdown", "id": "0b46b395-3f68-4dbd-b0c5-d67a0e623726", "metadata": { "tags": [] }, "source": [ "### We'll split each simulation in three segments (apparent from the graphs of the runs, further down): \n", "Time [0-0.03] fast changes \n", "Time [0.03-5.] medium changes \n", "Time [5.-8.] slow changes, as we approach equilibrium " ] }, { "cell_type": "markdown", "id": "5d0e1750-2a2d-4922-96cc-abc076667ac5", "metadata": {}, "source": [ "## RUN 1 - Initially, very large time steps are taken - at the edge of pushing some concentrations into negative values.\n", "The system automatic detects those problems, intercepts the problematic steps and re-runs them with 1/2 the time step. \n", "Negative concentrations are automatically avoided, but nonetheless the plots are ragged... and the solutions are eventually unstable." ] }, { "cell_type": "code", "execution_count": 6, "id": "e80645d6-eb5b-4c78-8b46-ae126d2cb2cf", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SYSTEM STATE at Time t = 0:\n", "5 species:\n", " Species 0 (A). Conc: 100.0\n", " Species 1 (B). Conc: 0.0\n", " Species 2 (C). Conc: 0.0\n", " Species 3 (E_high). Conc: 1000.0\n", " Species 4 (E_low). Conc: 0.0\n" ] } ], "source": [ "dynamics = ReactionDynamics(chem_data=chem_data)\n", "dynamics.set_conc(conc=initial_conc, snapshot=True)\n", "dynamics.describe_state()" ] }, { "cell_type": "code", "execution_count": 19, "id": "18acc0e2-b813-4a30-98a9-82aaed16cb19", "metadata": {}, "outputs": [], "source": [ "dynamics.set_diagnostics() # To save diagnostic information about the call to single_compartment_react()" ] }, { "cell_type": "code", "execution_count": 20, "id": "b7f58e02-569b-45c6-90b7-1c317519bbd2", "metadata": {}, "outputs": [], "source": [ "#dynamics.verbose_list = [\"substeps\"] # Uncomment for debugging information" ] }, { "cell_type": "code", "execution_count": 9, "id": "50ddd8e3-58c6-41f8-b874-ddc9a1d64d30", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The computation took 36 extra step(s) - automatically added to prevent negative concentrations\n", "44 total step(s) taken\n" ] } ], "source": [ "dynamics.single_compartment_react(time_step=0.0012, stop_time=0.03)\n", "#dynamics.get_history()" ] }, { "cell_type": "code", "execution_count": 22, "id": "39db6bc4-f7a4-49ff-ad02-3dc0262bc163", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "From time 0 to 0.0012, in 2 substeps of 0.0006 (each 1/2 of full step)\n", "From time 0.0012 to 0.0162, in 50 substeps of 0.0003 (each 1/4 of full step)\n", "From time 0.0162 to 0.024, in 13 substeps of 0.0006 (each 1/2 of full step)\n", "From time 0.024 to 0.0312, in 6 FULL steps of 0.0012\n" ] } ], "source": [ "dynamics.explain_time_advance()" ] }, { "cell_type": "markdown", "id": "229982f1-485f-4f84-8619-dd4c6ab343ef", "metadata": {}, "source": [ "### Notice how the system automatically intervened and reduced some steps in 1/2, to prevent negative concentrations" ] }, { "cell_type": "code", "execution_count": 11, "id": "fbf24a85-b152-4b24-b350-00bed12e7bb9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The computation took 130 extra step(s) - automatically added to prevent negative concentrations\n", "2053 total step(s) taken\n" ] } ], "source": [ "dynamics.single_compartment_react(time_step=0.0025, stop_time=5.)\n", "#dynamics.get_history()" ] }, { "cell_type": "code", "execution_count": 24, "id": "4c76da69-a77a-410d-bcd0-e255aa574a8f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "From time 0 to 0.0012, in 2 substeps of 0.0006 (each 1/2 of full step)\n", "From time 0.0012 to 0.0162, in 50 substeps of 0.0003 (each 1/4 of full step)\n", "From time 0.0162 to 0.024, in 13 substeps of 0.0006 (each 1/2 of full step)\n", "From time 0.024 to 0.0312, in 6 FULL steps of 0.0012\n", "From time 0.0312 to 0.0337, in 2 substeps of 0.00125 (each 1/2 of full step)\n", "From time 0.0337 to 0.17, in 218 substeps of 0.000625 (each 1/4 of full step)\n", "From time 0.17 to 0.172, in 1 FULL step of 0.0025\n", "From time 0.172 to 2.28, in 1686 substeps of 0.00125 (each 1/2 of full step)\n", "From time 2.28 to 5, in 1089 FULL steps of 0.0025\n" ] } ], "source": [ "dynamics.explain_time_advance()" ] }, { "cell_type": "markdown", "id": "c3fd34b4-3d87-4ae6-b69e-f8ddc2020a7e", "metadata": {}, "source": [ "### Notice how the system automatically periodically intervened and reduced some steps in 1/2, to prevent negative concentrations.\n", "The time step that was originally requested hovers near the cusp of being so large as to lead to negative concentrations" ] }, { "cell_type": "code", "execution_count": 13, "id": "938b9678-900a-4d6b-8cac-68d62545b09e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The computation took 49 extra step(s) - automatically added to prevent negative concentrations\n", "400 total step(s) taken\n" ] } ], "source": [ "dynamics.single_compartment_react(time_step=0.008, stop_time=8.)" ] }, { "cell_type": "code", "execution_count": 26, "id": "5f2ad190-f858-4227-a370-04e8b5d854be", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| \n", " | SYSTEM TIME | \n", "A | \n", "B | \n", "C | \n", "E_high | \n", "E_low | \n", "caption | \n", "
|---|---|---|---|---|---|---|---|
| 0 | \n", "0.00000 | \n", "100.000000 | \n", "0.000000 | \n", "0.000000 | \n", "1000.000000 | \n", "0.000000 | \n", "Initial state | \n", "
| 1 | \n", "0.00060 | \n", "99.460000 | \n", "0.540000 | \n", "0.000000 | \n", "1000.000000 | \n", "0.000000 | \n", "Interm. step, due to the fast rxns: [0, 1, 2] | \n", "
| 2 | \n", "0.00120 | \n", "98.923888 | \n", "1.073520 | \n", "0.002592 | \n", "1000.000000 | \n", "0.000000 | \n", "\n", " |
| 3 | \n", "0.00150 | \n", "98.658537 | \n", "1.337075 | \n", "0.004388 | \n", "999.999222 | \n", "0.000778 | \n", "Interm. step, due to the fast rxns: [0] | \n", "
| 4 | \n", "0.00180 | \n", "98.394140 | \n", "1.599676 | \n", "0.006183 | \n", "999.998445 | \n", "0.001555 | \n", "\n", " |
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 3623 | \n", "7.97045 | \n", "10.018391 | \n", "30.006875 | \n", "59.974735 | \n", "32.330220 | \n", "967.669780 | \n", "\n", " |
| 3624 | \n", "7.97845 | \n", "10.018012 | \n", "30.006785 | \n", "59.975202 | \n", "32.329439 | \n", "967.670561 | \n", "\n", " |
| 3625 | \n", "7.98645 | \n", "10.017978 | \n", "30.006691 | \n", "59.975330 | \n", "32.328339 | \n", "967.671661 | \n", "\n", " |
| 3626 | \n", "7.99445 | \n", "10.017485 | \n", "30.006608 | \n", "59.975908 | \n", "32.327699 | \n", "967.672301 | \n", "\n", " |
| 3627 | \n", "8.00245 | \n", "10.017621 | \n", "30.006514 | \n", "59.975865 | \n", "32.326462 | \n", "967.673538 | \n", "\n", " |
3628 rows × 7 columns
\n", "| \n", " | SYSTEM TIME | \n", "A | \n", "B | \n", "C | \n", "E_high | \n", "E_low | \n", "caption | \n", "
|---|---|---|---|---|---|---|---|
| 0 | \n", "0.0000 | \n", "100.000000 | \n", "0.000000 | \n", "0.000000 | \n", "1000.000000 | \n", "0.000000 | \n", "Initial state | \n", "
| 1 | \n", "0.0012 | \n", "98.920000 | \n", "1.080000 | \n", "0.000000 | \n", "1000.000000 | \n", "0.000000 | \n", "\n", " |
| 2 | \n", "0.0024 | \n", "97.855552 | \n", "2.134080 | \n", "0.010368 | \n", "1000.000000 | \n", "0.000000 | \n", "\n", " |
| 3 | \n", "0.0030 | \n", "97.337194 | \n", "2.648440 | \n", "0.014366 | \n", "999.993779 | \n", "0.006221 | \n", "\n", " |
| 4 | \n", "0.0036 | \n", "96.824887 | \n", "3.156616 | \n", "0.018497 | \n", "999.985232 | \n", "0.014768 | \n", "\n", " |
| ... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
| 2493 | \n", "7.9732 | \n", "8.663600 | \n", "30.034016 | \n", "61.302385 | \n", "33.266197 | \n", "966.733803 | \n", "\n", " |
| 2494 | \n", "7.9812 | \n", "11.674383 | \n", "29.976478 | \n", "58.349140 | \n", "30.352451 | \n", "969.647549 | \n", "\n", " |
| 2495 | \n", "7.9892 | \n", "7.609519 | \n", "30.046276 | \n", "62.344205 | \n", "34.296194 | \n", "965.703806 | \n", "\n", " |
| 2496 | \n", "7.9972 | \n", "13.130430 | \n", "29.945103 | \n", "56.924467 | \n", "28.948509 | \n", "971.051491 | \n", "\n", " |
| 2497 | \n", "8.0052 | \n", "5.686231 | \n", "30.076908 | \n", "64.236860 | \n", "36.165999 | \n", "963.834001 | \n", "\n", " |
2498 rows × 7 columns
\n", "