diff --git a/notebooks/Schlogl_Model.ipynb b/notebooks/Schlogl_Model.ipynb new file mode 100644 index 0000000..ec248bd --- /dev/null +++ b/notebooks/Schlogl_Model.ipynb @@ -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 +}