{ "cells": [ { "cell_type": "code", "execution_count": 2, "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AB_system! (generic function with 1 method)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n_upperbound = 30\n", "\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 = -κ₂*Apoints_matrix - κ₃*Apoints_matrix - κ₄*Bpoints_matrix - κ₁*V*ones(l1, l2)\n", " B = κ₁*V*ones(l1,l2) \n", " C = κ₂.*(Apoints_matrix.+1)\n", " D = κ₃.*(Apoints_matrix.+1)\n", " E = κ₄.*(Bpoints_matrix.+1)\n", "\n", " du .= A .* u + B .* u_padded[1:end-2, 2:end-1] + C .* u_padded[3:end, 2:end-1] + D .* u_padded[3:end, 1:end-2] + E .* u_padded[2:end-1, 3:end]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the initial condition, the time interval of interest and the parameters." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "u₀ = zeros(n_upperbound+1, n_upperbound+1)\n", "u₀[1, 1] = 1\n", "\n", "T = 20.0\n", "tspan = (0, T)\n", "\n", "κ₁, κ₂, κ₃, κ₄, V = [2, 0.2, 0.2, 0.1, 1]\n", "p = [κ₁, κ₂, κ₃, κ₄, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[36mODEProblem\u001b[0m with uType \u001b[36mMatrix{Float64}\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n", "timespan: (0.0, 20.0)\n", "u0: 31×31 Matrix{Float64}:\n", " 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " ⋮ ⋮ ⋱ ⋮ ⋮\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prob = ODEProblem(AB_system!, u₀, tspan, p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solve the ODE problem" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sol = solve(prob, Tsit5(); dt=0.001);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot your numerical solution " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "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,\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", "\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": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "#3 (generic function with 1 method)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "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": 10, "metadata": {}, "outputs": [], "source": [ "# Compute the stoichiometric matrix \n", "Products = [1 0; \n", " 0 0; \n", " 0 1; \n", " 0 0]\n", "Reactants = [0 0; 1 0; 1 0; 0 1]\n", "SM = (Products-Reactants)'\n", "\n", "# Compute the propensity function for each of the reactions\n", "ν₁ = function(νₛ, x) return νₛ*V end \n", "ν₂ = function(κ₂, x) return κ₂*x[1] end \n", "ν₃ = function(κ₃, x) return κ₃*x[1] end \n", "ν₄ = function(κ₄, x) return κ₄*x[2] end \n", "\n", "x₀ = [0, 0]\n", "\n", "N = 1e4 #number of realisations\n", "n_upperbound = 30\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": 11, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Plots\n", "\n", "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": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "AB_system! (generic function with 1 method)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n_upperbound = 20\n", "\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 = (-κ₁*Apoints_matrix\n", " -κ₂*Bpoints_matrix \n", " -κ₃*V*ones(l1,l2) \n", " -κ₄*Apoints_matrix.*Bpoints_matrix./2V \n", " -κ₅*Apoints_matrix.*Bpoints_matrix./2V \n", " -κ₆*Apoints_matrix\n", " -κ₇*Bpoints_matrix)\n", " B = κ₁*(Apoints_matrix.+1) \n", " C = κ₂*(Bpoints_matrix.+1)\n", " D = κ₃*V*ones(l1, l2)\n", " E = κ₄*(Bpoints_matrix.+1).*(Apoints_matrix.-1)/2V\n", " F = κ₅*(Bpoints_matrix.-1).*(Apoints_matrix.+1)/2V\n", " G = κ₆*(Apoints_matrix.+1)\n", " H = κ₇*(Bpoints_matrix.+1)\n", "\n", " du .= (A .* u \n", " +B .* u_padded[3:end, 1:end-2] \n", " +C .* u_padded[1:end-2, 3:end] \n", " +D .* u_padded[2:end-1, 1:end-2] \n", " +E .* u_padded[1:end-2, 3:end]\n", " +F .* u_padded[3:end, 1:end-2]\n", " +G .* u_padded[3:end, 2:end-1]\n", " +H .* u_padded[2:end-1, 3:end])\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Set the initial condition, the time interval of interest and the parameters." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "u₀ = zeros(n_upperbound+1, n_upperbound+1)\n", "u₀[6, 1] = 1\n", "\n", "T = 20.0\n", "tspan = (0, T)\n", "\n", "κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇, V = [0.001, 0.05, 1, 0.5, 0.5, 0.1, 0.1, 1]\n", "p = [κ₁, κ₂, κ₃, κ₄, κ₅, κ₆, κ₇, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[36mODEProblem\u001b[0m with uType \u001b[36mMatrix{Float64}\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n", "timespan: (0.0, 20.0)\n", "u0: 21×21 Matrix{Float64}:\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0 0.0" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prob = ODEProblem(AB_system!, u₀, tspan, p)" ] }, { "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": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "retcode: Success\n", "Interpolation: specialized 4th order \"free\" interpolation\n", "t: 1719-element Vector{Float64}:\n", " 0.0\n", " 0.0008935239142316075\n", " 0.009828763056547682\n", " 0.03321049437397584\n", " 0.06535569275601097\n", " 0.10599044662756726\n", " 0.15816166787338787\n", " 0.22430127209338013\n", " 0.3092836721569964\n", " 0.41741061245211064\n", " 0.5535437715196642\n", " 0.7230672238750465\n", " 0.9321321849319046\n", " ⋮\n", " 19.88548773617738\n", " 19.89661095274366\n", " 19.907734168867943\n", " 19.918857384550225\n", " 19.929980602000498\n", " 19.941103819008774\n", " 19.952227035575053\n", " 19.963350251699335\n", " 19.974473468854942\n", " 19.985596685568552\n", " 19.996719901840166\n", " 20.0\n", "u: 1719-element Vector{Matrix{Float64}}:\n", " [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]\n", " [5.68892774752621e-21 2.8967436979652444e-22 … 0.0 0.0; 3.1833903654196313e-16 1.3014206600950795e-17 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]\n", " [9.058327907898058e-16 5.437318577253168e-17 … 0.0 0.0; 4.60696569263814e-12 2.2881089715464695e-13 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]\n", " [3.87498730448105e-13 3.225982824192406e-14 … 0.0 0.0; 5.825508798628553e-10 4.218503469290466e-11 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]\n", " [1.0990215616238068e-11 1.2623116236598568e-12 … 2.2406039813618214e-35 1.5629105519431998e-38; 8.380570227468697e-9 8.643969440441407e-10 … 3.5581862049322706e-36 2.2023276729809498e-39; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]\n", " [1.1722946802919428e-10 1.8144793700989867e-11 … 7.877680626004728e-30 3.9449014967960807e-32; 5.5005090360056405e-8 7.75976015665223e-9 … 1.0960130924727072e-30 4.914089310657494e-33; … ; 1.4330349394473182e-44 7.148998025535565e-47 … 0.0 0.0; 4.944153829995199e-50 0.0 … 0.0 0.0]\n", " [8.132076112285438e-10 1.6729444522334766e-10 … 2.467458628342865e-26 2.5574162917129033e-28; 2.550182320243015e-7 4.802610544803399e-8 … 2.9556895568818128e-27 2.710477659537235e-29; … ; 9.234571500568983e-36 3.857597241362522e-37 … 0.0 0.0; 3.2294221447050595e-39 6.667004693446639e-41 … 0.0 0.0]\n", " [4.300487810603634e-9 1.1600114314467182e-9 … 1.3528881213813374e-23 2.15218621887971e-25; 9.477361672909733e-7 2.3311641137667708e-7 … 1.4110855320986711e-24 1.968456498297167e-26; … ; 7.524656460907391e-31 7.220194343337442e-32 … 0.0 0.0; 1.660423899194648e-33 1.4589294269122813e-34 … 0.0 0.0]\n", " [1.9320711766899166e-8 6.782936665629485e-9 … 3.148736545264884e-21 7.03149646853206e-23; 3.0745455646415475e-6 9.731324319399189e-7 … 2.8512878333053612e-22 5.548087366967353e-24; … ; 3.574840370445302e-27 4.328529412381049e-28 … 4.5488609740129714e-66 1.0892740139768127e-69; 2.1134554033539345e-29 2.0623829337880095e-30 … 0.0 0.0]\n", " [7.587594229767366e-8 3.4376511991059396e-8 … 4.011570388450536e-19 1.2047791909188304e-20; 8.897252734415552e-6 3.5665292584854826e-6 … 3.1422241695024396e-20 8.190246908586611e-22; … ; 4.5428060943126976e-24 3.98991331358139e-25 … -1.5698782821578377e-54 -4.3727365483426176e-58; 5.783997789688429e-26 -8.530131329696206e-27 … 2.2785190205844754e-58 1.6647975973927066e-61]\n", " [2.642809566306407e-7 1.530544865564229e-7 … 3.1303773134191336e-17 1.2359906261676859e-18; 2.320659812457713e-5 1.1581563476661848e-5 … 2.0499987216671185e-18 8.424334445791986e-20; … ; 3.908467382361656e-21 -4.936739536862531e-21 … -2.0117391436648936e-44 4.1649912568676165e-47; 2.6506151317701944e-22 -6.036659958205531e-22 … -1.830624803066249e-47 2.3971374788094056e-50]\n", " [8.220897522218591e-7 6.022370466529499e-7 … 1.365824072622452e-15 1.5320915922144799e-16; 5.4788515148385576e-5 3.34231337686609e-5 … -1.2416212873125891e-15 3.764140563994684e-16; … ; 4.625430973632303e-17 -2.566229055885219e-16 … -1.7995949638314365e-35 7.015499659569123e-38; 1.1732515078509493e-17 -5.737314640890776e-17 … -4.711928796496514e-38 1.2814369044790266e-40]\n", " [2.293207322695447e-6 2.1016933459087046e-6 … -8.327890841439829e-12 4.85388176414393e-12; 0.0001173043968190918 8.581193692193065e-5 … -9.286689393806823e-11 4.6304680844936715e-11; … ; 3.2408941488356846e-12 -3.454916481394463e-11 … -1.67701059521488e-26 9.734592097444134e-29; 1.6447964398890694e-12 -1.4817457064959472e-11 … -7.59278967114399e-29 2.902039664048834e-31]\n", " ⋮\n", " [8.565487238929868e-5 0.0005438598534505712 … 0.0006040321285447826 0.00020184160515099417; 0.0002633593185828176 0.0006438785083821997 … 4.578798752666225e-5 8.500011870054305e-6; … ; 0.00012823095422950898 2.6610094399973347e-5 … -2.5247844942761012e-6 9.198780583343948e-8; 4.138995809696703e-5 8.386113186153318e-6 … -8.898111098782845e-8 8.145288237316326e-10]\n", " [8.560036265742185e-5 0.0005435617924523537 … 0.0006045100688759206 0.00020200974955493167; 0.0002631865236412175 0.0006435096674504114 … 4.5826406866921506e-5 8.507474851714702e-6; … ; 0.00012834869078081173 2.6634605781278113e-5 … -2.524784301236347e-6 9.198779880024969e-8; 4.1429805396562647e-5 8.394229762890113e-6 … -8.898110418452326e-8 8.145287614545261e-10]\n", " [8.554593551780708e-5 0.0005432641149955555 … 0.0006049876366337407 0.00020217776884221694; 0.000263014052150995 0.0006431414341918903 … 4.586479832289413e-5 8.514932751601975e-6; … ; 0.00012846636367581704 2.6659103836910668e-5 … -2.5247835166338452e-6 9.198777021411034e-8; 4.146963213673437e-5 8.402342143760684e-6 … -8.89810765327403e-8 8.145285083314287e-10]\n", " [8.54915908380168e-5 0.0005429668206629497 … 0.0006054648317115505 0.00020234566295667967; 0.0002628419033714616 0.0006427738074074723 … 4.5903161880786564e-5 8.522385566248327e-6; … ; 0.00012858397285093041 2.6683588553347892e-5 … -2.524782140329637e-6 9.198772006994617e-8; 4.1509438292127545e-5 8.410450323460498e-6 … -8.898102802758239e-8 8.145280643175266e-10]\n", " [8.543732847491836e-5 0.0005426699089775325 … 0.0006059416540981171 0.00020251343187578966; 0.00026267007652950745 0.000642406785827133 … 4.59414975344913e-5 8.529833293680686e-6; … ; 0.00012870151826607605 2.670805992196475e-5 … -2.5247833418834664e-6 9.198776384726185e-8; 4.154922384536067e-5 8.418554298309764e-6 … -8.89810703740002e-8 8.145284519546181e-10]\n", " [8.538314830735917e-5 0.000542373379581598 … 0.0006064181035957137 0.00020268107551141222; 0.00026249857092304057 0.0006420403683312406 … 4.5979805262908066e-5 8.53727592901326e-6; … ; 0.0001288189998351149 2.6732517924547905e-5 … -2.524783952049868e-6 9.198778607801406e-8; 4.158898876346643e-5 8.426654061454513e-6 … -8.898109187813187e-8 8.145286488023949e-10]\n", " [8.53290502034588e-5 0.0005420772320575197 … 0.000606894180102212 0.0002028485938090993; 0.00026232738581740707 0.0006416745537290474 … 4.601808505263301e-5 8.544713468856855e-6; … ; 0.00012893641749545052 2.6756962547787057e-5 … -2.5247839706889233e-6 9.198778675710462e-8; 4.162873302145827e-5 8.43474960766707e-6 … -8.898109253501884e-8 8.145286548155083e-10]\n", " [8.527503403152272e-5 0.0005417814659875471 … 0.0006073698835169968 0.00020301598671497388; 0.00026215652047995275 0.0006413093408325978 … 4.605633689039155e-5 8.552145909848354e-6; … ; 0.0001290537711848176 2.6781393778442545e-5 … -2.5247833976612445e-6 9.198776587946245e-8; 4.1668456594471814e-5 8.442840931745267e-6 … -8.898107233977638e-8 8.145284699491695e-10]\n", " [8.522109965283186e-5 0.0005414860809142537 … 0.0006078452138033834 0.00020318325419770394; 0.0002619859741572274 0.0006409447284079246 … 4.60945607680601e-5 8.559573249626865e-6; … ; 0.00012917106085670355 2.6805811606555916e-5 … -2.5247843459604683e-6 9.198780042972792e-8; 4.170815946298682e-5 8.450928029576148e-6 … -8.898110576074166e-8 8.145287758831309e-10]\n", " [8.516724694329516e-5 0.0005411910764593093 … 0.000608320170801243 0.0002033503961825392; 0.00026181574614343094 0.0006405807153215859 … 4.613275666759151e-5 8.566995483903564e-6; … ; 0.00012928828643406283 2.6830216015813736e-5 … -2.524784702709804e-6 9.198781342750498e-8; 4.1747841597154303e-5 8.45901089494344e-6 … -8.898111833367628e-8 8.145288909751699e-10]\n", " [8.511347577177548e-5 0.000540896452204615 … 0.0006087947544144679 0.00020351741261730456; 0.00026164583571187077 0.0006402173003939546 … 4.617092457609637e-5 8.574412609393068e-6; … ; 0.000129405447855618 2.685460699318702e-5 … -2.5247844677687006e-6 9.198780486767983e-8; 4.1787502972474364e-5 8.467089522721085e-6 … -8.898111005363613e-8 8.145288151800624e-10]\n", " [8.509769005895239e-5 0.0005408099944098588 … 0.0006089350263452843 0.00020356677156094133; 0.0002615959614225145 0.0006401106629907257 … 4.618220438412285e-5 8.57660440652505e-6; … ; 0.00012944006887733253 2.686181442003586e-5 … -8.983795658957117e-7 3.2731492632224116e-8; 4.179922172259884e-5 8.469476484900958e-6 … -3.1661637673789745e-8 2.8982911323028054e-10]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = solve(prob, Tsit5())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot your numerical solution " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "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,\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", "\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": 17, "metadata": {}, "outputs": [], "source": [ "# Compute the stoichiometric matrix \n", "Products = [0 1; 1 0; 0 1; 2 0; 0 2; 0 0; 0 0]\n", "Reactants = [1 0; 0 1; 0 0; 1 1; 1 1; 1 0; 0 1]\n", "SM = (Products-Reactants)'\n", "\n", "# Compute the propensity function for each of the reactions\n", "ν₁ = function(κ₁, x) return κ₁*x[1] end \n", "ν₂ = function(κ₂, x) return κ₂*x[2] end \n", "ν₃ = function(κ₃, x) return κ₃*V end \n", "ν₄ = function(κ₄, x) return κ₄*x[1]*x[2] end\n", "ν₅ = function(κ₅, x) return κ₅*x[1]*x[2] end\n", "ν₆ = function(κ₆, x) return κ₆*x[1] end\n", "ν₇ = function(κ₇, x) return κ₇*x[2] end\n", " \n", "x₀ = [5, 0]\n", "\n", "N = 1e4 #number of realisations\n", "\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": 18, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "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": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "lotka_volterra! (generic function with 1 method)" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function lotka_volterra!(du, u, p, t)\n", " κ₁, κ₂, κ₃, V = p\n", " du .= [κ₁*u[1] - κ₂*u[1]*u[2]/V,\n", " κ₂*u[1]*u[2]/V-κ₃*u[2]]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set the initial condition, the time interval of interest and the parameters" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "u₀ = [100, 100]\n", "\n", "T = 15.0\n", "tspan = (0, T)\n", "\n", "κ₁, κ₂, κ₃, V = 1, 0.005, 0.6, 1.0\n", "p = [κ₁, κ₂, κ₃, V];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create an ordinary differential equation problem with the variables defined above." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[36mODEProblem\u001b[0m with uType \u001b[36mVector{Int64}\u001b[0m and tType \u001b[36mFloat64\u001b[0m. In-place: \u001b[36mtrue\u001b[0m\n", "timespan: (0.0, 15.0)\n", "u0: 2-element Vector{Int64}:\n", " 100\n", " 100" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prob = ODEProblem(lotka_volterra!, u₀, tspan, p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the ODE problem" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "sol = solve(prob)\n", "# OR sol2 = solve(prob2, #### (algorithm name); #### (parameters))" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Int64}:\n", " 100\n", " 100" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Compute the stoichiometric matrix \n", "Products = [2 0; 0 2; 0 0]\n", "Reactants = [1 0; 1 1; 0 1]\n", "SM = (Products-Reactants)'\n", "\n", "# Compute the propensity function for each of the reactions\n", "ν₁ = function(κ₁, y) return κ₁*y[1] end \n", "ν₂ = function(κ₂, y) return κ₂*y[1]*y[2]/V end \n", "ν₃ = function(κ₃, y) return κ₃*y[2] end \n", "\n", "T = 15.0\n", "y₀ = u₀" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "([0.0, 0.0012044985931863469, 0.00134141238889642, 0.0016697405556296407, 0.015746572237407157, 0.01751603173329227, 0.024299751591754112, 0.032463126450146707, 0.039558235432356155, 0.039828653767716146 … 14.963676998671465, 14.971657771198252, 14.977623837088712, 14.982827443875525, 14.983067962790312, 14.983983834267862, 14.98443805519641, 14.991874422596583, 14.994913360293769, 15.009874846727058], [100 100; 99 101; … ; 9 196; 9 195])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "tt, yy = gillespie_alg(SM, [κ₁, κ₂, κ₃], [ν₁, ν₂, ν₃], y₀, T)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "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": 26, "metadata": {}, "outputs": [], "source": [ "left_bin_edges = 0.0:0.001:T\n", "mm = zeros(length(left_bin_edges), 2)\n", "\n", "Nₛ = 200\n", "\n", "for nₛ in 1:Nₛ \n", " tt, yy = gillespie_alg(SM, [κ₁, κ₂, κ₃], [ν₁, ν₂, ν₃], y₀, T)\n", " \n", " \n", " ## Computing stochastic mean (code not optimized yet)\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", "\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": 27, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "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])" ] } ], "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 }