{ "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 - 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 - 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" ] }, { "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": [ "# Extend the sys.path variable, to contain the project's root directory\n", "import set_path\n", "set_path.add_ancestor_dir_to_syspath(2) # The number of levels to go up \n", " # to reach the project's home, from the folder containing this notebook" ] }, { "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 = [1] # 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", "