{
"cells": [
{
"cell_type": "markdown",
"id": "3bbe8002-bdf3-490c-bde0-80dd3713a3d0",
"metadata": {},
"source": [
"## A simple `A <-> B` 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",
"In PART 1, a time evolution of [A] and [B], with known rate constants, 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": "900caaa3-7883-466c-a2bd-5f61de65b151",
"metadata": {},
"outputs": [],
"source": [
"LAST_REVISED = \"Nov. 12, 2024\"\n",
"LIFE123_VERSION = \"1.0.0.rc.0\" # Library 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!"
]
},
{
"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, ReactionKinetics, PlotlyHelper, Numerical"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "1276f744-1b85-46e1-8d9d-f6fa90e73b49",
"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": "code",
"execution_count": null,
"id": "1190f092-c2ff-46f0-bce9-86e38ef13248",
"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 (kF = 12 / kR = 2 / delta_G = -4,441.7 / K = 6) | 1st order in all reactants & products\n",
"Set of chemicals involved in the above reactions: {'B', 'A'}\n"
]
}
],
"source": [
"# Instantiate the simulator and specify the accuracy preset\n",
"uc = UniformCompartment(preset=\"mid\", enable_diagnostics=True)\n",
"\n",
"# Reaction A <-> B (mostly in the forward direction)\n",
"uc.add_reaction(reactants=\"A\", products=\"B\",\n",
" forward_rate=12., reverse_rate=2.) \n",
" \n",
"uc.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",
"2 species:\n",
" Species 0 (A). Conc: 40.0\n",
" Species 1 (B). Conc: 10.0\n",
"Set of chemicals involved in reactions: {'B', 'A'}\n"
]
}
],
"source": [
"uc.set_conc({\"A\": 40., \"B\": 10.}) # Set the initial concentrations\n",
"uc.describe_state()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "2502cd11-0df9-4303-8895-98401a1df7b8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"39 total step(s) taken in 0.260 sec\n",
"Number of step re-do's because of elective soft aborts: 2\n",
"Norm usage: {'norm_A': 24, 'norm_B': 22, 'norm_C': 22, 'norm_D': 22}\n",
"System Time is now: 0.51497\n"
]
}
],
"source": [
"uc.single_compartment_react(initial_step=0.01, duration=0.5,\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}
| \n", " | SYSTEM TIME | \n", "A | \n", "B | \n", "caption | \n", "
|---|---|---|---|---|
| 0 | \n", "0.000000 | \n", "40.000000 | \n", "10.000000 | \n", "Set concentration | \n", "
| 1 | \n", "0.001600 | \n", "39.264000 | \n", "10.736000 | \n", "1st reaction step | \n", "
| 2 | \n", "0.003520 | \n", "38.400584 | \n", "11.599416 | \n", "\n", " |
| 3 | \n", "0.005440 | \n", "37.560376 | \n", "12.439624 | \n", "\n", " |
| 4 | \n", "0.007744 | \n", "36.579229 | \n", "13.420771 | \n", "\n", " |
| 5 | \n", "0.010509 | \n", "35.439829 | \n", "14.560171 | \n", "\n", " |
| 6 | \n", "0.013274 | \n", "34.344532 | \n", "15.655468 | \n", "\n", " |
| 7 | \n", "0.016038 | \n", "33.291632 | \n", "16.708368 | \n", "\n", " |
| 8 | \n", "0.018803 | \n", "32.279486 | \n", "17.720514 | \n", "\n", " |
| 9 | \n", "0.021568 | \n", "31.306517 | \n", "18.693483 | \n", "\n", " |
| 10 | \n", "0.024886 | \n", "30.184148 | \n", "19.815852 | \n", "\n", " |
| 11 | \n", "0.028204 | \n", "29.113912 | \n", "20.886088 | \n", "\n", " |
| 12 | \n", "0.031521 | \n", "28.093386 | \n", "21.906614 | \n", "\n", " |
| 13 | \n", "0.034839 | \n", "27.120262 | \n", "22.879738 | \n", "\n", " |
| 14 | \n", "0.038820 | \n", "26.006754 | \n", "23.993246 | \n", "\n", " |
| 15 | \n", "0.042802 | \n", "24.955312 | \n", "25.044688 | \n", "\n", " |
| 16 | \n", "0.046783 | \n", "23.962474 | \n", "26.037526 | \n", "\n", " |
| 17 | \n", "0.051561 | \n", "22.837477 | \n", "27.162523 | \n", "\n", " |
| 18 | \n", "0.056338 | \n", "21.787726 | \n", "28.212274 | \n", "\n", " |
| 19 | \n", "0.061116 | \n", "20.808189 | \n", "29.191811 | \n", "\n", " |
| 20 | \n", "0.066849 | \n", "19.711365 | \n", "30.288635 | \n", "\n", " |
| 21 | \n", "0.072582 | \n", "18.702575 | \n", "31.297425 | \n", "\n", " |
| 22 | \n", "0.078315 | \n", "17.774755 | \n", "32.225245 | \n", "\n", " |
| 23 | \n", "0.085195 | \n", "16.750734 | \n", "33.249266 | \n", "\n", " |
| 24 | \n", "0.092074 | \n", "15.825343 | \n", "34.174657 | \n", "\n", " |
| 25 | \n", "0.100330 | \n", "14.821829 | \n", "35.178171 | \n", "\n", " |
| 26 | \n", "0.108586 | \n", "13.934301 | \n", "36.065699 | \n", "\n", " |
| 27 | \n", "0.118492 | \n", "12.992362 | \n", "37.007638 | \n", "\n", " |
| 28 | \n", "0.130381 | \n", "12.018806 | \n", "37.981194 | \n", "\n", " |
| 29 | \n", "0.144646 | \n", "11.044979 | \n", "38.955021 | \n", "\n", " |
| 30 | \n", "0.158912 | \n", "10.265644 | \n", "39.734356 | \n", "\n", " |
| 31 | \n", "0.176031 | \n", "9.517222 | \n", "40.482778 | \n", "\n", " |
| 32 | \n", "0.196574 | \n", "8.834360 | \n", "41.165640 | \n", "\n", " |
| 33 | \n", "0.221225 | \n", "8.250593 | \n", "41.749407 | \n", "\n", " |
| 34 | \n", "0.250806 | \n", "7.791835 | \n", "42.208165 | \n", "\n", " |
| 35 | \n", "0.286304 | \n", "7.469313 | \n", "42.530687 | \n", "\n", " |
| 36 | \n", "0.328902 | \n", "7.274627 | \n", "42.725373 | \n", "\n", " |
| 37 | \n", "0.380018 | \n", "7.180328 | \n", "42.819672 | \n", "\n", " |
| 38 | \n", "0.441359 | \n", "7.148149 | \n", "42.851851 | \n", "\n", " |
| 39 | \n", "0.514967 | \n", "7.142696 | \n", "42.857304 | \n", "last reaction step | \n", "