{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#using Pkg\n", "#Pkg.add(\"DifferentialEquations\");\n", "#Pkg.add(\"StatsBase\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 1: systems with more than one species" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. System 1\n", "\n", "$\\emptyset \\xrightarrow{\\kappa_1}\\mathcal{A}\\xrightarrow{\\kappa_2}\\emptyset$\n", "\n", "$\\mathcal{A}\\xrightarrow{\\kappa_3}\\mathcal{B}\\xrightarrow{\\kappa_4}\\emptyset $" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Define a function that computes the RHS of the chemical master equation.\n", "#### You are free to choose whether to implement it in a vector or in a matrix form (e.g. rows for number of $\\mathcal{A}$ molecules, columns for number of $\\mathcal{B}$ molecules)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_upperbound = 30\n", "\n", "## If im matrix form:\n", "Avec = collect(0:n_upperbound)\n", "Bvec = collect(0:n_upperbound)\n", "\n", "p = Iterators.product(Avec, Bvec)\n", "points = collect.(p)\n", "points_reshaped = reduce(hcat,points)\n", "Apoints_matrix = reshape(points_reshaped[1,:], n_upperbound+1, n_upperbound+1)\n", "Bpoints_matrix = reshape(points_reshaped[2,:], n_upperbound+1, n_upperbound+1)\n", "\n", "function AB_system!(du, u, p, t)\n", " \n", " κ₁, κ₂, κ₃, κ₄, V = p\n", "\n", " u ./= sum(u) \n", " \n", " ####\n", " #### Padding here (if necessary)\n", " ####\n", " \n", " A = (-κ₂*####\n", " -κ₃*####\n", " -κ₄*####\n", " -κ₁*####)\n", " B = κ₁*#### \n", " C = κ₂*####\n", " D = κ₃*####\n", " E = κ₄*####\n", " \n", " du .= #### RHS here\n", " \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the initial condition, the time interval of interest and the parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u₀ = zeros(####, ####)\n", "u₀[1, 1] = ####\n", "\n", "T = 20.0\n", "tspan = ####\n", "\n", "κ₁, κ₂, κ₃, κ₄, V = ####\n", "p = [κ₁, κ₂, κ₃, κ₄, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prob = ####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solve the ODE problem" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol = ####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot your numerical solution " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots\n", "\n", "a = heatmap(0:n_upperbound, 0:n_upperbound, sol(0),\n", " c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = 0 sec\")\n", "\n", "b = heatmap(0:n_upperbound, 0:n_upperbound, sol(T/16), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/16) .*\" sec\")\n", "\n", "c = heatmap(0:n_upperbound, 0:n_upperbound, sol(T/8), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/8) .*\" sec\")\n", "\n", "d = heatmap(0:n_upperbound, 0:n_upperbound, sol(T/4), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/4) .*\" sec\")\n", "\n", "e = heatmap(0:n_upperbound, 0:n_upperbound, sol(T/2), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/2) .*\" sec\")\n", "\n", "f = heatmap(0:n_upperbound, 0:n_upperbound, sol(T), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T) .*\" sec\")\n", "\n", "layout = @layout [a b c \n", " d e f]\n", "Plots.plot(a, b, c, d, e, f; layout = layout, colorbar = false, size = (1400, 900), margin = 30Plots.px)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### How about draws with the Gillespie algorithm?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Distributions, StatsBase \n", "\n", "propensities = function(κ, x, ν, Nᵣ)\n", " ν_eval = zeros(Nᵣ)\n", " for r in 1:Nᵣ\n", " ν_eval[r] = ν[r](κ[r], x) \n", " end\n", " return ν_eval\n", "end\n", "\n", "gillespie_alg = function(SM, κ, ν, x₀, stoptime)\n", " t = 0.0\n", " x = x₀\n", " \n", " Nₛ, Nᵣ = size(SM)\n", " \n", " tt = [t]\n", " xx = copy(x) \n", " \n", " k=1\n", " \n", " while (t <= stoptime && k<1e6)\n", " \n", " # step 1\n", " ν_eval = propensities(κ, x, ν, Nᵣ) \n", " α = sum(ν_eval) \n", "\n", " # step 2\n", " τ = rand(Exponential(1/α)) \n", " t += τ\n", "\n", " # step 3\n", " index_j = StatsBase.sample(1:Nᵣ, Weights(ν_eval./α))\n", " \n", " x = x + SM[:, index_j]\n", " append!(tt, t)\n", " xx = cat(xx, x, dims=2)\n", " \n", " k+=1\n", " end\n", " \n", " return tt, xx'\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compute the stoichiometric matrix \n", "Products = ####\n", "Reactants = ####\n", "SM = (Products-Reactants)'\n", "\n", "# Compute the propensity function for each of the reactions\n", "ν₁ = function(κ₁, x) return #### end \n", "ν₂ = function(κ₂, x) return #### end \n", "ν₃ = function(κ₃, x) return #### end \n", "ν₄ = function(κ₄, x) return #### end \n", "\n", "x₀ = [0, 0]\n", "\n", "N = 1e4 # number of realisations\n", "pp = zeros(n_upperbound+1, n_upperbound+1) # matrix of probabilities\n", "\n", "for n in 1:N \n", " tt, xx = gillespie_alg(SM, [κ₁, κ₂, κ₃, κ₄], [ν₁, ν₂, ν₃, ν₄], x₀, T)\n", " \n", " nA, nB = xx[end, :] \n", "\n", " if (nA <= n_upperbound && nB <= n_upperbound) \n", " pp[nA+1, nB+1] += 1\n", " end\n", "end\n", "\n", "pp /= N;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = heatmap(0:n_upperbound, 0:n_upperbound, sol(T),\n", " c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Solution to the \\nChemical Master equation @ t = 20sec\")\n", "\n", "b = heatmap(0:n_upperbound, 0:n_upperbound, pp, c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Empirical density with the \\nGillespie algorithm @ t = 20sec\")\n", "\n", "layout = @layout [a b]\n", "Plots.plot(a, b; layout = layout, colorbar = false, size = (800, 400), margin = 30Plots.px)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. System 2\n", "\n", "$\\mathcal{A}\\xrightarrow{\\kappa_1}\\mathcal{0}$\n", "\n", "$\\mathcal{B}\\xrightarrow{\\kappa_2}\\mathcal{A}$\n", "\n", "$\\emptyset \\xrightarrow{\\kappa_3}\\mathcal{B}$\n", "\n", "$\\mathcal{A}+\\mathcal{B}\\xrightarrow{\\kappa_4}2\\mathcal{A}$\n", "\n", "$\\mathcal{A}+\\mathcal{B}\\xrightarrow{\\kappa_5}2\\mathcal{B}$\n", "\n", "$\\mathcal{A}\\xrightarrow{\\kappa_6}\\emptyset$\n", "\n", "$\\mathcal{B}\\xrightarrow{\\kappa_7}\\emptyset$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Define a function that computes the RHS of the chemical master equation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_upperbound = 20\n", "\n", "## Again, if im matrix form:\n", "Avec = collect(0:n_upperbound)\n", "Bvec = collect(0:n_upperbound)\n", "\n", "p = Iterators.product(Avec, Bvec)\n", "points = collect.(p)\n", "points_reshaped = reduce(hcat,points)\n", "Apoints_matrix = reshape(points_reshaped[1,:], n_upperbound+1, n_upperbound+1)\n", "Bpoints_matrix = reshape(points_reshaped[2,:], n_upperbound+1, n_upperbound+1)\n", "\n", "function AB_system!(du, u, p, t)\n", " \n", " κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇, V = p\n", "\n", " u ./= sum(u) \n", "\n", " l1, l2 = size(u)\n", " u_padded = [zeros(l1) u zeros(l1)]\n", " u_padded = ([zeros(l2+2) u_padded' zeros(l2+2)])'\n", " \n", " A = (-κ₁*####\n", " -κ₂*####\n", " -κ₃*#### \n", " -κ₄*#### \n", " -κ₅*#### \n", " -κ₆*####\n", " -κ₇*####)\n", " B = κ₁*#### \n", " C = κ₂*####\n", " D = κ₃*####\n", " E = κ₄*####\n", " F = κ₅*####\n", " G = κ₆*####\n", " H = κ₇*####\n", "\n", " du .= (A .* #### \n", " +B .* #### \n", " +C .* #### \n", " +D .* #### \n", " +E .* ####\n", " +F .* ####\n", " +G .* ####\n", " +H .* ####)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the initial condition, the time interval of interest and the parameters." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u₀ = zeros(####, ####)\n", "u₀[6, 1] = ####\n", "\n", "T = 20.0\n", "tspan = ####\n", "\n", "κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇, V = ####\n", "p = [κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prob = ####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solve the ODE problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can for example solve the problem using Tsit5()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol = ####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot your numerical solution " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = heatmap(0:n_upperbound, 0:n_upperbound, sol(0),\n", " c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = 0 sec\")\n", "\n", "b = heatmap(0:n_upperbound,\n", " 0:n_upperbound, sol(T/16), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/16) .*\" sec\")\n", "\n", "c = heatmap(0:n_upperbound,\n", " 0:n_upperbound, sol(T/8), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/8) .*\" sec\")\n", "\n", "d = heatmap(0:n_upperbound,\n", " 0:n_upperbound, sol(T/4), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/4) .*\" sec\")\n", "\n", "e = heatmap(0:n_upperbound,\n", " 0:n_upperbound, sol(T/2), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T/2) .*\" sec\")\n", "\n", "f = heatmap(0:n_upperbound,\n", " 0:n_upperbound, sol(T), c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"t = \".* string(T) .*\" sec\")\n", "\n", "layout = @layout [a b c \n", " d e f]\n", "Plots.plot(a, b, c, d, e, f;layout = layout, colorbar = false, size = (1400, 900), margin = 30Plots.px)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### How about draws with the Gillespie algorithm?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Compute the stoichiometric matrix \n", "Products = #### \n", "Reactants = #### \n", "SM = (Products-Reactants)'\n", "\n", "# Compute the propensity function for each of the reactions\n", "ν₁ = function(κ₁, x) return #### end \n", "ν₂ = function(κ₂, x) return #### end \n", "ν₃ = function(κ₃, x) return #### end \n", "ν₄ = function(κ₄, x) return #### end\n", "ν₅ = function(κ₅, x) return #### end\n", "ν₆ = function(κ₆, x) return #### end\n", "ν₇ = function(κ₇, x) return #### end\n", " \n", "x₀ = ####\n", "\n", "N = 1e4 # number of realisations\n", "pp = zeros(n_upperbound+1, n_upperbound+1) # vector of probabilities\n", "\n", "for n in 1:N \n", " tt, xx = gillespie_alg(SM, [κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇], [ν₁, ν₂, ν₃, ν₄, ν₅,ν₆, ν₇], x₀, T)\n", " \n", " nA, nB = xx[end, :] \n", "\n", " if (nA <= n_upperbound && nB <= n_upperbound) \n", " pp[nA+1, nB+1] += 1\n", " end\n", "end\n", "\n", "pp /= N;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "a = heatmap(0:n_upperbound, 0:n_upperbound, sol(T),\n", " c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Solution to the \\nChemical Master equation @ t = 20sec\")\n", "\n", "b = heatmap(0:n_upperbound, 0:n_upperbound, pp, c=:coolwarm,\n", " xlabel=\"# molecules of B\", ylabel=\"# molecules of A\", title=\"Empirical density with the \\nGillespie algorithm @ t = 20sec\")\n", "\n", "layout = @layout [a b]\n", "Plots.plot(a, b; layout = layout, colorbar = false, size = (800, 400), margin = 30Plots.px)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part 2: deterministic versus stochastic modelling" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lotka-Volterra\n", "\n", "Deterministic model:\n", "\n", "$$d\\begin{bmatrix}y_1(t) \\\\ y_2(t) \\end{bmatrix} = \\begin{bmatrix} \\kappa_1 y_1(t)- \\frac{\\kappa_2}{V} y_1(t)y_2(t) \\\\ \\frac{\\kappa_2}{V} y_1(t)y_2(t)-\\kappa_3y_2(t) \\end{bmatrix}dt, \\quad \\begin{bmatrix}y_1(0) \\\\ y_2(0) \\end{bmatrix} = \\begin{bmatrix}100 \\\\ 100 \\end{bmatrix}.$$\n", "\n", "Here, $y = (y_1, y_2)^\\top$ and $\\kappa_1 = 1$, $\\kappa_2 = 5\\times 10^{-3}$, and $\\kappa_3 = 0.6$, and $V = 1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "function lotka_volterra!(du, u, p, t)\n", " κ₁, κ₂, κ₃, V = p\n", " du .= [####,\n", " ####]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the initial condition, the time interval of interest and the parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "u₀ = ####\n", "\n", "T = 15.0\n", "tspan = ####\n", "\n", "κ₁, κ₂, κ₃, V = ####\n", "p = [κ₁, κ₂, κ₃, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prob = ####" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solve the ODE problem" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol = ####\n", "# OR sol2 = solve(prob2, #### (algorithm name); #### (parameters))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### How about draws with the Gillespie algorithm?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Products = ####\n", "Reactants = ####\n", "SM = (Products-Reactants)'\n", "\n", "ν₁ = function(κ₁, y) return #### end \n", "ν₂ = function(κ₂, y) return #### end \n", "ν₃ = function(κ₃, y) return #### end \n", "\n", "y₀ = u₀" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tt, yy = gillespie_alg(SM, [κ₁, κ₂, κ₃], [ν₁, ν₂, ν₃], y₀, T)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(sol, linewidth= 6.0, alpha=0.5, colour=[:darkblue :darkorange], size=(1200,400), linestyle=:dash, labels = [\"y1(t)\" \"y2(t)\"],\n", " xlabel=\"time\", ylabel=\"# individuals\", title=\"Lotka-Volterra model: 1 realisation vs ODE solution\") \n", "plot!(tt, yy, linetype=:steppre, linewidth = 2.0, alpha=1.0, legend=:topleft, labels = [\"realisation of Y1(t)\" \"realisation of Y2(t)\"], colour=[:darkblue :darkorange])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Stochastic mean vs ODE solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "left_bin_edges = 0.0:0.001:T # collection of times \n", "mm = zeros(length(left_bin_edges), 2) # vector to store stochastic mean\n", "\n", "Nₛ = 200 # number of realisations\n", "\n", "for nₛ in 1:Nₛ \n", " tt, yy = gillespie_alg(SM, [κ₁, κ₂, κ₃], [ν₁, ν₂, ν₃], y₀, T)\n", " \n", " ## Computing stochastic mean (code not optimized yet)\n", " ## inputs: left_bin_edges, tt, yy\n", " ## output: yy_extended (binned version of yy)\n", " idxs = searchsortedfirst.(Ref(left_bin_edges[1:end-1]), tt)\n", "\n", " for j in 1:length(idxs)\n", " if j == 1\n", " yy_extended = repeat(yy[j,:]', idxs[j])\n", " else\n", " updates = repeat(yy[j,:]', idxs[j]-idxs[j-1])\n", " if size(updates, 1) > 0\n", " yy_extended = [yy_extended; updates] \n", " else\n", " yy_extended[end, :] += yy[j,:]\n", " end\n", " end\n", " end\n", " ## Filling the end section\n", " if idxs[end] < length(left_bin_edges) \n", " yy_extended = [yy_extended; repeat(yy[end, :]', length(left_bin_edges) - size(yy_extended, 1))]\n", " end\n", "\n", " mm += yy_extended./Nₛ\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot(sol, linewidth= 6.0, alpha=0.5, colour=[:darkblue :darkorange], size=(1200,400), linestyle=:dash, labels = [\"y1(t)\" \"y2(t)\"],\n", " xlabel=\"time\", ylabel=\"# individuals\", title=\"Lotka-Volterra model: Estimate for the stochastic mean vs ODE solution\") \n", "plot!(left_bin_edges, mm, linetype=:steppre, linewidth = 2.0, alpha=1.0, legend=:bottomright, labels = [\"E(Y1(t)] estimate\" \"E[Y2(t)] estimate\"], colour=[:darkblue :darkorange])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.6.2", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.2" } }, "nbformat": 4, "nbformat_minor": 4 }