Gillespie Algorithm¶
using Plots
using DataInterpolations: LinearInterpolation
using StatsBase: Weights, sample
using Statistics: mean
using Random
Random.seed!(2024)Random.TaskLocalRNG()Do-it-yourself¶
Stochastic chemical reaction: Gillespie Algorithm (direct and first reaction method) Adapted from: Chemical and Biomedical Engineering Calculations Using Python Ch.4-3
function ssa_alg(model, u0::AbstractVector, tend, p, stoich; tstart=zero(tend), method=:direct)
t = tstart ## Current time
ts = [t] ## Time points
u = copy(u0) ## Current state
us = copy(u') ## States over time
while t < tend
a = model(u, p, t) ## propensities
if method == :direct
dt = randexp() / sum(a) ## Time step for the direct method
du = sample(stoich, Weights(a)) ## Choose the stoichiometry for the next reaction
elseif method == :first
dts = randexp(length(a)) ./ a ## time scales of all reactions
i = argmin(dts) ## Choose the most recent reaction to occur
dt = dts[i]
du = stoich[i]
else
error("Method should be either :direct or :first")
end
u .+= du ## Update time
t += dt ## Update time
us = [us; u'] ## Append state
push!(ts, t) ## Append time point
end
return (t=ts, u=us)
endssa_alg (generic function with 1 method)Propensity model for this example reaction. Reaction of A <-> B with rate constants k1 & k2
model(u, p, t) = [p.k1 * u[1], p.k2 * u[2]]model (generic function with 1 method)parameters = (k1=1.0, k2=0.5)(k1 = 1.0, k2 = 0.5)Stoichiometry for each reaction
stoich = [[-1, 1], [1, -1]]2-element Vector{Vector{Int64}}:
[-1, 1]
[1, -1]Initial conditions (Usually discrete values)
u0 = [200, 0]2-element Vector{Int64}:
200
0Simulation time
tend = 10.010.0Solve the problem using both direct and first reaction method
@time soldirect = ssa_alg(model, u0, tend, parameters, stoich; method=:direct)
@time solfirst = ssa_alg(model, u0, tend, parameters, stoich; method=:first) 0.436968 seconds (681.02 k allocations: 51.816 MiB, 98.73% compilation time)
0.006292 seconds (16.34 k allocations: 14.937 MiB)
(t = [0.0, 0.006874672879606227, 0.0070050915014185, 0.01329115182676454, 0.015335530116774018, 0.01942319769041779, 0.021346508038833215, 0.02345030020299408, 0.025497753122539372, 0.02734483697039726 … 9.937773197389458, 9.938360255993572, 9.964450740930666, 9.96564917153764, 9.975440668199402, 9.978594329657025, 9.982318425245289, 9.987213207313044, 9.995833233887383, 10.024199143896238], u = [200 0; 199 1; … ; 75 125; 74 126])Plot the solution from the direct method
plot(soldirect.t, soldirect.u, labels=["A" "B"], xlabel="time", ylabel="# of molecules", title="SSA (direct method)")
Plot the solution by the first reaction method
plot(solfirst.t, solfirst.u, labels=["A" "B"], xlabel="time", ylabel="# of molecules", title="SSA (first reaction method)")
Running 50 simulations serially
TODO: parallelize this loop using Threads.@threads or Distributed for better performance
@time sols = map(1:50) do _
ssa_alg(model, u0, tend, parameters, stoich; method=:direct)
end; 0.302136 seconds (744.82 k allocations: 752.091 MiB, 31.13% gc time, 22.88% compilation time)
Average values and interpolation
ts = range(0, tend, 101)
a_interp = [LinearInterpolation(sol.u[:, 1], sol.t; cache_parameters=true) for sol in sols]
b_interp = [LinearInterpolation(sol.u[:, 2], sol.t; cache_parameters=true) for sol in sols]
a_avg(t) = mean(i-> i(t), a_interp)
b_avg(t) = mean(i-> i(t), b_interp)b_avg (generic function with 1 method)Plot the solution
plot(ts, a_avg.(ts), label="Average [A]", line=(:black, :solid))
plot!(ts, b_avg.(ts), label="Average [B]", line=(:black, :dash))
for sol in sols
plot!(sol.t, sol.u[:, 1], line = (:blue, 0.05), label=false)
plot!(sol.t, sol.u[:, 2], line = (:red, 0.05), label=false)
end
plot!(xlabel="time", ylabel="# of molecules", title="SSA (first method) ensemble")
This notebook was generated using Literate.jl.