Ordinary Differential Equations and their application in Ecology

Introduction & objectives

Welcome to you, young Master 1 student!

It’s now been a few chapters since you plunged into this adventure full of numbers, formulas and other mathematical oddities (🤢). At the moment, you’re probably feeling a strong need for a change of orientation. I mean, we’re doing ecology, not maths! And raising goats in Auvergne sounds pretty cool…kinds of…🐐

But don’t give up! The first (really interesting 😜) chapter is right there in front of you!

Are you interested in population and epidemiological trends and dynamics? Would you like to be able to predict the future (and the past, no kidding) from the present?

Then this is the chapter for you! Work it as much as you can! You’ll end up shining in modelling (and maybe even in society).

At the end, you’ll:
- know how to define an ordinary differential equation
- be able to detect a problem involving the need to use differential equations
- be able to carry out modelling involving an ODE system in an ecological context

I. Dynamic systems & ODE definition

Let’s start simply!

First of all, let’s define the framework of the study.

We consider a dynamic system. Or in other words, a physical system whose state changes as a function of time. By “state”, we mean “physical variables” such as a number of individuals, a position, a temperature, etc., which are functions of time (= which evolve over time).

For example, suppose you wanted to study the evolution over time of a population of aliens newly landed on Earth (yes, it’s possible, I saw that in Men in Black 😎).

Your system is: how the Alien population change.
The state of this population is: the number of individuals, noted ‘👽’.
This state is a function of time (the number of individuals changes over time). We therefore note it: ‘👽(t)’.

So your aim is to study this system. But as you can imagine, you can’t start from scratch! To study the evolution of such a system, you need to know:
- its initial state: values of its characteristic variables at the initial time of the study (t=0).
- the processes underway in the biological system: the way of gaining (i.e. birth, migration) and/or losing (i.e. mortality) individuals for example. Don’t hesitate to draw the system on a piece of paper, to avoid getting lost (Figure 1)!
- the parameters associated with these processes

Let’s make it more concrete:

To keep things simple, we’ll take a discrete-time simulation, with a time step of 1.

Let define t like that: \[t \in \mathbb{N}^+\] This means that t belongs to the set of natural integer numbers. t can therefore take the values 0, 1, 2, 3, 4, 5, etc..

Here, the variable of interest (= the state of the system) is the number of individuals in the population. This number changes over time, and is therefore denoted 👽(t). The biological processes responsible for this change over time are:
- mortality: noted d (for death). d👽 corresponds to the proportion of individuals that will die in 1 unit of time
- birth: denoted b. This means that 👽 individuals will produce b👽 new individuals in 1 unit of time
- migration: denoted m. This means that 👽 increases by m individuals in 1 unit of time.

Let’s try to predict step by step the number of individuals at time (t+1), based on the information known at time (t).

What will the number of individuals at time (t+1) depend on?

Firstly, it depends on the number of individuals at time (t)!

We then note: 👽(t+1) = 👽(t)

But defined in this way, our population is constant. It doesn’t change at all over time.

How can the population increase? - Through migration - Through births

We then note: 👽(t+1) = 👽(t) + migrations + births

But at each time step, the population is also subject to mortality (which obviously reduces the number of individuals).

We then note: 👽(t+1) = 👽(t) + migrations + births - mortality

Let’s replace these words with their notations: 👽(t+1) = 👽(t) + m + b👽(t) - d👽(t)

Now, with parameter values and initial values, you can predict what will happen in the population at time (t+1). Great!

Now a question: is time discrete or continuous? Continuous, yes! In reality, the values of t are defined as follows: \[t \in \mathbb{R}^+\]

We could take the problem from the beginning… or start again from what we have already done, this time integrating the continuous dimension of time.

So we start from: 👽(t+1) = 👽(t) + m + b👽(t) - d👽(t)

We know what happens at each time step of 1. What we don’t know is what happens between each of these time steps! Before working with continuous time, we’ll try to work with a time step of less than 1.

The time step between each time step of 1 is noted: \[\Delta t < 1\]

We then incorporate this new time step into our formula: \[\begin{aligned}[t] & \phantom{=}\left(👽(t+1) = 👽(t) + m + b👽(t) - d👽(t) \right)\\ & \Leftrightarrow 👽(t+\Delta t)= 👽(t) + m\Delta t + b\Delta t👽(t) - d\Delta t👽(t) \\ & \Leftrightarrow 👽(t+\Delta t) -👽(t) = m\Delta t + b\Delta t👽(t) - d\Delta t👽(t) \\ & \Leftrightarrow \frac{👽(t+\Delta t) -👽(t)}{\Delta t} = m + b👽(t) - d👽(t) \\ \end{aligned}\]

Great! Now, what do you have to do to really work in continuous time, and not just be interested in a very small time step?

We need to be interested in a Δt of zero (or very close to 0).

This gives us: \[\lim_{\Delta t \to 0} \frac{👽(t+\Delta t) -👽(t)}{\Delta t} = \frac{d👽}{dt} = 👽'(t) \]

Graphically, this gives (Figure 2):

The blue line is characterized by the slope: \[\frac{👽(t+\Delta t) -👽(t)}{dt}\]

The red line is characterized by the slope: \[\frac{d👽}{dt}\]

There are 3 behaviors to understand: \[\begin{aligned}[t] & \frac{d👽}{dt}>0\Leftrightarrow\text{👽is increasing around t} \\ & \frac{d👽}{dt}=0\Leftrightarrow\text{👽 does not vary around t}\\ & \frac{d👽}{dt}<0\Leftrightarrow\text{👽 is decreasing around t}\\ \end{aligned}\]

This N’(t) thus represents an ‘instantaneous’ increase over time in the number of individuals in the Alien population.

And this is what we call a differential equation: \[N'(t) = \frac{d👽}{dt} = m + b👽(t) - d👽(t)\]

In other words, it is an equation that relates a variable (here, 👽(t)) to its derivatives (here, only d👽/dt). It describes the speed at which a variable (in this case 👽: the number of individuals in the Alien population) varies over time.

II. Single ODE system: examples of exponential & logistic models

Let’s explore two well-known models to better understand the ODE concept.

II-A. The exponential model or Malthus model

The biological hypothesis behind this first model is as follows: The population grows in proportion to its present size.

Let r be defined as follows: \[r \in \mathbb{N}^+\]

This means that a population of N individuals at time t will gain r*N individuals at the next time.

It could be illustrated as follows (Figure 3):

We’re going to define a discrete and a continuous model here to make things very clear!

II-A-1. Discrete model

Let
\[n \in \mathbb{N}^+\]

As mentioned above, the population grows by r*N(n) per unit time.

In other words, between N(n+1) and N(n), the population increases/decreases by r*N(n).

We therefore have: \[N(n+1)-N(n) = r*N(n)\]

This gives the following model: \[N(n+1) = r*N(n)+N(n)\]

Which gives us the following solution: \[N(n) = (1+r)^n * N(0)\]

We’re going to graphically represent this model using r studio. Here are the lines of code to achieve this (Figure 4):

r = 0.7 # r corresponds to the population growth rate; it ranges from -1 to 1; here, it's positive, so the population is growing.
N0 = 10 # N0 corresponds to the initial number of individuals. There are 10 here.
t = seq(0,10,1) # We simply define the duration of the simulation. Here, we start from time 0, up to time 10, with a step of 1.
pop = c() # We create a list to store the values of N(t)
for(n in t){
  N = ((1+r)^n) * N0 # Model solution
  pop = append(pop,N) # Storing results in the 'pop' list
}

plot(t, pop, main = "",
     xlab = "t",
     ylab = "N(t)") # Display it all in a beautiful graph (feel free to make prettier ones using ggplot2 ^^) 

title(main = "Exponential evolution of population as a function of discrete time", adj=0.75)
title(sub = "Figure 4")

II-A-2. Continuous model

Let \[t \in \mathbb{R}^+\]

We had: \[N(n+1)-N(n) = r*N(n)\]

Which gives:
\[\begin{aligned}[t] & \phantom{=}\left(N(n+1)-N(n) = r*N(n)\right)\\ & \Leftrightarrow N(t+\Delta t)-N(t) = r * \Delta t * N(t)\\ & \Leftrightarrow \frac{N(t+\Delta t)-N(t)}{\Delta t} = r * N(t) \\ \end{aligned}\]

This gives us: \[\lim_{\Delta t \to 0} \frac{N(t+\Delta t)-N(t)}{\Delta t} = \frac{dN(t)}{dt} = N'(t) = r*N(t)\]

This gives the following model: \[\frac{dN(t)}{dt} = r * N(t)\]

The solution of this model is: \[N(t) = \exp(r * t) * N_0\]

We’re going to try to simulate this model using r studio. We’re going to use a package you’ve seen or will no doubt be seeing very soon: the ‘deSolve’ package.

Install and load it:

#install.packages("deSolve")
#library(deSolve)

There are two important points to bear in mind when using this package:

• The main function you’ll be using is the ‘ode()’ function. It follows a very specific syntax that you need to remember:

ode(y=y0, times=time, func=function, parms=parameters)

with:

• ‘y0’ = initial value of the variable of interest y (for the first time value)

• ‘times’ = the times you wish to work on

• ‘function’ = the function defining the system; this function follows a very specific syntax that you’ll also need to remember (see below).

• ‘parms’ = parameter values

• concerning the ‘function’ argument of the ‘ode()’ function. This corresponds to a function which, although specific to each model, follows a very particular syntax, as explained above. The syntax is as follows:

name_of_function<- function(t,y,parameters){ instructions list(flux)} with:

• ‘t’ = corresponds to the ‘times’ argument of the ‘ode()’ function
  
• ‘y’ = corresponds to the variable of interest
  
• ‘parameters’ = corresponds to the parameters of the ‘parms’ argument of the ‘ode()’ function

Knowing all this, we can simulate the exponential model with the following script (Figure 5):

expo <- function(t, N, r){
  list(r*N)}
N0 = 10
r = 0.7
temps <- seq(from = 0, to = 10, by = 0.01)
sol <- ode(y = N0, times = temps, func = expo, parms = r)
plot(sol,main = "")
title(main = "Exponential evolution of population as a function of continuous time", adj=0.65)
title(sub = "Figure 5")

II-A-2. Summary

We can see that, whether using the discrete or continuous model, the population is growing. And it’s growing exponentially!

Try changing the parameters to vary the trajectories! Try, for example, r = 0; r < 0; N0 = 0, etc. This will enable you to fully understand the model!

II-B. The logistic model or Verhulst model

The biological assumptions behind this second model are:

• The population grows in proportion to its present size.

• A density-dependency factor (noted K) is introduced to take account of the limiting effect of resources (which are not infinite): this is intraspecific competition!

Let defined: \[r \in \mathbb{N}^+ ; K \in \mathbb{R}^+ ; t \in \mathbb{R}^+ \]

This means that a population of N individuals at time t will gain r*(1-(N/K))*N individuals at the next time.

This gives the following model: \[\frac{dN(t)}{dt} = r * (1-\frac{N(t)}{K}) * N(t)\]

The solution of this model is:

\[N(t) = \frac{K}{1 + \frac{K - N_0}{N_0}\exp(-rt)}\]

Using the explanations about ‘deSolve’, you can illustrate this as follows (Figure 6):

logi <- function(t, N, param){
  r = param[1]
  K = param[2]
  list(r*N*(1-N/K))}
N0 = 100
r = 0.7
K = 100
temps <- seq(from = 0, to = 10, by = 0.01)
sol <- ode(y = N0, times = temps, func = logi, parms = c(r, K))
plot(sol, main = "",
     xlim = c(0,10),
     ylim = c(0,200),
     col = 'green')
title(main = "Logistic evolution of population as a function of continuous time", adj=0.65)
title(sub = "Figure 6")
N0 = 200
sol <- ode(y = N0, times = temps, func = logi, parms = c(r,K))
lines(sol, col = 'red')
N0 = 2
sol <- ode(y = N0, times = temps, func = logi, parms = c(r,K))
lines(sol, col = 'blue')

There are a number of different situations here, which can be summarized as follows:

• If K > N0 > 0, then ∀t > 0, N(t) tends towards K (increase; blue trajectory on Figure 6)
• If N0 = K, then for ∀t >0, N(t) = N0 (green trajectory on Figure 6)
• If N0 > K, then for ∀t >0, N(t) tends towards K (decrease; red trajectory on Figure 6)

Once again, try changing the parameters to vary the trajectories! This will enable you to fully understand the model!

III. 2-equations ODE systems: Lotka Volterra model

An ordinary differential equation system with 2 equations represents the simultaneous rates of change of two related functions, where each equation typically involves the other one, expressing their mutual dependence.

ODE systems with 2 equations are widely used in different disciplines. We will only concentrate on the ecology field to not be overwhelmed by all the possibilities… and also because it’s the most interesting one duh.

In ecology, they are mainly used to model the dynamics of populations interacting.
Generally, if we want to quantify the effect of a species on another species… What easily comes to mind is that it all depends on the population size and the “level” of competitivity of this species.

You’ve surely already heard about Lotka-Volterra (LV) prey-predator model. You may not know this but it IS an ODE system with 2 equations! LV is the most known way to model the dynamic of predator-prey populations.

But Lotka-Volterra also has applications in microbiology (see https://journals.asm.org/doi/10.1128/aem.36.1.11-17.1978)

Lotka-Volterra

A LV prey-predator model represents both equations over time. Although we will mostly focus on the cyclic patterns between those equations, keep in mind that there can be a certain point in time where both equations are in equilibrium and the values that are depending on the time remain approximately constant.

Why LV is useful: - To see that the interactions between populations can create periodic dynamics over time;
- To study how biodiversity interacts with the parameters put into the model, such as the interaction coefficients between populations.

Assumptions for LV:
- The prey and predator live in a limited space;
- The encounters are random;
- The prey is only limited by the predator;
- A predator can consume an infinite number of preys

Biological intuition:
In this model, the predator population thrives when prey are plentiful, but eventually exhaust their resources, and the population declines. When the predator population has declined sufficiently, the prey can reproduce more, and the population increases again. This dynamic continues in a cycle of growth and decline.

Typically, here’s what these equations can look like:

\[ PP = \begin{cases} \frac{d_x}{d_t}= a x - b x y &\text{Prey} \\ \frac{d_y}{d_t}= -c y + h x y &\text{Predator} \end{cases}\] ****

It’s important to really understand how equations are constructed and what they do before getting started.
Don’t be afraid by the looks of all those letters, we’ll gonna break it down to you quite simply:

The main distinction to make in those equation is with “x” and “y”: make sure to recognize that “x” is referring to the prey population (as written in the \[ \frac{d_x}{d_t} \] at the beginning of the line) and that “y” is referring to the predator population (as written in the \[ \frac{d_y}{d_t} \] ). ax: Growth term for the prey population: constant a * current prey population. bxy: Death term for the prey population. Depends on the size of the prey population and the size of the predator’s population. The death of the prey population depends partly on the predator eating them. It’s the same logic the other way around for the second equation because it includes the growth rate of the predator population with cy and the death rate of the predator population with hxy.

In other terms: The prey grows at a linear rate (a) and gets eaten by the predator at the rate of (b). Indeed, the rate of predation on prey is assumed to be proportional to the frequency of encounters between predators and prey; it is equal to bx(t)y(t).
The predator gains a certain amount of vitality by eating the prey at a rate (c), while dying off at another rate (d).

More precisely, x(t) and y(t) functions follow a malthusian model:
The prey population follows an exponential growth in the absence of predators. This growth is represented in the equation with the term ax(t).
The term cy(t) is the natural death of the predators, it is an exponential decrease of the population, in absence of preys.

Demonstration

Let’s dive into R, to simulate and plot the evolution of the two populations, using the Lotka-Volterra prey-predator model.
To simulate this population dynamic, we will use given coefficients, initial levels of populations and also a given simulation time and given time step size. Indeed, We will use the so-called Euler method, as we will be moving through time, and incrementing new values of population levels at each time steps, to keep track of the evolution of the two populations.

What you will have to keep in mind: In this Euler method, we will think of dt as \[\Delta t\]

Thus, we can make this approximation: \[ \frac{d_x}{d_t} = ax - bxy \\ => dx = (ax - byx) * \Delta t \\ and \\ \frac{d_y}{d_t} = -cy + hxy \\ => dy = (-cy + hxy) * \Delta t \]

Defining the terms

Disclaimer: terms’ values are not intended to be ecologically relevant, but to serve as an example.

To make it more telling, let’s say that our prey population is a zebra population, whereas the predator population are a group of lions. Obviously, we would think that the zebra population starts at a higher population size than a simple group of lions. That’s what we tell to our model:

rm(list=ls()) #to clear all the environment

X_init<- 5     #starting value for x
Y_init<- 3     #starting value for y

The initial values for the populations’ sizes have been defined. Now, since LV is a dynamic model, evolving over time, and since we do not want our simulation to run endlessly, it’d be a good idea to set a time limit, wouldn’t it be? Let’s say we want to follow the population dynamics throughout 100 months.

For our example, we will talk about numbers of individuals (Xinit=5 individuals, Yinit=3 individuals), and we will consider the time unit to be months. But we could have chosen days or years, and then chose consistent coefficient values:

tend<-100    #how long the simulation is going to run for
delta_t<-0.01    #the size of the time step, as we will be moving through time with discrete steps

# Model coefficients
a<-0.8 #prey
b<-0.4 #prey
c<-0.6 #predator
h<-0.2 #predator

Defining the vectors

We create some vectors to keep track of the variables, as we are moving through time. Then those vectors will prove useful when we will want to plot the dynamics:

X<-c(X_init)    #we will append to this vector, it is going to be filled up with the successive levels of X, during all the time of the simulation
Y<-c(Y_init)    #idem with Y

t<-c(0)   #this vector will be filled up with all the time steps

Simulation

while (t[length(t)]<tend) 
{
  current_X<-X[length(X)]   #we define current_X as the most recent value of X
  current_Y<-Y[length(Y)]    #idem for Y
  current_t<-t[length(t)]
  
  delta_X<-(a*current_X-b*current_X*current_Y)*delta_t
  next_X<-current_X+delta_X
  
  delta_Y<-(-c*current_Y + h*current_X*current_Y)*delta_t
  next_Y<-current_Y+delta_Y
  
  next_t<-current_t+delta_t
  
  #We append the successive values to the vectors
  X<-append(X,next_X)
  Y<-append(Y,next_Y)
  t<-append(t,next_t)

}

Plotting

library(ggplot2)

LV <- data.frame(time = t, prey = X, predator = Y)

ggplot(LV, aes(x = time)) +
  geom_line(aes(y = prey, color = "Prey"), size = 1.5) +
  geom_line(aes(y = predator, color = "Predator"), size = 1.5) +
  labs(x = "Time (months)", y = "Number of Individuals") +
  scale_color_manual(values = c(Prey = "blue", Predator = "red")) +
  theme_classic()

You can see the cyclic nature of this dynamic as explained earlier: the prey thrives when the predator declines, the predator thrives when there’s a lot of prey…

This dynamic can be seen in a phase space plot i.e. the level of predator as a function of the prey level:

ggplot(LV, aes(x = prey, y = predator)) +
  geom_point() +
  labs(x = "Prey (nb of individuals)", y = "Predator (nb of individuals)") +
  theme_classic()

Adding a parameter to make in more realistic

At that point, the model doesn’t contain a lot of terms. It is a great mathematical setting, we can understand how it works and what is its point. However, in nature, there are a lot more parameters controlling our populations. Adding some of these parameters would allow the model to be more realistic.
To cite only one parameter, we can talk about the carrying capacity of the environment. It defines how the environment can tolerate the species, according to the resources it contains. The more individuals from a species there is, the less the environment can keep up, until a certain point where there is too much individuals for a limited resource so the population declines…

Actually, this notion of carrying capacity (K) can bring us to an unaddressed point: reaching an equilibrium for both prey and predator populations, instead of endlessly going in circles.

The equations become: \[ PP = \begin{cases} \frac{d_x}{d_t}= a x \ (1 - \frac{x}{K}) - b x y &\text{Prey} \\ \frac{d_y}{d_t}= -c y + h x y &\text{Predator} \end{cases}\]

Exercises to understand:

1. Supposing that you want to depict the possible seasonal variation of the coefficient “a”, rewrite the code to simulate a model that includes the variability of “a” among time. We can consider a variation of “a” over a time period of one year, with an amplitude of 0.6. Moreover, we consider a mean value of “a” of 0.8.

Answer:

rm(list=ls())

####DEFINING THE COEFFICIENTS AND VARIABLES####
X_init<-5
Y_init<-3
tend<-100
delta_t<-0.01


####DEFINING THE VECTORS####
X<-c(X_init)   
Y<-c(Y_init)

t<-c(0)   #this vector will be filled up with all the time steps

###coefficients of the model
#function to describe the variation of a among time
a_seasonal_variation <- function(t) {
  amplitude <- 0.6  #amplitude of the seaonal variation
  period <- 12      #seasonal period
  mean <- 0.8     #mean value of the growth rate
  return(mean + amplitude * sin(2 * pi * t / period))
}

b<-0.4
c<-0.6
h<-0.2

####SIMULATION####

while (t[length(t)]<tend) 
{
  current_X<-X[length(X)]   #we define current_X as the most recent value of X
  current_Y<-Y[length(Y)]    #idem for Y
  current_t<-t[length(t)]
  
  delta_X<-(a_seasonal_variation(current_t)*current_X-b*current_X*current_Y)*delta_t
  next_X<-current_X+delta_X
  
  delta_Y<-(-c*current_Y + h*current_X*current_Y)*delta_t
  next_Y<-current_Y+delta_Y
  
  next_t<-current_t+delta_t
  
  #We append the successive values to the vectors
  X<-append(X,next_X)
  Y<-append(Y,next_Y)
  t<-append(t,next_t)
  
}

####PLOTTING#####

library(ggplot2)

LV_exercice <- data.frame(time = t, prey = X, predator = Y)

ggplot(LV_exercice, aes(x = time)) +
  geom_line(aes(y = prey, color = "Prey"), size = 1.5) +
  geom_line(aes(y = predator, color = "Predator"), size = 1.5) +
  labs(x = "Time (months)", y = "Number of Individuals") +
  scale_color_manual(values = c(Prey = "blue", Predator = "red")) +
  theme_classic()

2. What factors, other than the carrying capacity of the environment, could be added that would have ecological significance?

Answer:

• sǝʇɐɹ uoᴉʇɔnpoɹdǝɹ ɟo ʎʇᴉlɐuosɐǝS
• ǝʇɐɹ uoᴉʇɐpǝɹd ɟo ʎʇᴉlɐuosɐǝS
• uoᴉʇɐɹɓᴉꟽ
• uoᴉʇᴉʇǝdɯoɔ ɔᴉɟᴉɔǝdsɐɹʇuI
• .ɔʇǝ

You can use a online generator to flip the text back up, or just take a picture and turn it!

IV. 3-equations ODE system

We will now make the system more complex by adding an additional equation. To solve these systems, two methods are possible. If we are in the case of a closed population, then we can remove an equation: see example SIR model. If we are in the case of an open population: see example of the Lotka-Volterra model with 3 species.

1) SIR model and resolution by hand

The most used models in epidemiology are mathematical models of infectious diseases such as the SIR model. These models are based on compartmentalization that divide the population into various possible disease states. The SIR model is composed of 3 compartments: initially healthy or susceptible individuals (S), infected individuals (I) and cured individuals or resistent (R).

\[SIR=\begin{cases}\frac{d_S}{d_t}=- \alpha SI+\gamma R &\text{Susceptible} \\ \frac{d_I}{d_t}=\alpha SI-\beta I&\text{Infected} \\ \frac{d_R}{d_t}=\beta I- \gamma R &\text{Resistent}\end{cases} \\ \\ \begin{align} With : & - \alpha \text{: disease transmission rate} \\ & - \beta \text{: cure rate} \\ & - \gamma \text{: coefficient of loss of immunity} \end{align}\]

We assume a closed population, so: \[ S+I+R = 1 \Leftrightarrow R=1-S-I \]

Thus, we can express R from S and I, so we can remove the equation from our system and we end up with a system with 2 equations.

Balance points

We look for balance points such as:

\[ SIR=\begin{cases}\frac{d_S}{d_t}=0 &\text{Susceptible} \\ \frac{d_I}{d_t}=0&\text{Infected} \end{cases} \]

To do this, we solve:

\[ \begin{align} & \begin{cases} - \alpha SI+\gamma R = 0 &\text{Susceptible} \\ \alpha SI-\beta I =0 &\text{Infected} \end{cases} \\ \Leftrightarrow & \begin{cases} -\alpha SI+\gamma (1-S-I) = 0 \\ \alpha SI-\beta I =0 \end{cases} \end{align} \]

We first focus on solving:

\[\begin{align} & \alpha SI - \beta I=0 \\ \Leftrightarrow \ & I(\alpha S- \beta) =0 \\ \Leftrightarrow \ & I=0\ or\ S= \frac{\beta}{\alpha} \end{align}\]

If I = 0, we integrate I=0 into our “susceptible” equation to obtain the solution for S:

\[\begin{align} & \ \ \ \ \ \begin{cases} -\alpha S *0 + \gamma (1-S-0)=0 \\ I=0 \end{cases} \\ & \Leftrightarrow \begin{cases} \gamma (1-S)=0 \\ I=O \end{cases} \\ &\text{We thus obtain our first solution:} \\ &\Leftrightarrow \begin{cases} S=1 \\ I=0 \end{cases} \end{align} \]

It’s the disease-free equilibrium (DFE).

If S = b/a, we integrate S into our “susceptible” equation to obtain the solution for I:

\[\begin{align} & \ \ \ \ \ \begin{cases} - \alpha \frac{\beta}{\alpha} I + \gamma (1- \frac{\beta}{\alpha} - I) =0\\ S = \frac{\beta}{\alpha} \end{cases}\\ & \Leftrightarrow \begin{cases} -\beta I + \gamma- \frac{\beta \gamma}{\alpha} -\gamma I = 0 \\ S = \frac{\beta}{\alpha} \end{cases} \\ & \Leftrightarrow \begin{cases} \beta I + \gamma I = \gamma (1 - \frac{\beta}{\alpha}) \\ S= \frac{\beta}{\alpha} \end{cases}\\ & \text{We thus obtain our first solution:}\\ & \Leftrightarrow \begin{cases} I= \frac{\gamma (1- \frac{\beta}{\alpha})}{\beta + \gamma} \\ S=\frac{\beta}{\alpha} \end{cases} \end{align} \]

It’s the endemic equilibrium (EE).

Jacobian matrix

Now, let’s perform linear stability analysis for each equilibrium point by computing the Jacobian matrix:

\[ J=\begin{bmatrix} \frac{\delta S}{\delta S} & \frac{\delta S}{\delta I} \\ \frac{\delta I}{\delta S} & \frac{\delta I}{\delta I} \end{bmatrix} \Leftrightarrow \begin{bmatrix} - \alpha I- \gamma & - \alpha S- \gamma \\ \alpha I & \alpha S- \beta \end{bmatrix} \]

Analyze the stability of each equilibrium point by examining the signs of the real parts of the eigenvalues. If the real parts are all negative, the equilibrium is stable; if any are positive, the equilibrium is unstable. The stability will depend on the specific values of the parameters α and β.

For the disease-free equilibrium:

\[ J(1,0)= \begin{bmatrix} -\gamma & - \alpha -\gamma \\ 0 & \alpha - \beta \end{bmatrix} \]

The Jacobian matrix for the pair (1,0) is an upper triangular matrix, so the eigenvalues are the diagonal elements: \[ \lambda _1 = - \gamma \ \text{and} \ \lambda _2 = \alpha - \beta \]

γ is a coefficient representing loss of immunity, so γ>0.

So: \[ \begin{align} & \lambda _1 <0 \\ & \text{If:} \ - \alpha > \beta , \ \lambda _2 > 0 \ \text{Unstable equilibrium point} \\& \ \ \ \ - \alpha < \beta , \ \lambda _2 < 0 \ \text{Asymptotically stable point} \end{align} \]

In order to better understand the behavior of the unstable point, we look for the sign of the determinant of the matrix:

\[ det(J(1,0))= - \gamma (\alpha - \beta ) \\ \Leftrightarrow \ \gamma (\beta - \alpha) \\ \gamma > 0 \ and \ \alpha > \beta \\ det(J(1,0)) < 0 \Rightarrow \text{Saddle point} \]

Saddle points exhibit a combination of stable and unstable behavior in different directions in the state space.

For the endemic equilibrium

\[ J(\frac{\beta}{\alpha},\frac{\gamma (1-\frac{\beta}{\alpha})}{\beta + \gamma})= \begin{bmatrix} \frac{\gamma (\beta- \alpha)}{\beta + \gamma} - \gamma & - \beta -\gamma \\ \frac{\gamma (\alpha - \beta)}{\beta + \gamma} & 0 \end{bmatrix} \]

Here, we do not have a characteristic matrix so to determine the eigenvalues, it will be necessary to solve the characteristic polynomial. We will first look at the sign of the determinant of the matrix to study the behavior of the system.

So, we calculate the determinant: \[\begin{align} & det(J(\frac{\beta}{\alpha},\frac{\gamma (1-\frac{\beta}{\alpha})}{\beta + \gamma})) = - (- \beta - \gamma) * \frac{\gamma (\alpha - \beta)}{\beta + \gamma} = \gamma (\alpha - \beta) \end{align}\]

If α < β, det(J) < 0: a saddle point.

If α > β, det(J) > 0: look the trace of the matrix.

\[Tr(J(S,I)) = \frac{\gamma (\beta - \alpha)}{\beta - \gamma} - \gamma \\ \\ \gamma > 0 \ \text{and} \ \alpha > \beta \\ \\ \text{So:} \ Tr(J(S,I)) <0 \]

In this case we have a stable point. The specific type of stability depends on the signs and nature of the eigenvalues.

Now, we solve the characteristic polynomial to obtain the eigenvalues:

$$ \[\begin{align} & det(J- \lambda Id) = 0\\ & \text{with Id : identity matrix}\\ &\Leftrightarrow \begin{bmatrix} \frac{\gamma (\beta - \alpha)}{\beta + \gamma}-\gamma -\lambda & - \beta - \gamma\\ \frac{\gamma (\alpha - \beta)}{\beta + \gamma} & - \lambda \end{bmatrix} = 0\\ &\Leftrightarrow - \lambda (\frac{\gamma (\beta - \alpha)}{\beta + \gamma}-\gamma -\lambda) - (- \beta - \gamma )(\frac{\gamma (\alpha - \beta)}{\beta + \gamma} ) = 0 \\ &\Leftrightarrow \lambda^2 + \lambda (\gamma + \frac{\gamma (\alpha- \beta)}{\beta + \gamma})+ \gamma (\alpha - \beta) = 0 \end{align}\]$$

The number of solutions is indicated by the value of the discriminant : \[ \begin{align} & \Delta = (\gamma + \frac{\gamma (\alpha - \beta)}{\beta + \gamma})^2 - 4 (\gamma (\alpha - \beta)) \\ &\text{With} \ \alpha > \beta \\ \\ \text{If} \ & \Delta < 0:\ \text{you will have two complex conjugate roots.} \\ &\Delta = 0:\ \text{you will have one real repeated root.} \\ &\Delta>0: \ \text{you will have two distinct real roots.} \end{align} \]

From now,we will assign numerical values: \[ \alpha = 4, \ \beta = 2 \ and \ \gamma = 1 \\ \Delta = (1 + \frac{1(4-2)}{1+2})^2 - 4(1(4-2))= -5.22 \\ So \ \Delta < 0 \]

\[ \begin{align} &\text{Solutions are given by the following formula:} \\ & \lambda= \frac{-(\gamma + \frac{ \gamma (\alpha - \beta)}{\beta + \gamma}) \pm \sqrt \Delta}{2}i \\ & \frac{\sqrt \Delta}{2}i: \ \text{imaginary part.} \\ & \frac{-(\gamma + \frac{\gamma (\alpha - \beta)}{\beta +\gamma})}{2}: \ \text{real part.} \\ \\ & \text{Real part}= -\frac{5}{3} \end{align} \]

The stability of a linear system is often determined by the real parts of the eigenvalues. In this case, the real part is negative which indicates a stable attractive node.

This sounds complicated, but don’t panic! In ecology we are lucky to have a great tool, R, which allows us to resolve this type of system. You can see how to do it on R at the end of Lotka Volterra 3 species.

Solving SIR model with R

Use R to study the behavior of our model as double interest.

• Verify if our theoretical calcul is okay

• Help you to have an idea of the behavior of the model

To show how you can use phaseR and deSolve to solve ODE we will use the SIR we just solved. We hope our calcul is good. In this chapter, we will show what’s the most common function you can use, but not all the parameters you can add. I know you are probably sobbing “wHyYyYyy?” (Who has never dreamed of knowing all the difference between the Euler and Runge Kunta order 4?) but we let you make your search if you need more information.

Step 1: Prepare you environnement

Let’s start by erasing all the data in our environment in order to don’t bother us when we work. Then we load the 3 packages we will use.

• The most famous library to solve and study ODE is deSolve you can also compute yourself if you want but it will be a harsh task. A great point of deSolve it’s its simplicity to use and we hope you will agree with us a the end of this chapter

• phaseR Performs a qualitative analysis of one- and two-dimensional autonomous ordinary differential equation systems, using phase plane methods

• We will also use ggplot to make a beautiful graph. It’s also a little more intuitive way to create and manipulate plot than using a basic ploting system in r

• rgl is a package which give you tools to creat 3d plot and objects

• tidyr is like ggplot but for data manipulation. It doesn’t bring very special tools but able us to transform easily our table

As we talk about programmation we have to talk about technical issues you can have. To prevent you from committing this, we put during the presentation some technical issues to think of before coding. If you want to do the same things as us, you will probably find yourself in some problems: don’t be afraid R is widely used and you will easily find people who have defeated the same problem. Don’t be ashamed to search for help!

#We erase all data use before
rm(list = ls())
#load library we need
library(ggplot2)
library(deSolve)
library(phaseR)
library(rgl)
library(tidyr)

Step 2: Define your model

It’s time to start modeling! The first step is to create our function which gives us the evolution of our population according to time. This function takes a step of time, the initial population value, and the value of different parameters. The output of this function needs to be a list of the new populations. To have the behavior of our function we have to use the same hypothesis as in the calculation part so we will not code how the Recover compartment evaluates (R). If we need the Recover evolution we will use: \[1 = S + I + R \Leftrightarrow R = 1 - S - I \] Technical consideration: When you build your function you are free to make some modification BUT the order of the input has to be time, population, and parameter value because the ODE function want you to do this!

mod_SIR <- function(t , dy , parameters){
#we take our population and give give to to different variable
  S <- dy[1]
  I <- dy[2]

#We give the value of parameters to the variable (there is greater way to do this but i let you find out)
  alpha <- parameters[1]
  beta <- parameters[2]
  gamma <- parameters[3]

#sol is the vector where we stock our solution
  sol <- numeric(2)

#We compute the equation need
  sol[1] <- - alpha * S * I + gamma * (1 - S - I)
  sol[2] <- alpha * S * I  - beta * I
  
#And list solution 
  list(sol)
}

Step 3: Study how the system evaluates among time

Now we have done the heart of our model let’s study it with ODE. It’s from the deSolve package to solve our model. To do this we define some quantities. For the time we use the seq function which creates a vector that goes from 0 to 600 by 0.01 step. For parameter values, we create a vector in the same order as in our function. And value we start to form

Then we use ODE to simulate how our population evolves at each time in our vectors. We take the output from ODE and manipulate it a little bit to make visualization easier.

#Definition of parameters we will need for the suite
param <- c(alpha = 0.2, beta = 0.1, gamma = 0.2)
times = seq(from = 0 , to = 600, by = 0.01)
depart = c(0.99 , 0.001 )

#Calculation of the evolution of our population 
resolnum <- ode(y = depart, times = times, func = mod_SIR , parms = param)

#Transfortmation of data in order to be use easely
resolnum <- as.data.frame(resolnum)
resolnum$R<- 1-  resolnum[,2] + resolnum[,3] 

colnames(resolnum) <- as.factor(c("Times" , "Susceptible" ,"Infectious" ,"Recovered"))



#Plotting time
ggplot(resolnum,aes(x=Times))+
  geom_line(aes(x = Times, y = Susceptible),color ="chartreuse" ,size =1.5)+
  geom_line(aes(x = Times,y = Infectious),color = "red" ,size =1.5)+
  geom_line(aes(x = Times,y = Recovered),color = "darkviolet" ,size =1.5)+
  labs(title = "Population evolution among times", y = "Population")+
  theme_classic()

Now we have an idea of how our population evaluates thought time we can try to find the behavior

Step 4: Equilibrium point and behavior

Let’s study the behavior of our model and its equilibrium points. We’re going to study S as a function of I. First function displays model behavior. Second draws the lines where the models cancel out.

Technical consideration: To have all nullclines you should extend the area you want to explore. For example, we only have interest in the area where S and I go to 0 at 1(we made the hypothesis the sum of all compartments should be equal to 1) but in our simulation, we go from -0.1 to 1.1.

#Plot the velocity field
  field<-flowField(mod_SIR ,xlim=c(-0.1,1.1),ylim=c(-0.1,1.1),parameters = param,add = F,xlab = "Susceptible" , ylab = "Infectious")
  #Draw the isoclines of our model
nullc<-nullclines(mod_SIR ,xlim=c(-0.1,1.1),ylim=c(-0.1,1.1),parameters = param,points = 100, add = T,add.legend = F,size =2)

As we can see at crossing of nullclines we have two equilibrium points. The first is at the center of our phase diagram, and the arrows seem to show that we’re tending towards this point. A second is present where I = 0 and S = 1 on the contrary, we seem to be moving away from it. If this graph allows us to observe 2 points of equilibrium, it is not easy to determine their behavior.

To look at the behavior near these points, we can use the trajectory function, which will plot the trajectories of our system as a function of a starting point. Here we repeat the trajectory functions in order to have a 51 scenario.

#we define our initials conditions
X0 <- sample(x = seq(0,1,0.01),size =1)
Y0 <- sample(x = seq(0,1-X0,0.01),size =1)
cond_init = c(X0, Y0)
#We simulate what's happened for one trajectory
  t<-trajectory(mod_SIR ,cond_init,tlim= c(0,50),parameters = param,col = "black",add = F,xlim=c(0,1),ylim = c(0,1),xlab = "Susceptible" , ylab = "Infectious")
  
#We repeat the same process for 50 other trajectory
for(i in 1:50){
  X0 <- sample(x = seq(0,1,0.01),size =1)
Y0 <- sample(x = seq(0,1-X0,0.01),size =1)

cond_init = c(X0, Y0)
  traj<-trajectory(mod_SIR ,cond_init,tlim= c(0,50),parameters = param,col = i,add = T)
}

What the vector field showed us seems to be confirmed by this graph. Let’s check our observations with new functions (promise they’re the last)

We’ll check the behavior of the center point using the stability function. We could also have used the findEquilibrum function to find all the equilibrium points, but we’ll let you have your fun.

# Checking stability of the point we see
stability(SIR,ystar = c(0.5,0.5),parameters = param )$classification[1]

tr = -0.1, Delta = 0.01, discriminant = -0.03, classification = Stable focus

[1] "Stable focus"

The two functions confirm our observations and the calculations made just before - what more could you ask for?

There are still other functions in the deSolve and pahse R packages that could be of interest, but which are useful in more specific cases, so we won’t look at them here.

2) Lotka-Volterra: A three-species model

Biotic interactions are complex, requiring the interaction of numerous parameters. If we take the Lotka-Volterra example, it only considers the interactions between two species, prey and predator. It assumes a constant mortality rate among predators, regardless of prey. However, predators are themselves the prey of other super-predators. Their mortality rate should therefore depend on the super-predators they encounter. Thus, the prey-predator-superpredator extension of the Lotka-Volterra model would be a good representation of natural trophic interactions.

In this case, we considerer a superpredator, a predator an a prey.

A superpredator is an organism that reaches adulthood at the top of the food chain. Thus, by definition, a superpredator is not the prey of any other species. The predator is consider to be in the middle of the food chain and the prey at the bottom.

The model can be written as the following system of differential equations:

\[ PPS=\begin{cases} \frac{d_x}{d_t}=ax-bxy &\text{Prey} \\ \frac{d_y}{d_t}=-cy+hxy-eyz&\text{Predator} \\ \frac{d_z}{d_t}=-fz+gyz &\text{Superpredator} \end{cases} \]

The model parameters are:

• a: x natural growth rate

• -f and -c: natural death rate of y and z

• The growth rate of y and z depends on the number of prey caught. It is modeled by hxy for y and by gyz for z.

• The death rate of x and y depends on the number of top predators y for x and z for y. It is modeled by -bxy for x and by -eyz for y.

This model is only useful for chains of predators. It doesn’t take into account the presence of several predators with different dynamics on the same prey.

Let’s try to solve this system: To do this, the first step is to find the balance of the model, that is to say the values (x,y,z) so that:

\[\frac{d_x}{d_t}=\frac{d_y}{d_t}=\frac{d_z}{d_t}=0\]

We take our previous systems that we factored:

\[ PPS=\begin{cases}\frac{d_x}{d_t}=x(a-by) &\text{Prey} \\ \frac{d_y}{d_t}=y(hxy-c-ez)&\text{Predator} \\ \frac{d_z}{d_t}=z(gy-f) &\text{Superpredator}\end{cases} \]

The obvious solution is (0,0,0) but is not the only one. Indeed, for the derivatives to be equal to 0, then we must find the following solution:

\[(a-by)=(gy-f)=(hx-c-ez)=0\]

We first solve: \[a-by=0\Leftrightarrow a=by \Leftrightarrow y=\frac{a}{b}\] \[gy-f=0\Leftrightarrow f=by \Leftrightarrow y=\frac{f}{g}\]

So, a solution for a yeq is: \[y_{eq}=\frac{a}{b}=\frac{f}{g}\] Therefore, there is an equilibrium when: \[ga=fb\]

In order to determine the other solutions, we must solve the following jacobian matrix: \[ J=\begin{bmatrix}a-by&-bx&0\\hy&-c+hx-ez&-ey\\0&gz&gy-f\end{bmatrix}\]

The steps are the same as for the SIR model but here the matrix is 3x3. We invite you to look at the chapter on matrix calculation to determine the solutions.

Solving Lotka-Volterra with R

Hey R is back! You loved last part with how to solve 2 equation ODE in R ? Now let’s find out how to solve a 3 equation system! Unfortunately, trying to solve this equation by hand leads to problems: it’s not easy to find all the equilibrium points. As mentioned in the previous section, phaseR is not suitable for solving three-equation systems. We can compute our function but unless you a nerd or MODE you don’t want to do this. the solution we choose in this paragraph it’s to only study numercally system. So let’s start again what we just do! We use the same library as in last part.

We create a function with the same limitation as in our first example.

rm(list = ls())
PPS <- function(t , dy , parameters){

  x <- dy[1]
  y <- dy[2]
  z <- dy[3]
  
  r1 <- parameters[1]
  r2 <- parameters[2]
  r3 <- parameters[3]
  b1 <- parameters[4]
  b2 <- parameters[5]
  c1 <- parameters[6]
  c2 <- parameters[7]
  
  sol <- numeric(3)
  
  sol[1] <- r1 * x - b1 * y * x
  sol[2] <- - r2 * y + b2 * x * y  - c1 * y * z 
  sol[3] <- - r3 * z + c2 * z * y 
 
  list(sol)
}

In order to determine the behavior of our system we should use the ode function to make many graph and try to understand numerically which parameter ar important. In the following example we pass this step and present you 2 cases we can have in our system.

This first case is a very special one, since we observe a runaway system, with the number of prey increasing exponentially without being limited by predators. While mathematically this model makes sense, it is not very common in nature, due to the absence of infinite resources.

param <- c(r1 = 0.01, r2 = 0.2 , r3 = 0.1, b1 = 0.04, b2 = 0.04, c1 = 0.09, c2 = 0.01)
temps = seq(from = 0 , to = 1000, by = 0.01)
depart = c(20 , 10 , 5)
resolnum <- ode(y = depart, times = temps, func = PPS, parms = param)

resolnum <- as.data.frame(resolnum)
colnames(resolnum) <- c("Times" , "Prey" ,"Predators" ,"SuperPredators")



ggplot(resolnum)+
  geom_line(aes(x = Times, y = Prey),color ="chartreuse" ,size =1.5)+
  geom_line(aes(x = Times,y = Predators),color = "red" ,size =1.5)+
  geom_line(aes(x = Times,y = SuperPredators),color = "darkviolet",size =1.5)+
  theme_classic()+
  labs(title = "Case 2: Balance system")

ggplot(resolnum)+
   geom_point(aes(x = Prey, y = Predators),color ="black" ,size =1.5)+
theme_classic()

ggplot(resolnum)+
   geom_point(aes(x = Prey, y = SuperPredators),color ="black" ,size =0.2)+
theme_classic()

ggplot(resolnum)+
   geom_point(aes(x =  SuperPredators, y = Predators),color ="black" ,size =1.5)+
theme_classic()

plot3d(resolnum$Prey,resolnum$Predators,resolnum$SuperPredators,xlab ="Prey",ylab ="Predators",zlab ="SuperPredators")
rglwidget()

In this case, our system moves towards a balance between prey and predators. This time, however, superpredators are excluded from our system.

For a more thorough analysis, we’d have to simulate the same model but draw our parameters randomly, but unfortunately this takes a lot of time, so we’d have to do it with very few simulations. To get a better idea of the model’s behavior, more advanced simulations should be carried out, but unfortunately the fun has to stop at some point.

If you eat graphics cards and want to delve deeper into the subject, the book Solving Differential Equations in R may be of some interest, even if it deals with much more complex and interesting problems than we do.

3) Generalizations of the three-species model and other examples

a) Generalizations of the three-species model

The three-species model presented previously is a simplified model that models a perfectly linear trophic chain. However, nature is more complex and a species can be preyed upon by multiple predators and predators can predate multiple prey items.

We invite you to go see the course of Lalith Devireddy (2016) entitled “Extending the Lotka-Volterra Equations” which develops this model and its resolutions.

Many models with 3-ODE systems exist to model ecological dynamics. Here we will present one additional ones to you but we will not develop their resolution.

b) NPZ model

The NPZ model makes it possible to model the lake production of a pelagic ecosystem. This model studies variations in nutrients, phytoplankton and zooplankton that are linked to each other.

Here are the equations of the system, but it will not be worked up here. \[ \begin{align} & NPZ = \begin{cases} \frac{dP}{dt}= \frac{V_m N}{k_s + N}P -mP -ZR_m(1-e^{-\Lambda P}) & \text{Pythoplankton} \\ \frac{dZ}{dt}= \gamma ZR_m(1-e^{- \Lambda P}-dZ) & \text{Zooplakton} \\ \frac{dN}{dt}=-\frac{V_m N}{k_s + N}P+mP+dZ+(1- \gamma ) ZR_m (1-e^{- \Lambda P}) & \text{Nutrients} \end{cases} \\ \\ & \text{Definitions of variable:} \\ & - V_m: \ \text{Maximum phytoplankton growth rate} \\ & - N: \ \text{Nutrient concentration} \\ & - k_s: \ \text{Half saturation constant for nutrients} \\ & - P: \ \text{Phytoplankton stock size} \\ & -m: \ \text{Phytoplankton mortality rate} \\ & -Z: \ \text{Zooplankton stock size} \\ & - \gamma : \ \text{Zooplankton growth efficiency} \\ & - R_m : \ \text{Maximum zooplankton ration} \\ & - \Lambda : \ \text{Ivlev constant} \\ & -d: \ \text{mortality rate} \end{align} \]

Conclusion & motivation words!

What can we conclude from all this?

An Ordinary Differential Equation (ODE): - has many applications: in mechanics, physics, geometry and biology - is an equation relating a variable to its derivative(s) - describes the rate at which a variable varies over time - has an infinite number of solutions, so it’s important to define the initial conditions

It is possible to solve ODE systems with one, two, three… equations that may contain ODEs of order 1, 2 or 3… But this is complex to solve!