Skip to content

Add Schlogl model benchmark notebook #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jorblancoa wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add Schlogl model benchmark notebook #14

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
332 changes: 332 additions & 0 deletions notebooks/Schlogl_Model.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,332 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Benchmark of well-mixed stochastic solvers for a Schlogl model\n",
"\n",
"STEPS simulation of the Schlogl model described in appendix A.4 of https://link.springer.com/978-3-319-63113-4, comparing the performance of the 'direct' and 'rssa' stochastic solvers for this model.\n",
"\n",
"## Model\n",
"\n",
"The model consists of three species, $\\{A,B,X\\}$, and four reactions implementing a trimolecular, autocatalytic reaction\n",
"scheme. The population of species A and B are highly abundant and assumed to be constant values (buffer molecules). The state of the system is thus represented by\n",
"the number of species X.\n",
"\n",
"This are the four reactions, each with a reaction rate $c_i$:\n",
"\\begin{align}\n",
" A + 2X &{\\overset{c_1}{\\rightarrow}} 3X \\\\\n",
" 3X &{\\overset{c_2}{\\rightarrow}} A + 2X \\\\\n",
" B &{\\overset{c_5}{\\rightarrow}} X \\\\\n",
" X &{\\overset{c_6}{\\rightarrow}} B \\\\\n",
"\\end{align}\n",
"\n",
"The initial condition of species is #X = 250, #A = 1.0e5 and #B = 2.0e5. The\n",
"rate constants of reactions are $c_1$ = 3.0e\u20137, $c_2$ = 1.0e\u20134, $c_3$ = 1.0e\u20133 and $c_4$ = 3.5.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"We have reaction rates $c_i$ in units of $s^{-1}$, need to convert these to units of $M^{(1-n)}s^{-1}$ for the parameter `kcst` in STEPS, where $n$ is the order of the reaction, i.e. the number of molecules on the LHS.\n",
"\n",
"The k is calculated diferently depending on the type of the reaction, as shown in the \"2.3: Reaction propensity with mass action kinetics\" of https://link.springer.com/978-3-319-63113-4"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# convert stochastic reaction rate in specified volume to STEPS kcst\n",
"# order is the number of molecules on the left hand side of the reaction\n",
"# maxSameMol is the maximum number of molecules of the same type on the left of the reaction\n",
"def get_kcst(stoch_rate, order, maxSameMol, volume):\n",
" import numpy as np\n",
"\n",
" A = 6.02214076e23 # Avogadro number: molecules in 1 mol\n",
" litres = 1e3 * volume\n",
" if maxSameMol == 1:\n",
" k = stoch_rate / np.power(A * litres, 1 - order)\n",
" elif maxSameMol == 2:\n",
" k = stoch_rate / np.power(A * litres, 1 - order) / 2\n",
" elif maxSameMol == 3:\n",
" k = stoch_rate / np.power(A * litres, 1 - order) / 6\n",
" return k"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy\n",
"\n",
"import steps.model as smodel\n",
"\n",
"mdl = smodel.Model()\n",
"A = smodel.Spec(\"A\", mdl)\n",
"B = smodel.Spec(\"B\", mdl)\n",
"X = smodel.Spec(\"X\", mdl)\n",
"vsys = smodel.Volsys(\"vsys\", mdl)\n",
"\n",
"volume = 1.6667e-21\n",
"\n",
"# Number of molecules\n",
"Na = 1.0e5\n",
"Nb = 2.0e5\n",
"Nx = 250\n",
"\n",
"# stochastic reaction rate\n",
"c1 = 3.0e-7\n",
"c2 = 1.0e-4\n",
"c3 = 1.0e-3\n",
"c4 = 3.5\n",
"\n",
"# Reactions\n",
"R1 = smodel.Reac(\n",
" \"R1\", vsys, lhs=[A, X, X], rhs=[X, X, X], kcst=get_kcst(c1, 3, 2, volume)\n",
")\n",
"R2 = smodel.Reac(\n",
" \"R2\", vsys, lhs=[X, X, X], rhs=[A, X, X], kcst=get_kcst(c2, 3, 3, volume)\n",
")\n",
"R3 = smodel.Reac(\"R3\", vsys, lhs=[B], rhs=[X], kcst=get_kcst(c3, 1, 1, volume))\n",
"R4 = smodel.Reac(\"R4\", vsys, lhs=[X], rhs=[B], kcst=get_kcst(c4, 1, 1, volume))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import steps.geom as swm\n",
"\n",
"wmgeom = swm.Geom()\n",
"\n",
"comp = swm.Comp(\"comp\", wmgeom)\n",
"comp.addVolsys(\"vsys\")\n",
"comp.setVol(volume)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import steps.rng as srng\n",
"import steps.solver as ssolver\n",
"\n",
"rng = srng.create(\"mt19937\", 256)\n",
"rng.initialize(23412)\n",
"\n",
"solver_direct = ssolver.Wmdirect(mdl, wmgeom, rng)\n",
"solver_rssa = ssolver.Wmrssa(mdl, wmgeom, rng)\n",
"\n",
"timeAxis = numpy.arange(0.0, 50.001, 0.001)\n",
"\n",
"\n",
"def simulate(solver, niter=1, initX=Nx):\n",
" solver.reset()\n",
" solver.setCompCount(\"comp\", \"A\", Na)\n",
" solver.setCompCount(\"comp\", \"B\", Nb)\n",
" solver.setCompCount(\"comp\", \"X\", Nx)\n",
"\n",
" import numpy\n",
"\n",
" # Allocate array for 50 sec\n",
" res = numpy.zeros([niter, 50001, 1])\n",
" # Posible X values\n",
" histogram = numpy.zeros([1000])\n",
"\n",
" for i in range(0, niter):\n",
" solver.reset()\n",
" solver.setCompCount(\"comp\", \"A\", Na)\n",
" solver.setCompCount(\"comp\", \"B\", Nb)\n",
" solver.setCompCount(\"comp\", \"X\", Nx)\n",
"\n",
" # A and B are constants\n",
" solver.setCompClamped(\"comp\", \"A\", True)\n",
" solver.setCompClamped(\"comp\", \"B\", True)\n",
" for t in range(0, 50001):\n",
" solver.run(timeAxis[t])\n",
" numberX = solver.getCompCount(\"comp\", \"X\")\n",
" res[i, t, 0] = numberX\n",
" histogram[int(numberX)] += 1\n",
"\n",
" return (res, histogram)\n",
"\n",
"\n",
"def plotres(res, res_=[]):\n",
" # Plot number of molecules of 'molA' over the time range:\n",
" plt.plot(timeAxis, res[:, 0], label=\"X\")\n",
" if len(res_) > 0:\n",
" plt.plot(timeAxis, res_[:, 0], label=\"X'\")\n",
"\n",
" plt.xlabel(\"Time (sec)\")\n",
" plt.ylabel(\"#molecules\")\n",
" plt.legend()\n",
" plt.show()\n",
"\n",
"\n",
"def plothistogram(histogram, histogram_=[]):\n",
" plt.plot(range(1000), histogram, label=\"appearances\")\n",
" if len(histogram_) > 0:\n",
" plt.plot(range(1000), histogram_, label=\"appearances'\")\n",
"\n",
" plt.xlabel(\"#molecules\")\n",
" plt.ylabel(\"appearances\")\n",
" plt.legend()\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Bistability\n",
"The Schlogl model is a remarkable example of a reaction network which exhibits bistability.\n",
"Most of the time just changing the initial population of X by one (250 and 249) we will be able to see 2 stable trends, one in blue around 600 and another on red around 100."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Direct"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res, histogram = simulate(solver_direct)\n",
"res_, histogram_ = simulate(solver_direct, initX=249)\n",
"\n",
"res_mean = numpy.mean(res, 0)\n",
"res_mean_ = numpy.mean(res_, 0)\n",
"\n",
"plotres(res_mean, res_mean_)\n",
"plothistogram(histogram, histogram_)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### RSSA"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"res, histogram = simulate(solver_rssa)\n",
"res_, histogram_ = simulate(solver_rssa, initX=249)\n",
"\n",
"res_mean = numpy.mean(res, 0)\n",
"res_mean_ = numpy.mean(res_, 0)\n",
"\n",
"plotres(res_mean, res_mean_)\n",
"plothistogram(histogram, histogram_)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Benchmark\n",
"Finally we check the execution time of both direct and rssa solvers. Both tend to an average of around 300 X molecules. As we saw before, this is due to the bistability of the model.\n",
"\n",
"Overal rssa performs slightly better."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import time\n",
"\n",
"start = time.time()\n",
"resDirect, histogram = simulate(solver_direct, niter=100)\n",
"end = time.time()\n",
"timeDirect = end - start\n",
"\n",
"start = time.time()\n",
"resRSSA, histogram = simulate(solver_rssa, niter=100)\n",
"end = time.time()\n",
"timeRSSA = end - start\n",
"\n",
"res_mean_Direct = numpy.mean(resDirect, 0)\n",
"res_mean_RSSA = numpy.mean(resRSSA, 0)\n",
"plotres(res_mean_Direct, res_mean_RSSA)\n",
"\n",
"print(\"Time on Direct solver (s)\", timeDirect)\n",
"print(\"Time on RSSA solver (s)\", timeRSSA)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x = numpy.arange(2)\n",
"plt.bar(x, height=[timeDirect, timeRSSA])\n",
"plt.xticks(x, [\"Direct\", \"RSSA\"])\n",
"\n",
"plt.xlabel(\"Solvers\")\n",
"plt.ylabel(\"time(s)\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}