{ "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: Feb. 5, 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 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": 1, "id": "d9efa3fd-e95d-4e1c-878a-81ae932b2709", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Added 'D:\\Docs\\- MY CODE\\BioSimulations\\life123-Win7' to sys.path\n" ] } ], "source": [ "import set_path # Importing this module will add the project's home directory to sys.path" ] }, { "cell_type": "code", "execution_count": 2, "id": "01bae555-3dcf-42c1-bddc-9477a37f49f8", "metadata": { "tags": [] }, "outputs": [], "source": [ "from experiments.get_notebook_info import get_notebook_basename\n", "\n", "from src.modules.reactions.reaction_data import ReactionData as chem\n", "from src.modules.reactions.reaction_dynamics import ReactionDynamics\n", "from src.modules.numerical.numerical import Numerical as num\n", "\n", "import numpy as np\n", "import plotly.express as px\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(reaction_data=chem_data)\n", "dynamics.set_conc(conc=initial_conc, snapshot=True)\n", "dynamics.describe_state()" ] }, { "cell_type": "code", "execution_count": 7, "id": "b31872ae-5571-4602-a552-10767b20c7af", "metadata": {}, "outputs": [], "source": [ "dynamics.set_diagnostics() # To save diagnostic information about the call to single_compartment_react()" ] }, { "cell_type": "code", "execution_count": 8, "id": "ed639e10-85a7-4068-868a-8f1fc530d8ce", "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": 10, "id": "b1e8bef4-fd94-485d-b141-40b2af02a97a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "From time 0 to 0.0024, in 2 FULL steps of 0.0012\n", "From time 0.0024 to 0.024, in 36 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": 12, "id": "ea96d051-2672-48a0-a9e3-d356cb979d6f", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "From time 0 to 0.0024, in 2 FULL steps of 0.0012\n", "From time 0.0024 to 0.024, in 36 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.17, in 111 substeps of 0.00125 (each 1/2 of full step)\n", "From time 0.17 to 0.257, in 35 FULL steps of 0.0025\n", "From time 0.257 to 0.259, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.259 to 0.276, in 7 FULL steps of 0.0025\n", "From time 0.276 to 0.277, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.277 to 0.29, in 5 FULL steps of 0.0025\n", "From time 0.29 to 0.291, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.291 to 0.304, in 5 FULL steps of 0.0025\n", "From time 0.304 to 0.305, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.305 to 0.322, in 7 FULL steps of 0.0025\n", "From time 0.322 to 0.324, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.324 to 0.341, in 7 FULL steps of 0.0025\n", "From time 0.341 to 0.342, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.342 to 0.355, in 5 FULL steps of 0.0025\n", "From time 0.355 to 0.356, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.356 to 0.374, in 7 FULL steps of 0.0025\n", "From time 0.374 to 0.375, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.375 to 0.392, in 7 FULL steps of 0.0025\n", "From time 0.392 to 0.394, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.394 to 0.416, in 9 FULL steps of 0.0025\n", "From time 0.416 to 0.417, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.417 to 0.435, in 7 FULL steps of 0.0025\n", "From time 0.435 to 0.436, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.436 to 0.459, in 9 FULL steps of 0.0025\n", "From time 0.459 to 0.46, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.46 to 0.482, in 9 FULL steps of 0.0025\n", "From time 0.482 to 0.484, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.484 to 0.511, in 11 FULL steps of 0.0025\n", "From time 0.511 to 0.512, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.512 to 0.545, in 13 FULL steps of 0.0025\n", "From time 0.545 to 0.546, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.546 to 0.579, in 13 FULL steps of 0.0025\n", "From time 0.579 to 0.58, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.58 to 0.622, in 17 FULL steps of 0.0025\n", "From time 0.622 to 0.624, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.624 to 0.676, in 21 FULL steps of 0.0025\n", "From time 0.676 to 0.677, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.677 to 0.78, in 41 FULL steps of 0.0025\n", "From time 0.78 to 0.781, in 1 substep of 0.00125 (1/2 of full step)\n", "From time 0.781 to 5, in 1688 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": 14, "id": "4933f334-e4f5-45e7-be91-3b11e667f244", "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.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", "| \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", "