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

Francois Rebaudo (IRD, laboratorio EGCE)

2018-04-30

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érmicos, 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 de vida. Se oponen a los organismos endotérmicos, como los mamíferos, que producen principalmente su propia fuente de calor para regular su temperatura.

En la mayoría de los ambientes, la temperatura varía con 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. La temperatura de los ectotermicos evoluciona en ambientes donde los cambios de temperatura varían, a veces considerablemente. Hablamos entonces de organismos poikilotérmicos, entre los que encontramos la mayoría de los insectos, que contrastamos con los organismos homotérmicos 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 la temperatura no es el único factor involucrado. Por ejemplo, se puede mencionar 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 un preferencio térmico, es decir un valor o conjunto de valores de temperatura, para los cuales los insectos pueden alcanzar su desarrollo óptimo. Este preferendo es variable dependiendo de la especie, con especies que se desarrollarán en un amplio rango de temperaturas, se denomina especies generalistas y, a diferencia de las especies que pueden desarrollarse solo en un rango muy estrecho de temperaturas, estas son las llamadas especies especializadas. Los insectos se encuentran a lo largo de este gradiente entre especialistas y generalistas, con variaciones que también se pueden observar dentro de una especie, especialmente entre las diferentes etapas de desarrollo. Para entender estas diferentes estrategias y caracterizarlas, se hace necesario estudiar la relación entre la temperatura y el desarrollo de los insectos.

Aqui nos interesarán los insectos poikilotérmicos, que, les recuerdo, corresponden a los insectos cuya temperatura cambiará como consecuencia de los cambios de la temperatura exterior. Estudiaremos más particularmente la tasa de desarrollo de insectos, generalmente expresada como la inversa de la cantidad de días que pasan desde una etapa de desarrollo a otra. Para este proposito, usaremos datos de la literatura.

En la mayoría de los ecosistemas, las condiciones ambientales están cambiando. Este es el caso de la temperatura. Para los insectos poikilotérmicos, esto dará como resultado un cambio en la temperatura de su cuerpo. Estos mismos insectos tendrán un rango de temperatura para el cual pueden desarrollarse, y más allá de lo cual 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 ya no pueden desarrollarse, que se encuentra en la literatura bajo los acrónimos CTmin y CTmax para Critical Thermal minima Critical Thermal maxima, respectivamente. Estos límites varían ampliamente de una especie a otra, con algunas especies como el escarabajo ártico capaz de soportar hasta -70 grados, mientras que otras especies 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 médicos generales a especialistas, también varía mucho de una especie a otra e incluso dentro de una especie entre las diferentes etapas fenológicas. Más allá de estos dos valores umbral, la respuesta de los insectos variará según el tiempo de exposición, la frecuencia de exposición, la velocidad a la que la temperatura se lleva más allá de los valores umbral, la tasa de exposición, la historia de vida de los insectos y, por supuesto, las especies consideradas. En este 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 ha estudiado y publicado sobre el efecto de la temperatura en los insectos es René-Antoine Ferchault de Réaumur1: sus observaciones sobre las mariposas se registrarán en sus “Mémoires pour servir à l’histoire des insectes” entre 1734 y 1742. Despues de Réaumur, muchos científicos buscaran 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 a colocar insectos en jaulas a diferentes temperaturas constantes, para medir el tiempo de desarrollo de cada etapa fenológica. Los científicos luego observarán varias características:

  1. La primera es que es muy difícil, si no imposible, 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
  2. La segunda característica es que hay variabilidad dentro de la misma población de insectos en cuanto a su respuesta a la temperatura
  1. 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

Más allá de estas observaciones, 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 analisar 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 2, un insecto del Orden de los Dipteros, que hemos copiados desde el articulo cientifico de Messenger and Flitters (1958)3. 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 tener 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. Por lo tanto 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 trbajar 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 analisar los datos de tasa de desarrollo de B. dorsalis en 
### funcion 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 almacenadas 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
## [15] 92.5 95.0 96.0 97.0 97.5

Los numeros usan el punto como separator de decimales y la coma como separator de numeros. 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
##  [8] 19.72222 21.11111 23.88889 26.66667 29.44444 30.83333 32.22222
## [15] 33.61111 35.00000 35.55556 36.11111 36.38889

Para guardar las temperaturas, podemos almacenarlas en un objeto. Para afectar las valores de temperaturas a un objeto, se usa los simbolos <-. Si queremos guardar las 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 las valores de tasa de desarrollo en un objeto llamado devRate. La tasa de desarrollo corresponde al opuesto del numero de días. Aquí las valores estan expresadas en horas, entonces dividimos el numero de horas por 24 y despues calcualmos 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 funcion 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 funcion 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 son listos para empezar el analisis 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 de 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. 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 estar diferente del CTmin. La temperatura base esta 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.

El uso del modelo lineal fue popularisado con el articulo cientifico de Campbell et al. (1974)4. Este modelo es valido dentro de la zona de linearidad de la respuesta de los insectos a la temperatura. Para seleccionar los datos que vamos a usar, puede ser util visualisar los datos de laboratorio en un grafico. Para esto vamos a usar la funcion 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 los ultimos cinco 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 funcion lm, vamos a usar la funcion generica del paquete devRate.

modLin <- devRateModel(eq = campbell_74, df = 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

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 modelado de grados-días se basa en la suposición de que los insectos necesitan una suma de temperaturas más allá de la temperatura base para completar su desarrollo. Este guía no discutirá los numerosas métodos disponibles para calcular los días de grado, pero 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 dividirlo el resultado de la suma por dos para luego eliminar de 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 por 2 es igual a 20; menos 3 es igual a 17 grados. La suma de los grados-días se obtiene sumando los grados de los valores de día para cada nuevo día. 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 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:

Podemos buscar en la base de datos del packete 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 hay 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 en epidemiología. 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 baja temperatura, en el área entre la temperatura base y la temperatura umbral d e 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 un modelo no hace 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ísticas 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 sustanciales 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 a los conjuntos de datos experimentales, y seleccionar el que mejor responda la problemática del estudio, y compararlo con los criterios de comparación estadísticas. Una posibilidad es de buscar en la base de datos del paquete cuales son los mas usados con la funcion devRateFind pero como hemos visto anterioramente no hay modelos no lineales que fueron usados al momento de escribir esta guía. Una alternativa es de usar los modelos mas usados como los modelos de Briere, Logan, y Lactin. Para esto, hay que definir valores inciales a los parametros de los modelos para que el aloritmo que va a hacer el ajuste encuentra una solucion. En caso de que no se sabe valores pertinentes para iniciar el algoritmo, se puede usar las valores promedias 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
  df = 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.693e-06
modNoLin_02 <- devRateModel(
  eq = briere2_99, # nombre del modelo
  df = 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.451e-06
modNoLin_03 <- devRateModel(
  eq = lactin2_95, # nombre del modelo
  df = 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.62e-06

Los resultados se peuden visualisar graficamente con la funcion devRatePlot.

par(mfrow = c(1, 3)) # para hacer tres graficos en la misma pagina
devRatePlot(eq = briere1_99, 
  nlsDR = modNoLin_01, 
  temp = datosLab$temp, 
  devRate = datosLab$devRate, 
  xlim = c(10, 40), ylim = c(0, 1.2))
devRatePlot(eq = briere2_99, 
  nlsDR = modNoLin_02, 
  temp = datosLab$temp, 
  devRate = datosLab$devRate, 
  xlim = c(10, 40), ylim = c(0, 1.2))
devRatePlot(eq = lactin2_95, 
  nlsDR = modNoLin_03, 
  temp = datosLab$temp, 
  devRate = datosLab$devRate, 
  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 selccionados.

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 prediccion de tasa de desarrollo para esas temperaturas en base al modelo de Briere-2. Despues especificamos que las valores de tasa de desarrollo negativas son de cero (porque no puede ser negativas), y que las valores que faltan tambien son de cero (en caso de que haya valores ausentes). El Topt es el valor de temperatura al cual la tasa de desarrollo esta maxima, el CTmin es el valor de temperatura minima al cual observamos desarrollo y el CTmax el valor de temperatura maxima al 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 fenologicos

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 teniendo en cuenta solo la temperatura, así como que 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. Luego, el segundo paso consiste en utilizar la relación entre la temperatura y la tasa de desarrollo cuantificada previamente.

De hecho, para cada dato datos de temperatura, se puede calcular una tasa de desarrollo que, al acumularse, conducirá al desarrollo completo de una especie determinada. Así obtenemos una distribución potencial basada en un modelo de fenológica, 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. El entorno teorico tendra una temperatura sacada de una ley Normal de parámetro mu igual a 15 grados para el promedio e de parametro sigma igual a 1 para la deviacion estandar, durante un período de 100 días.

Nuestro modelo de Bactrocera dorsalis solo es valido para los huevos, cuando nesesitamos un modelo para todos los estadios fenologicos. Vamos a usar otros datos sacados de la literatura. Vamos a usar los datos de T. solanivora, una polilla de la papa (Lepidoptera:Gelechiidae) desde el articulo cientifico de Crespo-Perez et al. 20115, usando Web Plot Digitizer6.

## 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, 
  df = datosLabTS_egg, 
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))
modTs01_larva <- devRateModel(
  eq = lactin1_95, 
  df = datosLabTS_larva, 
  startValues = list(aa = 0.177, Tmax = 36.586, deltaT = 5.631))
modTs01_pupa <- devRateModel(
  eq = lactin1_95, 
  df = 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, agregamos 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,]   15   48   73   89   NA
##  [2,]   13   42   67   79   NA
##  [3,]   12   43   66   82   NA
##  [4,]   11   39   58   73   96
##  [5,]   13   43   66   78   NA
##  [6,]   16   45   70   83   NA
##  [7,]   14   43   69   84   NA
##  [8,]   12   41   65   79   NA
##  [9,]   15   45   69   80   NA
## [10,]   13   41   66   81   NA
## [11,]   14   46   69   84   NA
## [12,]   16   42   66   80   NA
## [13,]   15   41   63   78   NA
## [14,]   15   43   68   83   NA
## [15,]   14   42   66   79   NA
## [16,]   15   41   66   81   NA
## [17,]   12   42   69   82   NA
## [18,]   14   44   69   86   NA
## [19,]   14   43   66   79   NA
## [20,]   14   42   67   80   NA
## [21,]   14   38   60   76   NA
## [22,]   12   37   66   77   NA
## [23,]   15   51   79   93   NA
## [24,]   17   50   76   90   NA
## [25,]   15   42   64   79   NA
## [26,]   12   37   56   70   99
## [27,]   15   44   68   82   NA
## [28,]   15   44   69   84   NA
## [29,]   15   47   66   78   NA
## [30,]   14   37   60   77   NA
## [31,]   14   42   60   75   NA
## [32,]   15   43   68   83   NA
## [33,]   15   39   63   75   NA
## [34,]   12   39   61   74   NA
## [35,]   15   44   68   81   NA
## [36,]   12   42   67   79   NA
## [37,]   12   53   76   90   NA
## [38,]   13   40   61   75   NA
## [39,]   11   36   59   74   NA
## [40,]   12   35   55   70   96
## [41,]   12   34   55   71  100
## [42,]   14   44   68   83   NA
## [43,]   14   44   70   88   NA
## [44,]   14   39   60   75   NA
## [45,]   13   42   65   82   NA
## [46,]   13   39   62   75   NA
## [47,]   15   48   73   90   NA
## [48,]   15   41   67   81   NA
## [49,]   14   38   61   74  100
## [50,]   15   43   65   79   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: 2.428e-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.817e-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.1461 
##  residual sum-of-squares: 0.008379
## 
## Number of iterations to convergence: 15 
## Achieved convergence tolerance: 2.684e-06
## 
## 
## [[3]]
##   [1] 15.18456 15.14433 16.19886 15.24167 16.29651 13.07006 16.28848
##   [8] 14.64102 13.13913 15.24720 13.93365 17.05948 13.17555 17.29715
##  [15] 15.42675 15.29738 15.23078 14.34221 14.70408 14.22697 12.74839
##  [22] 15.23978 13.33497 13.78270 15.09584 13.48274 15.70958 14.88342
##  [29] 15.89265 15.62743 13.33883 16.16068 14.14446 15.82007 14.93024
##  [36] 14.65255 14.98744 14.57426 15.49973 14.74779 16.28102 15.45959
##  [43] 16.34952 15.85349 15.58938 15.09143 14.81653 14.31061 15.88767
##  [50] 16.08188 16.02635 14.02566 15.40586 13.05218 15.80535 15.89373
##  [57] 15.78840 12.54592 15.94988 16.36732 15.37395 12.67563 14.37052
##  [64] 14.50978 14.60722 16.24739 17.06804 14.39566 16.94481 16.46284
##  [71] 15.41133 16.94884 14.49538 14.69488 14.02145 15.50982 14.06174
##  [78] 15.79245 14.43455 15.06353 16.40932 15.41470 16.40777 14.79806
##  [85] 16.49291 14.98735 14.33821 13.48703 14.35263 13.63938 13.81947
##  [92] 14.95415 15.27674 15.14415 14.65041 14.44791 14.81692 17.15782
##  [99] 16.00313 15.89464
par(mar = c(4, 4, 0, 0))
devRateIBMPlot(ibm = simul01)

## Threshold for visualization = 10% of individuals

En los resultados tenemos la fenología de cada uno de los 50 individuos, el modelo utilizado, luego 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 representado aquí. 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 presentado aquí y resaltar la importancia de la elección del modelo, y concluir que las incertidumbres son múltiples:

Están las incertidumbres relacionadas con el cambio climático: ¿cuáles serán las temperaturas que vendrán y cuánto variarán estas temperaturas? Hay 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? Y también las 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.

tempEspacio <- matrix(rnorm(100, mean = 15, sd = 1), ncol = 10) # mapa teorica
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 los insectos.

Limites

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

  1. 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
  2. las curvas de desarrollo dependientes de la temperatura generalmente se establecen a temperaturas constantes en el laboratorio, pero las condiciones externas son las de las temperaturas 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 la temperatura promedio en aumento, sino que también exacerbará los extremos, con temperaturas que fluctuarán aún más
  3. 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
  4. 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 hoy en día y que hay que conocer antes de interpretar los resultados de los modelos presentados en esta guia.

Conclusión

Esos modelos relacionando temperatura y tasa de desarrollo son la base de casi todos los modelos fenologicos predictivos que se pueden encontrar en la literaura. He escrito la presente guía con el proposito de ayudarles a analisar sus datos, interpretar los resultados, y construir sus modelos fenologicos conociendo las fortalezas y debilidades de este enfoque. Cualquier comentario o pregunta son bienvenidas, no duden en ponerse en contacto conmigo.


  1. https://en.wikipedia.org/wiki/Ren%C3%A9_Antoine_Ferchault_de_R%C3%A9aumur

  2. https://en.wikipedia.org/wiki/Bactrocera_dorsalis

  3. 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.

  4. 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.

  5. 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.

  6. Rohatgi, A. (2015) WebPlotDigitalizer: HTML5 Based Online Tool to Extract Numerical Data from Plot Images.