This package for Julia provides a tool for easy specification and simulation of stochastic agent-based models where internal continuous-time agent state dynamics explicitly modelled and affect the population level interactions between agents.
- Individual agent dynamics defined as ODEs, SDEs or jump-diffusion processes using ModelingToolkit.jl and Catalyst.jl.
- High-level symbolic specification of agent interaction rules.
- Stochastic sampling for rate-based interactions between agents.
The package can be installed through the Julia package manager.
pkg> add https://github.com/pihop/AgentBasedModeling.jl.gitExample providing cell population dynamics. Start by defining the continuous time dynamics of the agents.
using AgentBasedModeling
using ModelingToolkit
using Catalyst
using Distributions
@variables t τ(t) s(t)
@species C(τ)
@parameters a μ cv L Cs Cτ
der = Differential(t)
# Age increases linearly in time t while size grows exponentially.
@named CellDynamics = ODESystem([der(τ) ~ 1.0, der(s) ~ a*s], t)
Cell = AgentDynamics(CellDynamics, ())Define the division interaction dividing the cell and creating two daughter cells with half the size each.
# Hazard of gamma distribution.
gammahaz(μ,cv,x) = exp(logpdf(Gamma(1/cv, μ*cv), x) - logccdf(Gamma(1/μ, μ*cv), x))
# Adder rule for cell division.
γdiv(μ, cv, s, τ, a) = gammahaz(μ, cv, s - s*exp(-a*τ))
@register_symbolic γdiv(μ, cv, s, τ, a)
divide = @interaction begin
@channel γdiv($μ, $cv, $Cs, $Cτ, a), $C --> 2*$C
@sampler ExtrandeMethod(1/($μ * $cv), $L)
@connections (($Cτ => $τ, $Cs, => $s),)
@transition (($τ => 0.0, $s => 0.5*$Cs), ($τ => 0.0, $s => 0.5*$Cs))
endCombine the agent dynamics and interactions into a model.
cell_population_model = AgentsModel([divide,], Dict(C => Cell, ))Give some initial values and specify the model parameters. Population is initialised with a single cell with age 0 and size 0.1.
using OrdinaryDiffEq
init_pop = repeat([C => (τ => 0.0, s => 0.1)], 1)
# SimulationParameters expects vector of parameter value pairs, simulation timespan,
# saving timestep and a solver for the agent dynamics.
ps = [a => 0.5, L => 1.0, μ => 1.0, cv => 0.1]
tspan = (0, 10.0)
Δt = 1.0
simulation_params = SimulationParameters(ps, tspan, Δt;
snapshot=[PopulationSnapshot(C), ])Simulate the system and display the population size snapshots.
res = simulate(cell_population_model, init_pop, simulation_params)
res.snapshot[:C]The package is flexible --- see the following for more complex examples!