Run the World

run!(w, alg, tend, b, d; dt_saving=nothing, cb=(names = String[],agg =nothing))

Run w with algorithm alg, until tend is reached. User needs to provide b the birth function, which takes as arguments X,t, and provide d the death function, with arguments X,Y,t. The run stops if the number of agents reachesp["NMax"], where p is the parameter dictionary in the world w.

Returns a Simulation object. It is a container for snapshots of the world at every dt_saving time steps. It renders post processing easier, through dedicated methods to obtain time series of quantities.

Keyword arguments

• if dt_saving specified, world is saved every time steps.

If dt_saving not specified, sim contains an array of two elements, first corresponding to initial conditions and last corresponding to world in the last time step.

• if t_saving_cb specified, callbacks are computed at each steps time specified in the array.

This functionality is as of now only compatible with dt_saving not specified.

• cb correspond to callbacks function. Callbacks can be used

to extract properties of the world at each dt_saving time steps of your simulation.

Constructing the callbacks

A callback has to be of the form

cb = (names = String[], agg = Function[])

It is a tuple, with first value corresponding to the names of the aggregate properties of the world. The second correspond to the aggregation functions.

We provide here an example on how to extract the $\gamma$ diversity of a simulation biological population. $\gamma$` diversity can be calculated as the variance of the trait distribution of the population. Here is how we write the function

cb = (names = ["gamma_div"], agg = Function[w -> var((get_x(w,1)))])

Example

myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
b(X) = gaussian(X[1],0.,sigma_K)
d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
tend = 1.5
p = Dict{String,Any}();@pack! p = D,mu,NMax

myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents, myspace, p, 0.)
w1 = copy(w0)
sim = run!(w1,Gillepsie(),tend,b,d,cb=cb,dt_saving = .1)

Accessing the callbacks

You can easily access the properties, using sim as you would for a usual Dictionary.

using Plots
plot(get_tspan(sim),sim["gamma_div"])

struct Gillepsie <: EvoId.AbstractAlg

Gillespie algorithm.

Denote by $b_i$ and $d_i$ the birth and death rates of agent $i$. The total rate is given by the sum of all individual rates $R(t) = \left[ \sum_i b_i(t) + d_i(t) \right]$. A particular event, birth or death, is chosen at random with a probability equal to the rate of this event divided by the total rate $R$.

The original article by Gillsepie:

A general method for numerically simulating the stochastic time evolution of coupled chemical reactions

struct CFM <: EvoId.AbstractAlg

Champagnat Ferriere Méléard algorithm described in Champagnat and Ferriere founding article. This algorithm is similar to Gillepsie, excepts that it runs faster for higher number of individuals.

Indeed, at every time step, only the fitness of the individual picked at random is evaluated. In order to use it, you need to feed to the dictionary parameters p a constant Cbar<:Real that is the upperbound of the maximum of the sum of the birth and death rates (cf article).

Example

using EvoId,UnPack,Plots
myspace = (RealSpace{1,Float64}(),)
σ_b = .9;
σ_d = .7;
K0 = 1000
b(X,t) = 1.
d(X,Y,t) = gaussian(X[1],Y[1],σ_d)/K0 / gaussian(X[1],0.,σ_b)
Cbar = b([0],0.) + d([0],[0],0.)
D = (1e-2,)
mu = [.1]
NMax = 2000
tend = 1500
p = Dict{String,Any}();@pack! p = D,mu,NMax,Cbar
myagents = [Agent(myspace,(1e-2 * randn(),)) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@time sim = run!(w0,CFM(),tend, b, d, dt_saving = 4)

Algorithm described in article DOI : 10.1016/j.tpb.2005.10.004

struct WF <: EvoId.AbstractAlg

Wright Fisher algorithm.

In the Wright Fisher process the number of agents is constant through time. It is helpful to visualize it through marbles in jar At each time step, $N$ agents are picked up from previous generation to reproduce. Their number of offspring is proportional to their fitness, calculated as usual with birth and death rates. It takes thus only one time step to go trough one generation. Thus it is more suitable for numerical simulations than CFM or Gillespie.