Guía metodológica para el uso del paquete R devRate

Recordatorios sobre biología de insectos

Antes de comenzar, es necesario hacer algunos recordatorios sobre la biología de los insectos. Los insectos son organismos ectotérmos, lo que significa que la temperatura de su cuerpo es poco o no está controlada por un proceso biológico activo: en otras palabras, producen poco o ningún calor y son totalmente dependientes de la temperatura de su entorno. Se oponen a los organismos endotérmos, como los mamíferos, que producen en general su propia fuente de calor para regular su temperatura.

En la mayoría de los ambientes, la temperatura varía en el tiempo -por ejemplo, es más fría por la noche que durante el día- y en el espacio, por ejemplo, es más fría a la sombra que al sol. Los organismos ectotermos evolucionan en ambientes donde los cambios de temperatura varían, a veces considerablemente. Hablamos entonces de organismos poikilotérmos, entre los que encontramos la mayoría de los insectos, que contrastan con los organismos homotérmos que mantienen una temperatura relativamente constante independientemente de las condiciones externas.

La importancia de la temperatura ambiente en los procesos biológicos de los insectos es clave, aunque no es el único factor involucrado. Por ejemplo, otros factores pueden ser la humedad relativa, o el fotoperiodo (la duración del día y la noche). Aqui nos vamos a enfocar unicamente en la temperatura.

Los insectos tienen una preferencia térmica, es decir un valor o conjunto de valores de temperatura, dentro de los cuales pueden alcanzar un desarrollo óptimo. Esta preferencia es variable dependiendo de la especie, con especies que se desarrollarán en un amplio rango de temperaturas. A estas se denomina especies generalistas, a diferencia de las especies que pueden desarrollarse solo en un rango muy estrecho de temperaturas, que son las llamadas especialistas (los insectos se encuentran a lo largo de este gradiente entre especialistas y generalistas). Las variaciones de preferencia también se pueden observar dentro de una especie, especialmente entre las diferentes etapas de desarrollo. Para entender estas diferentes estrategias y caracterizarlas, es necesario estudiar la relación entre la temperatura y el desarrollo de los insectos.

Nos interesamos en esta occasión en los insectos poikilotérmos, que, recordamos, corresponden a los insectos cuya temperatura cambiará como consecuencia de los cambios de la temperatura exterior. Estudiaremos particularmente la tasa de desarrollo de insectos, generalmente expresada como la inversa de la cantidad de días que pasan de una etapa de desarrollo a otra. Para este proposito, usaremos datos disponibles en la literatura.

En la mayoría de los ecosistemas, las condiciones ambientales son variables. Este es el caso de la temperatura. Para los insectos poikilotérmos, esto dará como resultado un cambio en la temperatura de su cuerpo. Mas alla del rango en el cual pueden desarrollarse, sufrirán un estrés que puede llevarles a la muerte. Por lo tanto, existen umbrales críticos, un umbral mínimo y un umbral máximo, más allá de los cuales los insectos no pueden continuar su desarrollo. Esto límites se encuentran en la literatura bajo los acrónimos CTmin y CTmax para Critical Thermal minima Critical Thermal maxima, respectivamente. Los mismos varían ampliamente de una especie a otra, con algunas especies como el escarabajo ártico capaz de soportar hasta -70 grados, mientras que otras admiten hasta mas de 45 grados. Como se mencionó anteriormente, la sensibilidad a la temperatura, que clasifica los insectos de acuerdo con un gradiente de generalistas a especialistas varía mucho de una especie a otra e incluso dentro de una especie entre sus etapas fenológicas. Más allá de estos dos valores umbral, la respuesta de los insectos variará según el tiempo frecuencia y tasa de exposición, la velocidad a la que la temperatura se lleva más allá de los valores umbral, la historia de vida de los insectos y, por supuesto, las especies consideradas. En esta guía nos enfocaremos en lo que sucede entre estos umbrales.

Para estudiar qué ocurre entre estos dos valores umbral, mediremos el tiempo de desarrollo del insecto, expresado como la tasa de desarrollo. La tasa de desarrollo de los insectos comienza en el primer valor umbral, luego aumenta progresivamente hasta alcanzar un valor óptimo, luego desciende rápidamente al segundo valor umbral. La respuesta de los insectos a la temperatura es, por lo tanto, no lineal.

La temperatura a la cual el desarrollo es máximo se llama temperatura óptima de desarrollo u óptimo térmico (Topt). Es importante enfatizar que esta es la temperatura a la cual el desarrollo es óptimo, no la temperatura óptima para el insecto. De hecho, la temperatura óptima depende de otros criterios distintos del desarrollo. Un buen ejemplo es la supervivencia: a la temperatura óptima para el desarrollo, la supervivencia observada es diferente de la temperatura óptima para la supervivencia.

Los estudios sobre la respuesta de los insectos a la temperatura no son nuevos: el primer científico que estudio y publico sobre el efecto de la temperatura en los insectos es René-Antoine Ferchault de Réaumur¹: sus observaciones sobre las mariposas estan registradas en sus “Mémoires pour servir à l’histoire des insectes” entre 1734 y 1742. Despues de Réaumur, muchos científicos buscaron y siguen buscando cuantificar la relación entre la temperatura y la tasa de desarrollo.

Para hacer esto, el método más común consiste en un primer paso en colocar insectos en jaulas a diferentes temperaturas constantes, para medir el tiempo de desarrollo de cada etapa fenológica. Mientras haciendo este experimento, podemos observar algunas características de la respuesta de los insectos a la temperatura:

La primera que resulta muy difícil, si no imposible, de obtener datos de desarrollo cercanos a los mínimos y máximos, ya que a estas temperaturas la gran mayoría de los insectos de la misma especie y población no sobreviven
La segunda característica es que se observa variabilidad dentro de la misma población de insectos en cuanto a su respuesta a la temperatura

La tercera característica es que para el rango de temperaturas que se encuentra en el hábitat del insecto, la relación entre la temperatura y la tasa de desarrollo es casi lineal, pero fuera de este rango la relación es no lineal

El segundo paso es el ajuste de un modelo matemático a las mediciones realizadas en el laboratorio. El objetivo del ajuste es encontrar estimadores de parámetros de manera que el modelo esté lo más cerca posible de los puntos experimentales y al mismo tiempo retenga un significado biológico.

El objetivo de este documento es de proveer una guía para facilitar y analizar el ajuste de un modelo matemático a las mediciones realizadas en el laboratorio.

Datos del laboratorio

Los datos de laboratorio que vamos a usar son los de Bactrocera dorsalis ², un insecto del Orden de los Dipteros, que hemos obtenidos a partir del articulo de Messenger and Flitters (1958)³. Los datos coresponden al tiempo de desarrollo de los huevos expresado en horas a distintas temperaturas expresadas en Fahrenheit.

Para seguir la guía, es nesesario contar con el programa R instalado en su computadora. Para usar el paquete devRate, hay que instalarlo mediante:

install.packages("devRate")

Una vez el paquete instalado, no es necesario repetir a futuro la linea de instalacion. Sin embargo instalar el paquete no es suficiente para poder usarlo dentro de R. Hay que especificar que vamos a usarlo cada vez que usamos R. Para esto se puede usar:

library("devRate")

A partir de ahora vamos a trabajar desde un archivo que podriamos llamar B_dorsalis_devRate.R. Este archivo empeza con un comentario indicando de que se trata para acordarse de su contenido. Todo lo que viene despues del simbolo # son comentarios.

### script para analizar los datos de tasa de desarrollo de B. dorsalis en 
### función de la temperatura.
require("devRate") # para cargar el paquete devRate

Para que los datos sean reconocidos en R, tenemos que usar un formato especial. Los datos de temperatura son almacenados en un formato vector que se escribe con c() en R y que corresponde a una colección ordenada de elementos del mismo tipo:

c(55.0, 56.0, 57.0, 58.0, 60.0, 62.5, 65.0, 67.5, 70.0, 75.0, 80.0, 85.0, 87.5, 
  90.0, 92.5, 95.0, 96.0, 97.0, 97.5)

##  [1] 55.0 56.0 57.0 58.0 60.0 62.5 65.0 67.5 70.0 75.0 80.0 85.0 87.5 90.0 92.5
## [16] 95.0 96.0 97.0 97.5

Los numeros usan el punto como separator de decimales y la coma como separator de valores. Aqui las temperaturas estan en Fahrenheit, tenemos que convertirlas en Celsius con la formula T(°C) = (T(°F) - 32) / 1.8.

(c(55.0, 56.0, 57.0, 58.0, 60.0, 62.5, 65.0, 67.5, 70.0, 75.0, 80.0, 85.0, 87.5, 
  90.0, 92.5, 95.0, 96.0, 97.0, 97.5) - 32) / 1.8

##  [1] 12.77778 13.33333 13.88889 14.44444 15.55556 16.94444 18.33333 19.72222
##  [9] 21.11111 23.88889 26.66667 29.44444 30.83333 32.22222 33.61111 35.00000
## [17] 35.55556 36.11111 36.38889

Para guardar las temperaturas, podemos almacenarlas en un objeto. Para almacenar los valores de temperaturas en un objeto, se usa los simbolos <-. Si queremos guardar los valores de temperatura en un objeto llamado temp:

temp <- (c(55.0, 56.0, 57.0, 58.0, 60.0, 62.5, 65.0, 67.5, 70.0, 75.0, 80.0, 
           85.0, 87.5, 90.0, 92.5, 95.0, 96.0, 97.0, 97.5) - 32) / 1.8

De la misma manera, vamos a guardar los valores de tasa de desarrollo en un objeto llamado devRate. La tasa de desarrollo corresponde al opuesto del numero de días. Aquí los valores estan expresados en horas, entonces dividimos el numero de horas por 24 y despues calculamos el opuesto.

devRate <- 1/(c(263.0, 232.0, 170.5, 148.0, 121.3, 95.5, 74.0, 62.5, 51.5, 
                38.0, 30.5, 27.0, 25.0, 24.0, 23.5, 25.0, 26.5, 29.3, 34.3)/24)

Ahora podemos almacenar los datos de temperatura y de tasa de desarrollo en una tabla datosLab donde temperatura y tasa de desarrollo son las columnas. La función data.frame permite hacer la tabla.

datosLab <- data.frame(temp, devRate)

Si queremos hacer la tabla y llenarla, hay que especificar el nombre que queremos dar a las columnas con el simbolo = :

datosLab <- data.frame(
  temp = (c(55.0, 56.0, 57.0, 58.0, 60.0, 62.5, 65.0, 67.5, 70.0, 75.0, 80.0, 
            85.0, 87.5, 90.0, 92.5, 95.0, 96.0, 97.0, 97.5) - 32)/1.8, 
  devRate = 1/(c(263.0, 232.0, 170.5, 148.0, 121.3, 95.5, 74.0, 62.5, 51.5, 
            38.0, 30.5, 27.0, 25.0, 24.0, 23.5, 25.0, 26.5, 29.3, 34.3)/24))

Podemos verificar el contenido del objeto datosLab con la función print :

print(datosLab)

##        temp    devRate
## 1  12.77778 0.09125475
## 2  13.33333 0.10344828
## 3  13.88889 0.14076246
## 4  14.44444 0.16216216
## 5  15.55556 0.19785655
## 6  16.94444 0.25130890
## 7  18.33333 0.32432432
## 8  19.72222 0.38400000
## 9  21.11111 0.46601942
## 10 23.88889 0.63157895
## 11 26.66667 0.78688525
## 12 29.44444 0.88888889
## 13 30.83333 0.96000000
## 14 32.22222 1.00000000
## 15 33.61111 1.02127660
## 16 35.00000 0.96000000
## 17 35.55556 0.90566038
## 18 36.11111 0.81911263
## 19 36.38889 0.69970845

Nuestros datos estan listos para empezar el análisis con el paquete devRate.

Seleccionar un modelo y ajustarlo

El modelo lineal

Aunque tenga limitaciones, el método más común y simple es usar los resultados obtenidos en la parte casi lineal de la relación entre la temperatura y la tasa de desarrollo. El modelo matemático utilizado es un modelo lineal, la línea de ecuación y = bb * x + aa que todos conocen de sus clases de estadistica. El uso del modelo lineal fue popularizado por el articulo cientifico de Campbell et al. (1974)⁴. Este modelo se basa en la suposición de que la relación entre temperatura y desarrollo en insectos es lineal en un cierto rango de temperaturas y sirve de base para el concepto de “grado-dias”. La extensión de la línea de regresión lineal fuera de la zona de linealidad intersecta el eje x. El valor de la temperatura en esta intersección se denomina temperatura base, que suele ser diferente del CTmin. La temperatura base es también conocido como el “umbral de desarrollo” y puede calcularse simplemente dividiendo el opuesto de la intersección en y por el coeficiente “a”. Pero primero hay que ajustar el modelo a los datos de laboratorio.

Para seleccionar los datos que vamos a usar, puede ser util visualizar los datos de laboratorio en un grafico. Para esto vamos a usar la función plot.

plot(x = datosLab$temp, y = datosLab$devRate, 
     xlab = "Temperatura", ylab = "Tasa de desarrollo", 
     xlim = c(0, 40), ylim = c(0, 1.2))

Podemos observar que los ultimos cinco puntos salen de la zona de linearidad. Para no usar estos ultimos puntos, vamos a seleccinar una parte de los datos de laboratorio. Quitar los ultimos cinco elementos es igual a seleccionar unicamente las 14 primeras lineas de nuestra tabla. Se usan corchetes para seleccionar una parte de los datos, y una coma para especificar que se trata de lineas.

datosLab14 <- datosLab[1:14,]
plot(x = datosLab14$temp, y = datosLab14$devRate, 
     xlab = "Temperatura", ylab = "Tasa de desarrollo", 
     xlim = c(0, 40), ylim = c(0, 1.2))

Aun que se podria ajustar directamente el modelo lineal con la función lm, vamos a usar la función generica del paquete devRate.

modLin <- devRateModel(eq = campbell_74, dfData = datosLab14)
plot(x = datosLab$temp, y = datosLab$devRate, 
     xlab = "Temperatura", ylab = "Tasa de desarrollo", 
     xlim = c(0, 40), ylim = c(0, 1.2)) # datos del laboratorio
abline(modLin) # añadir modelo lineal en el grafico

print(modLin) # imprimir resultados del ajuste del modelo

## Nonlinear regression model
##   model: rT ~ aa + bb * T
##    data: data.frame(rT = devRate, T = temp)
##       aa       bb 
## -0.55176  0.04881 
##  residual sum-of-squares: 0.004761
## 
## Number of iterations to convergence: 1 
## Achieved convergence tolerance: 1.218e-06

El concepto de grados-días se basa en la hipótesis de que los insectos necesitan un mínimo de acumulación de temperatura para completar su desarrollo. El ciclo de vida de los insectos esta divido en fases. Cada fase del desarrrollo va a necesitar su mínimo de acumulación de temperatura para que el inscto pueda llegar a la fase siguiente y completar su desarrollo. La acumulación de temperatura no empeza desde cero grados, sino de una temperatura base que depende del insecto y de la fase de desarrollo.

Esta guía no discutirá los numerosos métodos disponibles para calcular los grados-días. Sin embargo, el método más simple es agregar el valor de temperatura mínima a la temperatura máxima de un día y luego dividir el resultado de la suma entre dos para luego restar a este resultado el valor de la temperatura base. Es decir, para un día determinado cuando la temperatura máxima es de 25 grados, y la temperatura mínima de 15 grados, y con una temperatura base de una especie de insecto dada de 3 grados, el cálculo de los grados los días acumulados en este día son 25 más 15 es igual a 40; dividido entre 2 es igual a 20; menos 3 es igual a 17 grados-días. El individuo completa su desarrollo cuando la suma de los grados-día alcanza el valor de su constante térmica, determinada a partir de la línea de regresión. La constante térmica denotada por “K” es igual a la inversa del coeficiente director de la línea de regresión, es decir, K es igual a 1 por bb con y = bbx + aa, donde y representa el desarrollo y x la temperatura.

Aqui bb = 0.04881 y aa = -0.55176, entonces:

la temperatura base es Tbase = -aa/bb = 11.30424
la constante térmica es K = 1/bb = 20.4876.

Podemos buscar en la base de datos del paquete devRate si hay otros modelos que fueron usados para la especie Bactrocera dorsalis.

devRateFind(species = "Bactrocera dorsalis")

##      equation occu paramNumb
## 1 campbell_74    2         2

campbell_74$startVal[campbell_74$startVal$genSp == "Bactrocera dorsalis",]

##     ordersp    familysp    genussp  species               genSp stage
## 679 Diptera Tephritidae Bactrocera dorsalis Bactrocera dorsalis   all
## 680 Diptera Tephritidae Bactrocera dorsalis Bactrocera dorsalis   all
##        param.aa    param.bb                 ref
## 679 -0.03256579 0.003289474 Jarosik et al. 2011
## 680 -0.01863014 0.002739726 Jarosik et al. 2011

Al momento de escribir esta guía, en la base de datos existen dos estudios que han trabajado con Bactrocera dorsalis usando el modelo lineal, pero ninguno ha trabajado unicamente con los huevos, así que no podemos comparar los resultados obtenidos.

Los modelos de grados-días tienen una aplicación directa para predecir, a veces con precisión, la aparición de un insecto basándose en las proyecciones de temperatura futuras, y así servir como un sistema experto para muchas aplicaciones agrícolas o epidemiológicas. El concepto de grados-días tiene la ventaja de ser simple y, por lo tanto, fácil de implementar a partir de un conjunto de datos restringidos, pero tiene sus límites. En primer lugar, tiende a subestimar el desarrollo de insectos a
temperaturas bajas, en el área entre la temperatura base y la temperatura umbral de desarrollo CTmin. Además, no permite determinar la temperatura umbral de desarrollo más baja (CTmin), ni superior (CTmax). Luego, sobreestima el desarrollo de los insectos tan pronto como uno se acerca a la temperatura óptima de desarrollo. Por lo tanto, solo puede usarse en un rango limitado de temperaturas y brinda poca información sobre la relación entre la temperatura y la tasa de desarrollo. La alternativa al concepto de grados-días es el uso de modelos no lineales.

Los modelos no lineales

Una de las principales motivaciones para el desarrollo de ecuaciones no lineales es el cálculo de CTmin, CTmax, Topt. Existen varios modelos no lineales. La lista de los que estan en el paquete devRate se puede encontrar en la documentacion mediante ?devRate o a traves del objeto devRateEqList.

names(devRateEqList)

##  [1] "janisch_32"         "davidson_44"        "campbell_74"       
##  [4] "stinner_74"         "logan6_76"          "logan10_76"        
##  [7] "sharpeDeMichele_77" "analytis_77"        "schoolfield_81"    
## [10] "schoolfieldHigh_81" "schoolfieldLow_81"  "taylor_81"         
## [13] "wang_82"            "poly2"              "harcourtYee_82"    
## [16] "poly4"              "ratkowsky_82"       "ratkowsky_83"      
## [19] "rootsq_82"          "bieri1_83"          "hilbertLogan_83"   
## [22] "wagner_88"          "lamb_92"            "lactin1_95"        
## [25] "lactin2_95"         "beta_95"            "beta_16"           
## [28] "wangengel_98"       "briere1_99"         "briere2_99"        
## [31] "bayoh_03"           "kontodimas_04"      "damos_08"          
## [34] "damos_11"           "shi_11"             "perf2_11"          
## [37] "regniere_12"

Dentro de todos estos modelos, la seleccion de uno en particular no tiene consenso en la comunidad cientifica. La mayoría de los modelos no lineales se basan en la descripción de la curva de respuesta a la temperatura de los insectos, con un enfoque matemático sin base biológica real, por lo que los diferentes modelos no pueden separarse en sus supuestos subyacentes, y los criterios estadísticos para decidir entre ellos no permiten que un modelo se afirme a sí mismo como mejor que los demás. Otros modelos no lineales, llamados modelos biofísicos, se basan en la activación enzimática, como el modelo Sharpe y DeMichèle de 1977, modificado en 1981 por Schoolfield et al. Lamentablemente, estos modelos son más complejos y requieren conjuntos de datos vastos para ser utilizados. Además, algunos estudios cuestionan su base teórica. Al final, se usan poco.

Para cuantificar la respuesta de los insectos a la temperatura, se trata de ajustar varios modelos al conjunto de datos experimentales, y seleccionar el que mejor responda a la problemática del estudio, y evaluarlo con los criterios de comparación estadísticos. Una posibilidad es de buscar en la base de datos del paquete cuales son los mas usados con la función devRateFind pero como hemos visto anteriormente no existen modelos no lineales que hayan sido usados al momento de escribir esta guía. Una alternativa es usar los modelos mas comunes como los de Briere, Logan, y Lactin. Para esto, hay que definir valores inciales para los parametros de los modelos para que el algoritmo que va a hacer el ajuste encuentre una solución. En caso de que no se conocen valores pertinentes para iniciar el algoritmo, se puede usar los valores promedios que se encuentran en la tabla startVal de cada modelo (por ejemplo briere1_99$startVal para el modelo de Briere-1).

devRateFind(species = "Bactrocera dorsalis")

##      equation occu paramNumb
## 1 campbell_74    2         2

modNoLin_01 <- devRateModel(
  eq = briere1_99, # nombre del modelo
  dfData = datosLab, # nombre de los datos de laboratorio
  startValues = list(aa = 0.01, Tmin = 10, Tmax = 40)) # valores iniciales
print(modNoLin_01)

## Nonlinear regression model
##   model: rT ~ aa * T * (T - Tmin) * (Tmax - T)^(1/2)
##    data: data.frame(rT = devRate, T = temp)
##        aa      Tmin      Tmax 
## 5.581e-04 1.103e+01 3.875e+01 
##  residual sum-of-squares: 0.02468
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 5.695e-06

modNoLin_02 <- devRateModel(
  eq = briere2_99, # nombre del modelo
  dfData = datosLab, # nombre de los datos de laboratorio
  startValues = list(
    aa = 0.01, Tmin = 10, Tmax = 40, bb = 2)) # valores iniciales
print(modNoLin_02)

## Nonlinear regression model
##   model: rT ~ aa * T * (T - Tmin) * (Tmax - T)^(1/bb)
##    data: data.frame(rT = devRate, T = temp)
##        aa      Tmin      Tmax        bb 
## 9.291e-04 9.101e+00 3.673e+01 4.049e+00 
##  residual sum-of-squares: 0.001665
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 3.416e-06

modNoLin_03 <- devRateModel(
  eq = lactin2_95, # nombre del modelo
  dfData = datosLab, # nombre de los datos de laboratorio
  startValues = list(
    aa = 0.03, Tmax = 30, deltaT = 5.0, bb = -1.5)) # valores iniciales
print(modNoLin_03)

## Nonlinear regression model
##   model: rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT) + bb
##    data: data.frame(rT = devRate, T = temp)
##       aa     Tmax   deltaT       bb 
##  0.02783 39.52991  2.12437 -1.33771 
##  residual sum-of-squares: 0.004803
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 3.668e-06

Los resultados se pueden visualizar graficamente con la función devRatePlot.

par(mfrow = c(1, 3)) # para hacer tres graficos en la misma pagina
devRatePlot(eq = briere1_99, 
  nlsDR = modNoLin_01, 
  xlim = c(10, 40), ylim = c(0, 1.2))
devRatePlot(eq = briere2_99, 
  nlsDR = modNoLin_02, 
  xlim = c(10, 40), ylim = c(0, 1.2))
devRatePlot(eq = lactin2_95, 
  nlsDR = modNoLin_03, 
  xlim = c(10, 40), ylim = c(0, 1.2))

Existen varios indices estadisticos para comparar los modelos, uno de esos es AIC. El modelo con menor AIC es el que se ajusta mejor a los datos. No significa que el modelo es valido, sino que tiene el mejor ajuste dentro de los modelos seleccionados.

c(AIC(modNoLin_01), AIC(modNoLin_02), AIC(modNoLin_03))

## [1]  -64.36146 -113.58638  -93.45456

En este caso, el modelo de Briere-2 se ajusta mejor a los datos que los otros modelos. Ahora podemos calcular los valores de CTmin, CTmax, y Topt. Primero vamos a simular datos de temperatura desde 0 hasta 45 grados. Despues vamos a hacer una predicción de la tasa de desarrollo para esas temperaturas en base al modelo de Briere-2. Despues especificamos que los valores de tasa de desarrollo negativos son de cero (porque no pueden ser negativos), y que los valores que faltan tambien son de cero (en caso de que haya valores ausentes). El Topt es el valor de temperatura a la cual la tasa de desarrollo esta máxima, el CTmin es el valor de temperatura mínima a la cual observamos desarrollo y el CTmax el valor de temperatura máxima a la cual observamos desarrollo.

tempS <- seq(from = 0, to = 45, by = 0.1) # temperaturas simuladas
devRateS <- predict(modNoLin_02, newdata = list(T = tempS)) # predicciones
devRateS[devRateS < 0] <- 0
devRateS[is.na(devRateS)] <- 0
c(AIC(modNoLin_01), AIC(modNoLin_02), AIC(modNoLin_03))

## [1]  -64.36146 -113.58638  -93.45456

Topt <- tempS[devRateS == max(devRateS)]
CTmin <- tempS[devRateS == min(devRateS[devRateS > 0 & 
  tempS < Topt])]
CTmax <- tempS[devRateS == min(devRateS[devRateS > 0 & 
  tempS > Topt])]
cat(paste0("Topt: ", Topt, "\nCTmin: ", CTmin, "\nTmax: ", CTmax))

## Topt: 33.3
## CTmin: 9.2
## Tmax: 36.7

Construir modelos fenológicos

Mediante el uso de la tasa de desarrollo, es posible modelar el área de distribución teórica de una especie, es decir, el área en la que el desarrollo sería posible considerando la temperatura, así como las diferentes etapas y el número de generaciones potenciales por año de una especie; el voltinismo. El primer paso es tener datos de temperatura. El segundo paso consiste en utilizar la relación entre la temperatura y la tasa de desarrollo cuantificada previamente.

De hecho, para cada dato de temperatura, se puede calcular una tasa de desarrollo que, al acumularse, conducirá al ciclo completo de desarrollo de una especie determinada. Así obtenemos una distribución potencial basada en un modelo de fenológia, que puede modificarse mediante la introducción de escenarios de cambio climático como un aumento general de las temperaturas o la introducción de eventos extremos más frecuentes.

Para esto vamos a imaginar un entorno teórico. Este entorno tendra una temperatura sacada de una ley Normal de parámetro mu igual a 15 grados para el promedio, de parametro sigma igual a 1 para la desviación estandar, durante un período de 100 días.

Nuestro modelo de Bactrocera dorsalis solo es valido para los huevos, mientras que necesitamos un modelo para todos los estadios fenologicos. Vamos a usar entonces otros datos sacados de la literatura. Los datos seran de T. solanivora, una polilla de la papa (Lepidoptera:Gelechiidae) del el articulo cientifico de Crespo-Perez et al. 2011⁵, usando Web Plot Digitizer⁶.

## cargar datos del laboratorio
datosLabTS_egg <- data.frame(
  temp = c(10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 16.0, 16.0, 17.0, 20.0, 20.0, 
    25.0, 25.0, 30.0, 30.0, 35.0), 
  devRate = c(0.031, 0.039, 0.072, 0.047, 0.059, 0.066, 0.083, 0.1, 0.1, 0.1, 0.143, 
    0.171, 0.2, 0.2, 0.18, 0.001))
datosLabTS_larva <- data.frame(
  temp = c(10.0, 10.0, 10.0, 13.0, 15.0, 15.5, 15.5, 15.5, 17.0, 20.0, 25.0, 
    25.0, 30.0, 35.0), 
  devRate = c(0.01, 0.014, 0.019, 0.034, 0.024, 0.029, 0.034, 0.039, 0.067, 0.05, 
    0.076, 0.056, 0.0003, 0.0002))
datosLabTS_pupa <- data.frame(
  temp = c(10.0, 10.0, 10.0, 13.0, 15.0, 15.0, 15.5, 15.5, 16.0, 16.0, 17.0, 
    20.0, 20.0, 25.0, 25.0, 30.0, 35.0), 
  devRate = c(0.001, 0.008, 0.012, 0.044, 0.017, 0.044, 0.039, 0.037, 0.034, 0.051, 
    0.051, 0.08, 0.092, 0.102, 0.073, 0.005, 0.0002))
## ajustar un modelo a los datos (Lactin-1)
## ver vignette quickUserGuide para mayor informacion
modTs01_egg <- devRateModel(
  eq = lactin1_95, 
  dfData = datosLabTS_egg, 
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))
modTs01_larva <- devRateModel(
  eq = lactin1_95, 
  dfData = datosLabTS_larva, 
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))
modTs01_pupa <- devRateModel(
  eq = lactin1_95, 
  dfData = datosLabTS_pupa, 
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))

A partir de los modelos ajustados a T. solanivora y una población de 50 individuos simulados, podemos representar la distribución de las
generaciones con el tiempo. Para agregar realismo a la simulación, añadimos variabilidad en la respuesta de los insectos a la temperatura que seguirá una distribución Normal de parametro mu igual al valor de la tasa de desarrollo y de parametro sigma = 0.025.

simul01 <- devRateIBM(
  tempTS = rnorm(n = 100, mean = 15, sd = 1),
  timeStepTS = 1,
  models = list(modTs01_egg, modTs01_larva, modTs01_pupa),
  numInd = 50,
  stocha = 0.025,
  timeLayEggs = 1)
print(simul01)

## [[1]]
##       g1s1 g1s2 g1s3 g2s1 g2s2
##  [1,]   14   37   56   70   93
##  [2,]   16   39   67   81   NA
##  [3,]   15   44   66   80   NA
##  [4,]   14   39   60   75  100
##  [5,]   14   37   64   77   NA
##  [6,]   14   45   68   84   NA
##  [7,]   16   43   71   85   NA
##  [8,]   16   41   62   77   NA
##  [9,]   17   45   71   85   NA
## [10,]   17   40   68   83   NA
## [11,]   13   42   66   81   NA
## [12,]   14   43   66   80   NA
## [13,]   13   42   69   81   NA
## [14,]   13   40   66   79   NA
## [15,]   15   38   69   85   NA
## [16,]   13   46   72   89   NA
## [17,]   14   44   74   89   NA
## [18,]   15   55   80   95   NA
## [19,]   15   40   67   80   NA
## [20,]   15   43   70   83   NA
## [21,]   15   42   69   83   NA
## [22,]   14   46   70   83   NA
## [23,]   14   42   72   85   NA
## [24,]   16   42   66   80   NA
## [25,]   14   38   63   77   98
## [26,]   15   47   68   82   NA
## [27,]   14   34   59   73   NA
## [28,]   14   43   71   84   NA
## [29,]   16   46   69   83   NA
## [30,]   15   35   58   72   99
## [31,]   14   47   69   84   NA
## [32,]   16   44   66   79   NA
## [33,]   14   43   65   79   NA
## [34,]   14   38   67   82   NA
## [35,]   15   44   68   83   NA
## [36,]   13   36   55   68   91
## [37,]   14   42   66   79   NA
## [38,]   14   40   63   78   NA
## [39,]   15   36   58   75   NA
## [40,]   15   46   73   85   NA
## [41,]   15   48   77   91   NA
## [42,]   14   41   64   79   NA
## [43,]   13   47   70   84   NA
## [44,]   14   41   63   78   NA
## [45,]   14   41   65   78   NA
## [46,]   13   42   62   75   NA
## [47,]   13   45   73   89   NA
## [48,]   12   35   54   71   98
## [49,]   15   43   67   82   NA
## [50,]   13   37   59   75   NA
## 
## [[2]]
## [[2]][[1]]
## Nonlinear regression model
##   model: rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT)
##    data: data.frame(rT = devRate, T = temp)
##      aa    Tmax  deltaT 
##  0.1583 34.9875  6.3033 
##  residual sum-of-squares: 0.00348
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 1.475e-06
## 
## [[2]][[2]]
## Nonlinear regression model
##   model: rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT)
##    data: data.frame(rT = devRate, T = temp)
##     aa   Tmax deltaT 
##  0.106 34.363  9.401 
##  residual sum-of-squares: 0.00382
## 
## Number of iterations to convergence: 16 
## Achieved convergence tolerance: 1.317e-06
## 
## [[2]][[3]]
## Nonlinear regression model
##   model: rT ~ exp(aa * T) - exp(aa * Tmax - (Tmax - T)/deltaT)
##    data: data.frame(rT = devRate, T = temp)
##      aa    Tmax  deltaT 
##  0.1224 34.4675  8.1460 
##  residual sum-of-squares: 0.008379
## 
## Number of iterations to convergence: 14 
## Achieved convergence tolerance: 9.492e-07
## 
## 
## [[3]]
##   [1] 15.16088 16.13961 15.10242 15.16444 12.69104 14.77750 15.13686 15.53528
##   [9] 12.70756 16.88487 14.46583 14.49209 15.06710 15.34523 14.62479 16.10230
##  [17] 13.26523 15.26766 14.91759 16.13879 16.80097 15.49237 16.40947 15.24969
##  [25] 16.07975 14.10906 14.63641 14.78356 16.14845 15.40372 14.65961 15.18095
##  [33] 14.90515 15.77565 16.70333 14.51308 14.99533 15.75146 16.70659 15.75039
##  [41] 15.08844 14.98185 15.29598 15.75244 14.56420 15.04183 16.03266 15.01618
##  [49] 15.82182 14.31113 15.43161 14.38182 13.94936 14.41584 13.89570 14.79407
##  [57] 13.67214 16.07015 13.84023 13.67259 15.13174 14.32570 14.21809 15.31591
##  [65] 16.02039 16.30514 16.06390 16.03364 14.36303 15.01091 14.28782 14.76062
##  [73] 16.53580 15.62920 15.34975 16.23871 15.37552 14.59470 15.68155 16.26867
##  [81] 14.76893 14.45899 15.65342 16.10102 14.40558 14.38435 14.79235 15.97309
##  [89] 14.52336 15.12639 15.04026 14.94276 12.34175 15.81280 15.17159 15.63882
##  [97] 14.04151 16.43450 15.33599 13.74792

par(mar = c(4, 4, 0, 0))
devRateIBMPlot(ibm = simul01)

## Threshold for visualization = 10% of individuals

En los resultados obtenemos la fenología de cada uno de los 50 individuos, el modelo utilizado, y la serie temporal de temperatura. En la visualización de los resultados, cada color representa una generación diferente y cada tipo de línea una etapa fenológica diferente.

Dado que estamos interesados en la respuesta de los insectos al cambio climático, es posible modificar el escenario inicial, es decir, modificar la serie temporal de temperaturas introduciendo un aumento en las temperaturas, por ejemplo pasando de un promedio de 15 a 17 grados.

simul02 <- devRateIBM(
  tempTS = rnorm(n = 100, mean = 17, sd = 1),
  timeStepTS = 1,
  models = list(modTs01_egg, modTs01_larva, modTs01_pupa),
  numInd = 50,
  stocha = 0.025,
  timeLayEggs = 1)
par(mar = c(4, 4, 0, 0))
devRateIBMPlot(ibm = simul02)

## Threshold for visualization = 10% of individuals

Tambien podemos aumentar la variabilidad de las temperaturas al cambiar la desviación estándar de la ley normal que representa las temperaturas, y pasar sigma de 1 a 2.

simul02 <- devRateIBM(
  tempTS = rnorm(n = 100, mean = 17, sd = 2),
  timeStepTS = 1,
  models = list(modTs01_egg, modTs01_larva, modTs01_pupa),
  numInd = 50,
  stocha = 0.025,
  timeLayEggs = 1)
par(mar = c(4, 4, 0, 0))
devRateIBMPlot(ibm = simul02)

## Threshold for visualization = 10% of individuals

Variando los parámetros incluso ligeramente, podemos ver que las consecuencias pueden ser importantes para la fenología de los insectos, a pesar de la simplicidad del modelo aquí representado. Podemos modificar el modelo para cambiar el tamaño de las poblaciones, aumentar el número de días de la simulación o modificar la variabilidad intrapoblacional en la respuesta a la temperatura para profundizar estos elementos y observar la complejidad de la respuesta de los insectos a los cambios globales.

También podemos usar un modelo diferente al utilisado y resaltar la importancia de la elección del mismo, y concluir que las incertidumbres son múltiples:

Incertidumbres relacionadas al cambio climático: ¿cuáles serán las temperaturas que vendrán y cuánto variarán estas temperaturas?
Incertidumbres sobre los insectos: ¿cuál es la variabilidad intrapoblacional de la respuesta de los insectos a la temperatura? ¿Cuál es el efecto de las temperaturas fluctuantes? ¿Cómo transferir estudios de laboratorio a condiciones externas?
Incertidumbres relacionadas con la caracterización de la respuesta de los insectos a la temperatura y en particular las incertidumbres relacionadas con la elección del modelo matemático.

Aquí el modelo presentado es bastante simple, podría enriquecerse teniendo en cuenta la variabilidad espacial en las temperaturas. En el ejemplo que sigue vamos a crear un entorno teórico para representar un mapa donde a cada punto corresponde un dato de temperatura que vamos a usar para calcular el numero de generaciones potenciales simulando 100 días a esta temperatura.

tempEspacio <- matrix(rnorm(100, mean = 15, sd = 1), ncol = 10) # mapa teorico
myDevRate <- 1/devRateMap(nlsDR = modTs01_egg, tempMap = tempEspacio) +
  1/devRateMap(nlsDR = modTs01_larva, tempMap = tempEspacio) +
  1/devRateMap(nlsDR = modTs01_pupa, tempMap = tempEspacio)
filled.contour(100 / myDevRate, # numero de generaciones en 100 dias
    col = rev(heat.colors(30)), 
    main = "#generaciones", plot.axes = FALSE)

Tambien se podria incorporar la respuesta de los insectos a otros factores como la humedad relativa, o la disponibilidad y calidad del recurso disponible para estos.

Límites

Aunque las predicciones de estos modelos son generalmente correctas a escala regional o continental, las predicciones a escalas más locales se enfrenten con las limitaciones de estos modelos:

solo se tiene en cuenta la temperatura, mientras que otros factores son importantes para el desarrollo de insectos, como la humedad relativa, el fotoperiodo o la fuente y la disponibilidad de comida
las curvas de desarrollo dependientes de la temperatura generalmente se establecen a temperaturas constantes de laboratorio, pero en condiciones externas las temperaturas son fluctuantes, lo que tiene un efecto en el desarrollo, incluso cuando las temperaturas no exceden los valores umbral de desarrollo. Este punto es aún más importante cuando sabemos que el cambio climático no solo influirá en el aumento de la temperatura promedia, sino que también exacerbará los extremos, con temperaturas que fluctuarán aún más
las tasas de desarrollo calculadas a temperaturas cercanas a Ctmin y Ctmax pueden variar significativamente de un modelo a otro, con la incertidumbre sobre el efecto de las temperaturas cercanas a estos valores umbral
dentro de la misma especie y en la misma etapa fenologica, hay una variabilidad en la respuesta de los insectos a la temperatura que no es tomada en cuenta por los modelos que cuantifican la relación entre la temperatura y la tasa de desarrollo. Los modelos mecanísticos de la fenología de insectos tenderán a sobreestimar o subestimar el desarrollo de una parte de la población.

Estos son los desafíos a los que se enfrentan los estudios de este tipo hoy en día y que hay que conocer antes de interpretar los resultados de los modelos presentados en esta guía.

Conclusión

Los modelos que relacionan temperatura y tasa de desarrollo son la base de casi todos los modelos fenologicos predictivos que se pueden encontrar en la literatura hoy en día. Se ha escrito la presente guía con el proposito de ayudar a analizar sus datos, interpretar los resultados, y construir sus modelos fenologicos conociendo las fortalezas y debilidades de este enfoque. Cualquier comentario o pregunta son bienvenida(o)s, no dude en ponerse en contacto conmigo.

https://en.wikipedia.org/wiki/Ren%C3%A9_Antoine_Ferchault_de_R%C3%A9aumur ↩︎
https://en.wikipedia.org/wiki/Bactrocera_dorsalis ↩︎
Messenger, P.S., and Flitters, N.E. (1958). Effect of constant temperature environments on the egg stage of three species of Hawaiian fruit flies. Annals of the Entomological Society of America, 51(2), 109-119.↩︎
Campbell, A., Frazer, B.D., Gilbert, N.G.A.P., Gutierrez, A.P., Mackauer, M. (1974). Temperature requirements of some aphids and their parasites. Journal of applied ecology, 431-438.↩︎
Crespo-Pérez, V., Rebaudo, F., Silvain, J.-F. Dangles, O. (2011) Modeling invasive species spread in complex landscapes: the case of potato moth in Ecuador. Landscape Ecology, 26, 1447–1461.↩︎
Rohatgi, A. (2015) WebPlotDigitalizer: HTML5 Based Online Tool to Extract Numerical Data from Plot Images.↩︎