We describe a toolkit to fit hydraulic vulnerability curves, such as the percent loss of xylem hydraulic conductivity ('PLC curves') as a function of the water potential. The toolkit is implemented as an R package, and is thus free to use and open source. The package fits the Weibull or sigmoidal function to measurements of PLC, conductance or conductivity, at corresponding leaf or stem water potentials. From the fitted curve, estimates of P _x (the water potential at which x% conductivity is lost, e.g. the P ₅₀ ), and slope parameter (S _x ) are provided together with confidence intervals (CI) around the fitted line. The CIs are estimated with the bootstrap. We also demonstrate the advantages of using mixed-effects models in situations where multiple individuals are measured on a species, as compared to the more traditional approach of fitting curves separately and averaging the parameters. We demonstrate the use of the new package with example data on seven species measured with two different techniques.

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 et al., 2015). During extreme droughts, embolism blockages can cause complete hydraulic failure in the xylem pathway and subsequent death of the plant (Kursar et al., 2009; Urli et al., 2013; Li et al., 2016a). Indeed, hydraulic failure is now considered to be a principal cause of tree mortality during drought (Brodribb and Cochard, 2009; Anderegg et al., 2016). This is particularly relevant since the hydraulic system of plants is finely tuned to their growing environment with the majority of woody plant species across major forest biomes maintaining narrow safety margins to hydraulic failure (Choat et al., 2012). Given the importance of vulnerability to embolism for plant survival during periods of environmental stress, it is essential that appropriate statistical methodology is developed for comparison of between and within species.

</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 et al., 2012). This variation in vulnerability results from differences in xylem anatomy, principally the physical characteristics and distribution of inter-conduit pit membranes (Tyree et al., 1994; Choat et al., 2008; Pittermann et al., 2010; Li et al., 2016b). Vulnerability to embolism is described by the relationship between loss of hydraulic conductivity in the xylem and xylem water potential, which produces a vulnerability curve (Fig. 1). Because of the importance of the trait to survival under conditions of drought stress, it is common to investigate differences in vulnerability to embolism among species, populations, plant components or experimental treatments. Vulnerability is compared using the xylem water potential at which some fixed proportion of hydraulic conductivity is lost, most commonly 50%, (known as the P₅₀) although other points are also used, e.g. P_e, P₂₅, P₈₈, P₉₅ (Tyree and Ewers, 1991; Tsuda and Tyree, 1997; Brodribb and Hill, 1999; Sparks and Black, 1999; Choat et al., 2003; Schuldt et al., 2015). These parameters are determined by fitting a response curve to the data for each species, population or organ.

</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 et al., 1988; Tyree et al., 1993). In techniques based on air injection, centrifugation, and imaging repeated measurements of percentage loss hydraulic conductivity or percent embolism can be made at different xylem pressures (Cochard et al., 1992; Cochard et al., 2005; Pockman et al., 1995; Choat et al., 2015). Vulnerability curves have been fitted to these data using a variety of functions including sigmoidal, Weilbull and polynomial (Pammenter and Van der Willigen , 1998; Ogle et al., 2009).

</p><p> Clearly, when making comparisons, statistical testing of the differences in P₅₀ 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₅₀ (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₅₀ (or of the slope parameter - see Methods). A hierarchical Bayesian method proposed by Ogle et al. (2009) estimates full uncertainties of the parameters but is difficult to implement and has not seen widespread uptake by the field. Apart from uncertainty of the parameters that describe the curve, it would also be very useful to be able to compute confidence intervals along the entire fitted curve, in order to judge statistical significance of differences between curves at any value of the water potential.

</p><p> Here we present an open source implementation of a fitting routine that can be used in estimating P₅₀ (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 fitplc package, a toolkit written in the widely used R language, is the first open source toolkit that fits PLC curves, and is easy to use and install. We have previously used and briefly described an old version of this package (Nolf et al., 2015). Here we describe the current functionality in detail, and demonstrate its use with several examples.

Hydraulic vulnerability parameters

The fitplc package fits non-linear curves to measurements of 'percent loss of conductivity' (PLC) or raw conductance/conductivity (K) at corresponding leaf or stem water pressure (P). From these curves, two parameters are estimated: the P_x (the value of P where x% of the conductivity is lost, for example P₅₀) and a slope parameter, i.e. the derivative (S_x, % MPa^-1). The package provides the choice of two models to fit the data: the Weibull (as reparameterized by Ogle et al. (2009)), and a sigmoidal model (Pammenter and Van der Willigen , 1998).

</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,

P L C = K m a x - K K m a x ⋅ 100

(1)

and

K = K m a x * ( 1 - P L C 100 )

Using the above relationships between the variables, it is straightforward to express any measurement as the relative conductance (K/K_max).

The Weibull model

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).

K / K m a x = ( 1 - x / 100 ) p

(2)

where

p = ( P / P x ) P x S x / V

where

V = ( x - 100 ) l o g ( 1 - x / 100 )

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₅₀ is the slope of the curve at P₅₀). Higher values of S_x thus indicate a steeper response to P. This reparameterization from Ogle et al. (2009) allows a straightforward implementation as a non-linear regression model, thus enabling computation of standard errors with textbook methods. We use ordinary non-linear least squares (with the 'nls' function in R) to fit Eq. 2 to observations of either PLC or K (converting either to K/K_max based on Eq. 1) measured at a range of water potentials (P). Confidence intervals (CIs) of the fitted parameters are estimated using standard profiling methods (Ritz and Streibig, 2008), and with the bootstrap (see below).

</p><p> Because we use non-linear 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 .

The sigmoidal model

We have also implemented a sigmoidal-exponential model proposed by Pammenter and Van der Willigen (1998), as follows. The equation is,

</p><p>

P L C = 100 1 + e x p ( a ( P - b ) )

(3)

Eq. 3 can be linearized to allow use of linear regression (Eq. 4),

</p><p>

l o g ( 100 / P L C - 1 ) = a P - a b

(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>

P x = l o g ( 1 / ( 1 - x / 100 ) - 1 ) / a + b

(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 non-linear parameter combinations, we report only the bootstrap confidence intervals for the sigmoidal model (these will be more appropriate anyway, given there is often non-constant error variance in hydraulic vulnerability curves, see discussion by Ogle et al. (2009)).

Bootstrap confidence intervals

We implemented the non-parametric 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 non-linear 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 non-parametric 'case resampling' method very similar to that described in the 'bootCase' function in the ‘ car ’ package (Fox and Weisberg, 2010).

Random effects

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 et al. (2002) when fitting non-linear curves of the response of photosynthesis to light.

</p><p> Using LME, we obtain a single population average curve (the fixed effect) as well as the branch-to-branch variance in P_x and S_x. Also estimated are the individual random effects (the so-called 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 non-linear mixed-effects model, and for the sigmoidal model a linear mixed-effects model. Both are fit using the ‘ nlme ’ R package (Pinheiro and Bates, 2006).

Weights

Following the discussion by Nolf et al. (2015), we also allow for the specification of a weighting function. Nolf et al. argued that the data around the P₅₀ should be weighted more heavily, when we are primarily interested in the fit of the curve in that region. We have made no a priori choice on the form of the weighting function, but allow for its specification when fitting the Weibull or sigmoidal models, and suggest more research is needed on the correct choice of the weighting function.

Implementation

We implemented the fitting routines as an R package (fitplc) (R Core Team, 2016), available on CRAN. The package also includes a few tools to visualize the fit along with the raw data. The package was implemented in native R (i.e., no compilation necessary), with minimal dependencies, allowing high portability and ease of installation. The implementation has a user-friendly command-line interface. When fitting a vulnerability curve, the data need to be read into in a 'dataframe' in R. Column names and order are arbitrary and can be set when fitting the curve. The required data are xylem pressure (in MPa) and hydraulic conductivity, either as raw measurements of scaled to some maximum value determined by the user. We return to this point in one of the examples below.

</p><p> The code is maintained in an online repository (

www.bitbucket.org/remkoduursma/fitplc

). 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 built-in help pages for the package in R.

Example data

To illustrate the functionality of the fitplc package, we use example data on 7 diverse species. All hydraulic vulnerability curves were measured with benchtop dehydration except Callitris and Pinus radiata, which were measured with the centrifuge technique (see e.g. Choat et al., 2010, and Cochard et al., 2005 for methods). The method for measurement is not relevant to the applicability of the fitplc package.

Example applications

In the following examples, we demonstrate the use of the fitplc package. All example R code is simplified for reasons of clarity, omitting formatting and other minor settings. In all examples except the last, we only demonstrate the Weibull model (the default model).

</p><p> The first example shows how we can estimate P₁₂, P₅₀ and P₈₈ 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 fitplc is stored in an object that can be readily plotted. In this case we decide to plot % embolism ('percent loss conductivity', PLC) as a function of xylem pressure ( what = "PLC" - a synonym is what = "embol" ); the default is to plot relative conductivity as a function of xylem pressure (see later figures). By default, the Weibull model is fit to the data, unless the argument model = "sigmoidal" is specified.

# 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₅₀ only, but note that any analysis can be quickly redone when some other parameter should be estimated. One of the strengths of the fitplc package is the use of the bootstrap to estimate confidence intervals of P₅₀ as well as around the entire fitted curve. In the following example, we fit the Weibull curve to data from two species differing in sensitivity to stem xylem pressure, and make a simple plot. The calculation of the bootstrap confidence intervals is done by default, but can be switched off by the user if it is not needed to save time (though it only takes a few seconds to fit a curve and perform the bootstrap).

</p><p> The example is shown in Fig. 2, where two species with very different vulnerability are compared in terms of their P₅₀ and fitted curve. In this case, PLC at any xylem pressure is different between the species, as indicated by the non-overlapping 95% confidence intervals. Clearly, the P₅₀'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 fitplc function as shown in the previous examples. It is straightforward to scale conductivity measurements to the maximum observed for the dataset, however this approach is not robust to outliers. Another option is to scale the data relative to the average conductivity calculated for high water potentials (e.g. > -1 MPa). An example of the consequences of these two options is shown in Fig. 3.

If the raw data have not been converted into PLC, but are available only as raw conductance (or conductivity) values, the fitcond function can be used. In this case, the user can either specify the maximum conductance for the dataset (perhaps measured independently or directly from the data, somehow), or a threshold water potential (in which case the maximum conductance will be calculated from the data where water potential is above this threshold). The user may also fit multiple curves at once (using fitconds); all other options to fitplc are supported in fitcond as well. In Fig. 4, we demonstrate the use of fitcond, and show a standard plot comparing three species.

</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 mixed-effects 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 individual-level predictions for each branch (i.e. the random effects).

Fig. 5 shows the results from using fitplc when a random effect is included. Simplified code for these examples is shown below.

# 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 )

Comparison of Weibull and sigmoidal models

We have so far demonstrated only the use of the Weibull model to estimate hydraulic parameters. The sigmoidal model is also implemented in the fitplc package, and can be fit with the following command:

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 log-transformed version as presented by Pammenter and Van der Willigen (1998). The first problem, as reported by Ogle et al. (2009), is that the sigmoidal model does not guarantee that PLC is zero at P = 0; the PLC is often higher and depends on the overall response of PLC to P. In contrast, in the Weibull model, PLC = 0 at P = 0, which is more biologically reasonable. The second problem we found is that curves with a very steep response of PLC to P are poorly fit by the sigmoidal model (see Fig. 6a), but for other curves the fits are very comparable (Fig. 6b). Finally, the log-transformed version (Eq. 4) is very sensitive to very low values of conductivity, because values close to zero are large negative values when log-transformed, and may have a large impact on estimated coefficients (as we found with one example dataset; not shown). For these reasons we recommend the Weibull model as the default choice, but encourage more comprehensive comparisons.

Statistics

The fitplc package computes confidence intervals (CIs) of the fitted parameters with standard profiling methods or the bootstrap. It is well known that profiling standard errors in non-linear regression are often biased downward, thus underestimating the uncertainty of the fitted parameters. When possible, it is recommended to use the bootstrap resampling approach, as this provides reliable estimates of the CI. For the example datasets we have presented so far, we compared the width of the CI as reported by both methods (Table 1). As expected, the profiling method frequently underestimates the width of the CI, particularly for the slope parameter S_x. It is also worthwhile to point out that S_x often has very wide confidence intervals, demonstrating that this parameter is probably poorly constrained by the data, making meaningful comparisons between species difficult.

</p><p> It is important to point out that in non-linear 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 t-test). For example, when 95% confidence intervals for P₅₀ do not overlap between two species, one may conclude that P₅₀ is significantly different between the two species (p < 0.05). For this reason, the fitplc package does not report the standard errors (although they can be extracted if really needed).

</p><p>

titre

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 mixed-effects 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.