Thermal performance curves in R

Fitting thermal performance curves in R with nls: the optimum, the critical thermal maximum and breadth, and why a symmetric curve misplaces the upper limit.

Give an insect, a lizard or a leaf a body temperature and it performs at some rate: it grows, sprints, develops or photosynthesises. Plot that rate against temperature and you get a thermal performance curve, one of the workhorses of thermal ecology (Huey and Stevenson 1979 American Zoologist 19(1):357-366; Kingsolver 2009 The American Naturalist 174(6):755-768). The curve rises slowly from a lower limit, peaks at a thermal optimum, then falls off a cliff to an upper limit. That last word, cliff, is the whole point: performance curves are asymmetric, and pretending otherwise misplaces the numbers ecologists most want.

The shape, and a model for it

Two features define the curve. It is bounded, between a critical thermal minimum (CTmin) and maximum (CTmax) where performance reaches zero, and it is left-skewed: the rise is gentle, the fall is steep. A convenient model that captures both is the Briere function (Briere et al. 1999 Environmental Entomology 28(1):22-29), which is exactly zero outside the thermal limits and asymmetric in between.

briere <- function(T, a, T0, Tm) {
  p <- a * T * (T - T0) * sqrt(pmax(Tm - T, 0))
  p[T < T0 | T > Tm] <- 0
  p
}

Here T0 is CTmin, Tm is CTmax, and a scales the height. For a curve with CTmin 8 and CTmax 42 degrees C, the optimum falls at 34.5 degrees, which is 78% of the way from the lower limit to the upper one. The optimum sits close to the upper limit, not in the middle: an organism performing near its best is only a few degrees from the temperature that shuts it down.

Fitting a curve to data

Real data are performance measurements at a set of temperatures, with noise. Fit the Briere model with nls.

set.seed(391)
Td <- rep(seq(10, 40, 2.5), each = 4)
y  <- pmax(briere(Td, 8e-4, 8, 42) + rnorm(length(Td), 0, 0.06), 0)
fitB <- nls(y ~ aa * Td * (Td - t0) * sqrt(pmax(tm - Td, 0)),
            start = list(aa = 6e-4, t0 = 6, tm = 44))
round(coef(fitB), 4)

     aa      t0      tm 
 0.0008  8.1980 41.9207 

The fit recovers the optimum at 34.5 degrees and CTmax at 41.9, both close to the truth. From the fitted curve you read off the three quantities that matter: the optimum, the upper limit, and the thermal breadth, the width of the temperature band over which performance stays above eighty per cent of its peak, here about 12.2 degrees.

dd  <- data.frame(Td = Td, y = y)
cur <- data.frame(T = TgB, P = PgB)
ggplot(dd, aes(Td, y)) +
  annotate("rect", xmin = b80[1], xmax = b80[2], ymin = -Inf, ymax = Inf, fill = col_bri, alpha = 0.10) +
  geom_line(data = cur, aes(T, P), colour = col_bri, linewidth = 1) +
  geom_point(colour = col_dat, size = 2, alpha = 0.7) +
  geom_vline(xintercept = Topt_B, colour = col_bri, linetype = "dashed") +
  annotate("text", x = Topt_B - 0.5, y = 0.2, label = "optimum", angle = 90, hjust = 0, colour = col_bri, size = 3.4) +
  labs(x = "Temperature (C)", y = "Performance", title = "An asymmetric thermal performance curve")

Figure 1: A thermal performance curve fitted with the Briere model. The optimum (dashed) sits close to the critical thermal maximum; the shaded band is the thermal breadth where performance exceeds eighty per cent of the peak.

Why the symmetry matters

It is tempting to reach for a Gaussian, the familiar bell, because it is symmetric and easy. Fit one to the same data and it will happily find a peak. The trouble is on the way down.

fitG <- nls(y ~ A * exp(-((Td - mo)^2) / (2 * s^2)),
            start = list(A = 2, mo = 34, s = 8))
round(coef(fitG), 3)

     A     mo      s 
 1.961 33.041 10.205 

The Gaussian places the optimum at 33.0 degrees, a little low, because a symmetric curve cannot lean toward the upper limit. Worse is what it implies for the thermal maximum. A Gaussian has no hard upper bound, but projected out to where performance would vanish it puts the limit near 53.5 degrees, about 11 degrees above the real CTmax of 42. The symmetric curve smears the steep descending limb across a temperature range the organism could never survive.

TgG <- seq(10, 58, 0.1); PgG <- cg["A"] * exp(-((TgG - cg["mo"])^2) / (2 * cg["s"]^2))
ggplot(dd, aes(Td, y)) +
  geom_point(colour = col_dat, size = 2, alpha = 0.7) +
  geom_line(data = cur, aes(T, P), colour = col_bri, linewidth = 1) +
  geom_line(data = data.frame(T = TgG, P = PgG), aes(T, P), colour = col_gau, linewidth = 1, linetype = "22") +
  geom_vline(xintercept = Tm, colour = col_ref, linetype = "dotted") +
  annotate("text", x = Tm + 0.4, y = 1.6, label = "true CTmax", hjust = 0, colour = col_ref, size = 3.4) +
  labs(x = "Temperature (C)", y = "Performance", title = "Symmetric curves overshoot the thermal limit")

Figure 2: Briere (green) and Gaussian (ochre) fits to the same data. They agree on the rising limb and disagree sharply on the fall: the symmetric curve pushes the apparent upper limit far too high.

What to trust, and what to check

The optimum and the breadth are interpolated from the bulk of the data and are reasonably safe. CTmax is different: it is an extrapolation to the edge of the curve, where there are few measurements and where the choice of model matters most. Treat it as the least certain number on the curve, and see the checking post for how far it can move. The two follow-ups take the curve somewhere useful: what it says about performance in a fluctuating climate, and which species a warming world actually threatens.

Related tutorials

• Jensen’s inequality and thermal variability
• Thermal safety margins and warming
• Checking a thermal performance analysis
• Nonlinear regression with nls

References

• Briere JF et al. 1999. Environmental Entomology 28(1):22-29 (10.1093/ee/28.1.22).
• Huey RB, Stevenson RD 1979. American Zoologist 19(1):357-366 (10.1093/icb/19.1.357).
• Kingsolver JG 2009. The American Naturalist 174(6):755-768 (10.1086/648310).
• Sinclair BJ et al. 2016. Ecology Letters 19(11):1372-1385 (10.1111/ele.12686).