Nratones <- 100
rratones <- 0.5
Kratones <- 1200
alpha <- 1/65
Ncoyotes <- 100
rcoyotes <- 0.14
Kcoyotes <- 15
beta <- 6
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.5.2
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(deSolve)
Lotka_volterra <- function(tiempo, estado, parametros){
  with(as.list(c(estado, parametros)),{
    # rate of change
    dx1 <- r1*x1*((k1 -x1 - alpha * x2)/k1) 
    dx2 <- r2*x2*((k2- x2 - beta * x1)/k2) 
    # return the rate of change
    list(c(dx1, dx2))
  }) #  xx = y[1];
  #  yy = y[2];
}
parametros <- c(r1 = 0.5, # Tasa de crecimiento SP1
                r2 = 0.14, # Tasa de Crecimiento SP2
                alpha = 0.1 , # Efecto de la SP2 sobre la SP1
                beta = -2, # Efecto de la SP1 sobre la SP2
                k1 = 120, # Capacidad de carga SP1
                k2 = 7) # Capacidad de carga SP2
estado <- c(x1 = 150, x2 = 10) # Poblaciones de SP1 y SP2
tiempo <- seq(0, 60, by = 1)
resultado <- ode(func = Lotka_volterra, y = estado, times = tiempo, parms = parametros)
resultadod <- as.data.frame(resultado)
ggplot(resultadod, aes(x=tiempo)) +
  geom_line(aes(y=x1), color="pink") +
  geom_line(aes(y=x2), color="purple") +
  labs(x= "Tiempo",y="Población")+
  geom_text(data = resultadod, aes(label = "x1", x = tail(tiempo,1), y = tail(x1,1)), color = "pink", hjust = -0.1, vjust = 0)+
  geom_text(data = resultadod, aes(label = "x2", x = tail(tiempo,1), y = tail(x2,1)), color = "purple", hjust = -0.1, vjust = 0)+
  theme_bw()
## Warning in geom_text(data = resultadod, aes(label = "x1", x = tail(tiempo, : All aesthetics have length 1, but the data has 61 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.
## Warning in geom_text(data = resultadod, aes(label = "x2", x = tail(tiempo, : All aesthetics have length 1, but the data has 61 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.

# Definir el modelo Lotka-Volterra para la dinámica de la plaga y las plantas
modelo_lotka_volterra <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    
    # Definir las ecuaciones diferenciales
    d_insectos <-  1-(tasa_crecimiento_insectos * insectos - tasa_interaccion * insectos * plantas - tasa_mortalidad_insectos * insectos)
    
    d_plantas <- (tasa_reproduccion_plantas * plantas - tasa_interaccion * insectos * plantas - tasa_mortalidad_plantas*plantas)
    
    # Devolver las derivadas
    return(list(c(d_insectos, d_plantas)))
  })
}
# Definir condiciones iniciales y parámetros del modelo
condiciones_iniciales <- c(insectos = 100, plantas = 1500)  # Población inicial de insectos y plantas
parametros <- c(tasa_crecimiento_insectos = 0.8,          # Tasa de crecimiento de los insectos
                tasa_mortalidad_insectos = 0.01,          # Tasa de mortalidad de los insectos
                tasa_reproduccion_plantas = 0.2,         # Tasa de reproducción de las plantas
                tasa_mortalidad_plantas = 0.3,           # Tasa de mortalidad de las plantas
                tasa_interaccion = 0.01)                 # Tasa de interacción entre insectos y plantas

# Definir el tiempo de simulación
tiempo <- seq(0, 100, by = 1)  # Simularemos 100 unidades de tiempo

# Resolver el modelo
resultado <- ode(y = condiciones_iniciales, times = tiempo, func = modelo_lotka_volterra, parms = parametros, maxsteps = 10000)
resultadod <- as.data.frame(resultado)
ggplot(resultadod, aes(x=tiempo)) +
  geom_line(aes(y=insectos), color="pink") +
  geom_line(aes(y=plantas), color="purple") +
  labs(x= "Tiempo",y="Población")+
  geom_text(data = resultadod, aes(label = "I", x = tail(tiempo,1), y = tail(insectos,1)), color = "pink", hjust = -0.1, vjust = 0)+
  geom_text(data = resultadod, aes(label = "P", x = tail(tiempo,1), y = tail(plantas,1)), color = "purple", hjust = -0.1, vjust = 0)+
  theme_bw()
## Warning in geom_text(data = resultadod, aes(label = "I", x = tail(tiempo, : All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.
## Warning in geom_text(data = resultadod, aes(label = "P", x = tail(tiempo, : All aesthetics have length 1, but the data has 101 rows.
## ℹ Please consider using `annotate()` or provide this layer with data containing
##   a single row.