Exserohilum turcicum (teleomorph Setosphaeria turcica) is an important pathogen of maize (Zea mays) in Africa. Like many other foliar plant pathogens, spores of S. turcica require moist, warm conditions to germinate and infect the host. Despite many decades of research on this pathogen, empirical data relating environmental conditions (temperature, moisture) to germination and infection rates, required to model the infection process, have rarely been published. Fortunately, rates of appressorium formation in relation to temperature and leaf wetness duration (dew period) were estimated under controlled conditions in a single study in Israel (1). These data can be used to parameterize an infection model based on the concept of variable failure time survival models.
In the following discussion, we introduce the concept of survival models for spore germination (and infection), the hazard function which is the instantaneous risk of transition from one state to another (for example, ungerminated to germinated spore), and the Weibull distribution which describes how the hazard could change over time. We then discuss how temperature affects the rate of biological processes, and how this rate can be interpreted as the effective amount of time passed in an accelerated failure time model. Finally, we parameterize a temperature-dependent accelerated failure time model using experimental data.
We are interested in the process of spore germination and subsequent infection of leaves. A population of spores on a leaf don’t all instantaneously germinate and infect once conditions are right (e.g. the leaf becomes wet or humidity hits a threshold). Empirical studies show that the fraction of spores that have germinated increases over time. In other words, the time of transition (i.e. the lifespan) \(L\) has a probability distribution known as the survival function \(S(t)\): \[ S(t) = \Pr(L>t) \] which is the probability at time \(t\) that a spore has not already germinated. So, \(S(0)=1\), and \(S\) declines as \(t\) increases. When considering a large population, \(S(t)\) can be thought of as the fraction of the population that has not germinated by \(t\). Analysis of how \(S\) declines with \(t\) is known as survival analysis.
In the case of leaf infection, we are more interested in fraction of spores that have germinated, rather than those which haven’t. The complement of \(S(t)\) is known as the lifetime distribution function \(F(t)\):
\[ F(t) = 1 - S(t) \]
We can differentiate \(F(t)\) to give the rate of transitions per unit time, also known as the event density \(f(t)\):
\[ f(t) = \frac{d}{dt}F(t) \]
In other words, a large \(f(t)\) means that many spores are germinating at time \(t\), while a large \(F(t)\) means that a lot of spores have already germinated by \(t\).
Survival processes are usually modelled with reference to the hazard function \(h(t)\), which is the instantaneous risk of transition (i.e. spore germination) at time \(t\), given that transition has not yet occurred:
\[ h(t) = \lim\limits_{\Delta t \rightarrow 0} \frac{P(t \leq L < t + \Delta t | L > t)}{\Delta t} = \frac{f(t)}{S(t)} \]
In other words, the hazard is the transition rate divided by the fraction still left to transition.
Integration of the hazard function over time gives us the cumulative hazard \(H(t)\), which is the accumulated risk of transition up to time \(t\):
\[ H(t) = \int_0^t{h(u)du} \]
As risk accumulates over time, the fraction of spores left to transition (\(S(t)\)) decreases and the fraction which have transitioned (\(F(t)\)) increases:
\[ S(t) = 1-F(t) = \exp(-H(t)) \]
Survival processes in nature can be modelled using a variety of different hazard functions. For example, if the hazard is constant (\(\theta\)) we get an exponential survival process:
\[ H(t) = \theta t \]
\[ F(t) = 1-\exp(-\theta t) \]
Here \(\theta > 0\) is the rate parameter, representing the constant risk of transition per unit time. An equivalent parameterisation uses the scale parameter \(\lambda = 1/\theta\), giving:
\[ F(t) = 1-\exp \left( - \frac{t}{\lambda} \right) \]
As \(t\) increases, the exponent declines toward zero and the fraction which have germinated tends towards 1.
The shapes of the \(F(t)\) curves produced in empirical studies of spore germination indicate that the exponential model, i.e. constant hazard, is not usually appropriate and that \(h(t)\) varies over time (2,3). For example, \(F(t)\) often has a lag phase, with germination rate (\(f(t)\)) increasing then decreasing, resulting in a sigmoid \(F(t)\) curve. This suggests that \(h(t)\) increases with \(t\) in some way. Various parametric equations can be used to flexibly model variation in \(h(t)\), for example the Gamma, Gompertz and Weibull models.
The Weibull model fits fungal spore germination data well (2), while maintaining simplicity with only two parameters, the shape \(\beta\) and the scale \(\lambda\):
\[ h(t) = \frac{\beta}{\lambda}\left(\frac{t}{\lambda}\right)^{\beta-1} \]
Both \(\beta\) and \(\lambda\) must be greater than zero. \(\beta\) affects the overall shape of \(h(t)\). When \(\beta = 1\) the Weibull reduces to the exponential case and \(h(t) = 1/\lambda\); when \(\beta < 1\), \(h(t)\) decreases with \(t\); and when \(\beta > 1\), \(h(t)\) increases with \(t\). The scale parameter \(\lambda\) stretches or compresses the time axis. Larger \(\lambda\) values correspond to slower germination, because more time must pass before the hazard reaches a given level.
The integral of \(h(t)\), the cumulative hazard \(H(t)\), has a conveniently simple form:
\[ H(t) = \left( \frac{t}{\lambda} \right)^\beta \]
The fraction (or probability) of spores transitioned is then:
\[ F(t) = 1 - S(t) = 1 - \exp \left[- \left( \frac{t}{\lambda} \right)^\beta \right] \]
Temperature influences the rate of nearly all biological processes (4). Rates are highest at an optimal temperature \(T_\text{opt}\), declining to (near) zero below a minimum \(T_\text{min}\) or above a maximum \(T_\text{max}\). These temperatures are usually termed the cardinal temperatures of a process. The cardinal temperatures of various processes relating to plant pathogenesis have been catalogued for hundreds of fungi and oomycetes (5).
Many different functions have been used to describe the thermal performance curve or temperature response curve that relates temperature to rate (6). The beta function has been used to model plant growth rates in relation to temperature \(T\) (7), and has been applied to spore germination rates (2):
\[ r(T)=\left( \frac{T_\text{max}-T}{T_\text{max}-T_{\text{opt}}} \right) \left( \frac{T-T_\text{min}}{T_{\text{opt}}-T_\text{min}} \right)^ \frac{T\text{opt}-T\text{min}}{T\text{max}-T\text{opt}} \]
for \(T_{\min} < T < T_{\max}\), 0 otherwise. In other words, the rate is 0 when conditions are colder than \(T_{\min}\) or warmer than \(T_{\max}\), and 1 when \(T = T_{\text{opt}}\).
There are several ways in which the survival process could incorporate the temperature-dependent rate \(r(T)\). One way would be to directly modify \(h(t)\), so that \(h(t,T) = h(t)r(T)\). The Weibull cumulative hazard is then:
\[ H(t,T) = r(T) \left( \frac{t}{\lambda} \right)^\beta \]
and \(F(t, T)\) would be calculated as before.
This would be relatively straightforward to model if temperature stayed constant over time, but in the real world temperature varies continuously. Because \(T\) varies arbitrarily with \(t\), \(H(t, T)\) must be calculated piecewise by numerically integrating over short time periods. Mean hourly temperatures are often used in physiological models, balancing the need to capture diurnal temperature variation with computational burden (2). However, this is a somewhat cumbersome approach, so instead let us consider accelerated failure time models which accelerate (or, in this case, decelerate) a process according to some external driver.
In accelerated failure time (AFT) models, external factors such as temperature alter the rate at which effective time accumulates. The underlying distribution of lifetimes (\(L\)) is assumed to be fixed, but the mapping from real time to effective time is rescaled by the environment. In biological terms, this is equivalent to the physiological clock ticking faster or slower depending on conditions. At the optimal temperature, effective time matches clock time; under suboptimal temperatures, effective time accrues more slowly, so it takes longer in real time for the same fraction of spores to germinate. This idea of temporal scaling is supported by empirical studies of organismal lifespans (8).
If temperature was constant, the amount of effective time \(\tau\) elapsed over real time duration \(t\) would be \(\tau = r(T)t\). Suboptimal temperatures would result effectively in less physiological time passing. As temperature varies in reality, again work piecewise with mean temperatures within a short enough time duration to capture important variation, normally hourly. For a period \(n\) time steps, with mean temperature \(T_i\) in the \(i\)th time step of length \(t_i\):
\[ \tau(t, T) = \sum_{i=1}^{n}r(T_i)t_i \]
When working with hourly temperature data, \(t_i=1\) and \(\tau\) is the sum of the hourly rates.
The fraction of spores transitioned by clock time \(t\) is then:
\[ F(t,T) = 1 - \exp \left[- \left( \frac{\tau(t,T)}{\lambda} \right)^\beta \right] \]
In experiments tracking spore germination or infection, it is possible that only a fraction \(s\) of spores are viable. These viable spores germinate according to the Weibull distribution, while the remaining fraction \((1-s)\) are non-viable (dead) and will never germinate. This type of model is known as a cure-fraction model, a term originating in the medical survival-analysis literature.
The presence of non-viable spores can be inferred if the fraction germinated plateaus below 1 after a sufficiently long period of observation, i.e.
\[ \lim_{t \to \infty} F(t) < 1. \]
If no plateau is observed, it can be difficult to distinguish whether non-viable spores are present or whether the germination process is simply very slow (corresponding to a large scale parameter \(\lambda\)).
Experiments that track spore germination over time at constant temperature allow us to parameterize the model (2,3). We’ll analyse data on E. turcicum spore germination and appressorium formation (1) in the R environment. Specifically, data abstracted (using the online tool WebPlotDigitizer (9)) from Figure 4 showing the percentage of spores (100 conidia were counted per temperature and time point) which have formed appressoria at various constant temperatures and dew periods (Figure 1).
The experiment was conducted at six temperatures (10–30°C at 5°C intervals) and appressorium formation noted hourly from 4–8h of wetness. As there was zero appressorium formation at the beginning of the experiment, we can include zero appressoria for all temperatures at \(t=0\).
# Load tidyverse libraries for data handling and graphics
library(tidyverse)
# Set theme for plotting
theme_set(theme_bw())
# Read data
dat <- read_csv("Levy1983.csv")
# Inspect data
glimpse(dat)Rows: 36
Columns: 3
$ lwd <dbl> 8, 8, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5,…
$ temp <dbl> 10, 15, 20, 25, 30, 35, 10, 15, 20, 25, 30, 35, 10, 15, 20, 25, 3…
$ germ <dbl> 2, 20, 70, 40, 21, 1, 2, 11, 50, 30, 16, 1, 2, 6, 25, 18, 9, 1, 2…There are 30 observations of three variables: lwd is the leaf wetness duration, which can be interpreted as clock time \(t\); temp is the temperature; and germ is the number of spores (out of 100) which have germinated and formed appressoria. Re-plotting the data by lwd, we see no evidence of a plateau (Figure 2). Therefore, we won’t be able to determine whether or not all spores were viable, and will not try to fit a cure fraction model.
ggplot(dat,
aes(lwd, germ, col = factor(temp))) +
geom_point() +
geom_line() +
scale_color_viridis_d(option = "B", begin = 0, end = 0.9) +
labs(col = "temp") +
theme(legend.position = c(0.2, 0.5))
The beta function describes the relationship between effective (or physiological) time and real (or clock) time. We will need to optimize the parameters (cardinal temperatures) of the model to fit the data. The beta function in R therefore includes some safeguards to help avoid unwanted behaviour during optimization:
# BETA FUNCTION
# Canonical beta temperature response (Yan & Hunt 1993)
# with safeguards for optimization
beta_safe <- function(x, tmin, topt, tmax){
# invalid cardinals → zero vector (keeps objective finite)
if (topt <= tmin || topt >= tmax) return(rep(0.0, length(x)))
a <- topt - tmin
b <- tmax - topt
k <- a / b
# Normalised coordinates on each side (in [0,1] inside the interval)
z1 <- (x - tmin) / a # 0 at tmin, 1 at topt
z2 <- (tmax - x) / b # 1 at topt, 0 at tmax
# Optional: log form for numeric stability in extreme k
r <- exp(log(pmax(z2, 0)) + k * log(pmax(z1, 0))) # r = z2 * z1^k
# Zero outside [tmin, tmax], clip to [0,1], set exact peak
r[(x < tmin) | (x > tmax)] <- 0
r <- pmin(pmax(r, 0), 1)
r[abs(x - topt) <= 1e-12 * (1 + abs(topt))] <- 1
r[!is.finite(r)] <- 0
r
}Based on the original figure, we can visualize the rate \(r(T)\) for a set of reasonable starting parameters: \(T_\text{min}=10\), \(T_\text{topt}=20\) and \(T_\text{max}=35\) (Figure 3).
tibble(
x = seq(0,40,by = 0.1),
r = beta_safe(x, 10, 20,35)) %>%
ggplot(aes(x,r)) +
geom_line() +
labs(x = "temperature", y = "r(T)") +
theme_classic()
Optimization is a numerical method for finding the best parameters for a model, by searching parameter space for solutions that match the data most accurately. Parameter space essentially the range of different values all the variables that we are interested in could take. There are various ways to measure “best”, such as least-squares differences between observed and fitted values, or maximum-likelihood, which identifies the parameter set that makes the observed data most probable.
Successful optimisation depends on being able to move smoothly through parameter space, but problems can arise if there are constraints. For example, in this case we must have \(T_{\text{min}} < T_{\text{opt}} < T_{\text{max}}\). To ensure the cardinals remain ordered, we re-parameterise the problem: instead of optimising directly over \((T_{\text{min}}, T_{\text{opt}}, T_{\text{max}})\), we optimise over \(\theta = (T_{\text{opt}}, \log d_1, \log d_2)\), where \(d_1 = T_{\text{opt}} - T_{\text{min}} > 0\) and \(d_2 = T_{\text{max}} - T_{\text{opt}} > 0\). We then reconstruct \(T_{\text{min}} = T_{\text{opt}} - d_1\) and \(T_{\text{max}} = T_{\text{opt}} + d_2\). Optionally, we can enforce a minimum gap (e.g. \(1^{\circ}\)C) by adding 1 to each distance.
# CARDINAL TEMPERATURE ORDERING
# Ordered re-parameterisation: guarantees tmin < topt < tmax (and optional min gap)
order_cardinals <- function(theta, gap_min = 1){
topt <- theta[[1]]
d1 <- exp(theta[[2]]) + gap_min
d2 <- exp(theta[[3]]) + gap_min
setNames(c(topt - d1, topt, topt + d2), c("Tmin", "Topt", "Tmax"))
}Each observation in the experiment consists of the number of germinated conidia out of 100 counted. This outcome is naturally modelled as binomially distributed, with success probability \(p(t,T)\) given by the Weibull–beta model. For each observation of \(y\) conidia germinated out of \(n=100\), the binomial likelihood is:
\[ L(y \mid n, p) = \binom{n}{y} p^{y} (1-p)^{n-y}. \]
where the binomial coefficient is:
\[ \binom{n}{y}=\frac{n!}{y!(n-y)!} \]
If we have \(m\) independent observations, the joint likelihood is the product:
\[ L(\theta) = \prod_{i=1}^m \binom{n_i}{y_i} p_i^{y_i} (1-p_i)^{n_i - y_i}, \]
where \(p_i = p(t_i, T_i \mid \theta)\) depends on the parameters \(\theta=(\beta,\lambda,T_{\min},T_{\text{opt}},T_{\max})\).
For optimization we’ll work with the negative log-likelihood. The optimisation algorithm will look for a minimum (the best solution will have the highest likelihood, and thus the lowest negative log likelihood). The log likelihood is:
\[ \ell(\theta) = \log L(\theta) = \sum_{i=1}^m \left[ \log\binom{n_i}{y_i} + y_i \log(p_i) + (n_i-y_i)\log(1-p_i) \right]. \]
The binomial coefficient \(\binom{n_i}{y_i}\) varies with the data but does not depend on the model parameters, so it can be dropped without affecting optimisation. The negative log-likelihood is then:
\[ \text{NLL}(\theta) = -\ell(\theta) = -\sum_{i=1}^m \big[ y_i \log(p_i) + (n_i-y_i)\log(1-p_i) \big]. \]
Optimization algorithms are more likely to converge on optimal parameter values when given starting values close to those optima. We can either estimate some starting values by eye from the plots, or use some simple calculations. For the cardinal temperatures, we’ll estimate \(T_\text{opt}^0\) as the temperature of fastest germination (20°C), \(T_\text{min}^0\) as the lowest temperature in the experiment with near zero germination (10°C) and \(T_\text{max}^0\) as the highest temperature with near zero germination (35°C).
For the shape \(\beta\) and scale \(\lambda\) of the Weibull distribution, we can use a mathematical property of the distribution, namely that the complementary log-log (cloglog) transformation results in a linear regression with slope \(\beta\) and intercept \(-\beta \log(\lambda)\). Recall that the fraction germinated is:
\[ F(t) = 1 - \exp \left[- \left( \frac{t}{\lambda} \right)^\beta \right] \] If we rearrange the expression and take logs, then:
\[ \log\!\big(-\log(1 - F(t))\big) = \beta \log(t) - \beta \log(\lambda) \]
If we plot \(y=\log\!\big(-\log(1 - F(t))\big)\) against \(x=\log(t)\) then starting parameters \(\beta^0\) and \(\lambda^0\) can be obtained.
# GET STARTING VALUES FOR OPTIMISATION
# dat: data frame with columns temp (deg C), lwd (h), germ (%)
get_start <- function(dat, Tmin0 = 10, Topt0 = 20, Tmax0 = 35,
min_gap = 1, n_per_point = 100) {
# 1) Pick the exact or nearest available temperature to Topt0
T_levels <- sort(unique(dat$temp))
if (length(T_levels) == 0) stop("No temperatures found in 'dat$temp'.")
Topt_use <- T_levels[which.min(abs(T_levels - Topt0))]
# 2) Estimate beta0, lambda0 at Topt_use via cloglog
dat_topt <- dat |>
dplyr::filter(temp == Topt_use, lwd > 0) |>
dplyr::transmute(
t = lwd,
# continuity-corrected proportion; assumes n = n_per_point (100 in Levy)
pcc = (round(germ) + 0.5) / (n_per_point + 1),
x = log(t),
y = log(-log(1 - pcc))
) |>
dplyr::filter(is.finite(x), is.finite(y), pcc > 0 & pcc < 1)
# Fit the cloglog model
fit_lm <- lm(y ~ x, data = dat_topt)
beta0 <- max(1.05, as.numeric(coef(fit_lm)[2])) # slope ≈ beta; floor > 1
a0 <- as.numeric(coef(fit_lm)[1]) # intercept ≈ -beta*log(lambda)
lambda0 <- exp(-a0 / beta0)
# Guard against non-finite or non-positive lambda
if (!is.finite(lambda0) || lambda0 <= 0) lambda0 <- 8
# 3) Map cardinals to optimiser parameter vector (ordered re-param with min_gap)
d1_0 <- max(Topt0 - Tmin0, min_gap) # Topt - Tmin
d2_0 <- max(Tmax0 - Topt0, min_gap) # Tmax - Topt
# Return named vector on the transformed scale expected by your optimiser
c(
log_beta = log(beta0),
log_lambda = log(lambda0),
Topt = Topt0,
log_d1 = log(d1_0 - min_gap), # because order_cardinals() adds +min_gap back
log_d2 = log(d2_0 - min_gap)
)
}
# Run the function
par0 <- get_start(dat)
par0 log_beta log_lambda Topt log_d1 log_d2
1.425890 2.047232 20.000000 2.197225 2.639057 We can now plot the Weibull model with these starting parameters (Figure 4).
# Back-transform starting parameters
beta0 <- exp(par0["log_beta"])
lambda0 <- exp(par0["log_lambda"])
cards <- order_cardinals(par0[c("Topt","log_d1","log_d2")])
# Generate predictions over 0–8 h
expand.grid(
t = seq(0, 8, by = 0.1),
temp = seq(10, 35, by = 5)
) %>%
# Weibull–beta cumulative germination
mutate(
rT = beta_safe(temp, cards["Tmin"], cards["Topt"], cards["Tmax"]),
F = 1 - exp(-((rT * t) / lambda0)^beta0),
germ_pred = 100 * F # scale to %
) %>%
# Plot predicted
ggplot(aes(t, germ_pred, col = factor(temp))) +
geom_line() +
scale_color_viridis_d(option = "B", begin = 0, end = 0.9) +
labs(x = "lwd", y = "germ", col = "temp") +
theme(legend.position = c(0.2, 0.5))
We can now optimize the parameters using our starting estimates. The optimization will minimize the negative log-likelihood using the L-BFGS-B algorithm. There are five steps: 1) A function to calculate the negative log-likelihood from the parameter set, temperature, time, the observed total number and germinated number of spores; 2) A get starting estimates, if not done already; 3) Toptimize our parameters using maximum likelihood; 4) back-transform our logged parameter estimates:
# 1) Binomial negative log-likelihood (compact)
# par = c(log_beta, log_lambda, Topt, log_d1, log_d2)
# y = germinated, n = total number
nll_binom <- function(par, temp, time, y, n, gap_min = 1){
beta <- exp(par[1])
lambda <- exp(par[2])
cards <- order_cardinals(par[3:5], gap_min = gap_min)
rT <- beta_safe(temp, cards["Tmin"], cards["Topt"], cards["Tmax"])
eta <- ((rT * time) / lambda)^beta
p <- 1 - exp(-eta)
p <- pmin(pmax(p, 1e-12), 1 - 1e-12)
-sum(dbinom(y, size = n, prob = p, log = TRUE))
}
# 2) Starting estimates
par0 <- get_start(dat, Tmin0 = 10, Topt0 = 20, Tmax0 = 35, min_gap = 1)
# 3) Optimization
fit <- optim(
par = par0,
fn = nll_binom,
temp = dat$temp,
time = dat$lwd,
y = round(dat$germ), # convert % to counts
n = rep(100, nrow(dat)), # 100 per point
method = "L-BFGS-B",
control = list(maxit = 2000)
)
stopifnot(fit$convergence == 0)
# 4) Back-transform
par_hat <- fit$par
beta_hat <- exp(par_hat[1])
lam_hat <- exp(par_hat[2])
cards <- order_cardinals(par_hat[3:5], gap_min = 1)
par_opt <- tibble(
parameter = c("beta (shape)", "lambda (scale)", "Tmin", "Topt", "Tmax"),
estimate = c(beta_hat, lam_hat, cards)
)
par_opt# A tibble: 5 × 2
parameter estimate
<chr> <dbl>
1 beta (shape) 3.46
2 lambda (scale) 8.64
3 Tmin 6.22
4 Topt 22.3
5 Tmax 37.6 The optimization suggests that the cardinal temperature range is slightly wider than the conditions under which the experiments were conducted. We can now visualize how well the optimized Weibull model fits the observations. First we need a function to predict the number of germinated spores, then we can apply this to our experimental time and temperature combinations.
# Predict proportion germinated p(t, T | par) in [0,1]
# par = c(log_beta, log_lambda, Topt, log_d1, log_d2)
predict_weibull_beta <- function(par, temp, time, gap_min = 1){
beta <- exp(par[1])
lambda <- exp(par[2])
cards <- order_cardinals(par[3:5], gap_min = gap_min) # Tmin, Topt, Tmax
rT <- beta_safe(temp, cards["Tmin"], cards["Topt"], cards["Tmax"])
eta <- ((rT * time) / lambda)^beta
p <- 1 - exp(-eta)
pmin(pmax(p, 0), 1) # clamp
}
# Add fitted proportion and percent to the original data
dat_pred <- dat %>%
mutate(
p_hat = predict_weibull_beta(par_hat, temp = temp, time = lwd, gap_min = 1),
fit_pct = 100 * p_hat
)
# Plot observed vs. fitted values
ggplot(data = dat_pred,
aes(fit_pct, germ)) +
geom_abline(slope = 1, intercept = 1) +
geom_point() +
labs(x = "predicted", y = "observed")
The plot of observed vs. predicted germination events suggests that the model is slightly over-predicting germinations at lower values (Figure 5). Nevertheless, we now have an optimized set of parameters with which to model E. turcicum appressorium formation over time at different temperatures.