Gillespie Algorithm

"""

Stochastic chemical reaction: Gillespie Algorithm

Adapted from: Chemical and Biomedical Enginnering Calculations Using Python Ch.4-3

Reaction of A <-> B with rate constants k1 & k2

"""

from random import choices, expovariate

import numpy as np

import matplotlib.pyplot as plt

%matplotlib inline

class Gillespie():

    def __init__(self, propensityFuncs, actionFuncs, parameters=None):

        self.propensityFuncs = propensityFuncs

        self.actionFuncs = actionFuncs

        self.parameters = parameters

    def run(self, u0, tend, tstart=0):

        # Setup

        t = tstart

        p = self.parameters

        u = np.asarray(u0)

        us = [u.copy()]

        ts = [t]

        while t < tend:

            # propensities of reactions

            ps = [f(u, p, t) for f in self.propensityFuncs]

            pTotal = sum(ps)

            dt = expovariate(pTotal)

            # Choose an action by the weight of each propensities, and update the state variable

            act = choices(actionFuncs, weights=ps)[0]

            u = np.asarray(act(u, p, t))

            t += dt

            us.append(u.copy())

            ts.append(t)

        return np.array(ts), np.array(us)

parameters = {"k1": 1.0, "k2": 0.5}

propensityFuncs = (lambda u, p, t: p["k1"] * u[0], 

                   lambda u, p, t: p["k2"] * u[1])

actionFuncs = (lambda u, p, t: [u[0] - 1, u[1] + 1], 

               lambda u, p, t: [u[0] + 1, u[1] - 1])

ssa = Gillespie(propensityFuncs = propensityFuncs,

                actionFuncs = actionFuncs,

                parameters = parameters)

ts, us = ssa.run([175, 25], 10.0)

fig, ax = plt.subplots()

ax.plot(ts, us[:, 0], label="A")

ax.plot(ts, us[:, 1], label="B")

ax.set_xlabel('time')

ax.set_ylabel('number of molecules')

ax.legend(loc='best')

# Packages

using Pkg

pkg"up; add StatsBase Plots DifferentialEquations"

# Coded in the same manner as the Python counterpart

struct Gillespie

  propensityFuncs

  actionFuncs

  parameters

end

function (g::Gillespie)(u0, tend, tstart=zero(tend))

  t = tstart

  ts = [t] # Records time

  u = copy(u0)

  us = copy(u)

end

#=

Stochastic chemical reaction: Gillespie Algorithm

Adapted from: Chemical and Biomedical Enginnering Calculations Using Python Ch.4-3

Reaction of A <-> B with rate constants k1 & k2

=#

using StatsBase, Random

# Coded in the same manner as the Python counterpart

struct Gillespie

    propensityFuncs

    actionFuncs

    parameters

end

function (g::Gillespie)(u0, tend, tstart=zero(tend))

    t = tstart

    ts = [t]

    u = copy(u0)

    us = copy(transpose(u))

    p = g.parameters

    while t < tend

        # propensities of reactions

        ps = [f(u, p, t) for f in g.propensityFuncs]

        pTotal = sum(ps)

        dt = randexp() / pTotal

        f = sample(g.actionFuncs, Weights(ps))

        u = f(u, p, t)

        t += dt

        us = vcat(us, transpose(u))

        push!(ts, t) # Record t

    end

    return (ts = ts, us = us)

end

parameters = (k1=1.0, k2=0.5)

propensityFuncs = [(u, p, t) -> p.k1 * u[1], (u, p, t) -> p.k2 * u[2]]

actionFuncs = [(u, p, t) -> u .+ [-1, 1] , (u, p, t) -> u .+ [1, -1]]

ssa = Gillespie(propensityFuncs, actionFuncs, parameters)

ts, us = ssa([175, 25], 10.0)

using Plots

plot(ts, us, xlabel="time", ylabel="# of molecules", title = "SSA", label=["A" "B"])