Overview

In DEB parameter estimation it is critical to get a good estimate of the body temperature at which growth, reproduction, development and aging rates are observed. A measure of mean ‘ambient’ temperature will often be inaccurate, leading to poor parameter estimates. For example, often the observations of these rates are made for organisms experiencing fluctuating body temperatures, as typically occurs in nature or if animals are cultured under cycling temperature regimes. Moreover, the temperatures experienced may be above ambient air temperature if the organism is in the boundary layer near the ground, exposed to sunlight or high levels of infrared radiation. It may also be below air temperature if it is cooling by evaporation (e.g. a frog) or to cold radiative environments (e.g. clear night skies).

This document shows how to estimate the body temperatures of animals in nature by two methods of increasing sophistication. First, it shows how to extract 6-hourly air temperature data for a location of interest (e.g. from where field growth rates have been observed) from a global database of historical daily weather data via the RNCP package (Kemp et al. 2012). Second it shows how to use the NicheMapR package to estimate microclimates (e.g. near-ground air temperatures, soil temperatures) and to simulate heat exchange for an organism of given physical characteristics (e.g. solar reflectance) and behavioural thermoregulatory responses.

In each case the resulting temperature vector is converted to a ‘constant temperature equivalent’, CTE, which can be used as the temperature estimate for zero-variate observations of e.g. age at birth, age at puberty, longevity and reproduction rate.

For the example we will use the green anole, Anolis carolinensis, which is in the AmP collection and has field growth observations for a number of locations in Russia (Roitberg and Smirina 2006).

Package Installation

First install and load the necessary packages (if you don’t have them already).

install.packages("rvest")
install.packages("RNCEP")
install.packages("ncdf4")
install.packages("devtools")
install.packages("stringr")
devtools::install_github("mrke/NicheMapR")
# to install with devtools on Windows you'll need rtools:
# https://cran.r-project.org/bin/windows/Rtools/ you can
# also download a binary from
# https://github.com/mrke/NicheMapR/releases/

Load the libraries.

library(rvest)
library(RNCEP)
library(ncdf4)
library(stringr)
library(NicheMapR)

Extracting NCEP air temperature for a location and time

The RNCEP package interfaces with two gridded (~2.5 degree) global data set on various meteorological variables that go back to as far as 1949 and as recently as the past month. You can read more about these data sets in (Kemp et al. 2012).

Here we are quering the ‘air.2m’ variable, which is the estimate of air temperature at the typical weather station height. There are a number of other options you could choose.

We need to specify the location as longitude and latitude (decimal degrees), as well as the starting and ending year and the starting and ending months. Note that you ONLY get the months you ask for, so if you’re trying to get a continuous period that spans a year, you should do separate queries per year and combine them.

# define query inputs
lonlat <- c(-94.656, 31.604)  # the Nacogdoches site of Michael (1972)
y.start <- 1966
y.finish <- 1970
m.start <- 1
m.finish <- 12
var <- "air.2m"
# you could also choose e.g. 'skt.sfc' (surface
# temperature) or 'tmp.0-10cm' (shallow soil)
lev <- "gaussian"

# run the query
tair <- NCEP.gather(variable = var, level = lev, months.minmax = c(m.start,
    m.finish), years.minmax = c(y.start, y.finish), lat.southnorth = c(lonlat[2],
    lonlat[2]), lon.westeast = c(lonlat[1], lonlat[1]), reanalysis2 = FALSE,
    return.units = TRUE, status.bar = FALSE)
## [1] Units of variable 'air.2m' are degK
# extract the actual coordinates you got data from (given
# that the grid is coarse)
lonlat.returned <- c(as.numeric(dimnames(tair)[[2]][2]), as.numeric(dimnames(tair)[[1]][1]))

# dates
dates1 <- seq(ISOdate(y.start, m.start, 1, tz = "UTC") - 3600 *
    12, ISOdate(y.finish + 1, 1, 1, tz = "UTC") - 3600 * 13,
    by = "hours")
dates <- dates1[seq(1, length(dates1), 6)]  # cut down to 6-hourly

t.air <- cbind(as.data.frame(dates), tair[1, 1, ])
colnames(t.air) <- c("date", var)
plot(t.air - 273.15, type = "l", ylim = c(-30, 50), ylab = "2m air temperature")  # plot, converting to deg C

# write to disk if you want write.csv(tair[1, 1, ],
# 'tair.csv')

Compute constant temperature equivalent (CTE) based on air temperature

Now we want to compute a ‘constant temperature equivalent’, CTE, which is an estimate of the mean temperature accounting for the nonlinearity of the thermal response curve.

First, some parameters are defined for the 5-parameter Arrhenius response curve. These are a rough guess for Anolis carolinensis.

T_A <- 9625
T_AL <- 20000
T_AH <- 20000
T_L <- 10 + 273.15
T_H <- 40 + 273.15
T_REF <- 20 + 273.15

Now we define the temperature vector of interest, T, as the 2 m air temperature NCEP data, and use the 5 parameter Arrhenius formula to convert it to the temperature correction factor, c_T.

T_b <- t.air$air.2m  # temps from simulation to be used to estimate CTE
Tjuly <- mean(T[which(format(dates, "%m") == "07")] - 273.15)  # get mean July air temperature

# temperature correction function
get_c_T <- function(T_A, T_L, T_H, T_AL, T_AH, T_REF, T_b) {
    c_T <- exp(T_A/(T_REF) - T_A/T_b) * (1 + exp(T_AL/T_REF -
        T_AL/T_L) + exp(T_AH/T_H - T_AH/T_REF))/(1 + exp(T_AL/T_b -
        T_AL/T_L) + exp(T_AH/T_H - T_AH/T_b))
}

# convert T_b each hour to temperature correction factor
c_T <- get_c_T(T_A, T_L, T_H, T_AL, T_AH, T_REF, T_b)
plot(c_T ~ t.air$date, type = "l", xlab = "date", ylim = c(0,
    5))

Now we compute the mean c_T and define a function for use in iteratively guessing values of body temperature that result in this mean c_T value, i.e. we turn the mean temperature correction factor back into a temperature. The ‘uniroot’ function is used to do the guessing, with bracketing starting T_L and T_H.

c_T.mean <- mean(c_T)  # get mean temperature correction factor

getCTE <- function(T_b) {
    # function finding the difference between a temperature
    # correction factor, c_T.mean, for a specified T_b
    # compared to the mean calculated one (aim to make this
    # zero)
    x <- get_c_T(T_A, T_L, T_H, T_AL, T_AH, T_REF, T_b) - c_T.mean
}
# search for a T_b (CTE) that gives the same temperature
# correction factor as the mean of the simulated
# temperature corrections
CTE <- uniroot(f = getCTE, c(T_L, T_H), check.conv = TRUE)$root -
    273.15  # converting from K to deg C

We now have the constant temperature equivalent. As you can see from the values below, the naïve temperature correction factor we’d get from the mean of the raw input air temperatures (17.37 °C) is around 0.47 but transformation of the 6-hourly temperatures into c_T gives a mean of 0.5. This corresponds to a CTE of 21.3 °C, i.e. 4 °C higher than the mean air temperature.

{
    cat(paste0("mean T_air = ", round(mean(T_b - 273.15), 2)),
        "\n")  # report value to console
    cat(paste0("mean c_T = ", round(c_T.mean, 2)), "\n")  # report value to console
    cat(paste0("naïve c_T = ", round(as.numeric(exp(T_A/(T_REF) -
        T_A/mean(T_b)) * (1 + exp(T_AL/T_REF - T_AL/T_L) + exp(T_AH/T_H -
        T_AH/T_REF))/(1 + exp(T_AL/mean(T_b) - T_AL/T_L) + exp(T_AH/T_H -
        T_AH/mean(T_b)))), 2)), "\n")
    # report value to console
    cat(paste0("CTE = ", round(CTE, 2)), "\n")  # report value to console
}
## mean T_air = 17.37 
## mean c_T = 1.18 
## naïve c_T = 0.7 
## CTE = 21.33

Using NicheMapR to estimate body temperature including thermoregulatory behaviour

The use of fluctuating air temperature is clearly making a difference to the temperature we would associate with field observations for this lizard. The next section shows an example of how to do this using the NicheMapR package (Kearney and Porter, 2017; Kearney and Porter 2020). This package includes a microclimate model, which computes aboveground and belowground conditions relevant to solving heat and water budgets of organisms (radiation, wind, humidity, air temperature). It also includes an ectotherm model of heat and water exchange in the context of behavioural regulation under spatially implicit microclimatic variation (Kearney and Porter 2020). I.e. extremes of available shade, and the range of possible depths, are provided as well as body temperature thresholds for different behavioural responses, and the model computes an hourly temperature trajectory given the microclimate.

Here we run the microclimate function ‘micro_ncep’ which forces the microclimate model via the RNCEP package used above, but also extracting information on humidity, wind speed, cloud cover, air pressure and solar radiation from the NCEP data. Alternatively, we could run the microclimate function ‘micro_usa’ which links to a historical daily gridded data set for the USA (and is much faster to obtain from the web).

# dstart <- '01/01/1966' dfinish <- '31/12/1970' micro <-
# micro_ncep(loc = lonlat, dstart = dstart, dfinish =
# dfinish, reanalysis = FALSE) save(micro, file =
# 'micro_Nacogdoches_66_77.Rda')
load("micro_Nacogdoches_66_77.Rda")

# dstart <- '01/01/1980' dfinish <- '31/12/1981' micro <-
# micro_usa(loc = lonlat, dstart = dstart, dfinish =
# dfinish) save(micro, file =
# 'micro_Nacogdoches_80_81_gridmet.Rda')

Next we have to specify the parameters determining how the species behaviourally thermoregulates. Brief explanations are provided next to each one.

Ww_g <- 3.5  # wet weight, grams
CT_max <- 36  #42    # critical thermal maximum, Wilson and Echternacht 1990
CT_min <- 11  # critical thermal minimum, Wilson and Echternacht 1987
T_F_max <- 36  # maximum T_b at which activity occurs, Wilson and Echternacht 1990, Jenssen et al. 1996
T_F_min <- 25  # minimum T_b at which activity occurs
T_pref <- 33.8  # preferred T_b (will try and regulate to this)
T_B_min <- 20  # minimum T_b at which leaves retreat
T_RB_min <- 15  # minimum T_b at which they will attempt to leave retreat
shape <- 3  # shape (3 = lizard)
shade_seek <- 1  # shade seeking?
burrow <- 1  # can it move below ground?
climb <- 1  # can it climb?
mindepth <- 3  # minimum depth (node, 1-10) allowed
maxdepth <- 10  # maximum depth (node, 1-10) allowed
nocturn <- 0  # nocturnal activity?
crepus <- 0  # crepuscular activity?
diurn <- 1  # diurnal activity?

Finally we run the ectotherm model and plot the results. In the figure you can see in grey the air temperature, as extracted via RNCEP above, but also in black the predicted temperature incorporating the thermoregulatory response, which includes moving underground to escape extremes of heat and cold, and choosing shade to attempt to stay close to the preferred range. In brown is the depth selected (cm) and the orange indicates the hours that foraging was possible. You can see the lizard is warmer in the winter due to being underground, and can also get warmer in the active season due to exposure to solar radiation.

ecto <- ectotherm(nocturn = nocturn, diurn = diurn, shape = shape,
    T_F_max = T_F_max, T_F_min = T_F_min, T_pref = T_pref, T_B_min = T_B_min,
    T_RB_min = T_RB_min, CT_max = CT_max, CT_min = CT_min, crepus = crepus,
    burrow = burrow, shade_seek = shade_seek, mindepth = mindepth,
    maxdepth = maxdepth)
dates <- micro$dates
sub <- which(format(micro$dates, "%Y") == "1967")
# sub <- which(format(micro$dates, '%Y') == '1981')
environ <- as.data.frame(ecto$environ)
metout <- as.data.frame(micro$metout)
# plot some results
plot(dates[sub], metout$TAREF[sub], ylim = c(-30, 40), type = "l",
    col = "lightblue", ylab = "temperature, deg C", xlab = "time")
points(dates[sub], environ$TC[sub], type = "l")
points(dates[sub], environ$ACT[sub], type = "h", col = "orange")
points(dates[sub], environ$DEP[sub], type = "h", col = "brown")

Now we repeat the procedure above for getting the CTE, using estimated body temperature as our input vector.

T_b <- environ$TC + 273.15  # temps from simulation to be used to estimate CTE
# convert T_b each hour to temperature correction factor
c_T <- get_c_T(T_A, T_L, T_H, T_AL, T_AH, T_REF, T_b)
plot(c_T ~ dates, type = "l", xlab = "date")

c_T.mean <- mean(c_T)  # get mean temperature correction factor
# search for a T_b (CTE) that gives the same temperature
# correction factor as the mean of the simulated
# temperature corrections
CTE <- uniroot(f = getCTE, c(T_L, T_H), check.conv = TRUE)$root -
    273.15  # converting from K to deg C
{
    # report values to console
    cat(paste0("mean T_b = ", round(mean(T_b) - 273.15, 2)),
        "\n")
    cat(paste0("mean c_T = ", round(c_T.mean, 2)), "\n")
    cat(paste0("naïve c_T = ", round(as.numeric(exp(T_A/(T_REF) -
        T_A/mean(T_b)) * (1 + exp(T_AL/T_REF - T_AL/T_L) + exp(T_AH/T_H -
        T_AH/T_REF))/(1 + exp(T_AL/mean(T_b) - T_AL/T_L) + exp(T_AH/T_H -
        T_AH/mean(T_b)))), 2)), "\n")
    cat(paste0("CTE = ", round(CTE, 2)), "\n")  # report value to console
}
## mean T_b = 25.01 
## mean c_T = 2.15 
## naïve c_T = 1.79 
## CTE = 26.82

When we incorporate behaviour, we end up with a higher raw input temperature of 25.0 °C and a higher CTE of 26.8 °C. This is very close to the value of 27 °C used in the AMP parameter estimation.

References

Kearney, M. R., and W. P. Porter. 2017. NicheMapR - an R package for biophysical modelling: the microclimate model. Ecography 40:664–674.

Kearney, M. R., & Porter, W. P. (2020). NicheMapR – an R package for biophysical modelling: The ectotherm and Dynamic Energy Budget models. Ecography, 43(1), 85–96. https://doi.org/10.1111/ecog.04680

Kemp, M. U., E. Emiel van Loon, J. Shamoun-Baranes, and W. Bouten. 2012. RNCEP: global weather and climate data at your fingertips: RNCEP. Methods in Ecology and Evolution 3:65–70.

Michael, E. D. (1972). Growth Rates in Anolis carolinensis. Copeia, 1972(3), 575. https://doi.org/10.2307/1442932