Examples

Logistic Growth with noise

This is a simple example of a logistic growth model with noise, where reaction is generated by @cfunc and @ufunc.

The determine differential equation of this model is:

\[\dot{x} = rx \left(1 - \frac{x}{K} \right)\]

using RecordedArrays
using EvolutionaryModelingTools
using Random
using Plots

# parameters
const r = 0.5 # growth rate
const d = 0.1 # mortality rate
const K = 100 # carrying capacity

# clock
c = ContinuousClock(100.0) # define a clock, the population will growth for 100 time unit

# state
x = recorded(DynamicEntry, c, 10)  # define a scalar to record population size

# reactions

@cfunc growth_c(r, x) = r * x
@ufunc growth_u!(ind, x) = x[ind] += 1
growth = Reaction(growth_c, growth_u!)

@cfunc death_c(d, x) = d * x
@ufunc death_u!(ind, x) = x[ind] -= 1
death = Reaction(death_c, death_u!)

@cfunc comp_c(r, K, x) = r * x * x / K
@ufunc comp_u!(ind::CartesianIndex{0}, x) = x[ind] -=1
comp = Reaction(comp_c, comp_u!)

# build the model
ps = (; r, d, K, x) # define the parameters for the model
rs = (growth, death, comp) # define the reactions for the model

# simulate
ps_new, t, state = gillespie(MersenneTwister(1), c, ps, rs)

# plot
timeseries(ps_new.x;
    grid=false, frame=:box, legend=false,
    title="Population Dynamics", xlabel="Time", ylabel="Population Size",
)

Evolutionary competitive dynamics

A simple example model to simulate the evolution of a competition, where the reactions is defined by @reaction.

Unlike the previous example, this model consists of various of species introduced by mutation. And their competition is defined by a competition coefficient matrix $M$, where the element $M_{ij}$ is the competition force from species $j$ to $i$, and the element $M_{ii}$ is the competition force from species $i$ to itself.

The determine differential equation of this model is:

\[\dot{x_i} = rx \left(1 - \frac{\sum_j x_j M_{ij}}{K} \right)\]

using RecordedArrays
using EvolutionaryModelingTools
using Random
using Plots


# parameters
const r = 0.5 # growth rate
const d = 0.1 # mortality rate
const K = 100 # carrying capacity
const μ = .01 # mutation rate

# clock
c = ContinuousClock(200.0) # define a clock, the population will growth for 10 time unit

# state
x = recorded(DynamicEntry, c, [10]) # define a vector to record population size
m = recorded(StaticEntry, c, ones(1, 1)) # define a matrix to record competition coefficient

# reactions
@reaction growth begin
    r * (1 - μ) * x
    x[ind] += 1
end

@reaction mutation begin
    @. r * μ * x
    begin
        push!(x, 1) # add a new species
        n = length(x)
        resize!(m, (n, n)) # resize the competition matrix
        # new species compete with all other species and itself
        m[:, n] = rand(rng, n)
        m[n, 1:n-1] = rand(rng, n-1)
    end
end

@reaction competition begin
    begin
        cc = r / K # baseline competition coefficient
        @. cc * x * m * x'
    end
    begin
        i = ind[1]
        x[i] -= 1
        if x[i] == 0 # check extinction
            deleteat!(x, i)
            resize!(m, (not(i), not(i)))
        end
    end
end

# hook function to check if all species is extinct
@ufunc check_extinction(x) = isempty(x) ? :extinct : :finnish

# build the model
ps = (; r, d, K, μ, x, m) # define the parameters for the model
rs = (growth, mutation, competition) # define the reactions for the model

# simulate
ps_new, t, state = gillespie(check_extinction, MersenneTwister(1), c, ps, rs)

# plot
timeseries(ps_new.x;
    grid=false, frame=:box, legend=false,
    title="Population Dynamics", xlabel="Time", ylabel="Population Size",
)

SIR model with vital dynamics and logistic population

A simple SIR model with vital dynamics and logistic population, where the host population follows the logistic growth model.

This model may contains more reactions than the previous example, which can be generated easily by @reaction_eq.

The deterministic differential equation of this model is:

\[\begin{aligned} \dot{S} &= r (S + I + R) - S(d + c (S + I + R) + \beta I) \\ \dot{I} &= I (\beta S - d - c (S + I + R) - \nu) \\ \dot{R} &= R (\nu - d - c (S + I + R)) \end{aligned}\]

using RecordedArrays
using EvolutionaryModelingTools
using Random
using Plots

const T = 100.0
const β = 0.005
const ν = 0.01
const α = 0.2
const r = 0.5
const d = 0.1
const c = 0.001
const S = 100
const I = 10
const R = 0

# epidemic dynamics
@reaction_eq infection β S + I --> I + I
@reaction_eq recovery ν I --> R
@reaction_eq virulence α I --> 0 # death of infection host caused by virus
# demography dynamics, generate reactions with @eval
for sym in (:S, :I, :R)
    r_name = Symbol(:growth_, sym)
    d_name = Symbol(:death_, sym)
    @eval @reaction_eq $r_name r $sym --> S + $sym # growth
    @eval @reaction_eq $d_name d $sym --> 0 # death
    for sym′ in (:S, :I, :R)
        c_name = Symbol(:competition_, sym, sym′)
        @eval @reaction_eq $c_name c $sym + $sym′ --> $sym′ # competition
    end
end

const REACTIONS = (
    infection, recovery, virulence,
    growth_S, growth_I, growth_R,
    death_S, death_I, death_R,
    competition_SS, competition_SI, competition_SR,
    competition_IS, competition_II, competition_IR,
    competition_RS, competition_RI, competition_RR
)

clock = ContinuousClock(T)
ps = (;β, ν, α, r, d, c,
    S=recorded(DynamicEntry, clock, S),
    I=recorded(DynamicEntry, clock, I),
    R=recorded(DynamicEntry, clock, R),
)
ps_new, t, term = gillespie(MersenneTwister(1), clock, ps, REACTIONS)

plt = plot(grid=false, frame=:box,
    title="Population Dynamics", xlabel="Time", ylabel="Population Size"
)
timeseries!(plt, ps_new.S; label="S")
timeseries!(plt, ps_new.I; label="I")
timeseries!(plt, ps_new.R; label="R")
plt