Population ecology: using deSolve to solve ODEs in R

We will create functions and use the package deSolve to numerically integrate population dynamics. This package contains several solvers for diferential equations that perform numerical integration.

library(deSolve)
library(ggplot2) # because we will plot things
library(tidyr) # because we will manipulate some data

Logistic growth model

Let’s start with a simple example using the logistic growth model.

\[\frac{dN}{dt} = rN(1 − \alpha N)\]

# Creating a function for logistic growth
logGrowth <- function(t, y, p) {
  N <- y[1]
  with(as.list(p), {
    dN.dt <- r * N * (1 - a * N)
    return(list(dN.dt))
  })
}

We will specify the parameters and then solve the ODE.

# named vector with parameters
p <- c(r = 1, a = 0.001)
# initial condition
y0 <- c(N = 10)
# time steps
t <- 1:20

Then, we can use the ode() function to solve the model for our defined parameters.

# give the function and the parameters to the ode function
out_log <- ode(y = y0, times = t, func = logGrowth, parms = p)

Let’s explore the output. The object out is a matrix with a first column containing the time stamp and the following columns, the state variable(s).

head(out_log)

     time         N
[1,]    1  10.00000
[2,]    2  26.72369
[3,]    3  69.45310
[4,]    4 168.66426
[5,]    5 355.46079
[6,]    6 599.85982

Plotting the result using ggplot2. We will need to convert the out object into a data.frame.

df_log <- as.data.frame(out_log)
ggplot(df_log) +
  geom_line(aes(x = time, y = N)) +
  theme_classic()

Lotka-Volterra competition model

\[\frac{dN_1}{dt} = r_1 N_1 (1 − \alpha_{11} N_1 − \alpha_{12} N_2)\]

\[\frac{dN_2}{dt} = r_2 N_2 (1 − \alpha_{22} N_2 − \alpha_{21} N_1)\]

First, we creating a function to represent the ODEs.

LVComp <- function(t, y, p) {
  N <- y
  with(as.list(p), {
    dN1.dt <- r[1] * N[1] * (1 - a[1, 1] * N[1] - a[1, 2] * N[2])
    dN2.dt <- r[2] * N[2] * (1 - a[2, 1] * N[1] - a[2, 2] * N[2])
    return(list(c(dN1.dt, dN2.dt)))
  })
}

Then, we solve the system using defined parameters and the function ode(). Here, we have to define the \(\alpha\) matrix with the competition coefficients:

\[ \left(\begin{array}{cc} \alpha_{11} & \alpha_{12}\\ \alpha_{21} & \alpha_{22} \end{array}\right) \left(\begin{array}{cc} 0.02 & 0.01\\ 0.01 & 0.03 \end{array}\right) \]

# LV parameters
a <- matrix(c(0.02, 0.01, 0.01, 0.03), nrow = 2)
r <- c(1, 1)
p2 <- list(r, a)
N0 <- c(10, 10)
t2 <- c(1:100)

Solving the system of ODE.

out_lv <- ode(y = N0, times = t2, func = LVComp, parms = p2)
head(out_lv)

     time        1        2
[1,]    1 10.00000 10.00000
[2,]    2 18.06979 15.99489
[3,]    3 26.20962 20.10821
[4,]    4 32.07829 21.66437
[5,]    5 35.53914 21.79929
[6,]    6 37.43160 21.46853

Notice that now we have two columns because we have two state variables.

We can again plot the solution. But first, we have to convert out data in a format in which every variable is represented in a column and every observation is represented in a row. We will use the function pivot_longer() from the package tidyr.

df_lv <- pivot_longer(as.data.frame(out_lv), cols = 2:3)

ggplot(df_lv) +
  geom_line(aes(x = time, y = value, color = name)) +
  labs(x = "Time", y = "N", color = "Species") +
  theme_classic()