Stochastic simulations

Stochastic simulations#

Gillespie Algorithm#

using Plots
using DataInterpolations: LinearInterpolation
using StatsBase: Weights, sample
using Statistics: mean
using Random
using DisplayAs: PNG    ## For faster rendering
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)
end

ssa_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
   0

Simulation time

tend = 10.0

10.0

Solve 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.443216 seconds (783.81 k allocations: 54.753 MiB, 98.59% compilation time)
  0.006636 seconds (16.33 k allocations: 14.929 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,
    xlabel="time", ylabel="# of molecules",
    title="SSA (direct method)", label=["A" "B"]) |> PNG
_images/929003d6e40f7a1677513004831b00ccae06fb4363d97a0c29caa35b8d38f204.png

Plot the solution by the first reaction method

plot(solfirst.t, solfirst.u,
    xlabel="time", ylabel="# of molecules",
    title="SSA (1st reaction method)", label=["A" "B"]) |> PNG
_images/309e984f5796524430fcf1189c0aa111fa16853c7e7db6816108766a59821189.png

Running 50 simulations

numRuns = 50

@time sols = map(1:numRuns) do _
    ssa_alg(model, u0, tend, parameters, stoich; method=:first)
end;

  0.480305 seconds (941.19 k allocations: 766.588 MiB, 54.23% gc time, 13.67% compilation time)

Average values and interpolation

ts = range(0, tend, 101)
a_avg(t) = mean(sols) do sol
    A = LinearInterpolation(sol.u[:, 1], sol.t)
    A(t)
end

b_avg(t) = mean(sols) do sol
    A = LinearInterpolation(sol.u[:, 2], sol.t)
    A(t)
end

b_avg (generic function with 1 method)

Plot the solution

fig1 = plot(xlabel="Time", ylabel="# of molecules", title="SSA (first method) ensemble")

for sol in sols
    plot!(fig1, sol.t, sol.u, linecolor=[:blue :red], linealpha=0.05, label=false)
end

fig1 |> PNG
_images/f36956c5b4f40c2de4163336d19980e070133de201c8c80644cc255ce7e189d2.png

Plot averages

plot!(fig1, a_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle=:solid, label="Average [A]") |> PNG
plot!(fig1, b_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle=:dash, label="Average [B]") |> PNG
fig1 |> PNG
_images/fa92faa8304c5550be4c4e879b15b819f9504a97489e7ec13f3d80dbc5e7eec4.png
Contents