Show the code
# import packages
library(tidyverse)
library(deSolve)
library(phaseR)
library(ggquiver)
# set the plot theme
theme_set(theme_bw(base_size = 13))Differential Equations and Dynamical Systems in R
Examples for Biomathematics
Example: A Single Species Population Model
Consider the logistic growth model for a single population:
\[ \frac{dN}{dt} = rN\left(1-\frac{N}{K}\right) \]
where \(N\) is the population size, \(r\) is the intrinsic growth rate, and \(K\) is the carrying capacity. It can be shown that one can choose a scaling so that \(r=1\) and \(K=1\). Then the model becomes
\[ \frac{dN}{dt} = N\left(1-N\right) \]
which is a scalar autonomous differential equation. We will use R to obtain a numerical solution to an initial value problem for this equation and then plot the solution. Further, since the equation is autonomous, we can plot the phase line and use it to classify the equilibrium solutions. We will do this both from scratch and using the phaseR package.
Numerical Solution for IVP
We will use the deSolve package to obtain a numerical solution to the initial value problem
\[ \frac{dN}{dt} = N\left(1-N\right), \quad N(0)=N_0 \]
for \(N_0=0.1\). We first define the right-hand side of the differential equation as a function of \(t\) and \(N\).
# define the right-hand side of the differential equation
rhs <- function(t, N, parms) {
list(N * (1 - N))
}
# define initial value
N0 <- c(N=0.1)
# define time interval
t0 <- 0
tend <- 10
# define time points for output
times <- seq(t0, tend, by = 0.1)
# solve the IVP
sol <- ode(y = N0, times = times, func = rhs)
# print the first few solution value
sol[1:5, ] time N
[1,] 0.0 0.1000000
[2,] 0.1 0.1093670
[3,] 0.2 0.1194954
[4,] 0.3 0.1304233
[5,] 0.4 0.1421895Now, let’s plot the solution. To do so, we will use the ggplot2 package. Note that this requires us to format out solution as a data frame
# format the solution as a data frame
sol_df <- data.frame(sol)
# plot the solution
ggplot(sol_df, aes(x = time, y = N)) +
geom_line(linewidth=1) +
labs(x = "time", y = "N")
Phase Line: Manual Approach
We can plot the phase line for the differential equation by hand as follows:
Show the code
logistic_growth <- function(N) {
N * (1 - N)
}
phase_line_data <- data.frame(N = c(-0.5,1.5))
ggplot(phase_line_data, aes(x = N)) +
geom_function(fun=logistic_growth,color = "steelblue", linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
geom_vline(xintercept = 1, linetype = "dashed", color = "black") +
geom_segment(aes(x = 0.9, y = 0, xend = 1, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "purple") +
geom_segment(aes(x = 0.7, y = 0, xend = 0.8, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "purple") +
geom_segment(aes(x = 0.5, y = 0, xend = 0.6, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "purple") +
geom_segment(aes(x = 0.3, y = 0, xend = 0.4, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "purple") +
geom_segment(aes(x = 0.1, y = 0, xend = 0.2, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "purple") +
geom_segment(aes(x = 1.2, y = 0, xend = 1.1, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "orange") +
geom_segment(aes(x = 1.4, y = 0, xend = 1.3, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "orange") +
geom_segment(aes(x = -0.1, y = 0, xend = -0.2, yend = 0),
arrow = arrow(length = unit(0.3, "cm")), color = "orange") +
geom_point(aes(x = 0, y = 0), color = "darkred", size = 4, shape = 1) +
geom_point(aes(x = 1, y = 0), color = "darkred", size = 4) +
labs(title = "Logistic Growth Phase Line",
x = "Population",
y = "Population Rate of Change",
caption = "Phase line for logistic growth") 
It is also useful to plot the slope field for the differential equation in the \(t,y\) - plane. We can do this as follows:
Show the code
# define the right-hand side of the differential equation
rhs_vf <- function(t,N) {
return(N * (1 - N))
}
expand.grid(t=seq(0,5,0.5), N=seq(-0.5,1.5,0.1)) |>
ggplot(aes(x=t,y=N)) +
geom_quiver(aes(u=t,v=rhs_vf(t,N)),color="darkgrey", vecsize = 5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_hline(yintercept = 1, linetype = "dashed", color = "blue") +
labs(x = "t",
y = "N") 
Let’s add solutions corresponding to \(N(0) = 0.1\) and \(N(0) = 0.5\) to the vector field plot:
Show the code
# define initial value
N0_05 <- c(N=0.5)
# solve the IVP
sol_05 <- ode(y = N0_05, times = times, func = rhs)
# format the solution as a data frame
sol_df_05 <- data.frame(sol_05)
sol_df_f <- sol_df |>
filter(time <= 5)
sol_df_05_f <- sol_df_05 |>
filter(time <= 5.1)
expand.grid(t=seq(0,5,0.5), N=seq(-0.5,1.5,0.1)) |>
ggplot(aes(x=t,y=N)) +
geom_quiver(aes(u=t,v=rhs_vf(t,N)),color="darkgrey", vecsize = 5) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_hline(yintercept = 1, linetype = "dashed", color = "blue") +
geom_line(data=sol_df_f,aes(x=time,y=N),linewidth=1) +
geom_line(data=sol_df_05_f,aes(x=time,y=N),linewidth=1) +
labs(x = "t",
y = "N") 
Phase Line: Using phaseR
We can also use the phaseR package to plot the phase line:
Show the code
logistic_phasePortrait <- phasePortrait(rhs,
ylim = c(-0.5, 1.5),
points = 10,
frac = 0.5,
state.names = "N")
Further, we can use phaseR to plot the solutions to the IVP:
Show the code
logistic_flow <- flowField(rhs,
xlim = c(0, 10),
ylim = c(-0.5, 1.5),
parameters = NULL,
points = 15,
system = "one.dim",add=FALSE,
state.names = "N")
logistic_solutions <- trajectory(rhs,
y0 = c(0,0.5,1,1.5),
tlim=c(0,10),
system = "one.dim")
Example: Interacting Species Population Model
The Lotka-Volterra model is a model of predator-prey interactions. It is a system of two differential equations:
\[ \begin{aligned} \frac{dN}{dt} &= rN - \alpha N P \\ \frac{dP}{dt} &= \beta N P - \gamma P \end{aligned} \]
where \(N\) is the population of prey, \(P\) is the population of predators, \(r\) is the growth rate of the prey population, \(\alpha\) is the rate at which predators kill prey, \(\beta\) is the rate at which predators reproduce, and \(\gamma\) is the rate at which predators die. This system can be non-dimensionalized to:
\[ \begin{aligned} \frac{dx}{dt} &= x - xy = x(1 - y) \\ \frac{dy}{dt} &= \rho x y - \rho y = \rho y(x - 1) \end{aligned} \]
For this non-dimensionalized system, we will:
Plot the corresponding vector field including the nullclines for a particular value of \(\rho\).
Compute a numerical solution to the IVP for a particular initial condition and value of \(rho\). Then, we will plot this solution as a time series. Also, we will add the solution to the vector field plot.
Use
phaseRto plot a phase portrait for a particular value of \(\rho\).
Vector Field and Nullclines
Figure 7 shows the vector field for the Lotka-Volterra model with nullclines. Here, \(\rho = 0.5\).
colors <- c("dx/dt=0" = "#551A8BFF", "dy/dt=0" = "#FF7F00FF")
rho_val <- 0.5
lv_vf <- expand.grid(x=seq(0,3,0.15), y=seq(0,2,0.15)) |>
ggplot(aes(x=x,y=y)) +
geom_quiver(aes(u=x - x*y,v= rho_val*x*y - rho_val*y),color="darkgray") +
geom_hline(aes(color="dx/dt=0",yintercept = 0.0),linetype = "dashed") +
geom_vline(aes(color="dx/dt=0",xintercept = 1.0),linetype = "dashed") +
geom_vline(aes(color="dy/dt=0",xintercept = 0.0),linetype = "dashed") +
geom_hline(aes(color="dy/dt=0",yintercept = 1.0),linetype = "dashed") +
labs(x = "x",
y = "y",
color = "Null-cline") +
scale_color_manual(values = colors)
lv_vf
Numerical Solution of IVP
Now, to compute numerical solutions we need a function to define the right-hand side of the system of ODEs:
Show the code
lotka_volterra_rhs <- function(t, y, parms) {
x <- y[1]
y <- y[2]
rho <- parms[1]
dx <- x - x*y
dy <- rho*x*y - rho*y
return(list(c(dx, dy)))
}We can use the deSolve package to solve the IVP:
Show the code
parms <- c(rho = 0.5)
y0 <- c(x = 0.5, y = 0.5)
times <- seq(0, 25, by = 0.1)
lotka_volterra_sol <- ode(y = y0, times = times, func = lotka_volterra_rhs, parms = parms)
lotka_volterra_sol_df <- as.data.frame(lotka_volterra_sol)
lotka_volterra_sol_df |>
ggplot(aes(x = time)) +
geom_line(aes(y = x, color = "Prey"),linewidth=1) +
geom_line(aes(y = y, color = "Predator"),linewidth=1) +
labs(x = "Time",
y = "Population",
color = "Species") +
scale_color_manual(values = c("Prey" = "steelblue", "Predator" = "purple"))
Figure 8 shows the numerical solution to the Lotka-Volterra model for a particular initial condition and value of \(\rho=0.5\). Let’s add this solution to the vector field plot.
Show the code
lv_vf +
geom_path(data = lotka_volterra_sol_df, aes(x = x, y = y, color = "Solution"), linewidth = 1, color="black") 
Figure 9 shows the numerical solution we obtained for the Lotka-Volterra system plotted in the phase plane.
Phase Portrait
Finally, we can use the phaseR package to plot the phase portrait for the Lotka-Volterra model with \(\rho = 0.5\). Figure 10 shows the phase portrait for the Lotka-Volterra model obtained using the phaseR package.
Code
system_p_flowField <- flowField(lotka_volterra_rhs,
xlim = c(0, 3),
ylim = c(0, 2),
parameters = c(rho=rho_val),
points = 19,
add = FALSE)
system_p_nullclines <- nullclines(lotka_volterra_rhs,
xlim = c(0, 3),
ylim = c(0, 2),
parameters = c(rho=rho_val),
points = 500)
state <- matrix(c(0.5,0.5,0.4,0.7,0.4,0.9),
3, 2, byrow = TRUE)
system_p_trajectory <- trajectory(lotka_volterra_rhs,
y0 = state,
tlim = c(0, 25),
parameters = c(rho=rho_val),add=TRUE)
Expand for Session Info
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.3.2 (2023-10-31)
os macOS Sonoma 14.3.1
system aarch64, darwin20
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2024-02-23
pandoc 3.1.12.1 @ /opt/homebrew/bin/ (via rmarkdown)
quarto 1.4.550 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
deSolve * 1.40 2023-11-27 [1] CRAN (R 4.3.1)
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
ggplot2 * 3.4.4 2023-10-12 [1] CRAN (R 4.3.1)
ggquiver * 0.3.3 2023-11-17 [1] CRAN (R 4.3.1)
lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.1)
phaseR * 2.2.1 2022-09-02 [1] CRAN (R 4.3.0)
purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.0)
readr * 2.1.5 2024-01-10 [1] CRAN (R 4.3.1)
sessioninfo * 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.1)
tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.3.1)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
[1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
──────────────────────────────────────────────────────────────────────────────