EvolutionaryModelingTools
A simple package provides an easy way to build evolutionary biology models and simulate them by Gillespie's direct method algorithm.
Why?
DifferentialEquations.jl
is a brilliant suite solving deferential equations,
including simulating jump process with Gillespie algorithm, and it is a very good choice.
However, it is not suitable for solving differential equations
with "variable length" state,
which is the main reason why I created this package.
For example, in a SIR model, the host population is composed of "Susceptible", "Infected" and "Recovered".
In a normal case, there are only one type of "host" and one type "virus" in the system.
Thus the state of host population can be represent as a vector u = [S, I, R]
,
In a complex case, the host population can be composed of many types of "host" and infected by many types of "virus".
In a system with n
types of "hosts" and m
types of "viruses",
the state of host population can also be represented as a vector
by concatenating the state of each component of host u = vcat(S, vec(I), vec(R))
,
where S
is a vector of length n
and I
, R
are matrixes of size n × m
.
However, in evolutionary biology,
the "mutation" and "extinction" will change the types of hosts and viruses,
which means the n
and m
changes during the evolution,
and the length of the state vector u
will also change.
How to use?
The package Catalyst.jl
provides a simple way
to build biochemical reaction for DifferentialEquations.jl
.
Similarly, this package provides a macro @reaction_eq
which generate reaction(s),
with given equation.
For example, the infect reaction of above SIR model with multi-type of viruses can be defined as:
@reaction_eq infect β S + I[i] --> 2I[i]
where i
donates the virus type. The equation S + I[i] --> 2I[i]
means
the an host infected by the virus i
infect a susceptible host with rate β
,
then convert the the susceptible host to a infectious host.
This expression will not only generate one reaction but a group of reactions trough the index i
.
However, the mutation and extinction can not be defined easily by the macro @reaction_eq
currently.
Thus the alternative macro @reaction
provides a low-level way to build reaction(s)
@reaction mutation begin
@quickloop μ * I[i]
begin
i = ind[1]
I[i] -= 1 # the host individual is converted to another type
push!(I, 1) # add a new type of infectious host to the system
push!(α, randn() + α[i]) # the virulence of the new host type is generated randomly with mean `α[i]`
end
end
This expression defined mutation of virus which contains two parts:
- The
@quickloop μ * I[i]
defines how to calculate the rate of mutation, - The
begin ... end
block defines what happens when the host is mutated, whereind
is a preserved variable which is used to store the index of the mutated host.
Once you have defined all reactions, put them together as a tuple:
reactions = (infect, mutation, ...)
and define the initial state and parameters of the system as a named tuple:
params = (β = 0.1, μ = 0.001, ..., S = scalar(100), I = [1], R = [0], α = [0.5])
where the scalar
will create a type similar to Number
,
but it can be update in-place like Ref
like S[] += 1
.
Once reaction and parameters is defined, you can use gillespie
to simulate the system:
max_time = 100 # the maximum time of simulation
params′, t, term = gillespie(max_time, params, reactions)
where the gillespie
function returns an tuple t, ps′, term
where t
is the time when the simulation ends, params′
is an updated params
,
and term
is a flag indicating whether the simulation is finished after the maximum time
or break with given flag.
Note: Changes of the state will not be recorded by default,
but you can use my another package RecordedArrays
to record them, like this:
using RecordedArrays
c = DiscreteClock(max_time) # clock store information of max_time
S = recorded(DynamicEntry, c, 100) # create a recorded number as S with the clock c
I = recorded(DynamicEntry, c, [1]) # create a recorded vector as I with the clock c
R = recorded(DynamicEntry, c, [0]) # create a recorded vector as R with the clock c
α = recorded(StaticEntry, c, [0.5]) # create a recorded vector as α with the clock c
params = (; β = 0.1, μ = 0.001, ..., S, I, R, α) # create new params with recorded S, I, R, α
gillespie(c, params, reactions) # run the simulation with the clock and new params
More information about RecordedArrays
, see its
documentation.