Tarea 1

1. Modelo de Lotka Volterra

Considere el modelo \[ \begin{align} \frac{d x_1(t)}{dt}= &r_1x_1\left[1-\frac{x_1(t)+\alpha_{12}x_2(t)}{K_1}\right]\\ \frac{d x_2(t)}{dt}= &r_2x_2\left[1-\frac{x_2(t)+\alpha_{21}x_1(t)}{K_2}\right] \end{align} \] donde \(x_1\) y \(x_2\) representan dos poblaciones que comparten un ecosistema. Las variables \(x_1\) y \(x_2\) pueden representar los tamaños de las dos poblaciones en términos absolutos.

El modelo representado por las ecuaciones anteriores se denomina modelo de Volterra-Lotka y representa muchas relaciones posibles, según sean los signos de \(\alpha_{12}\) y \(\alpha_{21}\), como se muestra a continuación:

\(\alpha_{12}\)	\(\alpha_{21}\)	Tipo de relación
\(-\)	\(-\)	Mutualismo
\(-\)	\(0\)	Comensal
\(0\)	\(-\)	Comensal
\(+\)	\(-\)	Parasitismo
\(-\)	\(+\)	Parasitismo
\(+\)	\(+\)	Competencia

Cada uno debe seleccionar dos de estos casos y responder las preguntas a continuación para cada uno de estos casos.

1. Haga una interpretación biológica de cada uno de los parámetros del modelo.
2. Describa qué es lo que representan los diferentes términos de las ecuaciones diferenciales.
3. Considere un caso en el cual una de las poblaciones se extingue (por ejemplo \(x_2\)) y analice el comportamiento asintótico (en el límite cuando \(t\to \infty\)) de la otra. Cuál es el tamaño de la población en condiciones estacionarias \(\lim_{t\to\infty} x_1(t)=:\bar{x}_1\)? Analice cómo se afecta el valor de \(\bar{x}_1\) con cada uno de los parámetros. Para esto, haga simulaciones con valores iniciales \(x_1^0>\bar{x}_1\) y también \(x_1^0<\bar{x}_1\).
4. Ahora introduzca la segunda población. Observe que las soluciones se aproximan a valores constantes, llamados puntos de equilibrio (llámelos \(\bar{x}_1\) y \(\bar{x}_2\)). Encuentre los puntos de equilibrio para cuatro combinaciones de los parámetros \(\alpha_{12}\) y \(\alpha_{21}\). Observe el efecto de estos parámetros sobre los valores estacionarios \(\bar{x}_1\) y \(\bar{x}_2\).
5. A continuación, evalúe el efecto de los parámetros \(r_1\) y \(r_2\) sobre la respuesta del modelo.
6. Analice la estabilidad de los puntos de equilibrio; para ello suponga varias condiciones iniciales ligeramente diferentes del punto de equilibrio (por encima y por debajo de los valores de equilibrio) y observe el comportamiento de las soluciones. El equilibrio es estable si las soluciones ligeramente diferentes del equilibrio resultan en trayectorias que “regresan” a éste. En caso contrario es inestable.

Solución de los dos casos

Importación de librerías útiles. La primera, para resolver ecuaciones diferenciales, la segunda para graficar y ordenar. La librería tyidyverse importa otras librerías como ggplot2.

library(deSolve)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

Función que representa el modelo de Lotka-Volterra:

LotkaVolt = function (t, state, parametros){
  with(as.list(c(state, parametros)),{
    # rate of change
    dx1 <- r1*x1*(1-(x1 + a12 * x2)/k1) 
    dx2 <- r2*x2*(1-(x2 + a21 * x1)/k2) 
    # return the rate of change
    list(c(dx1, dx2))
  }) #  xx = y[1];
#  yy = y[2];
}

2. Describa qué es lo que representan los diferentes términos de las ecuaciones diferenciales.

Este punto aplica para ambos casos. Las derivadas representan la tasa instantánea de crecimiento y están coondicionadas a las tasas de crecimiento \(r\), el número de especies, \(x\), el efecto de la otra especie, \(\alpha\) y la capacidad de carga, \(K\). Este último valor es restado a 1 para simular una extinsión si se alzanzare la capacidad de carga.

Cuando tanto \(\alpha_{12}\) como \(\alpha_{21}\) son negativos: Mutualismo

Parámetros para este caso:

parametros <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = -0.9, # effect of species 2 on the growth of species 1
                a21 = -0.5, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

state <- c(x1 = 20, x2 = 20) # number of individuals from 1st and 2nd species
tiempos <- seq(0, 50, by = 0.1) # sequence for the passage of time
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)

Para resolver las ecuaciones diferenciales ordinales usamos deSolve::ode, como parámetros están y, estado inicail de las variables, times, secuecnias para la cual se desea el resultado, el primer valor de dicha secuencia es el tiempo inicial; los otros se sobreentienden. La solución de la ecuación se guarda en el objeto out a manera de matrix con valores para cada variable en cada tiempo de la secuencia.

matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

### To plot with ggplot2
#outd <- as.data.frame(out)
#ggplot(outd, aes(x=time)) +
#  geom_line(aes(y=x1), color="green") +
#  geom_line(aes(y=x2), color="red") +
#  theme_bw()
#tail(outd)

1. Haga una interpretación biológica de cada uno de los parámetros del modelo.

Dado que queremos evaluar el caso de mutualismo fijamos todos los parámetros excepto \(\alpha_{12}\) y \(\alpha_{21}\). En ese caso ambos \(\alpha\) son negativos y el efecto de una especie sobre la otra resulta beneficioso para la otra; pero como \(\alpha_{21}\) es manor que \(\alpha_{12}\), la especie dos se beneficia menos que la especie uno y su crecimiento no est tan alto, el beneficio no es simétrico. En general ambaas especies tienen problaciones altas, > 400, cuando se estabilizan (comprando por ejemplo con el caso del comensalismo [ver más abajo], en que la especie más abundante llega a un máximo de ~300).

3. Considere un caso en el cual una de las poblaciones se extingue (por ejemplo \(x_2\)) y analice el comportamiento asintótico (en el límite cuando \(t\to \infty\)) de la otra. Cuál es el tamaño de la población en condiciones estacionarias \(\lim_{t\to\infty} x_1(t)=:\bar{x}_1\)?

Para eso vamos a crear un conjunto de parámetros que reduzca \(x_2\) a la extinsión.

parametros2 <- c(r1 = 0.9, # 1st species growth rate
                r2 = -0.1, # 2nd species growth rate
                a12 = -0.9, # effect of species 2 on the growth of species 1
                a21 = -0.5, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

state2 <- c(x1 = 20, x2 = 200) # number of individuals from 1st and 2nd species
tiempos2 <- seq(0, 80, by = 0.1) # sequence for the passage of time
out2 <- deSolve::ode(func = LotkaVolt, y = state2, times = tiempos2, parms = parametros2)

matplot(out2[ , 1], out2[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

Al ajustar a menos 0 \(r_2\) se logra una caída en la población de la especie 2, con esto simulamos una extinsión. Con esto podemos ver que en ausencia virtual de \(x_2\), \(x_1\) alcanza un equilibrio, \(\bar{x_1} = 150\) cuando \(\lim_{t\to\infty} x_1(t)\). Esto indica que la supervivencia de \(x_1\) no depende en su totalidad de \(x_2\).

3.2 Analice cómo se afecta el valor de \(\bar{x}_1\) con cada uno de los parámetros. Para esto, haga simulaciones con valores iniciales \(x_1^0>\bar{x}_1\) y también \(x_1^0<\bar{x}_1\).

## Effect on bar x with different x0 
for(i in 1:4){
  x1 <- round(runif(1, 1, 900), digits=0)
  x2 = 100 # let's not vary x2
  param <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = -0.9, # effect of species 2 on the growth of species 1
                a21 = -0.5, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

  start <- c(x1 = x1, x2 = x2) # number of individuals from 1st and 2nd species
  time <- seq(0, 40, by = 0.1) # sequence for the passage of time
  result <- deSolve::ode(func = LotkaVolt, y = start, times = time, parms = param)
  
  matplot(result[ , 1], result[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
          main = paste("Poblaciones; x1 inicial =", x1, "y x2 inicial =", x2), lwd = 2)
  legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)
}

Vemos que a gran mayoría de las veces los equilibrios se mantienen en 500

4. Ahora introduzca la segunda población. Observe que las soluciones se aproximan a valores constantes, llamados puntos de equilibrio (llámelos \(\bar{x}_1\) y \(\bar{x}_2\)). Encuentre los puntos de equilibrio para cuatro combinaciones de los parámetros \(\alpha_{12}\) y \(\alpha_{21}\). Observe el efecto de estos parámetros sobre los valores estacionarios \(\bar{x}_1\) y \(\bar{x}_2\).

Para simular la introducción de la segunda población vamos a comenzar con valores de poblacion de 150 y 1 para \(x_1\) y \(x_2\) respectivamente.

for(i in 1:4){
  a12 <- round(runif(1, 0.1, 0.9), digits = 2) * -1 # random number from 0.1 to 0.9; -1 to keep negative values
  a21 <- round(runif(1, 0.1, 0.9), digits = 2) * -1
  param <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = a12, # effect of species 2 on the growth of species 1
                a21 = a21, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

  start <- c(x1 = 20, x2 = 20) # number of individuals from 1st and 2nd species
  time <- seq(0, 40, by = 0.1) # sequence for the passage of time
  result <- deSolve::ode(func = LotkaVolt, y = start, times = time, parms = param)
  
  matplot(result[ , 1], result[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
          main = paste("Poblaciones para dos especies relacionadas", "a12 =", a12, ", a21 =", a21), lwd = 2)
  legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)
}

Podemos ver que según sea \(\alpha\) más grande o más pequeño el beneficio disminuye o aumenta para la especie en cuestión. Esto se ve reflejado en los valores de \(\bar{x_1}\) y \(\bar{x_2}\)

5. A continuación, evalúe el efecto de los parámetros \(r_1\) y \(r_2\) sobre la respuesta del modelo.

Para ello vamos a cambiar la tasa de crecimiento de \(x_1\), pero dejando \(r_2\) como antes, y volvemos a graficar. Cambiémosla a una más pequeña:

parametros["r1"] = 0.2
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

Lo que vemos es que el punto de equilibrio se alzanza con valores similares antes de cambiar \(r_1\). Si bien, al principio, la población de la especie 1, que es a la que le hemos reducido la tasa de crecimiento, crece menos que la población de la especie 2 (tiempo 0 a 10) a la larga la sobrepasa. En este caso el hecho de que \(\alpha_{12}\) sea mayor que \(\alpha_{21}\) parece ser más preponderante, aunque si \(r_1\) es muy pequeña se espera que la población no logre aumentar aún con ayuda de la especie 2.

6. Analice la estabilidad de los puntos de equilibrio; para ello suponga varias condiciones iniciales ligeramente diferentes del punto de equilibrio (por encima y por debajo de los valores de equilibrio) y observe el comportamiento de las soluciones. El equilibrio es estable si las soluciones ligeramente diferentes del equilibrio resultan en trayectorias que «regresan» a éste. En caso contrario es inestable.

Hemos visto que al cambiar \(r_1\) el equilibrio es el mismo, cambiemos también \(r_2\) y el número de individuos, y grafiquemos:

parametros["r2"] = 0.1
state["x1"] = 300 # Twice the population when x2 is virtually non-existent
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

De nuevo, con condiciones diferente, al principio es las poblaciones tienen valores diferentes pero a la larga se llega a un mismo nivel. Así, se puede concluir que el equilibrio es estable, \(\bar{x_1} \approx 500\) y \(\bar{x_2} \approx 400\). Sin embargo, los parámetros que pueden cambiar dicho equilibrio son la capacidad, \(K\), de carga y el efecto de una especie sobre otra, \(\alpha\).

---

Cuando \(\alpha_{12} = 0\) y \(\alpha_{21}\) es negativo: Comensalismo

Parámetros para este caso:

parametros <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = 0, # effect of species 2 on the growth of species 1
                a21 = -0.9, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

state <- c(x1 = 20, x2 = 20) # number of individuals from 1st and 2nd species
tiempos <- seq(0, 50, by = 0.1) # sequence for the passage of time
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)

Gráfica:

matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

1. Haga una interpretación biológica de cada uno de los parámetros del modelo.

Dado que queremos evaluar el caso de mutualismo fijamos todos los parámetros excepto \(\alpha_{12}\) y \(\alpha_{21}\). En este caso vemos cómo la especie dos se beneficia de la especie 1 pero no así en el sentido inverso. Esto es representado con un \(\alpha_{21} > 0\) y un \(\alpha_{12} = 0\). Hay que notar que auque la especie dos se beneficie de la otra, su crecimeinto no es tan alto como el de la especie más abundante en el caso del mutualismo.

3. Considere un caso en el cual una de las poblaciones se extingue (por ejemplo \(x_2\)) y analice el comportamiento asintótico (en el límite cuando \(t\to \infty\)) de la otra. Cuál es el tamaño de la población en condiciones estacionarias \(\lim_{t\to\infty} x_1(t)=:\bar{x}_1\)?

Para eso vamos a crear un conjunto de parámetros que reduzca \(x_2\) a la extinsión.

parametros2 <- c(r1 = 0.9, # 1st species growth rate
                r2 = -0.1, # 2nd species growth rate
                a12 = 0, # effect of species 2 on the growth of species 1
                a21 = -0.9, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

state2 <- c(x1 = 20, x2 = 200) # number of individuals from 1st and 2nd species
tiempos2 <- seq(0, 80, by = 0.1) # sequence for the passage of time
out2 <- deSolve::ode(func = LotkaVolt, y = state2, times = tiempos2, parms = parametros2)

matplot(out2[ , 1], out2[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

Al ajustar a menos 0 \(r_2\) se logra una caída en la población de la especie 2, con esto simulamos una extinsión. Podemos ver que en ausencia virtual de \(x_2\), \(x_1\) alcanza un equilibrio, \(\bar{x_1} = 150\) cuando \(\lim_{t\to\infty} x_1(t)\). Esto indica que la supervivencia de \(x1\) no depende en su totalidad de \(x2\).

3.2 Analice cómo se afecta el valor de \(\bar{x}_1\) con cada uno de los parámetros. Para esto, haga simulaciones con valores iniciales \(x_1^0>\bar{x}_1\) y también \(x_1^0<\bar{x}_1\).

## Effect on bar x with different x0 
for(i in 1:4){
  x1 <- round(runif(1, 1, 900), digits=0)
  x2 = 100 # let's not vary x2
  param <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = 0, # effect of species 2 on the growth of species 1
                a21 = -0.9, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

  start <- c(x1 = x1, x2 = x2) # number of individuals from 1st and 2nd species
  time <- seq(0, 40, by = 0.1) # sequence for the passage of time
  result <- deSolve::ode(func = LotkaVolt, y = start, times = time, parms = param)
  
  matplot(result[ , 1], result[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
          main = paste("Poblaciones; x1 inicial =", x1, "y x2 inicial =", x2), lwd = 2)
  legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)
}

Vemos que la gran mayoría de las veces los equilibrios se mantienen en ~180

4. Ahora introduzca la segunda población. Observe que las soluciones se aproximan a valores constantes, llamados puntos de equilibrio (llámelos \(\bar{x}_1\) y \(\bar{x}_2\)). Encuentre los puntos de equilibrio para cuatro combinaciones de los parámetros \(\alpha_{12}\) y \(\alpha_{21}\). Observe el efecto de estos parámetros sobre los valores estacionarios \(\bar{x}_1\) y \(\bar{x}_2\).

Para simular la introducción de la segunda población vamos a comenzar con valores de poblacion de 180 y 1 para \(x_1\) y \(x_2\) respectivamente. Uno de los \(\alpha\) no lo vamos a variar porque estamos en el caso de comensalismo.

for(i in 1:4){
  a12 <- 0 # Because this is commercialism
  a21 <- round(runif(1, 0.1, 0.9), digits = 2) * -1
  param <- c(r1 = 0.9, # 1st species growth rate
                r2 = 0.9, # 2nd species growth rate
                a12 = a12, # effect of species 2 on the growth of species 1
                a21 = a21, # effect of species 1 on the growth of species 2
                k1 = 150, # 1st species carrying capacity
                k2 = 150) # 2nd species carrying capacity

  start <- c(x1 = 20, x2 = 20) # number of individuals from 1st and 2nd species
  time <- seq(0, 40, by = 0.1) # sequence for the passage of time
  result <- deSolve::ode(func = LotkaVolt, y = start, times = time, parms = param)
  
  matplot(result[ , 1], result[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
          main = paste("Poblaciones para dos especies relacionadas", "a12 =", a12, ", a21 =", a21), lwd = 2)
  legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)
}

Vemos que el equilibrio de la población de la especie no comensal, \(x_1\), se mantiene invariable. Mientras que al varial el efecto de la otra especie sobre sí, o bien tiene más éxito o bien lo disminuye. En caso de que ambos \(\alpha\) sean de cero las dos especies serían independeintes

5. A continuación, evalúe el efecto de los parámetros \(r_1\) y \(r_2\) sobre la respuesta del modelo.

Para ello vamos a cambiar la tasa de crecimiento de \(x_1\), pero dejando \(r_2\) como antes, y volvemos a graficar. Cambiémosla a una más pequeña:

parametros["r1"] = 0.2
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

Cuando la presencia de la especie 1 es escasa, la población de la especie dos se demora más en crecer en compración con la gráfica de arriba. El efecto de la especie uno sobre la especie dos afecta el crecimiento de la especie dos al menos hasta el tiempo 30, en adelante las poblaciones se equilibran como en la gráfica de arriba.

6. Analice la estabilidad de los puntos de equilibrio; para ello suponga varias condiciones iniciales ligeramente diferentes del punto de equilibrio (por encima y por debajo de los valores de equilibrio) y observe el comportamiento de las soluciones. El equilibrio es estable si las soluciones ligeramente diferentes del equilibrio resultan en trayectorias que «regresan» a éste. En caso contrario es inestable.

Hemos visto que al cambiar \(r_1\) el equilibrio es el mismo, cambiemos también \(r_2\) y el número de individuos, y grafiquemos:

parametros["r2"] = 0.1
state["x2"] = 400
out <- deSolve::ode(func = LotkaVolt, y = state, times = tiempos, parms = parametros)
matplot(out[ , 1], out[ , 2:3], type = "l", xlab = "Tiempo", ylab = "Población",
        main = "Poblaciones para dos especies relacionadas", lwd = 2)
legend("bottomright", c("Especie 1", "Especie 2"), col = 1:2, lty = 1:2)

De nuevo, con condiciones diferente, al principio es las poblaciones tienen valores diferentes pero a la larga se llega a un mismo nivel. Así, se puede concluir que el equilibrio es estable, \(\bar{x_1} \approx 250\) y \(\bar{x_2} \approx 150\). Sin embargo, los parámetros que pueden cambiar dicho equilibrio son la capacidad, \(K\), de carga y el efecto de una especie sobre otra, \(\alpha\).