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: 229-element Vector{Float64}:
    0.0
    0.11103570399195549
    0.478256922222615
    0.9297265883848238
    1.4009022415569263
    1.983534100723224
    2.6478984664670446
    3.459523296085478
    4.42232538304872
    5.493776491618197
    ⋮
  963.7089895796448
  968.2315680978994
  972.8714111765078
  977.6641628818189
  982.3027276479369
  986.8668156114379
  991.5405602285765
  996.3015982155536
 1000.0
u: 229-element Vector{Vector{Float64}}:
 [0.06940661404940796, 0.4569924776363492, 0.7406976984851047]
 [0.07373305415181056, 0.41614045920315246, 0.76103659719395]
 [0.09062510864979705, 0.30644119379523393, 0.8099083720981389]
 [0.11755193862348964, 0.21888917921104298, 0.8354919294515492]
 [0.15360576600979003, 0.16513664104208187, 0.8358228156950895]
 [0.20930836196207372, 0.12879923731002949, 0.8177407798934304]
 [0.2845364587805067, 0.11011022754715397, 0.7871599402793337]
 [0.38054518780334723, 0.1068874425330645, 0.7470298765227934]
 [0.4710093618654694, 0.12128569717401096, 0.7045491571596437]
 [0.5121579919234602, 0.15068353860344594, 0.6714209741800676]
 ⋮
 [0.41662673350423873, 0.18908897661444526, 0.6504426195176103]
 [0.4167132091679696, 0.1889982409564515, 0.6504598973040386]
 [0.41670622342283786, 0.18891846157295006, 0.6505070319901195]
 [0.4166092462649183, 0.18885075445575858, 0.6505795270118627]
 [0.41641119490528283, 0.18888466334983112, 0.650631636097843]
 [0.41632813776032646, 0.1889734844536669, 0.6506140074585344]
 [0.41633493562705826, 0.18904965751791714, 0.6505720148712764]
 [0.4164242564192656, 0.18911319419290135, 0.6505055495970493]
 [0.4165239971636296, 0.18898240776460867, 0.6505392873125507]

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.4165239971636296
 0.18898240776460867
 0.6505392873125507

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. Moreover, species extinctions can be printed to the console with the verbose keyword argument.

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.6897057785564756
 5.474642132648367e-175

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

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

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