Solving differential equations for ecological problems.

Up until a few months ago, when people asked what the advantages of Julia were, I usually mentioned its speed, maintainability, and how easy it was to run your code in a distributed way. Now, I usually add that Chris Rackauckas created the best differential equations package available (DifferentialEquations.jl). Check out the comparison with other packages. I thought it would be useful to give a brief overview of one very useful feature it has when solving ecological problems.

Update (Feb. 2018): There have been syntax changes in the package, discussed [here][updates].

[updates]: {% link _posts/2018-02-12-differential-equations-updated.md %}

Let’s use a model of the Allee effect as an example. The Allee effect describes a situation in which the growth rate becomes positive only when a minimum population density is reached . We can write it as

$$\frac{1}{N}\frac{d}{dt}N = r\left(\frac{N}{A}-1\right)\left(1-\frac{N}{K}\right)$$

This population will have a negative growth rate when $0 < N < A$, positive when $A < N < K$, and reaches an equilibrium for $N^\star = {0, A, K}$. We can write this equation as:

function dndt(N, A, K, r)
    return N*r*(N/A-1)*(1-N/K)
end

In the DifferentialEquations.jl package [@RackNie17], things proceed a little bit differently. The problem to integrate must be a function $f(t,u)$, where $t$ is the time, and $u$ is the state of the problem (here, population size).

using DifferentialEquations
using Plots # We'll need this

# Definition of the parameters
r, A, K = 1.2, 0.3, 2.0

# Function giving the derivative
f(t, n) = n*r*(n/A-1)*(1-n/K)

# T0 and final time
t = (0.0,3.0)

# Definition of the problem
n0 = 1.0
p = ODEProblem(f, 1.0, t)
solution = solve(p)

The package also interacts really well with Plots.jl, so that we can write

plot(solution)

And get the following output:

{: .center-image }

But what about multiple parameters? It is easy to wrap the above code within a function:

function allee_effect(n0;
 A::Float64 = 0.3,
 K::Float64 = 1.2,
 r::Float64 = 1.0,
 timesteps::Float64 = 5.0)

 @assert K &gt; 0.0
 @assert A &gt;= 0.0
 @assert K &gt; A
 @assert r &gt; 0.0
 @assert timesteps &gt; 0.0
 f(t, n) = n*r*(n/A-1)*(1-n/K)
 t = (0.0,timesteps)
 return solve(ODEProblem(f, n0, t))

end

We can now apply this function to get a sense of, for example, what the consequence of the initial population size is:

p1 = plot([0.0],[0.0], frame=:semi)
allee_effect.(linspace(0.1,1.0,50)) .|>
    (x) -> plot!(p1, x, leg=false, c=:grey, lw=0.5)
yaxis!(p1, "Population density", (0, 1.3))
xaxis!("Time")
p1

This gives the following output:

{: .center-image }

We can see the initial value under which the population is not going extinct. And speaking of going extinct…

Models of population dynamics represent the changes in population sizes over time. One direct consequence is that we should almost always expect to have values that are positive or null when the initial population densities are positive or null; this is the theory. In practice, the fact that our models have a mix of linear and non-linear effects, and the fact that we often simulate systems with a reasonably large number of species, means that the numerical integration can go awry. When the population sizes reach values that are very close to zero, it takes very little to bring the integrator to return a value that would make the population size negative at the next step.

This particular bug has been used many times as a feature: when a population size goes negative, we can assume that this represents an extinction. This assumption is very problematic: the system is in a state where it breaks the solver, but we decide that the breakage is not significant enough; for the sake of simplicity, a lot of code for population dynamics assumed that negative density = extinction. Ideally, we would like our numerical tools to only return solutions that make sense.

Enter callbacks. The DifferentialEquations.jl package implements the strategies in  , to ensure that things that should remain positive, actually do. In the few trials we ran with it on a variety of models, it works flawlessly. All that is needed to make sure the solutions remain as positive as possible is to add a PositiveDomain callback – we need to change the final line of the function:

return solve(ODEProblem(f, n0, t), callback=PositiveDomain())

When comparing multiple simulations, it is almost always easier / required to work with the same timesteps from a problem to the other. This too can be done with a solver option: we need to modify the final line so that it reads

 return solve(ODEProblem(f, n0, t),
             callback=PositiveDomain(),
             saveat=linspace(0.0,timesteps,100),
             dense=false)

This way, all outputs will have the same size.