Tarefa 2 - Gillespie e Lotka-Volterra

Atividade 2: Algoritmo de Gillespie implementado ao modelo de Lotka-Volterra

Definir a função do algoritmo e as taxas dos eventos:

library(ggplot2)
Warning: package 'ggplot2' was built under R version 4.3.3
#definir a função do algoritmo de Gillespie para o modelo de Lotka-Volterra
gil_lkv <- function(a, b, c, d, x, y, t_end) {
  
  #definir o estado inicial
  t <- 0
  results <- data.frame(Time = t, Prey = x, Predator = y)
  
  #loop principal do Algoritmo de Gillespie
  while (t < t_end && (x > 0 || y > 0)) {
    
    #Definir taxas de eventos
    nasc_presa <- a * x  #taxa de nascimento das presas
    predacao <- b * x * y  #taxa de predação (presa consumida, predador nasce)
    nasc_predad <- c * x * y  #raxa de nascimento dos predadores
    morte_pred <- d * y  #taxa de morte dos predadores
    
    #calcular a taxa total
    taxa_total <- nasc_presa + predacao + nasc_predad + morte_pred
    
    #parar, se não houver mais interações possíveis
    if (taxa_total == 0) break
    
    #tempo até o próximo evento (distribuição exponencial)
    tau <- rexp(1, taxa_total)
    t <- t + tau  # Atualiza o tempo
    
    #escolher qual evento ocorre (probabilidade proporcional às taxas)
    r <- runif(1, 0, taxa_total)
    
    if (r < nasc_presa) {
      x <- x + 1  #nascimento de uma presa
    } else if (r < nasc_presa + predacao) {
      x <- x - 1  #predação (presa consumida)
    } else if (r < nasc_presa + predacao + nasc_predad) {
      y <- y + 1  #nascimento de um novo predador
    } else {
      y <- y - 1  #morte de um predador
    }
    
    #armazenar os valores num data frame
    results <- rbind(results, data.frame(Time = t, Prey = x, Predator = y))
  }
  
  return(results)
}

Adicionar os parâmetros e realizar a simulação:

#definir parâmetros do modelo
a <- 0.025  #taxa de nascimento das presas
b <- 0.025  #taxa de predação
c <- 0.025  #taxa de nascimento dos predadores 
d <- 0.025  #taxa de morte dos predadores 

#definir condições iniciais
x <- 30  #pop inicial de presas
y <- 100  #pop inicial de predadores


t_end <- 50  #tempo máximo da simulação

set.seed(124)
#executar a simulação
sim_data <- gil_lkv(a, b, c, d, x, y, t_end)

#visualizar os primeiros valores da simulação
head(sim_data)
        Time Prey Predator
1 0.00000000   30      100
2 0.01571086   29      100
3 0.01591710   28      100
4 0.02701302   27      100
5 0.02841659   27      101
6 0.03211574   27      102
#grafico das populações ao longo do tempo
ggplot(sim_data, aes(x = Time)) +
  geom_line(aes(y = Prey, color = "Presas")) +
  geom_line(aes(y = Predator, color = "Predadores")) +
  labs(title = "Simulação de Lotka-Volterra com o Algoritmo de Gillespie",
       x = "Tempo", y = "População") +
  scale_color_manual(name = "Espécies", values = c("Presas" = "blue", "Predadores" = "red")) +
  theme_minimal()

#definir parâmetros do modelo
a <- 0.025  #taxa de nascimento das presas
b <- 0.025  #taxa de predação
c <- 0.025  #taxa de nascimento dos predadores 
d <- 0.025  #taxa de morte dos predadores 

#definir condições iniciais
x <- 100  #pop inicial de presas
y <- 30  #pop inicial de predadores


t_end <- 50  #tempo máximo da simulação

set.seed(124)
#executar a simulação
sim_data <- gil_lkv(a, b, c, d, x, y, t_end)

#visualizar os primeiros valores da simulação
head(sim_data)
        Time Prey Predator
1 0.00000000  100       30
2 0.01571086   99       30
3 0.01591234   98       30
4 0.02649114   97       30
5 0.02779559   97       31
6 0.03115746   97       32
#criar gráfico da evolução das populações ao longo do tempo
ggplot(sim_data, aes(x = Time)) +
  geom_line(aes(y = Prey, color = "Presas")) +
  geom_line(aes(y = Predator, color = "Predadores")) +
  labs(title = "Simulação de Lotka-Volterra com o Algoritmo de Gillespie",
       x = "Tempo", y = "População") +
  scale_color_manual(name = "Espécies", values = c("Presas" = "blue", "Predadores" = "red")) +
  theme_minimal()

  1. Adotando os mesmos parâmetros para as taxas e alterando apenas população inicial de presas e predadores, vemos que para 50 unidade de tempo ocorre a extinção das presas e queda acentuada da população de predadores - que ainda são mantidos devido a baixa taxa de mortalidade.

  • Se a baixa → pode levar a redução da população de presas;

  • Se b elevada → predação ocorre rápido;

  • Se c elevada → mais predadores consumindo presas;

  • Se d baixa → predadores vivem mais tempo, consomem mais presas.

#definir parâmetros do modelo
a <- 0.006  #taxa de nascimento das presas
b <- 0.0025  #taxa de predação
c <- 0.001  #taxa de nascimento dos predadores 
d <- 0.07  #taxa de morte dos predadores 

#definir condições iniciais
x <- 40  #pop inicial de presas
y <- 15  #pop inicial de predadore


t_end <- 50  #tempo máximo da simulação

set.seed(124)
#executar a simulação
sim_data <- gil_lkv(a, b, c, d, x, y, t_end)

#visualizar os primeiros valores da simulação
head(sim_data)
       Time Prey Predator
1 0.0000000   40       15
2 0.7102326   39       15
3 0.7194087   38       15
4 1.2048763   37       15
5 1.2652089   37       14
6 1.4364248   37       13
#criar gráfico da evolução das populações ao longo do tempo
ggplot(sim_data, aes(x = Time)) +
  geom_line(aes(y = Prey, color = "Presas")) +
  geom_line(aes(y = Predator, color = "Predadores")) +
  labs(title = "Simulação de Lotka-Volterra com o Algoritmo de Gillespie",
       x = "Tempo", y = "População") +
  scale_color_manual(name = "Espécies", values = c("Presas" = "blue", "Predadores" = "red")) +
  theme_minimal()

  1. A elevada taxa de morte dos predadores e a população inicial de presas maior que a de predadores promove a coexistência por alguns anos.

  2. A população de presas cai e com ela a de predadores (influenciado pela baixa taxa de nascimento), levando a população de predadores próximo da extinção.

#definir parâmetros do modelo
a <- 0.5  #taxa de nascimento das presas
b <- 0.02  #taxa de predação
c <- 0.008  #taxa de nascimento dos predadores 
d <- 0.5  #taxa de morte dos predadores 

#definir condições iniciais
x <- 20  #pop inicial de presas
y <- 18  #pop inicial de predadore


t_end <- 50

set.seed(124)
sim_data <- gil_lkv(a, b, c, d, x, y, t_end)

head(sim_data)
        Time Prey Predator
1 0.00000000   20       18
2 0.08279534   19       18
3 0.08388417   18       18
4 0.14257713   17       18
5 0.15001686   17       17
6 0.17058979   17       16
#criar gráfico da evolução das populações ao longo do tempo
ggplot(sim_data, aes(x = Time)) +
  geom_line(aes(y = Prey, color = "Presas")) +
  geom_line(aes(y = Predator, color = "Predadores")) +
  labs(title = "Simulação de Lotka-Volterra com o Algoritmo de Gillespie",
       x = "Tempo", y = "População") +
  scale_color_manual(name = "Espécies", values = c("Presas" = "blue", "Predadores" = "red")) +
  theme_minimal()

  1. A população de presas e predadores aumenta e diminui simultâneamente, devido a baixa taxa de nascimentos de predadores e taxas semelhantes de mortalidade de predadores e de nascimento de presas.

  2. Presas e predadores são extintos depois de uma oscilação entre o tamanho de ambas as populações.

#definir parâmetros do modelo

a <- 0.05  #taxa de nascimento das presas
b <- 0.02  # taxa de predação
c <- 0.008  #taxa de nascimento dos predadores 
d <- 0.5  #taxa de morte dos predadores 

#definir condições iniciais
x <- 20  #pop inicial de presas
y <- 18  #pop inicial de predadore

t_end <- 50  

set.seed(124)
sim_data <- gil_lkv(a, b, c, d, x, y, t_end)

head(sim_data)
       Time Prey Predator
1 0.0000000   20       18
2 0.1199048   20       19
3 0.1213509   19       19
4 0.1986413   19       20
5 0.2076241   19       19
6 0.2327343   19       18
#criar gráfico da evolução das populações ao longo do tempo
ggplot(sim_data, aes(x = Time)) +
  geom_line(aes(y = Prey, color = "Presas")) +
  geom_line(aes(y = Predator, color = "Predadores")) +
  labs(title = "Simulação de Lotka-Volterra com o Algoritmo de Gillespie",
       x = "Tempo", y = "População") +
  scale_color_manual(name = "Espécies", values = c("Presas" = "blue", "Predadores" = "red")) +
  theme_minimal()

  1. Com a redução na taxa de nascimento das presas, a extinção dos predadores ocorre mais rapidamente.