Integrating Lotka Volterra Model into the Gillespie Algorithm

First, we visualize the model Lotka-Volterra Model

\[ \begin{cases} \frac{dx}{dt} = ax - bxy \\ \frac{dy}{dt} = dxy - cy \end{cases} \]

We sum all the terms contained in the model: \[ A = ax + bxy + dxy + cy \]

We assign values that together form a whole, called reactions \[ Rtotal = R1+R2+R3+R4 \]

We determine the reactions

\[X (ax) x + 1\] \[X (bxy) x - 1 \] \[Y(dxy) y +1\] \[Y(cy) y - 1 \]

What are the events involved?

Prey birth \[ax/A \] Death of prey by predation \[ bxy/A \] Conversion of prey to predator \[dxy/A \] Natural death of predator \[ cy/A \]

Implementing the Gillespie algorithm for the Lotka-Volterra model.

The necessary library

library(ggplot2)
library(dplyr)
library(magrittr)
library(tidyr)

Creating the function, initialization state, sum of reactions and data frame to store the information

LotkaVolterraGillespie <- function(a, b, c,  d, X0, Y0, t_end) {
  
  # Initialize state
  X <- X0  
  Y <- Y0  
  t <- 0   
  
  all.reactions <- a+b+c+d
  
  dataPP <- data.frame(Time = c(t), 
                       Prey = c(X), 
                       Predator = c(Y))
  
 # Lotka-Volterra into Gillespie simulation loop 
  
  while (t < t_end && X > 0 && Y > 0) {
    R1 <- a * X                   # Prey birth      
    R2 <- b * X * Y               # Prey death by predation 
    R3 <- d * X * Y               # Predator birth by predation
    R4 <- c * Y                   # Predator natural death     
    R_total <- R1 + R2 + R3 + R4  # The sum of the reactions
    
    if (R_total == 0) break       # If there is no more interactions
    tau <- rexp(1, R_total)       # Generating time until the next event
    event <- runif(1, 0, R_total) # Choosing which reaction will occur
    
    # When which event will occur? 
    
    if (event < R1) { 
      X <- X + 1                   # Prey birth 
    } else 
      if (event < R1 + R2) {
      X <- X - 1                   # Prey death by predation
    } else 
      if (event < R1 + R2 + R3) {
      Y <- Y + 1                   # Predator birth by predation
    } else {
      Y <- Y - 1                   # Predator natural death
    }
    
    t <- t + tau                   # Updating time
    dataPP <- rbind(dataPP, c(t, X, Y)) # Storing results
  }
  return(dataPP)
}

Choosing the parameters

set.seed(124)     # random number sequence
a = 1.0           # maximum prey birth
b = 0.02          # low predation
c = 1.0           # maximum predator birth
d = 0.01          # low natural death
X0 = 50           # initial number of prey
Y0 = 20           # initial number of predator
t_end = 100       # when the time stops

simulation <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end)   # same parameters
simulation1 <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end)  # same parameters

Plotting simulation

plot1 <- ggplot(simulation, aes(x = Time)) +
  geom_point(aes(y = Prey, color = "Prey"), size = 2) +
  geom_point(aes(y = Predator, color = "Predator"), size = 2) +
  labs(title = "Lotka-Volterra in Gillespie Simulation by Amanda Loureiro", 
       x = "Time", y = "Population") +
  scale_color_manual(values = c("#caf085", "#FFC0CB"),
                     name = "Animals",  
                     labels = c("Prey", "Predator")) +
  theme_minimal()+
  theme(
      legend.position = "bottom")

So far, we have implemented the Lotka-Volterra model in the Gillespie algorithm, but let’s continue and explore the results.

Viewing the plot of the first simulation

In this graph, it is possible to see that populations become extinct before the end time. However, even before they become extinct, they fluctuate equally, but with the predator population always greater than the prey population.

Plotting the second simulation

plot2 <- ggplot(simulation1, aes(x = Time)) +
  geom_point(aes(y = Prey, color = "Prey"), size = 2) +
  geom_point(aes(y = Predator, color = "Predator"), size = 2) +
  labs(title = "Lotka-Volterra in Gillespie Simulation by Amanda Loureiro", 
       x = "Time", y = "Population") +
  scale_color_manual(values = c("#caf085", "#FFC0CB"),
                     name = "Animals",  
                     labels = c("Prey", "Predator")) +
  theme_minimal()+
   theme(
      legend.position = "bottom",
      legend.text = element_text(size = 12, face = "italic")+
  guides(fill = guide_legend(override.aes = list(size = 5))))

Viewing the plot of the second simulation

We have the same observations as the previous graph, except that extinction occurs much earlier.

Plotting the combination of two simulations

combination <- ggplot() +
  geom_point(data = simulation, aes(x = Time, y = Prey), color = "#66CDAA") +
  geom_point(data = simulation1, aes(x = Time, y = Prey), color = "#caf085")+
  geom_point(data = simulation1, aes(x = Time, y = Predator), color = "#DE3163")+
    geom_point(data = simulation, aes(x = Time, y = Predator), color = "#FFC0CB")+
  labs(title = "Stochastic Lotka-Volterra Model (Gillespie Algorithm)", x = "Time", y = "Population Size") +
  theme_minimal()

Prey in shades of green Predator in shades of pink

Viewing the plot of the two combined simulations

Just a way to see that even with the same parameters, two simulations can be totally different due to stochasticity.

Multiple replications (Monte Carlo)

So far, we have done 1 simulation. But what if we did countless simulations? To do this, we use the Monte Carlo method, with multiple replications.

# Defining parameters and initial conditions

a <- 0.3      # Low prey growth 
b <- 0.02     # Very low death by predation 
c <- 0.3      # Low predator reproduction rate
d <- 0.01     # Very low predator mortality rate
X0 <- 10      # Initial prey population size
Y0 <- 10      # Initial predator population size
t_end <- 200  # Simulation time

set.seed(123) # Random number sequence
  
ntimes <- 100  # Number of simulations
final_pop_prey <- c()
final_pop_pred <- c()
extinctions_prey <- 0
extinctions_predator <- 0
simulations_list <- list()

Running the simulations

for (i in 1:ntimes) {
  nsimulations <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end)
  nsimulations$Simulation <- i
  simulations_list[[i]] <- nsimulations
  
# Capturing final population sizes
  
  final_pop_prey[i] <- nsimulations$Prey[nrow(nsimulations)]  # Final prey population
  final_pop_pred[i] <- nsimulations$Predator[nrow(nsimulations)]  # Final predator population
  
# Checking for extinction (if population reaches zero)
  if (final_pop_prey[i] == 0) {
    extinctions_prey <- extinctions_prey + 1
  }
  if (final_pop_pred[i] == 0) {
    extinctions_predator <- extinctions_predator + 1
  }
}

Viewing the average of the simulations

all_simulations <- bind_rows(simulations_list)
time_grid <- seq(0, t_end, by = 1)
interpolated_data <- all_simulations %>%
  group_by(Simulation) %>%
  complete(Time = time_grid) %>%
  fill(Prey, Predator, .direction = "down")  

mean_data <- interpolated_data %>%
  group_by(Time) %>%
  summarize(MeanPrey = mean(Prey, na.rm = TRUE), 
            MeanPredator = mean(Predator, na.rm = TRUE))

Plotting average over time

ggplot(mean_data, aes(x = Time)) +
  geom_line(aes(y = MeanPrey, color = "Prey"), size=1) +
  geom_line(aes(y = MeanPredator, color = "Predator"), size=1) +
  labs(title = "Simulations Averages", x = "Time", y = "Average Population") +
  scale_color_manual(values = c("Prey" = "#caf085", "Predator" = "#DE3163"),
                     name = "Animals") +
  theme_minimal()

Here we can see the average of several simulations with different parameters than the previous ones. The joint oscillation continues, but with the number of prey greater than the number of predators. However, as in the previous simulations, coexistence is not eternal and extinction occurs before the end time.

Calculating extinction probabilities

extinction_prob_prey <- extinctions_prey / ntimes
extinction_prob_predator <- extinctions_predator / ntimes

Outputting the extinction probabilities

cat("Extinction probability for prey:", extinction_prob_prey, "\n")

Extinction probability for prey: 0.4

cat("Extinction probability for predator:", extinction_prob_predator, "\n")

Extinction probability for predator: 0.6

In this scenario, the probability of extinction of predators is higher than that of prey.

Histogram for final prey population

df_prey <- data.frame(Final_Population = final_pop_prey)

plot.prey <- ggplot(df_prey, aes(x = factor(Final_Population))) +
  geom_bar(fill = "#caf085", alpha = 0.7) +
  labs(title = "Distribution of Final Prey Populations", x = "Final Prey Population", y = "Frequency") +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  ) +
  scale_x_discrete(breaks = seq(1, length(unique(df_prey$Final_Population)), by = 2))

Histogram for final predator population

Only to realize that, in the end, most individuals in the final populations are zero, indicating the probability of extinction after all the simulations.

An image to illustrate the problem

The end.

--- title: "*Activity II - Mathematical Models Applied to Ecology*" author: "Amanda Beatriz Loureiro" date: "03 Mar 2025" format: html editor: visual cap-location: bottom title-block-banner: true code-tools: true theme: custom.scss toc: true toc-title: "Sections" editor_options: chunk_output_type: console --- ### Integrating Lotka Volterra Model into the Gillespie Algorithm *First, we visualize the model Lotka-Volterra Model* $$ \begin{cases} \frac{dx}{dt} = ax - bxy \\ \frac{dy}{dt} = dxy - cy \end{cases} $$ *We sum all the terms contained in the model:* $$ A = ax + bxy + dxy + cy $$ *We assign values that together form a whole, called reactions* $$ Rtotal = R1+R2+R3+R4 $$ *We determine the reactions* $$X (ax) x + 1$$ $$X (bxy) x - 1 $$ $$Y(dxy) y +1$$ $$Y(cy) y - 1 $$ ::: callout-tip ## What are the events involved? Prey birth $$ax/A $$ Death of prey by predation $$ bxy/A $$ Conversion of prey to predator $$dxy/A $$ Natural death of predator $$ cy/A $$ ::: ### Implementing the Gillespie algorithm for the Lotka-Volterra model. *The necessary library* ```{r} #| warning: false library(ggplot2) library(dplyr) library(magrittr) library(tidyr) ``` ### Creating the function, initialization state, sum of reactions and data frame to store the information ```{r} LotkaVolterraGillespie <- function(a, b, c, d, X0, Y0, t_end) { # Initialize state X <- X0 Y <- Y0 t <- 0 all.reactions <- a+b+c+d dataPP <- data.frame(Time = c(t), Prey = c(X), Predator = c(Y)) # Lotka-Volterra into Gillespie simulation loop while (t < t_end && X > 0 && Y > 0) { R1 <- a * X # Prey birth R2 <- b * X * Y # Prey death by predation R3 <- d * X * Y # Predator birth by predation R4 <- c * Y # Predator natural death R_total <- R1 + R2 + R3 + R4 # The sum of the reactions if (R_total == 0) break # If there is no more interactions tau <- rexp(1, R_total) # Generating time until the next event event <- runif(1, 0, R_total) # Choosing which reaction will occur # When which event will occur? if (event < R1) { X <- X + 1 # Prey birth } else if (event < R1 + R2) { X <- X - 1 # Prey death by predation } else if (event < R1 + R2 + R3) { Y <- Y + 1 # Predator birth by predation } else { Y <- Y - 1 # Predator natural death } t <- t + tau # Updating time dataPP <- rbind(dataPP, c(t, X, Y)) # Storing results } return(dataPP) } ``` ### Choosing the parameters ```{r} set.seed(124) # random number sequence a = 1.0 # maximum prey birth b = 0.02 # low predation c = 1.0 # maximum predator birth d = 0.01 # low natural death X0 = 50 # initial number of prey Y0 = 20 # initial number of predator t_end = 100 # when the time stops simulation <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end) # same parameters simulation1 <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end) # same parameters ``` ## Plotting simulation ```{r} plot1 <- ggplot(simulation, aes(x = Time)) + geom_point(aes(y = Prey, color = "Prey"), size = 2) + geom_point(aes(y = Predator, color = "Predator"), size = 2) + labs(title = "Lotka-Volterra in Gillespie Simulation by Amanda Loureiro", x = "Time", y = "Population") + scale_color_manual(values = c("#caf085", "#FFC0CB"), name = "Animals", labels = c("Prey", "Predator")) + theme_minimal()+ theme( legend.position = "bottom") ``` *So far, we have implemented the Lotka-Volterra model in the Gillespie algorithm, but let's continue and explore the results.* ### Viewing the plot of the first simulation ```{r, warning=FALSE, echo=FALSE} plot1 ``` *In this graph, it is possible to see that populations become extinct before the end time. However, even before they become extinct, they fluctuate equally, but with the predator population always greater than the prey population.* ### Plotting the second simulation ```{r} plot2 <- ggplot(simulation1, aes(x = Time)) + geom_point(aes(y = Prey, color = "Prey"), size = 2) + geom_point(aes(y = Predator, color = "Predator"), size = 2) + labs(title = "Lotka-Volterra in Gillespie Simulation by Amanda Loureiro", x = "Time", y = "Population") + scale_color_manual(values = c("#caf085", "#FFC0CB"), name = "Animals", labels = c("Prey", "Predator")) + theme_minimal()+ theme( legend.position = "bottom", legend.text = element_text(size = 12, face = "italic")+ guides(fill = guide_legend(override.aes = list(size = 5)))) ``` ### Viewing the plot of the second simulation ```{r, warning=FALSE, echo=FALSE} plot2 ``` *We have the same observations as the previous graph, except that extinction occurs much earlier.* # Plotting the combination of two simulations ```{r} combination <- ggplot() + geom_point(data = simulation, aes(x = Time, y = Prey), color = "#66CDAA") + geom_point(data = simulation1, aes(x = Time, y = Prey), color = "#caf085")+ geom_point(data = simulation1, aes(x = Time, y = Predator), color = "#DE3163")+ geom_point(data = simulation, aes(x = Time, y = Predator), color = "#FFC0CB")+ labs(title = "Stochastic Lotka-Volterra Model (Gillespie Algorithm)", x = "Time", y = "Population Size") + theme_minimal() ``` *Prey in shades of green* *Predator in shades of pink* ### Viewing the plot of the two combined simulations ```{r} #| warning: false #| echo: false combination ``` *Just a way to see that even with the same parameters, two simulations can be totally different due to stochasticity.* # Multiple replications (Monte Carlo) *So far, we have done 1 simulation. But what if we did countless simulations? To do this, we use the Monte Carlo method, with multiple replications.* ```{r} # Defining parameters and initial conditions a <- 0.3 # Low prey growth b <- 0.02 # Very low death by predation c <- 0.3 # Low predator reproduction rate d <- 0.01 # Very low predator mortality rate X0 <- 10 # Initial prey population size Y0 <- 10 # Initial predator population size t_end <- 200 # Simulation time set.seed(123) # Random number sequence ntimes <- 100 # Number of simulations final_pop_prey <- c() final_pop_pred <- c() extinctions_prey <- 0 extinctions_predator <- 0 simulations_list <- list() ``` # Running the simulations ```{r} for (i in 1:ntimes) { nsimulations <- LotkaVolterraGillespie(a, b, c, d, X0, Y0, t_end) nsimulations$Simulation <- i simulations_list[[i]] <- nsimulations # Capturing final population sizes final_pop_prey[i] <- nsimulations$Prey[nrow(nsimulations)] # Final prey population final_pop_pred[i] <- nsimulations$Predator[nrow(nsimulations)] # Final predator population # Checking for extinction (if population reaches zero) if (final_pop_prey[i] == 0) { extinctions_prey <- extinctions_prey + 1 } if (final_pop_pred[i] == 0) { extinctions_predator <- extinctions_predator + 1 } } ``` # Viewing the average of the simulations ```{r} all_simulations <- bind_rows(simulations_list) time_grid <- seq(0, t_end, by = 1) interpolated_data <- all_simulations %>% group_by(Simulation) %>% complete(Time = time_grid) %>% fill(Prey, Predator, .direction = "down") mean_data <- interpolated_data %>% group_by(Time) %>% summarize(MeanPrey = mean(Prey, na.rm = TRUE), MeanPredator = mean(Predator, na.rm = TRUE)) ``` # Plotting average over time ```{r, warning=FALSE} ggplot(mean_data, aes(x = Time)) + geom_line(aes(y = MeanPrey, color = "Prey"), size=1) + geom_line(aes(y = MeanPredator, color = "Predator"), size=1) + labs(title = "Simulations Averages", x = "Time", y = "Average Population") + scale_color_manual(values = c("Prey" = "#caf085", "Predator" = "#DE3163"), name = "Animals") + theme_minimal() ``` *Here we can see the average of several simulations with different parameters than the previous ones. The joint oscillation continues, but with the number of prey greater than the number of predators. However, as in the previous simulations, coexistence is not eternal and extinction occurs before the end time.* ### Calculating extinction probabilities ```{r} extinction_prob_prey <- extinctions_prey / ntimes extinction_prob_predator <- extinctions_predator / ntimes ``` ### Outputting the extinction probabilities ```{r} cat("Extinction probability for prey:", extinction_prob_prey, "\n") cat("Extinction probability for predator:", extinction_prob_predator, "\n") ``` *In this scenario, the probability of extinction of predators is higher than that of prey.* ### Histogram for final prey population ```{r} df_prey <- data.frame(Final_Population = final_pop_prey) plot.prey <- ggplot(df_prey, aes(x = factor(Final_Population))) + geom_bar(fill = "#caf085", alpha = 0.7) + labs(title = "Distribution of Final Prey Populations", x = "Final Prey Population", y = "Frequency") + theme_minimal() + theme( panel.grid.major = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1) ) + scale_x_discrete(breaks = seq(1, length(unique(df_prey$Final_Population)), by = 2)) ``` ```{r, echo=FALSE} #| warning: false #| echo: false plot.prey ``` # Histogram for final predator population ```{r, echo=FALSE} df_pred <- data.frame(Final_Population = final_pop_pred) plot.predator <- ggplot(df_pred, aes(x = factor(Final_Population))) + geom_bar(fill = "#FFC0CB", alpha = 0.7) + labs(title = "Distribution of Final Predator Populations", x = "Final Predator Population", y = "Frequency") + theme_minimal() + theme( panel.grid.major = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text(angle = 45, hjust = 1) ) + scale_x_discrete(breaks = seq(1, length(unique(df_pred$Final_Population)), by = 2)) ``` ```{r} #| warning: false #| echo: false plot.predator ``` *Only to realize that, in the end, most individuals in the final populations are zero, indicating the probability of extinction after all the simulations.* ## An image to illustrate the problem ![](images/DALL·E%202025-02-25%2011.28.16%20-%20A%20natural%20scene%20depicting%20a%20predator%20and%20prey%20in%20the%20wild.%20The%20predator%20is%20a%20large,%20muscular%20wolf%20with%20a%20sleek,%20dark%20fur%20coat,.webp) ### The end.