Stochastic simulations

Stochastic simulations#

Gillespie Algorithm#

using StatsBase         ## Weights() and sample()
using Plots
using Interpolations
using Statistics        ## mean()
using Random            ## randexp()
using DisplayAs: PNG
Random.seed!(2024)

Random.TaskLocalRNG()

Stochastic chemical reaction: Gillespie Algorithm (direct method) Adapted from: Chemical and Biomedical Enginnering Calculations Using Python Ch.4-3

function ssa_direct(model, u0::AbstractArray, tend, p, stoich; tstart=zero(tend))
    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
        dt = randexp() / sum(a)          ## Time step for the direct method
        du = sample(stoich, Weights(a))  ## Choose the stoichiometry for the next reaction
        u .+= du   ## Update time
        t += dt    ## Update time
        us = [us u]  ## Append state
        push!(ts, t) ## Append time point
    end
    us = collect(us') ## Transpose to make columns as state variables, rows as observations
    return (t = ts, u = us)
end

ssa_direct (generic function with 1 method)

Stochastic chemical reaction: Gillespie Algorithm (first reaction method) Adapted from: Chemical and Biomedical Enginnering Calculations Using Python Ch.4-3

function ssa_first(model, u0, tend, p, stoich; tstart=zero(tend))
    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
        dts = randexp(length(a)) ./ a   ## time scales of all reactions
        i = argmin(dts) ## Choose which reaction to occur
        dt = dts[i]
        du = stoich[i]
        u .+= du  ## Update state
        t += dt   ## Update time
        us = [us u]  ## Append state variable to record
        push!(ts, t) ## Append time point to record
    end
    us = collect(us') ## Make column as variables, rows as observations
    return (t = ts, u = us)
end

ssa_first (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_direct(model, u0, tend, parameters, stoich)
@time solfirst = ssa_first(model, u0, tend, parameters, stoich)

  0.436588 seconds (518.98 k allocations: 42.560 MiB, 98.96% compilation time)
  0.257400 seconds (368.71 k allocations: 32.579 MiB, 20.32% gc time, 99.33% compilation time)

(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

sols = map(1:numRuns) do _
    ssa_direct(model, u0, tend, parameters, stoich)
end;

A interpolation function for each solution

itpsA = map(sols) do sol
    linear_interpolation(sol.t, sol.u[:, 1], extrapolation_bc = Line())
end;

itpsB = map(sols) do sol
    linear_interpolation(sol.t, sol.u[:, 2], extrapolation_bc = Line())
end;

Calculate average A and B levels in the ensemble.

a_avg(t) = mean(i->i(t), itpsA)
b_avg(t) = mean(i->i(t), itpsB)

b_avg (generic function with 1 method)

Plot the soluton

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

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

fig1 |> PNG
_images/6c6bee57b1551e3bfe7df3f8ef097f0e95ead4adcf5cc09da9fa7d566b2ec7ae.png

Plot averages

plot!(fig1, a_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle = :solid, label="Avarage [A]") |> PNG
plot!(fig1, b_avg, 0.0, tend, linecolor=:black, linewidth=3, linestyle = :dash, label="Avarage [B]") |> PNG

fig1 |> PNG
_images/a1ad2769f469499cc11f09293c1523331ace8594d3f97fba7521cc2f020eb765.png