{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Classroom problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Comparment-based approach to Diffusion" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Distributions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\"\"\"\n", " L - length in mm of the one-dimensional domain where molecules of chemical species A diffuse \n", " K - number of compartments \n", " h - length of a single compartment in μm (volume: h³)\n", " D - Diffusion constant in mm²sec⁻¹\n", " d - rate constant (d = D/h^2)\n", "\"\"\"\n", "\n", "L = 1 \n", "K = 40\n", "h = #### \n", "D = 0.0001\n", "d = ####;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gillespie_alg_pure_diffusion = function(A₀, Nt)\n", "\n", " \"\"\"\n", " Inputs:\n", " A₀ - a K-dimensional vector with the initial amount of molecules insided each compartment. \n", " Nₜ - maximum number of jumps between compartments.\n", " Outputs:\n", " tt - a partitioned time interval [0, T] with times sampled from an exponential distribution\n", " AA - the molecule count for each time t ∈ tt\n", " \"\"\"\n", " \n", " t = 0.0\n", "\n", " tt = [t]\n", " \n", " AA = zeros(Nt, K)\n", " AA[1, :] = A₀\n", "\n", " \n", " for t_idx in 1:Nt-1\n", " \n", " # Step 1\n", " # Generate two random numbers r₁, r₂ (r₁ used to sample from the Exponential distribution, r₂ to sample the compartment)\n", " \n", " r₁, r₂ = ####\n", "\n", " # Step 2\n", " # Compute propensity functions\n", " α = ####\n", " α₀ = ####\n", "\n", " # Step 3\n", " # Compute the time at which the next chemical reaction takes place as t+τ, where τ is sampled from an Exponential distribution with intensity α₀\n", " τ = (1/α₀)*log(1/r₁) \n", " t += τ\n", "\n", " # Step 4\n", " # Find j ∈ [1, 2, ..., K-1] for which one molecule diffuses to the right or j ∈ [2, 3, ..., K] for which one molecule diffuses to the left\n", " prob_right = sum(α[1:K-1]/α₀)\n", " if r₂ < prob_right\n", " index_j = find_j(1:K-1, α[1:K-1]/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j] = ####AA[t_idx, index_j]-1\n", " AA[t_idx+1, index_j+1] = ####AA[t_idx, index_j+1]+1\n", " else\n", " cte = α₀*prob_right\n", " index_j = find_j(2:K,([α[2]+cte; α[3:end]])/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j] = ####\n", " AA[t_idx+1, index_j-1] = ####\n", " end\n", "\n", " append!(tt, t)\n", "\n", " end \n", " return tt, AA\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "find_j = function(idx, ww, r)\n", " \"\"\"\n", " Inputs:\n", " idx - indexing vector. \n", " ww - vector of probabilities for each i ∈ idx.\n", " r - random value in the interval (0, 1)\n", " Output:\n", " a single sampled i, where i ∈ idx.\n", " \"\"\"\n", " ww_cumsum = cumsum(ww)\n", " return idx[findall(a -> a > r, ww_cumsum)[1]]\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A₀ = ####\n", "\n", "Nt = #### # maximum number of jumps between compartments.\n", "\n", "tt, AA = gillespie_alg_pure_diffusion(A₀, Nt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots\n", "anim = @animate for t_idx in 1:500:Nt\n", " bar(0.5*h/1000:h/1000:L, AA[t_idx,:], colour=:darkgrey,\n", " title=\"t: \".*\"$(round(tt[t_idx]/60, digits = 1))\".*\" minutes\",\n", " label = \"Compartment-based diffusion\",\n", " xlabel = \"x [mm]\",\n", " ylabel = \"number of molecules\",\n", " ylims = (0,500)) # 2:end since end = 1, periodic condition\n", "end\n", "gif(anim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Compartment-based approach to Reaction-Diffusion with Production and Degradation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathcal{A}_1 \\overset{d}{\\underset{d}\\rightleftarrows}\\mathcal{A}_2\\overset{d}{\\underset{d}\\rightleftarrows}\\mathcal{A}_3\\overset{d}{\\underset{d}\\rightleftarrows}\\dots\\overset{d}{\\underset{d}\\rightleftarrows}\\mathcal{A}_K,$ \n", "\n", "$\\mathcal{A}_i\\xrightarrow{\\kappa_1}\\emptyset\\quad$ for $i=1,2,\\dots,K,$\n", "\n", "$\\emptyset \\xrightarrow{\\kappa_2}\\mathcal{A}_i\\quad$ for $i=1,2,\\dots,K/5.$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "κ₁, κ₂ = ####)\n", "\n", "prepatternfactor = [0, 1/5];" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "propensities_degra_prod_diffusion = function(A, d, κ, prepatternfactor)\n", " κ₁, κ₂ = κ\n", " K = length(A)\n", " α = zeros(K)\n", " for i in 1:K\n", " α[i] = ####\n", " end\n", " return #####\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gillespie_alg_degra_prod_diff = function(A₀, Nt)\n", " \n", " \"\"\"\n", " Inputs:\n", " A₀ - a K-dimensional vector with the initial amount of molecules inside each compartment. \n", " Nₜ - the maximum number of jumps for one simulation.\n", " Outputs:\n", " tt - a partitioned time interval [0, T] with times sampled from an exponential distribution.\n", " AA - the history of counts, counts saved at each time t ∈ tt.\n", " \"\"\"\n", " \n", " K = length(A₀)\n", " \n", " t = 0.0\n", "\n", " tt = [t]\n", " AA = Integer.(zeros(Nt, K))\n", " AA[1,: ] = A₀\n", "\n", " κ = [κ₁, κ₂]\n", " \n", " for t_idx in 1:Nt-1\n", " \n", " # Step 1\n", " # Generate two random numbers r₁, r₂ (r₁ used to sample from the Exponential distribution, r₂ to sample the compartment index)\n", " r₁, r₂ = ####\n", " \n", " # Step 2\n", " # Compute propensity functions\n", " \n", " α = propensities_degra_prod_diffusion(AA[t_idx,:], d, κ, prepatternfactor)\n", " α₀ = ####\n", "\n", " # Step 3\n", " # Compute the time at which the next chemical reaction takes place as t+τ, where τ is sampled from an Exponential distribution with intensity α₀\n", " τ = #### \n", " t += τ\n", "\n", " # Step 4\n", " # Find j ∈ [1, 2, ..., K-1] for which one molecule diffuses to the right or j ∈ [2, 3, ..., K] for which one molecule diffuses to the left or\n", " # j ∈ [1, 2, ..., K] for which one molecule degrades or j ∈ [1, 2, ..., K/5] for which one molecule is produced.\n", " \n", " prob_right = sum(α[1:K-1]/α₀)\n", " prob_left = sum(α[2:K]/α₀)\n", " prob_degra = sum(α[(K+1):2K]/α₀)\n", " #prob_prod = sum(α[(2K+1):end]/α₀)\n", "\n", " #Diffusion to the right\n", " if r₂ < prob_right \n", " index_j = find_j(1:K-1, α[1:K-1]/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j] = ####\n", " AA[t_idx+1, index_j+1] = ####\n", " \n", " #Diffusion to the left\n", " elseif r₂ < prob_right + prob_left \n", " cte = α₀*prob_right\n", " index_j = find_j(2:K, ([α[2]+cte; α[3:K]])/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j] = ####\n", " AA[t_idx+1, index_j-1] = ####\n", " \n", " #Degradation of each volume\n", " elseif r₂ < prob_right + prob_left + prob_degra\n", " cte = α₀*(prob_right + prob_left)\n", " index_j = find_j(1:K, ([α[K+1]+cte; α[K+2:2K]])/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j] = ####\n", " \n", " #Production at the first Kfactor volumes\n", " else \n", " cte = α₀*(prob_right + prob_left + prob_degra)\n", " index_j = find_j(1:Integer(K*prepatternfactor[end]), [α[2K+1]+cte; α[(2K+2):end]]/α₀, r₂)\n", " AA[t_idx+1, :] = copy(AA[t_idx, :])\n", " AA[t_idx+1, index_j+1] += ####\n", " end\n", " append!(tt, t)\n", "\n", " end \n", " return tt, AA\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A₀ = ####\n", "\n", "Nt = 200000 # maximum number of jumps between compartments.\n", "\n", "tt, AA = gillespie_alg_degra_prod_diff(A₀, Nt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots\n", "anim = @animate for t_idx in 1:500:Nt\n", " bar(h/1000:h/1000:L, AA[t_idx,:], colour=:darkgrey,\n", " title=\"t: \".*\"$(round(tt[t_idx]/60, digits = 1))\".*\" minutes\",\n", " label = \"Compartment-based approach to reaction-diffusion\",\n", " xlabel = \"x [mm]\",\n", " ylabel = \"number of molecules\",\n", " ylims = (0,150)) # 2:end since end = 1, periodic condition\n", "end\n", "gif(anim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The French flag model " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bound₁ = 80\n", "bound₂ = 30;" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = h/1000:h/1000:L\n", "y = AA[end,:]\n", "\n", "plot(x[y .> bound₁], y[y .> bound₁]; \n", " seriestype=:bar, colour=:darkblue, label = \"N > 80\")\n", "plot!(x[bound₁ .>= y .> bound₂], y[bound₁ .>= y .> bound₂]; \n", " seriestype=:bar, colour=:white, label = \"80 >= N > 30\")\n", "plot!(x[y .<= bound₂], y[y .<= bound₂]; \n", " seriestype=:bar, colour=:red, label = \"30 >= N\")\n", "plot!(x, bound₁*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", "plot!(x, bound₂*ones(length(x)), linestyle=:dash, colour=:black, label=false)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. PDE solution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xx = range(.0, stop = L, length = K+1) # bin limits\n", "Δx = xx[2]-xx[1]\n", "xx_centre = xx[1:end-1] .+ Δx/2 # bin centroids\n", "\n", "A₀ = zeros(K)\n", "p = [κ₁, κ₂]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#import Pkg; Pkg.add([\"DiffEqOperators\", \"DifferentialEquations\"])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using DiffEqOperators\n", "# Second order approximation to the second derivative\n", "order = 2\n", "deriv = 2\n", "∂² = CenteredDifference(deriv, order, h/1000, K)\n", "nc = Neumann0BC(h/1000, 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "function discretized_pde!(du,u,p,t)\n", " du .= ####\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations\n", "\n", "T = tt[end]\n", "sol = solve(ODEProblem(discretized_pde!, A₀, (0, T)), Tsit5())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using Plots\n", "\n", "x = (h/1000)/2:h/1000:L\n", "y = AA[end,:]\n", "\n", "p₁ = plot(x[y .> bound₁], y[y .> bound₁],\n", " ylims=(0,150),\n", " seriestype=:bar, colour=:darkblue, label = \"N > 80\",\n", " xlabel = \"x [mm]\", ylabel = \"number of molecules\",\n", " title = \"compartment-based\")\n", " plot!(x[bound₁ .>= y .> bound₂], y[bound₁ .>= y .> bound₂]; \n", " seriestype=:bar, colour=:white, label = \"80 >= N > 30\")\n", " plot!(x[y .<= bound₂], y[y .<= bound₂]; \n", " seriestype=:bar, colour=:red, label = \"30 >= N\")\n", " plot!(x, bound₁*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", " plot!(x, bound₂*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", "\n", "\n", "y = sol(tt[end])*h^3\n", "\n", "yblue = y.*(y.>bound₁)\n", "ywhite = y.*(bound₂.>y.>=bound₁)\n", "yred = y.*(y.<=bound₂)\n", "\n", "p₂ = plot(x, y, fillrange=y.-yblue, \n", " ylims=(0,150),\n", " linetype=:steppre, colour = :black, fill=:darkblue, label = \"N > 80\",\n", " xlabel = \"x [mm]\", ylabel = \"number of molecules\",\n", " title = \"PDE solution\")\n", " plot!(x, y, fillrange=y.-ywhite, \n", " linetype=:steppre, colour = :black,fill=:white, label = \"80 >= N > 30\")\n", " plot!(x, y, fillrange=y.-yred, \n", " linetype=:steppre, colour = :black,　fill=:red,　label = \"30 >= N\")\n", " plot!(x, bound₁*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", " plot!(x, bound₂*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", "\n", "layout = @layout [p₁ p₂]\n", "Plots.plot(p₁, p₂; layout = layout, colorbar = false, size = (1400, 800), margin = 30Plots.px)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Velocity-jump process" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gillespie_alg_degra_prod_velocity = function(Δt, T)\n", "\n", " \n", " \"\"\"\n", " Inputs:\n", " Δt - length of time step. \n", " T - stopping time.\n", " Outputs:\n", " tt - a partitioned time interval [0, T] with times sampled from an exponential distribution.\n", " xx - the history of counts, counts saved at each time t ∈ tt.\n", " \"\"\"\n", " \n", " \n", " Nt = Integer(round(T/Δt))\n", "\n", " t = 0.0\n", "\n", " tt = [t]\n", "\n", " κ = [κ₁, κ₂]\n", " ensemble = [[x₀], [v₀]]\n", "\n", " for t_idx in 1:Nt-1\n", " \n", " xx = ensemble[1][]\n", " vv = ensemble[2][]\n", " \n", " # Generate a random numer r ∼ U(0,1) \n", " # Assume that a particle moves along the x-axis at a constant speed s and compute the position of the molecule at the future time t+Δt \n", " # Remember to apply reflective boundaries\n", " \n", " for i in 1:length(xx)\n", " xx[i] = abs.(xx[i]+vv[i]*Δt)\n", " if xx[i]>L\n", " xx[i] = ####\n", " end\n", " r = rand()\n", " # check if the particle turns in the time interval t+Δt\n", " if r<λ*Δt\n", " vv[i] = ####\n", " end\n", " end\n", "\n", " # for each molecule, check whether they have to be removed from the system\n", " id_to_remove = []\n", " for i in 1:length(xx)\n", " r₁ = rand()\n", " if r₁ < κ₁*Δt \n", " append!(id_to_remove, i)\n", " end\n", " end\n", "\n", " valid_id = setdiff(collect(1:length(xx)), id_to_remove)\n", " xx = xx[valid_id]\n", " vv = vv[valid_id]\n", "\n", " # check whether a new molecule should be introduced to the system\n", " r₂ = rand()\n", " if r₂ < (κ₂*(h)^2*L*1000/5)*Δt\n", " r₃ = rand()\n", " xx = [xx; ####] \n", " vv = [vv; v₀] \n", " end\n", "\n", " ensemble = [[xx], [vv]]\n", " t = t+Δt\n", " append!(tt, t)\n", " end \n", " return tt, ensemble[1]\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = 10^(-2)\n", "λ = s^2/(2D);\n", "\n", "Δt = 0.01\n", "x₀ = [0.0]\n", "v₀ = [s]\n", "\n", "tt, xx = gillespie_alg_degra_prod_velocity(Δt, T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#import Pkg; Pkg.add(\"StatsBase\")\n", "using StatsBase" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#x = h/1000:h/1000:L\n", "x = (h/1000)/2:h/1000:L\n", "y = copy(xx[end])\n", "\n", "\n", "firstEdge = 0.0\n", "lastEdge = L\n", "binSize = h/1000\n", "EdgeRange = (firstEdge:binSize:lastEdge)\n", "\n", "histogram_fit = fit(Histogram, y, h/1000:h/1000:L)\n", "weights = histogram_fit.weights\n", "edges = histogram_fit.edges[1]\n", "\n", "idx = findall(>(0), weights .> bound₁ )\n", "xblue = edges[idx]\n", "yblue = weights[idx]\n", "\n", "idx = findall(>(0), bound₁ .>= weights .> bound₂)\n", "xwhite = edges[idx]\n", "ywhite = weights[idx]\n", "\n", "idx = findall(>(0), weights .<= bound₂)\n", "xred = edges[idx]\n", "yred = weights[idx]\n", "\n", "p₃ = plot(xblue, yblue,\n", " ylims=(0,150),\n", " seriestype=:bar, colour=:darkblue, label = \"N > 80\",\n", " xlabel = \"x [mm]\", ylabel = \"number of molecules\",\n", " title = \"velocity-jump process\")\n", " plot!(xwhite, ywhite,\n", " seriestype=:bar, colour=:white, label = \"80 >= N > 30\")\n", " plot!(xred, yred,\n", " seriestype=:bar, colour=:red, label = \"30 >= N\")\n", " plot!(x, bound₁*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", " plot!(x, bound₂*ones(length(x)), linestyle=:dash, colour=:black, label=false)\n", "\n", "\n", "layout = @layout [p₁ p₂ p₃]\n", "ptotal = Plots.plot(p₁, p₂, p₃; layout = layout, colorbar = false, size = (2400, 800), m" ] } ], "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 }