Seed germination: fitting hydro-time models with R

Andrea Onofri ·  Added on March 23, 2020 ·  16 min read


THE CODE IN THIS POST WAS UPDATED ON JANUARY 2022

I am locked at home, due to the COVID-19 emergency in Italy. Luckily I am healthy, but there is not much to do, inside. I thought it might be nice to spend some time to talk about seed germination models and the connections with survival analysis.

We all know that seeds need water to germinate. Indeed, the absorption of water activates the hydrolytic enzymes, which break down food resources stored in seeds and provide energy for germination. As the consequence, there is a very close relationship between water content in the substrate and germination velocity: the higher the water content the quickest the germination, as long as the availability of oxygen does not become a problem (well, water and oxygen in soil may compete for space and a high water content may result in oxygen shortage).

Indeed, it is relevant to build germination models, linking the proportion of germinated seeds to water availability in the substrate; these models are usually known as hydro-time (HT) models. The starting point is the famous equation of Bradford (1992), where the germination rate (GR) for the \(i-th\) seed in the lot is expressed as a linear function of water potential in the substrate (\(\Psi\)):

\[ GR_i = \textrm{max} \left( \frac{\Psi - \Psi_{b(i)}}{\theta_H}; 0 \right) \quad \quad \quad \quad (1)\]

In that equation, \(\Psi_{b(i)}\) is the base water potential for the \(i-th\) seed and \(\theta_H\) is the hydro-time constant, expressed as MPa day or MPa hour. The concept is relatively simple: we just need to remember that the water can only move from a position with a higher water potential to a position with a lower water potential. Therefore, a seed cannot germinate when its base water potential is higher than the water potential in the substrate.

When \(\Psi > \Psi_b(i)\), the germination rate of the \(i-th\) seed is linearly related to \(\Psi\): the higher this latter value, the higher the germination rate. Now we should consider that the germination rate is the inverse of the germination time (\(GR = 1/t\)); thus, the higher the GR, the shortest the germination time. Germination is achieved at the time \(t\) when:

\[ t \, \left( \Psi - \Psi_{b(i)} \right) = \theta_H \quad \quad \quad (2)\]

Elsewhere in this website, I show that Equation 1 can be fitted to germination data in a two-steps fashion. In this page we will see how we can embed Equation 1 into a time-to-event model, to predict the proportion of germinated seeds, depending on time and water content in the substrate. As usual, let’s start from a practical example.


The dataset

The germination of rapeseed (Brassica napus L. var. oleifera, cv. Excalibur) was tested at fourteen different water potentials (0, -0.03, -0.15, -0.3, -0.4, -0.5, -0.6, -0.7, -0.8, -0.9, -1, -1.1, -1.2, -1.5 MPa), which were created by using a polyethylene glycol solution (PEG 6000). For each water potential level, three replicated Petri dishes with 50 seeds were incubated at 20°C. Germinated seeds were counted and removed every 2-3 days for 14 days.

The dataset was published by Pace et al. (2012) and it is available as rape in the drcSeedGerm package, which needs to be installed from github (see below). Furthermore, the package drcte is necessary to fit time-to-event models and it should also be installed from gitHub. Please, make sure you have the most updated version for both packages.

The following code loads the necessary packages, loads the dataset rape and shows the first six lines.

# library(devtools)
# install_github("OnofriAndreaPG/drcSeedGerm")
# install_github("OnofriAndreaPG/drcte")
library(drcSeedGerm)
library(drcte)
data(rape)
head(rape)
##   Psi Dish timeBef timeAf nSeeds nCum propCum
## 1   0    1       0      3     49   49    0.98
## 2   0    1       3      4      0   49    0.98
## 3   0    1       4      5      0   49    0.98
## 4   0    1       5      7      0   49    0.98
## 5   0    1       7     10      0   49    0.98
## 6   0    1      10     14      0   49    0.98

We can see that the data are grouped by assessment interval: ‘timeAf’ represents the moment when germinated seeds were counted, while ’timeBef’ represents the previous inspection time (or the beginning of the assay). The column ’nSeeds’ is the number of seeds that germinated during the time interval between ‘timeBef’ and ‘timeAf. The ’propCum’ column contains the cumulative proportions of germinated seeds and it is not necessary for time-to-event models.


Building hydro-time models

Models based on the distribution of germination time

How can we rework Equation 1 to predict the proportion of germinated seeds, as a function of time and water potential? One line of attack follows the proposal we made in a relatively recent paper (Onofri at al., 2018). We start from the idea that the time course of the proportion of germinated seeds (\(P\)) is expected to increase over time, according to a S-shaped curve, such as the usual log-logistic cumulative probability function (other cumulative distribution functions can be used; see our original paper):

\[P(t) = \frac{ d }{1 + \exp \left\{ b \left[ \log(t) - \log( e ) \right] \right\} } \quad \quad \quad (3)\]

where \(e\) is the median germination time, \(b\) is the slope at the inflection point and \(d\) is the maximum germinated proportion. Considering that the germination rate is the inverse of germination time, we can write:

\[P(t) = \frac{ d }{1 + \exp \left\{ b \left[ \log(t) - \log(1 / GR_{50} ) \right] \right\} } \quad \quad \quad (4)\]

where \(GR_{50}\) is the median germination rate in the population. We can now express \(GR_{50}\), \(b\) and \(d\) as linear/nonlinear functions of \(\Psi\) (temperature and other environmental variables can be included as well. See our original paper):

\[P(t, \Psi) = \frac{ h_1(\Psi) }{1 + \exp \left\{ b \left[ \log(t) - \log(1 / \left[ GR_{50}(\Psi) \right] ) \right] \right\} } \quad \quad \quad (5)\]

In our paper, for \(GR_{50}\), we slightly reparameterised Equation 1 above:

\[ GR_{50}(\Psi) = \textrm{max} \left( \frac{\Psi - \Psi_{b(50)}}{\theta_H}; 0 \right) \quad \quad \quad \quad (6)\]

where \(\Psi_{b(50)}\) is the median base water potential in the population.

For \(d\), we used a shifted exponential distribution, which implies that germination capability is fully determined by the distribution of base water potential within the population and no germination occurs at \(\Psi \leq \Psi_b\):

\[d = h_1(\Psi ) = \textrm{max} \left\{ G \, \left[ 1 - \exp \left( \frac{ \Psi - \Psi_{b(50)} }{\sigma_{\Psi_b}} \right) \right]; 0 \right\} \quad \quad \quad (7)\]

In the above equation, \(\sigma_{\Psi_b}\) represents the variability of \(\Psi_b\) within the population, which determines the steepness of the increase in \(d\) as \(\Psi\) increases. \(G\) is the germinable fraction, accounting for the fact that \(d\) may not reach 1, regardless of time and water potential.

The parameter \(b\) was assumed to be constant and independent on \(\Psi\). In the end, our hydro-time model is composed by four sub-models:

  1. a cumulative probability function (log-logistic, in our example), based on the three parameters \(d\), \(b\) and \(e = 1/GR_{50}\);
  2. a sub-model expressing \(d\) as a function of \(\Psi\);
  3. a sub-model expressing \(GR_{50}\) as a function of \(\Psi\);
  4. a sub-model expressing \(b\) as a function of \(\Psi\), although, this was indeed a simple identity model \(b(\Psi) = b\).

This hydro-time model was implemented in R as the HTE1() function, and it is available within the drcSeedGerm package, together with the appropriate self-starting routine. It can be fitting by using the drmte() function in the drcte package.

modHTE <- drmte(nSeeds ~ timeBef + timeAf + Psi, 
                data = rape, fct = HTE1())
summary(modHTE)
## 
## Model fitted: Hydro-time model with shifted exponential for Pmax and linear model for GR50
## 
## Robust estimation: no 
## 
## Parameter estimates:
## 
##                         Estimate Std. Error  t-value   p-value    
## G:(Intercept)          0.9577918  0.0063667  150.438 < 2.2e-16 ***
## Psib:(Intercept)      -1.0397239  0.0047017 -221.138 < 2.2e-16 ***
## sigmaPsib:(Intercept)  0.1108891  0.0087598   12.659 < 2.2e-16 ***
## thetaH:(Intercept)     0.9061385  0.0301594   30.045 < 2.2e-16 ***
## b:(Intercept)          4.0273963  0.1960845   20.539 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As seeds are clustered in Petri dishes, in order not to violate the independence assumption, it is preferable to get cluster robust standard errors. One possibility is to use the grouped version of the sandwich estimator, as available in the ‘sandwich’ package (Berger, 2017), which is implemented in the ‘drcte’ package.

summary(modHTE, robust = T, units = Dish)
## 
## Model fitted: Hydro-time model with shifted exponential for Pmax and linear model for GR50
## 
## Robust estimation: Cluster robust sandwich SEs 
## 
## Parameter estimates:
## 
##                         Estimate Std. Error   t value  Pr(>|t|)    
## G:(Intercept)          0.9577918  0.0080923  118.3591 < 2.2e-16 ***
## Psib:(Intercept)      -1.0397239  0.0047063 -220.9206 < 2.2e-16 ***
## sigmaPsib:(Intercept)  0.1108891  0.0121879    9.0983 < 2.2e-16 ***
## thetaH:(Intercept)     0.9061385  0.0410425   22.0781 < 2.2e-16 ***
## b:(Intercept)          4.0273963  0.1934513   20.8187 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Once the model is fitted, we may be interested in using the resulting curve to retrieve some biologically relevant information. For example, it may be interesting to retrieve the germination rates for some selected percentiles (e.g., the 30th, 20th and 10th percentiles). This is possible using the quantile() function in ‘drcte’ package, that is a wrapper for the original ED() function in the package drc. It works very much like the ED() function, although the environmental covariates, such as the selected \(\Psi\) level, can be specified as a numeric value.

#Naive standard errors
quantile(modHTE, probs = c(0.30, 0.20, 0.10), rate = T, Psi = -1)
## 
## Estimated quantiles
## 
##       Estimate        SE
## 1:30% 0.000000 0.0000000
## 1:20% 0.035792 0.0058117
## 1:10% 0.051304 0.0062243

In this example, we see that the \(GR_{30}\) cannot be calculated, as the germination capacity did not reach 30% at the selected water potential level (\(-1 \,\, MPa\)).

As we said, cluster robust standard errors are recommended, which can be requested by setting the ‘robust = TRUE’ argument and passing a vector for the grouping factor. If necessary, germination times can be obtained in a similar way, by removing the ‘rate’ argument

#Cluster robust standard errors
quantile(modHTE, probs = c(0.30, 0.20, 0.10), rate = T, Psi = -1, robust = T,
         units = Dish)
## 
## Estimated quantiles
## 
##       Estimate        SE
## 1:30% 0.000000 0.0000000
## 1:20% 0.035792 0.0054538
## 1:10% 0.051304 0.0058715
 r
#Germination times
quantile(modHTE, probs = c(0.30, 0.20, 0.10), Psi = -1, robust = T,
         units = Dish)
## 
## Estimated quantiles
## 
##       Estimate     SE
## 1:30%      Inf    NaN
## 1:20%   27.939 4.2572
## 1:10%   19.492 2.2308

Last, but not least, we can predict the proportion of germinated seeds at given time and water potential level.

predict(modHTE, se.fit=T, robust = T,
         units = Dish,
        newdata = data.frame(Time=c(10, 10, 10), 
                             Psi=c(-1.5, -0.75, 0))
        )
##   Time   Psi Prediction          SE
## 1   10 -1.50  0.0000000 0.000000000
## 2   10 -0.75  0.8794025 0.017562676
## 3   10  0.00  0.9576590 0.008070333

Models based on the distribution of \(\Psi_b\)

Another type of hydro-time model was proposed by Bradford (2002) and later extended by Mesgaran et al (2013). This approach starts always from Equation 1; from that equation, considering that the germination time is the inverse of the GR, we can easily get to the following equation:

\[\Psi_b = \Psi - \frac{\theta_H}{t} \quad \quad \quad (8)\]

where \(t\) is the germination time. What does this equation tell us? Let’s assume that the hydro-time to germination is 10 \(MPa \, d\) and the environmental water potential is -1 \(MPa\). A single seed germinates in exactly one day, if its base water potential is \(-1 - 10/1 = -11\). If the base water potential is higher, germination will take more than one day; if it is lower, germination will take less than one day. But now, the following questions come: how many seeds in a population will be able to germinate in one day? And in two days? And in \(t\) days?

We know that the seeds within a population do not germinate altogether in the same moment, as they have different individual values of base water potential. If the population is big enough, we can describe the variation of \(\Psi_b\) within the population by using some density function, possibly parameterised by way of a location (\(\mu\)) and a scale (\(\sigma\)) parameter:

\[ \Psi_b \sim \phi \left( \frac{\Psi_b - \mu}{\sigma} \right) \quad \quad \quad (9)\]

This is easier to understand if we make a specific example. Let’s assume that the distribution of \(\Psi_b\) values within the population is gaussian, with mean \(\mu = -9\) and standard deviation \(\sigma = 1\). Let’s also assume that the hydro-time parameter (\(\theta_H\)) is constant within the population. We have the situation depicted in the figure below.

The red left tail represents the proportion of seeds that germinate during the first day, as they have base water potentials equal to or lower than -11. By using the gaussian cumulative distribution function we can easily see that that proportion is 0.228:

pnorm(-1 - 10/1, mean = -9, sd = 1)
## [1] 0.02275013

More generally, we can write:

\[ P(t, \Psi) = \Phi \left\{ \frac{\Psi - (\theta_H / t) -\mu }{\sigma} \right\} \quad \quad \quad (10)\]

where \(\Phi\) is the selected cumulative distribution function. The above model returns the proportion of germinated seeds (G), as a function of time and water potential in the substrate. According to Bradford (2002), \(\Phi\) is cumulative gaussian.

Let’s think more deeply about Equation 9 (Bradford, 2002). This function was built to represent the cumulative distribution function of base water potential within the population. However, it can be as well taken to represent the cumulative distribution function of germination time within the population. Obviously, while the first distribution is gaussian, the second one is not: indeed, the germination time appears at the denominator of the expression \(\theta_H / t\). It doesn’t matter: every cumulative distribution function for germination time can be fit by using time-to-event methods!

We implemented this model in R as the function HTnorm() that is available within the drcSeedGerm package and it is meant to be used with the drm() function, in the drc package.

Mesgaran et al (2013) suggested that the distribution of base water potential within the population may not be gaussian and proposed several alternatives, which we have all implemented within the package. In all, drcSeedGerm contains six possible distributions:

  1. gaussian distribution (function HTnorm())
  2. logistic distribution (function HTL())
  3. Gumbel (function HTG())
  4. log-logistic (function HTLL())
  5. Weibull (Type I) (function HTW1())
  6. Weibull (Type II) (function HTW2())

The equations are given at the end of this page; for gaussian, logistic and log-logistic distributions, \(\Psi_{b(50)}\) is the median base water potential within the population. For the gaussian distribution, \(\sigma_{\Psi b}\) corresponds to the standard deviation of \(\Psi_b\) within the population.

Distributions based on logarithms (the log-logistic and all other distributions thereafter) are only defined for positive amounts. On the contrary, we know that base water potential is mostly negative. Therefore, shifted distributions need to be used, by introducing a shifting parameter \(\delta\) which ‘moves’ the distribution to the left, along the x-axis, so that negative values are possible (see Mesgaran et al., 2013).

Let’s fit the above functions to the ‘rape’ dataset. But, before, let me highlight that providing starting values is not necessary, as self-starting routines are already implemented for all models.

mod1 <- drmte(nSeeds ~ timeBef + timeAf + Psi, 
              data = rape, fct = HTnorm())
mod2 <- drmte(nSeeds ~ timeBef + timeAf + Psi,
              data = rape, fct = HTL())
mod3 <- drmte(nSeeds ~ timeBef + timeAf + Psi,
              data = rape, fct = HTG())
mod4 <- drmte(nSeeds ~ timeBef + timeAf + Psi,
            data = rape, fct = HTLL())
mod5 <- drmte(nSeeds ~ timeBef + timeAf + Psi,
            data = rape, fct = HTW1())
mod6 <- drmte(nSeeds ~ timeBef + timeAf + Psi,
            data = rape, fct = HTW2())

What is the best model for this dataset? Let’s use the Akaike’s Information Criterion (AIC) to decide:

AIC(mod1, mod2, mod3, mod4, mod5, mod6, modHTE)
##         df      AIC
## mod1   291 3516.914
## mod2   291 3300.824
## mod3   291 3097.775
## mod4   290 2886.608
## mod5   290 2889.306
## mod6   290 3009.023
## modHTE 289 2832.481

The first model modHTE considers explicitly the distribution of germination times and it is the best fitting of all. The other models consider explicitly the distribution of base water potential, while the distribution of germination times is indirectly included. Among these models, the gaussian is the worse fitting, while the log-logistic is the best one (mod4).

For this latter model, we take a look at the value of estimated parameters, with cluster robust standard errors.

summary(mod4, robust = T, units = Dish)
## 
## Model fitted: Hydrotime model with log-logistic distribution of Psib (Mesgaran et al., 2013)
## 
## Robust estimation: Cluster robust sandwich SEs 
## 
## Parameter estimates:
## 
##                     Estimate Std. Error  t value Pr(>|t|)    
## thetaH:(Intercept)  0.676804   0.074160   9.1263   <2e-16 ***
## delta:(Intercept)   1.136671   0.101804  11.1652   <2e-16 ***
## Psib50:(Intercept) -0.947798   0.020575 -46.0647   <2e-16 ***
## sigma:(Intercept)   0.371971   0.172825   2.1523   0.0322 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Germination rates and times for a certain percentile (e.g. GR50, GR30), can be obtained as shown above.

We can use the predict() function to plot the ‘modHTE’ and ‘mod4’ objects together in the same graph.

Thanks for reading!

Prof. Andrea Onofri
Department of Agricultural, Food and Environmental Sciences
University of Perugia (Italy)
Borgo XX Giugno 74
I-06121 - PERUGIA


References

  1. Berger, S., Graham, N., Zeileis, A., 2017. Various versatile variances: An object-oriented implementation of clustered covariances in R. Faculty of Economics and Statistics, University of Innsbruck, Innsbruck.
  2. Bradford, K.J., 2002. Applications of hydrothermal time to quantifying and modeling seed germination and dormancy. Weed Science 50, 248–260.
  3. Mesgaran, M.B., Mashhadi, H.R., Alizadeh, H., Hunt, J., Young, K.R., Cousens, R.D., 2013. Importance of distribution function selection for hydrothermal time models of seed germination. Weed Research 53, 89–101. https://doi.org/10.1111/wre.12008
  4. Onofri, A., Benincasa, P., Mesgaran, M.B., Ritz, C., 2018. Hydrothermal-time-to-event models for seed germination. European Journal of Agronomy 101, 129–139.
  5. Onofri, A., Mesgaran, M.B., Neve, P., Cousens, R.D., 2014. Experimental design and parameter estimation for threshold models in seed germination. Weed Research 54, 425–435. https://doi.org/10.1111/wre.12095
  6. Pace, R., Benincasa, P., Ghanem, M.E., Quinet, M., Lutts, S., 2012. Germination of untreated and primed seeds in rapeseed (brassica napus var oleifera del.) under salinity and low matric potential. Experimental Agriculture 48, 238–251.
  7. Ritz, C., Jensen, S. M., Gerhard, D., Streibig, J. C. (2019) Dose-Response Analysis Using R CRC Press. Achim Zeileis, Torsten Hothorn (2002). Diagnostic Checking in Regression Relationships. R News 2(3), 7-10. URL: https://CRAN.R-project.org/doc/Rnews/

Some further detail

Let us conclude this page by giving some detail on all equations.

The equation for the model HTnorm(). Here, we show all other equations, as implemented in our package.

HTL()

\[ P(t, \Psi) = \frac{1}{1 + exp \left[ - \frac{ \Psi - \left( \theta _H/t \right) - \Psi_{b(50)} } {\sigma} \right] }\]

HTG()

\[ P(t, \Psi) = \exp \left\{ { - \exp \left[ { - \left( {\frac{{\Psi - (\theta _H / t ) - \mu }}{\sigma }} \right)} \right]} \right\} \]

HTLL()

\[ P(t, \Psi) = \frac{1}{1 + \exp \left\{ \frac{ \log \left[ \Psi - \left( \frac{\theta _H}{t} \right) + \delta \right] - \log(\Psi_{b50} + \delta) }{\sigma}\right\} }\]

HTW1()

\[ P(t, \Psi) = exp \left\{ - \exp \left[ - \frac{ \log \left[ \Psi - \left( \frac{\theta _H}{t} \right) + \delta \right] - \log(\Psi_{b50} + \delta) }{\sigma}\right] \right\}\]

HTW2()

\[ P(t, \Psi) = 1 - exp \left\{ - \exp \left[ \frac{ \log \left[ \Psi - \left( \frac{\theta _H}{t} \right) + \delta \right] - \log(\Psi_{b50} + \delta) }{\sigma}\right] \right\}\]