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.

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 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 , things proceed a little bit differently. The problem to integrate must be a function f(t,u), where t is 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:

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 > 0.0
 @assert A >= 0.0
 @assert K > A
 @assert r > 0.0
 @assert timesteps > 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:

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.

More reading.

4 thoughts on “Solving differential equations for ecological problems.

  1. This looks fantastic! The interface looks clean, and straightforward to work with.

    Do you know if the developers plan to develop support for bifurcation / continuation methods for ODEs? I haven’t dove into Julia yet, but I’ll say I would make the plunge the day that someone provides a decent, rigorous bifurcation tool that doesn’t have the user interface pains of AUTO-07p or the the rigid interface of XPPAUT.

    1. Timothée says:

      Good question! I’m not using these methods, but the documentation has a page on bifurcation analysis: http://docs.juliadiffeq.org/stable/analysis/bifurcation.html

      1. Nice, thanks! That’s really exciting. Interesting that they’re using PyDSTool as the backend; I’ve tried to use it for continuation before, and found it pretty arcane and really poorly documented, but it looks like this interface to it is a lot more straightforward to use.

        Guess I need to add “learn Julia” to the todo list for the next few weeks!

      2. Also, after digging into the help files: it has an SDE solver, with decent options for different noise types and higher-order solvers? And it actually uses readable code and could be taught without needing a degree in comp sci? Sign me up for this!

Leave a Reply

Your email address will not be published. Required fields are marked *