Animation demonstrating the mapping between time series and phase plane representations of the Lotka-Volterra predator prey equations

phase_planes.R

# from https://github.com/ha0ye/personal-website/blob/master/figure%20code/causality-figures.R

library(tidyverse)
library(gganimate) # devtools::install_github("thomasp85/gganimate")

## generate model data ----

# params
t_max <- 100
h <- 0.01 # step size
num_points <- t_max / h
a <- 0.02 # search efficiency
r <- 0.5 # prey growth rate
m <- 0.35 # predator death rate
c <- 0.1 # conversion rate

# motion equations
dx <- function(x1, x2)
{
  r*x1 - a*x1*x2
}
dy <- function(x1, x2)
{
  a*c*x1*x2 - m*x2
}

# initialize data 
x <- vector("numeric", length = num_points)
y <- vector("numeric", length = num_points)
x[1] = 100
y[1] = 10

# generate data using Runge-Kutta
for (i in 1:(num_points - 1))
{
  xx = x[i];
  yy = y[i];
  kx1 = dx(xx, yy);
  ky1 = dy(xx, yy);
  kx2 = dx(xx + h/2*kx1, yy + h/2*ky1);
  ky2 = dy(xx + h/2*kx1, yy + h/2*ky1);
  kx3 = dx(xx + h/2*kx2, yy + h/2*ky2);
  ky3 = dy(xx + h/2*kx2, yy + h/2*ky2);
  kx4 = dx(xx + h*kx3, yy + h*ky3);
  ky4 = dy(xx + h*kx3, yy + h*ky3);

  x[i + 1] = xx + h / 6.0 * (kx1 + 2*kx2 + 2*kx3 + kx4);
  y[i + 1] = yy + h / 6.0 * (ky1 + 2*ky2 + 2*ky3 + ky4);
}

# rescale data to [0, 1]
x <- (x - min(x)) / (max(x) - min(x))
y <- (y - min(y)) / (max(y) - min(y))

## preview animation ----
if (FALSE) 
{
  num_frames <- 100
  t_vec <- 1:10000
  frame_idx <- round(0.00001 + seq(from = 0, to = num_frames, length.out = length(t_vec)))
  
  # time series
  to_plot <- data.frame(x = x[t_vec] * 0.3 + 0.65, 
                        y = y[t_vec] * 0.03 + 0.65, 
                        t = t_vec, 
                        frame = frame_idx)
  
  # arrows
  arrows <- data.frame(x = to_plot$x, 
                       t = t_vec,
                       xend = (to_plot$x-0.3)*10000, 
                       y = to_plot$y,
                       yend = (to_plot$y-0.6)*5, 
                       seg_idx = 1:50)

  xlims <- range(to_plot$t)
  ylims <- c(0, 1)
  
  p <- ggplot(to_plot, 
              aes(x = t, y = x)) + 
    scale_y_continuous(limits = ylims, expand = c(0.01, 0)) + 
    scale_x_continuous(limits = xlims, expand = c(0.01, 0)) +
    
    # arrow 
    geom_segment(aes(x = t, xend = xend, y = x, yend = 0.2),
                 data = arrows, size = 0.5, color = "red", linetype = 2) +
    geom_segment(aes(x = t, xend = 3400, y = y, yend = yend),
                 data = arrows, size = 0.5, color = "blue", linetype =2) +
    # geom_segment(aes(x = x, xend = xend, y = y, yend = yend), 
    #              data = filter(segment_data, seg_idx == 2), size = 2) + 
    # geom_segment(aes(x = x, xend = xend, y = y, yend = yend), 
    #              data = filter(segment_data, seg_idx == 3), size = 2) + 
    # geom_segment(aes(x = x, xend = xend, y = y, yend = yend), 
    #              data = filter(segment_data, seg_idx == 4), size = 2) + 
    # geom_segment(aes(x = x, xend = xend, y = y, yend = yend), 
    #              data = filter(segment_data, seg_idx == 5), size = 2) + 
    # geom_segment(aes(x = x, xend = xend, y = y, yend = yend), 
    #              data = filter(segment_data, seg_idx == 6), size = 2) + 
    
    # time series for x and z
    geom_line(size = 1, color = "red") + 
    geom_line(aes(x = t, y = y), size = 1, color = "blue") + 
    geom_line(aes(x = (x-0.3)*10000, y = (y-0.6)*5)) + 
    geom_point(aes(x = (x-0.3)*10000, y = (y-0.6)*5), size = 2) + 
    theme_void() + 
    annotate("segment", x = 1, y = 0.6, xend = 9999, yend = 0.6, size = 2) + 
    annotate("segment", x = 1, y = 0.6, xend = 1, yend = 1, size = 2)  +
    annotate("segment", x = 3400, y = 0.2, xend = 6600, yend = 0.2, size = 2, color = "red") + 
    annotate("segment", x = 3400, y = 0.2, xend = 3400, yend = 0.45, size = 2, color = "blue") 
  
  animate(p + transition_reveal(1, t), 
          nframes = num_frames, 
          width = 400, height = 400)
  
  anim_save(here::here("R/phase_planes.gif"))
}