Water is transported through the xylem under tension and is thus prone to cavitation, the rapid phase change from liquid to water vapour (Dixon and Joly, 1895; Tyree and Sperry, 1989). This results in the formation of gas bubbles (emboli) that block xylem conduits and reduce the hydraulic conductivity of the xylem. The probability of embolism formation is greatly increased by environmental stresses such as drought. During drought, the tension necessary to draw water from the soil increases as soil water content declines. At some critical threshold, embolism begins to spread through the xylem network via air seeding of pit membranes (Tyree and Zimmermann, 2013; Choat
</p><p> Plants vary widely in their vulnerability to embolism with vulnerability strongly related to the minimum seasonal water potential experienced in the field (Choat
</p><p> A range of methods are now employed to generate vulnerability curves and can generally be divided into two classes; those that expose a sample (branch, root, leaf) to repeated measurements and those in which each sample is only measured once. For instance, the bench dehydration method uses a collection of individual samples to construct a single composite curve (Sperry
</p><p> Clearly, when making comparisons, statistical testing of the differences in P_{50} is not possible unless the confidence interval (CI) is reported. The widely used method by Pammenter and Van der Willigen (1998) does not discuss the CI of P_{50} (or of the slope parameter), but instead recommends that statistical comparisons are made based on replicate curves, not based on the uncertainty of the parameter estimates from individual curves. In many cases, however, it is not possible to collect many replicate curves because of the time consuming nature of measurements or limitations in plant material available for destructive harvests. It is theref ore also useful to report the CI of the fitted parameters for each individual curve, which allows statistical comparison across different groups based on composite curves.
</p><p> In many cases the uncertainty of vulnerability curve parameter estimates has not been reported in the plant hydraulics literature, and yet conclusions have been drawn on apparent differences between groups. Indeed, the widely used method by Pammenter and Van der Willigen (1998) does not discuss confidence intervals or standard errors of P_{50} (or of the slope parameter  see Methods). A hierarchical Bayesian method proposed by Ogle
</p><p> Here we present an open source implementation of a fitting routine that can be used in estimating P_{50} (or any other point on the curve) and its statistical uncertainty. There is a clear need for such a utility to allow for fully reproducible results in plant hydraulics research. The
The
</p><p> The user either provides the data as PLC, or as K. Either variable is first converted to the conductance relative to the maximum value (K/K_{max}) as,
(1)
and
Using the above relationships between the variables, it is straightforward to express any measurement as the relative conductance (K/K_{max}).
The Weibull model is given by Eq. 2, where K/K_{max} is expressed as a function of P, and the two parameters to be estimated (P_{x} and S_{x}).
(2)
where
where
where P_{x} is the xylem pressure (P) where x% of the conductivity is lost, S_{x} is the derivative (% MPa^{1}) at x (e.g. S_{50} is the slope of the curve at P_{50}). Higher values of S_{x} thus indicate a steeper response to P. This reparameterization from Ogle
</p><p> Because we use nonlinear regression for the Weibull model, suitable starting values have to be estimated in order for the solution to converge. We estimate starting values for P_{x} and S_{x} from the fit of the sigmoidal model, which can be fit using linear regression and thus always yields a solution .
We have also implemented a sigmoidalexponential model proposed by Pammenter and Van der Willigen (1998), as follows. The equation is,
</p><p>
(3)
Eq. 3 can be linearized to allow use of linear regression (Eq. 4),
</p><p>
(4)
Eq. 4 can be directly used in a linear regression, yielding ab (as the intercept), and a (as the slope), from which P_{x} can be estimated as
</p><p>
(5)
Finally, the slope parameter S_{x} is found as the derivative of Eq. 3 evaluated at P_{x} (equation not shown for brevity). Because the parameters P_{x} and S_{x} are nonlinear parameter combinations, we report only the bootstrap confidence intervals for the sigmoidal model (these will be more appropriate anyway, given there is often nonconstant error variance in hydraulic vulnerability curves, see discussion by Ogle
We implemented the nonparametric bootstrap (Efron and Gong, 1983) to estimate confidence intervals for the estimated parameters, both for the Weibull and sigmoidal models. A common problem with bootstrapping nonlinear regression models is that they do not always converge, which limits application to smaller sample sizes (Fox, 2002). We avoid this problem by estimating P_{x} and S_{x} using the sigmoidal model first, and use these as accurate starting values, thus permitting bootstrap confidence intervals even with small sample sizes.
</p><p> With the bootstrap, the original data are resampled with replacement N times (typically, N = 1000), and each time the regression is applied to the resampled dataset to provide estimates of the fitted parameters as well as the fitted curve at a range of values of P. The bootstrap was implemented as a simple nonparametric 'case resampling' method very similar to that described in the 'bootCase' function in the ‘ car ’ package (Fox and Weisberg, 2010).
When fitting an ensemble of curves, for example when measurements were made on many branches for a single individual or species, we use a (non)linear mixed effects (LME) model approach (Pinheiro and Bates, 2006). In this situation, it would not be appropriate to combine all individual data points and fit a single curve, as measurements are correlated within branches (and thus, the SE will typically be underestimated). Instead, the common approach is to fit curves to each branch separately, and average the coefficients. However, in some cases the data do not lend themselves well to allow curve fitting on each branch separately (in particular, when some individual branches show a very poor fit). In addition, LME will be in general more robust to individuals with poor data (Pinheiro and Bates, 2006). A similar approach was taken by Peek
</p><p> Using LME, we obtain a single population average curve (the fixed effect) as well as the branchtobranch variance in P_{x} and S_{x}. Also estimated are the individual random effects (the socalled BLUPs, see Robinson (1991)), which can be used to visualize the differences between individual branches within the individual. The BLUPs are more robust than fitting curves to each branch separately particularly when some individuals show a poor fit. For the Weibull model, we use a nonlinear mixedeffects model, and for the sigmoidal model a linear mixedeffects model. Both are fit using the ‘ nlme ’ R package (Pinheiro and Bates, 2006).
Following the discussion by Nolf
We implemented the fitting routines as an R package (
</p><p> The code is maintained in an online repository (
). Issues or feature requests can be submitted there, and installation instructions are also available. This article is not intended as a full technical reference, for that we refer to the builtin help pages for the package in R.
To illustrate the functionality of the
In the following examples, we demonstrate the use of the
</p><p> The first example shows how we can estimate P_{12}, P_{50} and P_{88} from a hydraulic vulnerability curve, and make a simple plot of the fitted curve and the raw data (Fig. 1). Note how the result from
# Fit three curves </p><p> fit12 < fitplc (mydata, x= 12 ) </p><p> fit50 < fitplc (mydata, x= 50 ) </p><p> fit88 < fitplc (mydata, x= 88 ) </p><p> </p><p> # Plot </p><p> plot (fit12, what= "PLC" ) </p><p> plot (fit50, add= TRUE , what= "PLC" ) </p><p> plot (fit88, add= TRUE , what= "PLC" )
In the following examples we focus on estimating P_{50} only, but note that any analysis can be quickly redone when some other parameter should be estimated. One of the strengths of the
</p><p> The example is shown in Fig. 2, where two species with very different vulnerability are compared in terms of their P_{50} and fitted curve. In this case, PLC at any xylem pressure is different between the species, as indicated by the nonoverlapping 95% confidence intervals. Clearly, the P_{50}'s are also very different. Plots such as these are easy to construct with the default plotting routine, and refer to the online repository for the complete code to the examples.
When the conductivity has been scaled to a maximum value, we can use the
If the raw data have not been converted into PLC, but are available only as raw conductance (or conductivity) values, the
</p><p>
The basic code to fit the example is shown below: here we fit the Weibull curve to conductivity directly, and scale the data by maximum conductivity calculated from the data where the water potential > 0.3 MPa (an arbitrary choice for this example).
fitc < fitconds (mydata, group= "Species" , WP_Kmax =  0.3 )
</p><p>
Finally we demonstrate the inclusion of a random effect. In both example datasets, multiple branches were measured for a species. The data points are not independent since the measurements on a single branch will be correlated. It is thus appropriate to fit a mixedeffects model. The fitted model is visualized in Fig. 5, including the mean response (i.e. the 'fixed effect') estimated by the model as well as the individuallevel predictions for each branch (i.e. the random effects).
Fig. 5 shows the results from using
# Fit curves with 'individual' as a random effect (name of the variable in the dataframe ) </p><p> fitran < fitplc (mydata, random= individual) </p><p> </p><p> # Standard plot, add predictions for random effect. </p><p> plot ( fitran, plotrandom= TRUE )
We have so far demonstrated only the use of the Weibull model to estimate hydraulic parameters. The sigmoidal model is also implemented in the
fitsig < fitplc (mydata, model = "sigmoidal" )
A comprehensive comparison of the Weibull and sigmoidal models is well outside the scope of this paper. However, we found a few problems with the sigmoidal model, in particular the logtransformed version as presented by Pammenter and Van der Willigen (1998). The first problem, as reported by Ogle
The
</p><p> It is important to point out that in nonlinear regression, the confidence interval cannot be simply computed from the standard error (as the usual approximation of 2 * SE). Instead, the confidence interval is computed based on profiling the likelihood functions of the parameters (Ritz and Streibig, 2008), and is frequently asymmetric (see Table 1). It is thus important to use confidence intervals for inferences, not the standard errors directly (e.g. as input to a ttest). For example, when 95% confidence intervals for P_{50} do not overlap between two species, one may conclude that P_{50} is significantly different between the two species (p < 0.05). For this reason, the
</p><p>
titre
titre du tableau  








S_{50}  20.81  12.04  34.59  5.75  34.44 
P_{50}  4.31  3.82  4.8  3.74  6.33  

S_{50}  61.93  32.9 

36.48  350.91 
P_{50}  3.43  3.21  3.64  3.15  3.69  

S_{50}  124.92  67.86 

54.84  615.46 
P_{50}  1.42  1.31  1.52  1.27  1.54  

S_{50}  12.5  7.91  19.23  7.73  21.05 
P_{50}  11.61  10.83  12.35  10.6  12.65  

S_{50}  30.63  20.71  43.84  20.56  45.32 
P_{50}  2.87  2.54  3.21  2.53  3.25  

S_{50}  39.32  24.85  66.71  26.02  92.96 
P_{50}  3.07  2.78  3.38  2.66  3.41  

S_{50}  64.26  36.35  132.53  47.44  165.69 
P_{50}  2.34  2.11  2.55  2.14  2.65 
We presented a new package that includes a range of statistical methods to fit hydraulic vulnerability curves, and report appropriate uncertainty measures. The main novelty of the approach is the use of the bootstrap to compute confidence intervals of the parameters, as well as around the entire fitted curve, and the use of (non)linear mixedeffects models in cases where an ensemble of curves is measured. We discuss the difference in approach taken with the more common methods reported in the literature below, and provide some ideas for future work.
Our implementation of a fitting routine for vulnerability curves is novel because it reports uncertainty of not just the fitted parameters, but also around the entire fitted curve, using the bootstrap. We do not claim that this is an entirely new method, indeed, this approach has been used frequently to report uncertainty of fitted nonlinear curves. The bootstrap is widely used in ecology (Crowley, 1992), for example to report uncertainties in phylogenies ( e.g., over 30 thousand citations to Felsenstein, 1985). We hope that the userfriendly implementation presented here will allow wider uptake in the field of plant hydraulics.
</p><p> The use of the bootstrap to approximate the uncertainties differs from the fully Bayesian approach presented by Ogle
The use of (non)linear mixedeffects models (LME) is still quite rare in the analysis of response curves, although the technique is far from new (Peek
We leave the door open for future improvements and contributions to the
We thank Markus Nolf, Danielle Creek, Chris Blackman and Rosana Lopéz for testing and suggesting features. This manuscript is entirely reproducible, including all figures and analyses. The code to generate this manuscript, which was written with the rmarkdown package, as well as all examples can be downloaded from :
.