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",
)