intro_to_det_stoch_systems_part_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to deterministic and stochastic modelling of simple biological systems: part 1\n",
    "## Dr Marc Sturrock"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DiffEqBiological.jl is a domain specific language (DSL) for writing chemical\n",
    "reaction networks in Julia. The generated chemical reaction network model can\n",
    "then be translated into a variety of mathematical models which can be solved\n",
    "using components of the broader\n",
    "[DifferentialEquations.jl](http://juliadiffeq.org/) ecosystem. Biomodelling.jl is a package for simulating gene expression networks, cell growth, division and partitioning. \n",
    "\n",
    "In this tutorial we'll provide an introduction to using DiffEqBiological to\n",
    "specify chemical reaction networks, and then to solve deterministic and stochastic models generated from them. In the second part, we will then use Biomodelling.jl to understand the concepts of intrinsic and extrinsic noise.\n",
    "\n",
    "\n",
    "Let's start by using the DiffEqBiological\n",
    "`@reaction_network` macro to specify a simply chemical reaction network; the\n",
    "reactions associated with the central dogma of molecular biology.\n",
    "\n",
    "We first add and import the basic packages we'll need, and use Plots.jl for making\n",
    "figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First load packages\n",
    "using DifferentialEquations, DiffEqBiological, Plots, Latexify\n",
    "gr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now construct the reaction network. The basic types of arrows and predefined rate laws one can use are discussed in detail within the DiffEqBiological [Chemical Reaction Models documentation](http://docs.juliadiffeq.org/latest/models/biological.html). Here we use a mix of zero order and first order rate laws. Note, $\\varnothing$ corresponds to the empty state, and is used for zeroth order production and first order degradation reactions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geneexpression = @reaction_network begin\n",
    "    b_1, 0 --> mRNA\n",
    "    d_1, mRNA --> 0\n",
    "    b_2, mRNA --> mRNA + protein\n",
    "    d_2, protein --> 0\n",
    "end b_1 d_1 b_2 d_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use Latexify to look at the corresponding reactions and understand the generated rate laws for each reaction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(geneexpression; env=:chemical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use Latexify to look at the corresponding ODE model for the chemical system which is created using the law of mass action"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(geneexpression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To solve the ODEs we need to specify the values of the parameters in the model, the initial condition, and the time interval to solve the model on. To do this it helps to know the orderings of the parameters and the species. Parameters are ordered in the same order they appear after the end statement in the `@reaction_network` macro. Species are ordered in the order they first appear within the @reaction_network macro. We can see these orderings using the speciesmap and paramsmap functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speciesmap(geneexpression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "paramsmap(geneexpression)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving a deterministic version of the system: ODEs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Knowing these orderings, we can create parameter and initial condition vectors, and setup the ODEProblem we want to solve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters [b_1,d_1,b_2,d_2]\n",
    "p = (10.0,0.1,1.0,0.001)\n",
    "\n",
    "# initial condition [mRNA0,protein0]\n",
    "u₀ = [0.0,0.0]\n",
    "\n",
    "# time interval to solve on\n",
    "tspan = (0., 5000.0)\n",
    "\n",
    "# create the ODEProblem we want to solve\n",
    "prob = ODEProblem(geneexpression, u₀, tspan, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we are all set to solve the ODEs. We can now use any ODE solver from within the DiffEq package. We'll just use the default DifferentialEquations solver for now, and then plot the solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol1 = solve(prob, saveat=10.)\n",
    "plot(sol1,layout=(1,2),legend=:bottomright,ylabel=[\"mRNA\" \"protein\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more on choices of ODE solvers, see the JuliaDiffEq [documentation](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving a stochastic version of the system: Stochastic simulation algorithm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now look at a stochastic chemical kinetics model of the repressilator, modeling it with jump processes. Here we will construct a DiffEqJump JumpProblem that uses Gillespie's Direct method, and then solve it to generate one realization of the jump process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we redefine the initial condition to be integer valued\n",
    "u₀ = [10,1000]\n",
    "\n",
    "# next we create a discrete problem to encode that our species are integer valued:\n",
    "dprob = DiscreteProblem(geneexpression, u₀, tspan, p)\n",
    "\n",
    "# now we create a JumpProblem, and specify Gillespie's Direct Method as the solver:\n",
    "jprob = JumpProblem(dprob, Direct(), geneexpression, save_positions=(false,false))\n",
    "\n",
    "# now let's solve and plot the jump process:\n",
    "sol2 = solve(jprob, SSAStepper(), saveat=10.)\n",
    "plot(sol2,layout=(1,2),legend=:bottomright,ylabel=[\"mRNA\" \"protein\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While SSAs generate exact realisations for stochastic chemical kinetics jump process models, a problem with the SSA is that it can be prohibitively slow, particularly when large numbers of molecules are simulated.  $\\tau$-leaping methods offer a fast alternative by discretising in time the underlying time-change representation of the stochastic process. The DiffEqJump package has limited support for $\\tau$-leaping methods in the form of the basic Euler's method type approximation proposed by Gillespie (Biomodelling.jl contains more tau leaping algorithms). We can simulate a $\\tau$-leap approximation to the repressilator by using the  RegularJump representation of the network to construct a JumpProblem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rjs = regularjumps(geneexpression)\n",
    "lprob = JumpProblem(dprob, Direct(), rjs)\n",
    "lsol = solve(lprob, SimpleTauLeaping(), dt=.1)\n",
    "plot(lsol,layout=(1,2),legend=:bottomright,ylabel=[\"mRNA\" \"protein\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By taking the mean of the stochastic solutions across 100 different simulations (representing gene expression in 100 different cells), we can see that the ordinary differential equation model gives an accurate mean-field representation of gene expression for a cell population. This is not parameter set dependent and you can try changing p defined above to whatever you like and rerunning cells to confirm this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = zeros(2,50001,10)\n",
    "for i = 1:10\n",
    "    lsol = solve(lprob, SimpleTauLeaping(), dt=.1)\n",
    "    output[:,:,i] = lsol[1:2,:]\n",
    "end\n",
    "using Statistics\n",
    "mean_stochastic_sol = mean(output,dims=3)[:,:,1]\n",
    "plot(sol1)\n",
    "plot!(lsol.t,mean_stochastic_sol',linestyle=:dash,linewidth=4,labels=[\"mRNA stochastic\" \"protein stochastic\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we say this in general? Is the ODE model solution always consistent with the mean of the stochastic model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## mRNA-miRNA-protein system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The paper of [Paulsson et al.](https://www.pnas.org/content/97/13/7148) provides a neat example of how even a simple system can yield different results for the stochatic and deterministic analogues. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we define the system, which consists of zero and first order reactions as well as a single second order reaciton. This system can be interpreted in a number of different ways, for example, $I$ can be interpreted as an mRNA molecule, $P$ as a protein molecule, and $S$ as a microRNA molecule  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = @reaction_network begin\n",
    "   k, 0 --> I\n",
    "   ka*S, I --> 0\n",
    "   kp, I --> P\n",
    "   1, P --> 0\n",
    "   b, 0 --> S\n",
    "   d, S --> 0\n",
    "end k ka kp b d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system can be expressed as stochastic chemical reactions in the following way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(ps; env=:chemical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or as deterministic ODEs as the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we can solve the deterministic ODE system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters [k,ka,kp,b,d]\n",
    "params = (1000.0,100.0,100.0,1.0,2.0)\n",
    "\n",
    "#initial conditions\n",
    "u0 = [1.0, 1.0, 1.0]\n",
    "\n",
    "# time interval to solve on\n",
    "tspan = (0., 100.)\n",
    "\n",
    "# create the ODEProblem we want to solve\n",
    "oprob = ODEProblem(ps, u0, tspan, params)\n",
    "\n",
    "# time points to output solution at\n",
    "timestep = 1.0\n",
    "ts = collect(0:timestep:maximum(tspan));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the ODE system leads to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_sol  = solve(oprob, alg_hint=:stiff,saveat=ts);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the corresponding stochastic system for 10 different cells, averaging and plotting alongside the ODE results leads to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dprob = DiscreteProblem(ps, u0, tspan, params)\n",
    "jprob = JumpProblem(dprob, Direct(), ps, save_positions = (false, false))\n",
    "output = zeros(Int64,100, length(ts));\n",
    "for i = 1:100\n",
    "   jsol = solve(jprob, SSAStepper(), saveat = ts)\n",
    "   output[i,:] = jsol[2,:]\n",
    "end\n",
    "\n",
    "plot(ode_sol.t,output',title=\"\",labels=\"\",xlabel=\"\",ylabel=\"number of protein molecules\",grid=false,linewidth=2.0)\n",
    "plot!(ode_sol.t,mean(output,dims=1)',label=\"mean SSA\",linewidth=6.0,c=:blue)\n",
    "plot!(ode_sol.t,ode_sol[2,:],label=\"ODE\",linestyle=:dash,linewidth=6.0,c=:black,xlabel=\"time\",title=\"b = $(params[4]), d = $(params[5])\",legend=:bottomright)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see there is a discrepancy between the mean number of protein molecules between the stochastic system and solution of the deterministic system. However, we find that by increasing the reaction rates (while keeping the number of miRNA molecules constant), we get good agreement. Hence for some stochastic systems "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = (1000.0,100.0,100.0,1000.0,2000.0)\n",
    "oprob = ODEProblem(ps, u0, tspan, params)\n",
    "dprob = DiscreteProblem(ps, u0, tspan, params)\n",
    "jprob = JumpProblem(dprob, Direct(), ps, save_positions = (false, false))\n",
    "output = zeros(Int64,100, length(ts));\n",
    "for i = 1:100\n",
    "   jsol = solve(jprob, SSAStepper(), saveat = ts)\n",
    "   output[i,:] = jsol[2,:]\n",
    "end\n",
    "\n",
    "plot(ode_sol.t,output',title=\"\",labels=\"\",xlabel=\"\",ylabel=\"number of protein molecules\",grid=false,linewidth=2.0)\n",
    "plot!(ode_sol.t,mean(output,dims=1)',label=\"mean SSA\",linewidth=6.0,c=:blue)\n",
    "plot!(ode_sol.t,ode_sol[2,:],label=\"ODE\",linestyle=:dash,linewidth=6.0,c=:black,xlabel=\"time\",title=\"b = $(params[4]), d = $(params[5])\",legend=:bottomright,ylims=(0,1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As soon as a stochastic system contains a second order reaction (e.g. A + B --> C) it is possible that the mean of the stochastic system will deviate from the corresponding ODE system.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repressilator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The repressilator is a genetic regulatory network consisting of at least one feedback loop with at least three genes, each expressing a protein that represses the next gene in the loop. In biological research, repressilators have been used to build cellular models and understand cell function. There are both artificial and naturally-occurring repressilators. Recently, the naturally-occurring repressilator clock gene circuit in Arabidopsis thaliana (A. thaliana) and mammalian systems have been studied. The repressilator can be defined as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repressilator = @reaction_network begin\n",
    "    hillr(P₃,α,K,n), ∅ --> m₁\n",
    "    hillr(P₁,α,K,n), ∅ --> m₂\n",
    "    hillr(P₂,α,K,n), ∅ --> m₃\n",
    "    (δ,γ), m₁ ↔ ∅\n",
    "    (δ,γ), m₂ ↔ ∅\n",
    "    (δ,γ), m₃ ↔ ∅\n",
    "    β, m₁ --> m₁ + P₁\n",
    "    β, m₂ --> m₂ + P₂\n",
    "    β, m₃ --> m₃ + P₃\n",
    "    μ, P₁ --> ∅\n",
    "    μ, P₂ --> ∅\n",
    "    μ, P₃ --> ∅\n",
    "end α K n δ γ β μ;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again we can use Latexify to look at the corresponding reactions and understand the generated rate laws for each reaction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(repressilator; env=:chemical)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use Latexify to look at the corresponding ODE model for the chemical system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latexify(repressilator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we can solve the deterministic ODE system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters [α,K,n,δ,γ,β,μ]\n",
    "p = (.5, 40, 2, log(2)/120, 5e-3, 20*log(2)/120, log(2)/60)\n",
    "\n",
    "# initial condition [m₁,m₂,m₃,P₁,P₂,P₃]\n",
    "u₀ = [0.,0.,0.,20.,0.,0.]\n",
    "\n",
    "# time interval to solve on\n",
    "tspan = (0., 10000.)\n",
    "timestep = 1.0\n",
    "ts = collect(0:timestep:maximum(tspan));\n",
    "\n",
    "# create the ODEProblem we want to solve\n",
    "oprob = ODEProblem(repressilator, u₀, tspan, p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve(oprob, saveat=ts)\n",
    "plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution looks oscillatory in nature, we can see if the same is true for the mean of trajectories from the stochastic model..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve JumpProblem using Gillespie’s Direct Method\n",
    "u0 =  [0,0,0,20,0,0] # change to int type for stochastic system..\n",
    "maxcell = 1000 #number of cells to simulate for\n",
    "dprob = DiscreteProblem(repressilator, u0, tspan, p)\n",
    "jprob = JumpProblem(dprob, Direct(), repressilator, save_positions = (false, false))\n",
    "output = zeros(Int64,maxcell, length(ts),6);\n",
    "for i = 1:maxcell\n",
    "   jsol = solve(jprob, SSAStepper(), saveat = ts)\n",
    "    for j = 1:6\n",
    "       output[i,:,j] = jsol[j,:]\n",
    "    end\n",
    "end\n",
    "\n",
    "gr(size=(1200,600))\n",
    "plot(sol.t,mean(output,dims=1)[1,:,:],title=\"\",label=[\"\" \"\" \"\" \"\" \"\" \"mean SSA\"],xlabel=\"\",ylabel=[\"number of molecules\" \"\" \"\" \"number of molecules\" \"\" \"\"],grid=false,linewidth=6.0,layout=(2,3))\n",
    "plot!(sol.t,output[1,:,:],title=\"\",label=[\"\" \"\" \"\" \"\" \"\" \"single SSA\"],grid=false,linewidth=2.0,layout=(2,3))\n",
    "plot!(sol,label=[\"\" \"\" \"\" \"\" \"\" \"ODE\"],linestyle=:dash,linewidth=6.0,c=:black,legend=:topright,xlabel=[\"\" \"\" \"\" \"time\" \"time\" \"time\"],title=[\"m1\" \"m2\" \"m3\" \"P1\" \"P2\" \"P3\"],layout=(2,3),grid=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see qualitative (but not quantitative) agreement between the deterministic ODE solution and a single trajectory of the stochastic system. However, the mean of several stochastic simulations differs both qualitatively and quantitatively from the ODE solutions. "
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}