{
"cells": [
{
"cell_type": "markdown",
"id": "3bbe8002-bdf3-490c-bde0-80dd3713a3d0",
"metadata": {},
"source": [
"## A simple `A + B <-> C` reaction whose rate constants are to be estimated from a given time evolution of [A] and [B] \n",
"### (values given on a *variable-time* grid.)\n",
"\n",
"Assume the reaction is known to be 1st order (won't verify that.) \n",
"\n",
"This is the counterpart of experiment `mystery_reaction_1` for a more complex reaction; we'll also proceed faster. \n",
"\n",
"In PART 1, a time evolution of [A] and [B] is obtained by simulation \n",
"\n",
"In PART 2, the time evolutions generated in Part 1 are taken as a _starting point,_ to estimate the rate constants of `A <-> B` \n",
"\n",
"In PART 3, we'll repeat what we did in Part 2, but this time showing the full details of how the answer is arrived at"
]
},
{
"cell_type": "markdown",
"id": "4118b513-1c98-455d-b58f-37898bb03bdd",
"metadata": {},
"source": [
"### TAGS : \"numerical\", \"uniform compartment\", \"under-the-hood\""
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e4c76b22-243f-4c4b-b6ca-7d9bf89ec74c",
"metadata": {},
"outputs": [],
"source": [
"LAST_REVISED = \"Sep. 8, 2024\"\n",
"LIFE123_VERSION = \"1.0.0.beta.38\" # Version this experiment is based on"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "97b7f9e2-99b8-42a1-bede-b9e551b9e024",
"metadata": {},
"outputs": [],
"source": [
"#import set_path # Using MyBinder? Uncomment this before running the next cell!\n",
" # Importing this local file will add the project's home directory to sys.path"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "3924c013",
"metadata": {
"lines_to_next_cell": 2,
"tags": []
},
"outputs": [],
"source": [
"#import sys\n",
"#sys.path.append(\"C:/some_path/my_env_or_install\") # CHANGE to the folder containing your venv or libraries installation!\n",
"# NOTE: If any of the imports below can't find a module, uncomment the lines above, or try: import set_path\n",
"\n",
"import ipynbname\n",
"import numpy as np\n",
"\n",
"from life123 import check_version, UniformCompartment, PlotlyHelper, Numerical"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8966a716-74bd-4c5d-8f4a-5733b3806973",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OK\n"
]
}
],
"source": [
"check_version(LIFE123_VERSION)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b5208e04-e3b8-4769-8937-d61751e55a47",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"id": "9329208b-070f-4902-8f37-0f11ddf75ed6",
"metadata": {},
"source": [
"# PART 1 - We'll generate the time evolution of [A] and [B] by simulating a reaction of KNOWN rate constants...\n",
"## but pretend you don't see this section! (because we later want to estimate those rate constants)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "72b4245c-de4e-480d-a501-3495b7ed8bc4",
"metadata": {
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of reactions: 1 (at temp. 25 C)\n",
"0: A + B <-> C (kF = 5 / kR = 1 / delta_G = -3,989.7 / K = 5) | 1st order in all reactants & products\n",
"Set of chemicals involved in the above reactions: {'C', 'B', 'A'}\n"
]
}
],
"source": [
"# Instantiate the simulator and specify the accuracy preset\n",
"dynamics = UniformCompartment(preset=\"mid\")\n",
"\n",
"# Reaction A <-> B (mostly in the forward direction)\n",
"dynamics.add_reaction(reactants=[\"A\", \"B\"], products=\"C\",\n",
" forward_rate=5., reverse_rate=1.) \n",
" \n",
"dynamics.describe_reactions()"
]
},
{
"cell_type": "markdown",
"id": "98a9fbe5-2090-4d38-9c5f-94fbf7c3eae2",
"metadata": {},
"source": [
"### Run the simulation"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ae304704-c8d9-4cef-9e0b-2587bb3909ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SYSTEM STATE at Time t = 0:\n",
"3 species:\n",
" Species 0 (A). Conc: 40.0\n",
" Species 1 (B). Conc: 20.0\n",
" Species 2 (C). Conc: 5.0\n",
"Set of chemicals involved in reactions: {'C', 'B', 'A'}\n"
]
}
],
"source": [
"dynamics.set_conc({\"A\": 40., \"B\": 20., \"C\": 5.}, snapshot=True) # Set the initial concentrations\n",
"dynamics.describe_state()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "2502cd11-0df9-4303-8895-98401a1df7b8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Some steps were backtracked and re-done, to prevent negative concentrations or excessively large concentration changes\n",
"53 total step(s) taken\n",
"Number of step re-do's because of elective soft aborts: 1\n",
"Norm usage: {'norm_A': 19, 'norm_B': 17, 'norm_C': 17, 'norm_D': 17}\n"
]
}
],
"source": [
"dynamics.enable_diagnostics() # To save diagnostic information for the simulation run, below\n",
"\n",
"dynamics.single_compartment_react(initial_step=0.001, duration=0.05,\n",
" snapshots={\"initial_caption\": \"1st reaction step\",\n",
" \"final_caption\": \"last reaction step\"},\n",
" variable_steps=True)"
]
},
{
"cell_type": "markdown",
"id": "199b6238-f6ed-44a8-8130-a003444d7658",
"metadata": {},
"source": [
"### Plots of changes of concentration with time"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "a2c0e793-5457-46a5-9150-388c9f562cf0",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "Chemical=A
SYSTEM TIME=%{x}
Concentration=%{y}
SYSTEM TIME=%{x}
Concentration=%{y}
SYSTEM TIME=%{x}
Concentration=%{y}
| \n", " | SYSTEM TIME | \n", "A | \n", "B | \n", "C | \n", "caption | \n", "
|---|---|---|---|---|---|
| 0 | \n", "0.000000 | \n", "40.000000 | \n", "20.000000 | \n", "5.000000 | \n", "Initialized state | \n", "
| 1 | \n", "0.000400 | \n", "38.402000 | \n", "18.402000 | \n", "6.598000 | \n", "1st reaction step | \n", "
| 2 | \n", "0.000600 | \n", "37.696646 | \n", "17.696646 | \n", "7.303354 | \n", "\n", " |
| 3 | \n", "0.000800 | \n", "37.031002 | \n", "17.031002 | \n", "7.968998 | \n", "\n", " |
| 4 | \n", "0.001000 | \n", "36.401921 | \n", "16.401921 | \n", "8.598079 | \n", "\n", " |
| 5 | \n", "0.001240 | \n", "35.687511 | \n", "15.687511 | \n", "9.312489 | \n", "\n", " |
| 6 | \n", "0.001480 | \n", "35.017928 | \n", "15.017928 | \n", "9.982072 | \n", "\n", " |
| 7 | \n", "0.001768 | \n", "34.263512 | \n", "14.263512 | \n", "10.736488 | \n", "\n", " |
| 8 | \n", "0.002114 | \n", "33.422717 | \n", "13.422717 | \n", "11.577283 | \n", "\n", " |
| 9 | \n", "0.002528 | \n", "32.497253 | \n", "12.497253 | \n", "12.502747 | \n", "\n", " |
| 10 | \n", "0.003026 | \n", "31.492903 | \n", "11.492903 | \n", "13.507097 | \n", "\n", " |
| 11 | \n", "0.003524 | \n", "30.598990 | \n", "10.598990 | \n", "14.401010 | \n", "\n", " |
| 12 | \n", "0.004121 | \n", "29.639181 | \n", "9.639181 | \n", "15.360819 | \n", "\n", " |
| 13 | \n", "0.004718 | \n", "28.795266 | \n", "8.795266 | \n", "16.204734 | \n", "\n", " |
| 14 | \n", "0.005315 | \n", "28.048707 | \n", "8.048707 | \n", "16.951293 | \n", "\n", " |
| 15 | \n", "0.005912 | \n", "27.384727 | \n", "7.384727 | \n", "17.615273 | \n", "\n", " |
| 16 | \n", "0.006510 | \n", "26.791395 | \n", "6.791395 | \n", "18.208605 | \n", "\n", " |
| 17 | \n", "0.007107 | \n", "26.258967 | \n", "6.258967 | \n", "18.741033 | \n", "\n", " |
| 18 | \n", "0.007823 | \n", "25.683487 | \n", "5.683487 | \n", "19.316513 | \n", "\n", " |
| 19 | \n", "0.008540 | \n", "25.174287 | \n", "5.174287 | \n", "19.825713 | \n", "\n", " |
| 20 | \n", "0.009257 | \n", "24.721753 | \n", "4.721753 | \n", "20.278247 | \n", "\n", " |
| 21 | \n", "0.009973 | \n", "24.318020 | \n", "4.318020 | \n", "20.681980 | \n", "\n", " |
| 22 | \n", "0.010690 | \n", "23.956587 | \n", "3.956587 | \n", "21.043413 | \n", "\n", " |
| 23 | \n", "0.011407 | \n", "23.632031 | \n", "3.632031 | \n", "21.367969 | \n", "\n", " |
| 24 | \n", "0.012123 | \n", "23.339792 | \n", "3.339792 | \n", "21.660208 | \n", "\n", " |
| 25 | \n", "0.012840 | \n", "23.076005 | \n", "3.076005 | \n", "21.923995 | \n", "\n", " |
| 26 | \n", "0.013700 | \n", "22.789650 | \n", "2.789650 | \n", "22.210350 | \n", "\n", " |
| 27 | \n", "0.014560 | \n", "22.535388 | \n", "2.535388 | \n", "22.464612 | \n", "\n", " |
| 28 | \n", "0.015420 | \n", "22.309033 | \n", "2.309033 | \n", "22.690967 | \n", "\n", " |
| 29 | \n", "0.016280 | \n", "22.107053 | \n", "2.107053 | \n", "22.892947 | \n", "\n", " |
| 30 | \n", "0.017140 | \n", "21.926451 | \n", "1.926451 | \n", "23.073549 | \n", "\n", " |
| 31 | \n", "0.018000 | \n", "21.764669 | \n", "1.764669 | \n", "23.235331 | \n", "\n", " |
| 32 | \n", "0.018860 | \n", "21.619505 | \n", "1.619505 | \n", "23.380495 | \n", "\n", " |
| 33 | \n", "0.019720 | \n", "21.489062 | \n", "1.489062 | \n", "23.510938 | \n", "\n", " |
| 34 | \n", "0.020580 | \n", "21.371693 | \n", "1.371693 | \n", "23.628307 | \n", "\n", " |
| 35 | \n", "0.021612 | \n", "21.244815 | \n", "1.244815 | \n", "23.755185 | \n", "\n", " |
| 36 | \n", "0.022644 | \n", "21.132875 | \n", "1.132875 | \n", "23.867125 | \n", "\n", " |
| 37 | \n", "0.023675 | \n", "21.033975 | \n", "1.033975 | \n", "23.966025 | \n", "\n", " |
| 38 | \n", "0.024707 | \n", "20.946489 | \n", "0.946489 | \n", "24.053511 | \n", "\n", " |
| 39 | \n", "0.025739 | \n", "20.869015 | \n", "0.869015 | \n", "24.130985 | \n", "\n", " |
| 40 | \n", "0.026771 | \n", "20.800342 | \n", "0.800342 | \n", "24.199658 | \n", "\n", " |
| 41 | \n", "0.028010 | \n", "20.727233 | \n", "0.727233 | \n", "24.272767 | \n", "\n", " |
| 42 | \n", "0.029248 | \n", "20.663960 | \n", "0.663960 | \n", "24.336040 | \n", "\n", " |
| 43 | \n", "0.030486 | \n", "20.609146 | \n", "0.609146 | \n", "24.390854 | \n", "\n", " |
| 44 | \n", "0.031725 | \n", "20.561619 | \n", "0.561619 | \n", "24.438381 | \n", "\n", " |
| 45 | \n", "0.033211 | \n", "20.512134 | \n", "0.512134 | \n", "24.487866 | \n", "\n", " |
| 46 | \n", "0.034697 | \n", "20.470471 | \n", "0.470471 | \n", "24.529529 | \n", "\n", " |
| 47 | \n", "0.036183 | \n", "20.435365 | \n", "0.435365 | \n", "24.564635 | \n", "\n", " |
| 48 | \n", "0.037966 | \n", "20.399844 | \n", "0.399844 | \n", "24.600156 | \n", "\n", " |
| 49 | \n", "0.039749 | \n", "20.370985 | \n", "0.370985 | \n", "24.629015 | \n", "\n", " |
| 50 | \n", "0.041889 | \n", "20.342829 | \n", "0.342829 | \n", "24.657171 | \n", "\n", " |
| 51 | \n", "0.044457 | \n", "20.316603 | \n", "0.316603 | \n", "24.683397 | \n", "\n", " |
| 52 | \n", "0.047538 | \n", "20.293560 | \n", "0.293560 | \n", "24.706440 | \n", "\n", " |
| 53 | \n", "0.051236 | \n", "20.274774 | \n", "0.274774 | \n", "24.725226 | \n", "last reaction step | \n", "