Tutorial 7: Introducing Differential Equations

[ Info: Precompiling PlotsExt [ee28d484-87b8-5b84-8c5c-a74b800bbc50] (cache misses: wrong dep version loaded (2))

This document illustrates how to construct and solve differential equations in Julia using the DifferentialEquations.jl package.

In particular, we are interested in modelling a two species Lotka-Volterra like (predator-prey/consumer-resource) system. Such systems are fundamental in ecology and form the building blocks of complex networks and the models that represent them.

For this tutorial you’ll need the following two packages:

DifferentialEquations.jl to solve the differential equations (same ‘engine’ as the one used by the BioEnergeticFoodWebs.jl model)
Plots.jl to visualise the results

While we already have Plots.jl in this project, we just need to ] add DifferentialEquations.

The DifferentialEquations.jl can be a tad clunky so the below code might take a while to compile. If you hit errors, we recommend removing (] rm DifferentialEquations) and reinstalling (] add DifferentialEquations) the package.

You probably want to start a new script for this exercise. File -> New File -> julia

An introduction to Differential Equations.

Differential equations are frequently used to model the change in variables of interest through time. These changes are often referred to as derivatives (or du/dt). In this case, we are interested in modelling changes in the abundance of a consumer and its resource as a function of the system’s key processes (growth, ingestion of food/foraging and mortality) and its parameters.

This type of model can be formalised as a simple Lotka-Volterra predator prey model, consisting of a set of differential equations:

Resource dynamics: \(\frac{dR}{dt} = r R (1-\frac{R}{K}) - \alpha R C\)
Consumer dynamics: \(\frac{dC}{dt} = e \alpha R C - m C\)

where \(R\) and \(C\) are the abundances of the resource and consumer respectively, \(r\) is the resource’s growth rate, \(K\) is the producer’s carrying capacity, \(\alpha\) is the consumer’s ingestion rate, \(e\) is the assimilation efficiency and \(m\) is the consumer’s mortality rate.

To recall from your ecology modules, this is a system of equations with logistic growth for the producers (plants) and a Type I functional response describing how consumer foraging (the loss of prey caused by predators eating them) varies with prey density, a simple conversion efficiency that translates the functional response (eating the prey) into babies of the consumer (predator reproduction) and then a mortality rate for the consumer/predator.

A few details to note.

there is only one source of density dependence in this model/system - logistic growth of the producer. Density dependence is required to generate coexistence and stability.
the Type I functional response means that the consumer foraging rate increases with prey density and does not saturate. It means that the mortality rate that consumers impose on the prey is constant and independent of the prey density. It is thus density independent and does not contribute to stability. There are two other functional responses. The Type II, which has saturating consumption rate and generates inverse density dependence, and the Type III, which has accelerating consumption rate with prey density, and thus generates increasing mortality in the prey with increasing prey density, and thus density dependence. It is only Type III functional responses that are linked to consumers/predators stabilising their prey populations. It is also a pretty rare functional response!
the consumer mortality rate is also constant and thus density independent.

How does this relate to the Bioenergetic Food Web Model?

The complex model of consumer - resource interactions that is represented by the BEFW model has at it’s heart two equations: one for producers and one for consumers. So at it’s heart, it is not too different in structure to the model above! However, the differences, spelled out in Delmas et al. (2017), are

the parameters are functions of body size
because of this, we can expand the number of resource and consumer equations for multiple values of body size
this creates a potentially big system of equations, instead of just 2.

However, and this is important, the process of solving these equations, whether it is two of them, or 100 of them, is essentially what we are introducing below. If you were to look inside of the BEFW model, you would see the same functions and structure as we introduce here.

The Three Steps For Simulating Food Web Dynamics

There are 3 major steps involved in constructing and solving this model in Julia (these happen to be the same three steps in all programming languages, including R where the deSolve package accomplishes much of this same functionality):

Define a function for your model (i.e. transform the above differential equations into a function that can be read by the solver). This function tells the solver how the variables of interest (here \(R\) and \(C\)) change over time.
Define the problem. Here, the problem is defined by the function, the parameters (\(r\), \(\alpha\), \(e\) and \(m\)), the initial conditions and the timespan of the simulation. In this step you provide the solver with all the details it needs to find the solution.
Solve!

Step 1. Define the function

Here we construct a function for our model. The function needs to accept the following:

du (derivatives) - a vector of changes in abundance for each species
u (values) - a vector of abundance for each species
p (parameters) - a list of parameter values
t (time) - timespan

function LV_model(du,u,p,t)
   # growth rate of the resource (modelled as a logistic growth function)
   GrowthR = p.growthrate * u[1] * (1 - u[1]/p.K) 
   # rate of resource ingestion by consumer (modelled as a type I functional response)
   IngestC = p.ingestrate * u[1] * u[2]
   # mortality of consumer (modelled as density independent)
   MortC = p.mortrate * u[2]
   # calculate and store changes in abundance (du/dt):
   # change in resource abundance
   du[1] = GrowthR - IngestC
   # change in consumer abundance
   du[2] = p.assimeff * IngestC - MortC
end

LV_model (generic function with 1 method)

You’ll notice that in the above function (LV_model), we’ve specified specific parameters using the p.name notation. This is because we’ve opted to store our parameters in a named tuple called p. p is created below but it’s worth noting that when this notation is used e.g., p.growthrate, we are telling Julia that we want to use the value of growthrate that is stored as a named part of our tuple (a fancy name for an object in Julia) p.

Step 2. Define the problem

To define the problem we first have to fix the system’s parameters, the initial values and the timespan of the simulation:

p = (
    growthrate = 1.0, # growth rate of resource (per day)
    ingestrate = 0.2, # rate of ingestion (per day)
    mortrate = 0.2,   # mortality rate of consumer (per day)
    assimeff = 0.5,   # assimilation efficiency
    K = 10            # carrying capacity of the system (mmol/m3)
    )

(growthrate = 1.0, ingestrate = 0.2, mortrate = 0.2, assimeff = 0.5, K = 10)

Here, we have chosen to define p as a named tuple (similar to a list in R). A vector or dictionary would also work, however, named tuples are advantageous because they allow us to use explicit names and are unmutable meaning that once it’s created you can’t change it.

Initial values: For simplicity, we start with \(R = C = 1\):

u0 = [1.0; 1.0]

2-element Vector{Float64}:
 1.0
 1.0

Timespan:

tspan = (0.0,100.0) # you have to use a Pair (tuple with 2 values) of floating point numbers.

(0.0, 100.0)

We then formally define the problem by passing the function (LV_model), the parameters (listed in our named tuple p), the initial values (u0) and the timespan (tspan) to ODEProblem():

prob = ODEProblem(LV_model, u0, tspan, p)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

Step 3. Solve

To solve the problem, we pass the ODEProblem object to the solver.

Here we have chosen to use the default algorithm because it’s a simple problem, however there are several available - see here for more information. These two final steps (define and solve the problem) are analogous to using the deSolve package in R.”

sol = solve(prob)

The solver produces 2 objects: sol.t and sol.u that respectively store the time steps and the variables of interest through time. Let’s have a look.

Step 4. Visualise the outputs”

Once the problem has been solved, the results can be explored and plotted. In fact, the DifferentialEquations.jl package has its own built in plotting recipe that provides a very fast and convenient way of visualising the abundance of the two species through time:

plot(sol, 
    ylabel = "Abundance", 
    xlabel = "Time", 
    title = "Lotka-Volterra", 
    label = ["prey" "predator"], 
    linestyle = [:dash :dot], 
    lw = 2)

One thing to note here, when plotting in Julia you don’t need to separate label names (label = [prey predator]) or linestyles (linestyle = [:dash :dot]) with a comma as you would in R. This will also be case for most plotting options in Julia.

Step 5. Play some games

Now that you see how this works, perhaps you can do a bit of an experiment by varying things in the p = (parameters) part of the code?

Can you increase the number of time steps?
What happens when you increase K?
What happens when you increase or decrease growthrate?
What happens when you change the ingestrate

EXTRA CREDIT

If you are really keen…. how about trying to run the model in a loop over 3 values of growthrate, collect the final population size of the consumer and resource in a data frame with three columns (growthrate, finalCons, finalRes), and plot these values?

References

Delmas, Eva, Ulrich Brose, Dominique Gravel, Daniel B. Stouffer, and Timothée Poisot. 2017. “Simulations of Biomass Dynamics in Community Food Webs.” Edited by Richard Fitzjohn. Methods in Ecology and Evolution 8 (7): 881–86. https://doi.org/10.1111/2041-210x.12713.