Intro

Anteriormente vimos cómo generar un modelo de supervivencia en función del tiempo y de la temperatura con la función Weibull. Ahora veremos cómo utilizar dicho modelo para obtener predicciones espaciales con base en las condiciones de temperatura ambiental. Los requerimientos son:

  1. Identificar qué datos de temperatura utilizar
  2. Desarrollar la función para simular la ecuación diferencial
  3. Diseñar escenarios para simular

Escenario simple

Para proyectar nuestro modelo Weibull en el espacio de manera sencilla, utilizaremos una capa raster de temperatura anual promedio. Como corte de tiempo, fijaremos 12 h = 720 min.

Importaremos primero la capa ya recortada de temperatura anual promedio para la región del mundo donde tiene sentido biológico la representación de los patrones de supervivencia:

library(terra)
## terra 1.7.29
bio1 <- rast("Nicho/Bio1.tif")
plot(bio1)

Ahora importamos los coeficientes del modelo y la función que escribimos anteriormente para hacer los cálculos.

coef.pk <- read.csv("Nicho/Coeficientes-pk.csv") #Importando coeficientes
source("Nicho/weibull.R") #Importando la función

Para hacer los cálculos, tenemos que transformar la capa bio1 a una tabla:

bio1.df <- as.data.frame(bio1, xy = T)

Y continuamos tal y como lo hicimos anteriormente, usando los valores de bio1 como temperatura:

pars <- as.list(coef.pk$valor)
names(pars) <- coef.pk$par

Sup.bio1 <- weib(Tiempo = 12,
              Temp = bio1.df$Bio1,
              pars = pars)

Ahora tenemos que añadir las predicciones de supervivencia a bio1.df y transformar a raster:

bio1.df$Sup <- Sup.bio1
Sup.r <- rast(bio1.df[, -3])
plot(Sup.r)

Escenario realista

Comenzaremos creando la función de integración de la ecuación diferencial porque al sustituir la temperatura por el tiempo hemos eliminado la posibilidad de una forma cerrada como en la función Weibull del ejercicio anterior.

weib.ode <- function(t, y, params){
  with(params, {
      S <- y
      
      Temp <- Tmax - Tdif * cos(pi * t/24)^2
      
      rho <- exp(ap + Bp * Temp)
      
      kappa <- exp(ak + Bk * Temp)
      
      dS <- - rho * kappa * (-log(S))^(1 - 1/kappa) * S
      
      list(c(dS))
  })
}

El semestre pasado vimos cómo integrar sistemas de ecuaciones diferenciales ordinarias con el paquete deSolve, lo cual volveremos a hacer en esta ocasión, pero veremos cómo utilizar los datos de temperatura de una capa ráster. Para ello importaremos las capas de temperatura min y max en dicho formato, y las transformaremos a data.frame:

Tmax <- rast("Nicho/Tmax-01.tif")
Tmin <- rast("Nicho/Tmin-01.tif")

Tdif <- Tmax - Tmin

Tmax.df <- as.data.frame(Tmax, xy = T)
Tdif.df <- as.data.frame(Tdif, xy = F)

Temps.df <- data.frame(Tmax.df, Tdif = Tdif.df$`Tmax-01`)
names(Temps.df) <- c("x", "y", "Tmax", "Tdif")

Ahora, resolveremos el modelo para cada par de datos de temperatura máxima y diferencia dando como resultado una trayectoria de supervivencia a lo largo de 12h (720 min) para cada píxel. Para ahorrarnos datos, sólo utilizaremos el valor final de supervivencia. Primero preparamos el valor inicial y el vector que contiene los puntos en el tiempo para evaluar la supervivencia

library(deSolve)

y <- 0.9999
t <- seq(0, 12, length(100))

Estos los reciclaremos en todas las simulaciones por píxel por medio de una iteración con for. Yo utilizo doParallel para aumentar la velocidad:

library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:terra':
## 
##     extract
library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
registerDoParallel(cores = 8)

sims <- foreach(i = 1:nrow(Temps.df), .combine = c) %dopar% {
  params <- pars
  params$Tmax <- Temps.df$Tmax[i]
  params$Tdif <- Temps.df$Tdif[i]
  out <- lsoda(y = y, times = t,
               func = weib.ode,
               parms = params)
  return(out[nrow(out), 2])
}

Ahora transformaremos las predicciones de supervivencia con temperatura fluctuante en una capa raster:

Sup.fluc <- data.frame(Temps.df[, c("x", "y")], Sup = sims)
Sup.fluc.r <- rast(Sup.fluc)
plot(Sup.fluc.r)

Hemos aprendido a implementar una serie de análisis para utilizar conocimiento adquirido en laboratorio para generar predicciones sobre la posible susceptibilidad de Hendra virus a las características ambientales. Para continuar trataremos de responder las siguientes preguntas con base en la biología del organismo de que se trata:

  1. ¿Qué posible utilidad tendría este modelo tan simple?
  2. ¿Bajo qué supuestos se le daría ese uso?
  3. ¿Con qué propósito utilizarían la información ggenerada?
  4. ¿Cuáles son las limitaciones?
  5. ¿Cómo incorporarían nuevos factores en este modelo?
  6. ¿Es factible generar un modelo más complejo que este?