Simulate the Community Dynamics

Previous sections tackled how to create a model representing the desired ecological community. We now explain how to simulate the dynamics of this community. In short, we provide a function simulate that takes a model and a time interval as input and returns the temporal trajectory of the community. This function uses the DifferentialEquations package to solve the system of ordinary differential equations.

Basic Usage

Let's first illustrate how to simulate a simple community of three species.

using EcologicalNetworksDynamics, Plots

foodweb = Foodweb([3 => 2, 2 => 1])
m = default_model(foodweb)
B0 = rand(3) # Vector of initial biomasses.
t = 1_000
sol = simulate(m, B0, t)

retcode: Success
Interpolation: 3rd order Hermite
t: 37-element Vector{Float64}:
    0.0
    0.0893046973306578
    0.2543704634170648
    0.4476808040190955
    0.7042126611131245
    1.0676040579400243
    1.4702212250092703
    2.0118295421563563
    2.6446883425854306
    3.4216185585320815
    ⋮
   62.559735580999835
   71.76384003148556
   84.684292420364
  104.61272123657251
  131.22029689272821
  176.4814202424337
  276.2957519171365
  608.6933300828676
 1000.0
u: 37-element Vector{Vector{Float64}}:
 [0.4118674060453822, 0.6756183990877278, 0.12941553786801596]
 [0.36155935086621527, 0.6898338370808441, 0.13433746541862487]
 [0.29082195770590513, 0.703376194332363, 0.14411474947493375]
 [0.2362938711828531, 0.702604703103685, 0.1565709976746376]
 [0.1933498872912498, 0.6836733226857434, 0.17455486060402628]
 [0.16255064384872228, 0.6380610846414618, 0.202340328476311]
 [0.14823091416187184, 0.5759580771496242, 0.23523155681755573]
 [0.14494722427585843, 0.4878561277676562, 0.2800314183904344]
 [0.15432081074357984, 0.39263174530766715, 0.3271974873966053]
 [0.17903534213015768, 0.3012586878446336, 0.3693261586090291]
 ⋮
 [0.41372697591971, 0.18991973008133906, 0.6413924482499089]
 [0.41503658077116484, 0.18953805814992727, 0.6454071940707571]
 [0.4158659192159939, 0.18921037738806587, 0.6483310838636664]
 [0.4163561894380674, 0.18904360436960285, 0.6499707380502244]
 [0.4165018070715581, 0.18899020171135178, 0.6504605824459283]
 [0.41652552230295314, 0.18898179562890235, 0.6505424873271627]
 [0.41652409796000167, 0.1889823264028421, 0.650537414781329]
 [0.41652438089029553, 0.18898222316008392, 0.6505383797589629]
 [0.4165243384772354, 0.18898223831677236, 0.6505382364339435]

We can access the solution of the simulation with the output of the simulate function. We list below some useful properties of the solution:

sol.t # Time steps.
sol.u # Biomasses at each time step.
sol.u[1] # Biomasses of the first time step.
sol.u[end] # Biomasses of the last time step.

3-element Vector{Float64}:
 0.4165243384772354
 0.18898223831677236
 0.6505382364339435

The solution can be plotted with the plot function from the Plots package.

plot(sol)

The duration of the simulation can be changed with, for instance to reduce the simulation time to 100 time units:

smaller_t = 100
sol = simulate(m, B0, smaller_t)
sol.t[end] # The last time step.

100.0

Callbacks

We will now go through some advanced features of the simulate function. First, the callback keyword argument allows specifying a function that will be called at each time step of the simulation. We provide a built-in callback extinction_callback which extinguishes the species whose biomass falls below a given threshold. This threshold is set by default to 1e-12, but can be changed.

foodweb = Foodweb([3 => 1, 2 => 1]) # Two predators feeding on one prey.
m = default_model(foodweb, Metabolism([0, 0.1, 100.0])) # Predator (3) has a too high metabolic rate to survive.
sol = simulate(m, [1, 1, 1], 100_000; callback = nothing) # No callback.
sol[end]

3-element Vector{Float64}:
 0.18898223650461363
 0.6897057785564755
 5.474600390482506e-175

callback = extinction_callback(m, 1e-6)
sol = simulate(m, [1, 1, 1], 100_000; callback) # High extinction threshold.
sol[end]

3-element Vector{Float64}:
 0.18894484615604795
 0.6897503358890235
 0.0

callback = extinction_callback(m, 1e-12)
sol = simulate(m, [1, 1, 1], 100_000; callback) # Low extinction threshold.
sol[end]

3-element Vector{Float64}:
 0.18898223650461363
 0.6897057785564753
 0.0

Other callback functions are available in the DiffEqCallbacks package, and can be used in the same way.

Choose a Specific Solver

Depending on your needs, you may want to choose a specific solver for the simulation. As we use the solve function of the DifferentialEquations package, we can pass any solver available in this package (see the list of available solvers). Indeed, we allow the user to pass any keyword argument of the solve function to the simulate function.

import DifferentialEquations: Tsit5

sol = simulate(m, [1, 1, 1], 1_000; alg = Tsit5())
sol.alg

Tsit5(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),)