Plantecophys - An R Package for Analysing and Modelling Leaf Gas Exchange Data

Here I present the R package 'plantecophys', a toolkit to analyse and model leaf gas exchange data. Measurements of leaf photosynthesis and transpiration are routinely collected with portable gas exchange instruments, and analysed with a few key models. These models include the Farquhar-von Caemmerer-Berry (FvCB) model of leaf photosynthesis, the Ball-Berry models of stomatal conductance, and the coupled leaf gas exchange model which combines the supply and demand functions for CO2 in the leaf. The 'plantecophys' R package includes functions for fitting these models to measurements, as well as simulating from the fitted models to aid in interpreting experimental data. Here I describe the functionality and implementation of the new package, and give some examples of its use. I briefly describe functions for fitting the FvCB model of photosynthesis to measurements of photosynthesis-CO2 response curves ('A-Ci curves'), fitting Ball-Berry type models, modelling C3 photosynthesis with the coupled photosynthesis-stomatal conductance model, modelling C4 photosynthesis, numerical solution of optimal stomatal behaviour, and energy balance calculations using the Penman-Monteith equation. This open-source package makes technically challenging calculations easily accessible for many users and is freely available on CRAN.


Introduction
Since the advent of portable gas exchange instruments [1,2], a wealth of data on leaf gas exchange of CO 2 and H 2 O has been collected [3]. These data play a central role in physiological plant ecology [4], to better understand and quantify inter-specific differences in photosynthesis and transpiration, and to quantify and model the rapid response to changes in environmental drivers such as light, humidity and temperature. Not only do leaf gas exchange data allow detailed studies of the underlying plant physiology, they are also used to parameterize an important component of process-based models of vegetation function used to predict global water and carbon cycling [5,6].
The photosynthesis model of Farquhar, von Caemmerer and Berry [7] (the 'FvCB model') is widely used in interpreting and modelling leaf gas exchange, by providing comparable metrics of the photosynthetic capacity, and predicting the response of photosynthesis to changes in the CO 2 concentration inside the leaf air space (C i ). This widely cited model is embedded in many process-based models of vegetation function [5,8]. The key prediction of the model is the response of photosynthesis to [CO 2 ] inside the leaf (either chloroplastic [CO 2 ], C c , or intercellular [CO 2 ], C i ). It can also account for changes in leaf temperature if the various temperature sensitivities are parameterized [9][10][11]. To employ the model, it is generally fit to observations of net photosynthesis along a range of [CO 2 ] concentrations, yielding well-known measures of photosynthetic capacity (V cmax and J max , and optionally R d ) [12].
I do not repeat a detailed description of the FvCB model here, as it has been described many times [10,11]. But generally it is of the form, where A n is the net rate of CO 2 assimulation, A c is the gross photosynthesis rate when Rubisco activity is limiting, A j when RuBP-regeneration is limiting, and R d the rate of dark respiration (see Fig 1A). A c and A j are non-linear functions of the chloroplastic CO 2 concentration (C c ), both of the form k 1 (C c Γ Ã )/(k 2 +C c ), where Γ Ã is the CO 2 compensation point without R d , and k 1 and k 2 are different parameter combinations for A c and A j . The details of these functions and the temperature dependence of the various parameters are described elsewhere [11].
In the practical application of the FvCB model, when leaf gas exchange is measured with a portable gas exchange instrument, estimates of C c are difficult to obtain because they require an estimate of the mesophyll conductance (g m ). In this case, it is customary to us the intercellular [CO 2 ] concentration (C i ) as the driver of photosynthesis. This approach is useful because C i can be estimated from concurrent measurements of CO 2 and H 2 O flux [13]. In the remainder of this article I will use C c as the driver of photosynthesis, but point out that this can be replaced by C i if the user does not have an estimate of g m . When g m is known, C c is calculated from C i with Eq 2.
where g m is the mesophyll conductance (mol m -2 s -1 ). Although this method assumes that g m is constant for a given leaf, it is well known that g m responds dynamically to fluctuations in environmental drivers [14] although some of the variation in g m may due to artefacts related to (photo-)respiratory effects on measured g m with standard methods [15]. Because no model has been developed to date that adequately captures the variation in g m , it is a constant parameter in the implementation presented here. However, it can still be used to study the effects of nonconstant g m on rates of photosynthesis and its response to environmental drivers, as the parameter can be varied in model simulations. Fitting the FvCB model to data requires some finesse because net photosynthesis is modelled as a minimum function of two non-linear equations that is sometimes difficult to fit. Moreover, sample sizes collected are often small due to time constraints in the field. A widely used published method requires the user to specify the transition of V cmax to J max limitation [16], a process that is both arbitrary and prevents batch analysis. Another method [17] requires online submission of data and fits the model without much control or knowledge of the fitting process (following [18]), and does not report standard errors of the estimated parameters. Undoubtedly many more implementations of the fitting process have been developed over the years, but few of these are made publicly available (but see available online tools [19,20]). What is missing is an open-source tool that can be used for reproducible and transparent analysis of A-C i curves.
Through Eq 1, we have a dependency of photosynthesis on the availability of the substrate, C c . To estimate C c itself, we need C i , which can be estimated when with stomatal conductance to CO 2 . From Fick's law, we can relate A n to g s and C i as, where g s is the conductance to H 2 O (the factor 1.6 converts to conductance to CO 2 ). We now have two equations for A n : the 'demand function' (Eq 1), and the 'supply function' (Eq 3). At steady state these two equations should be equal, which can be graphically shown as in Fig 1B  (cf. [21]). Because g s itself responds to environmental drivers, another expression is needed to end up with a fully coupled model of leaf gas exchange. The most widely-used, though empirical, g s model is the Ball-Berry [22] class of models. This model posits an entirely empirical equation that describes the response of g s to air humidity, CO 2 and A n . This way, effects of leaf temperature and PPFD-both of which are known to affect g s -are modelled through the dependency of A n on these drivers. A general form of the Ball-Berry model is, where D is the vapour pressure deficit (kPa), g 0 and g 1 are empirical parameters, and f(D) can be one of many functions that describe the response to the vapour pressure deficit (D, [23,24]) or relative humidity [22]. An alternative approach to modelling g s is through the hypothesis that stomata act optimally in the sense that they maximize photosynthesis while minimizing water loss. This hypothesis was first developed by Cowan and Farquhar [25] and has seen many applications. Medlyn et al. [24] showed that the optimality hypothesis, when coupled to the FvCB model, leads to an expression analogous to the Ball-Berry type models (Eq 4), but with a different D response function (f(D) in Eq 3) compared to the original Ball-Berry model. Finally we can combine the biochemical demand function of photosynthesis (Eq 1) with the supply function (Eq 2) and an expression for the dependency of g s on environmental drivers (Eq 3). This 'coupled' leaf gas exchange model [23,26,27] is implemented in many processbased ecosystem and global land surface models [5,8,28,29]. This model allows prediction of A n , g s and leaf transpiration rate in response to all major environmental drivers (except soil water limitation), and incorporates key leaf traits (g 1 , V cmax , J max , R d , and their temperature dependencies).
Despite the widespread use of the FvCB model and the coupled leaf gas exchange model, tools to analyse data and perform simulations are scattered and subject to little standardization. Fitting the FvCB model to CO 2 response curves is a standard procedure but different methods can yield different parameter values, making comparisons difficult. The coupled leaf gas exchange model is not straightforward to implement, and I do not know of any standalone open-source implementations. I here describe the plantecophys package, implemented in the R language [30]. The code is freely available (without restrictions), and managed with a version control system. The package is the result of our work on leaf and canopy modelling of photosynthesis and stomatal conductance [24,[31][32][33][34][35][36], with many additions based on user requests.

Design and Implementation
The main functions The main tools included in the plantecophys package are to a) fit A-C i curves to estimate V cmax , J max and R d , b) fit Ball-Berry type models, c) simulate from the coupled leaf gas exchange model and d) calculate the optimal stomatal conductance. The key functions in the package are summarized in Table 1.

Language
The 'plantecophys' package is implemented in R, has no dependencies on other packages, and does not require compilation (i.e. it is written in native R only). As such it builds easily, and is highly portable. The source code is maintained with git version control, and is hosted in an online repository (http://www.bitbucket.org/remkoduursma/plantecophys), from which a development version of the package can easily be installed. The repository includes an issue tracker, where users can suggest changes or report bugs. This paper describes version 0.6.6 (git SHA b9a18c9). All code used in this article (including the code to generate the article written in markdown, all figures and full example code), can be downloaded from Ref. [37]. The repository also includes code to demonstrates how to extract additional statistics from fitted A-C i curves.

Results and Discussion
Fitting A-C i curves The fitaci function fits the FvCB model, yielding estimates of V cmax , J max and R d and their standard errors. Instead of fitting the minimum function (Eq 1), fitaci fits the hyperbolic minimum of A c and A j , which avoids a discontinuity (Eq 5).
where θ is a shape parameter, set to 0.9999, and A m is the hyperbolic minimum of A c and A j .
The fit of the FvCB model to data is achieved with non-linear least squares, and standard errors of the parameters are estimated with standard methods (nls function in base R, see [38]). The fitaci function includes methods to estimate appropriate starting values from the data, and attempts the fits along a wide range of possible starting values. Optionally, R d can be provided as a known value, otherwise it is estimated from the A-C i curve. The user does not have to provide the transition point (see Fig 1), as this is estimated by fitaci automatically. It is however an option to fix the transition point (via the citransition argument), which may be helpful to check whether the best fit was achieved. Finally, the user can provide an estimate of mesophyll conductance (g m ) (following [39]), in which case the fitted values of V cmax and J max can be interpreted as chloroplastic rates.
Because the fitting uses non-linear least squares, standard methods can be employed to estimate standard errors (SE), confidence intervals, and correlation of the fitted parameters. The fitaci function returns by default the SE and confidence intervals, and the built-in help page for the fitaci function shows how the nlstools package can be used to provide a detailed overview of the statistics of the non-linear least squares fit.
Required inputs are measurements of A n and C i , and optionally leaf temperature (T leaf ), and photosynthetically active radiation (PAR). Also required are estimates of Michaelis-Menten constants (K c , K o or the combination K m ) and Γ Ã . In the FvCB model, J max , V cmax and leaf respiration (R d ) (and other parameters like Γ Ã , K c and K o ) all depend non-linearly on T leaf . The Photosyn function incorporates standard temperature sensitivities for all parameters of the FvCB model (following [11]). Optionally, measured (or otherwise modelled) K m and Γ Ã can be provided as input.
The function takes a dataframe as input, which includes measurements of A n , C i and optionally T leaf and PAR, and is easily used like this: # Fit FvCB model f <-fitaci(mydata) # Print a summary with coefficients and more f # Make standard plot lot(f) The output of the above example is shown in Fig 2. Additionally, the batch utility fitacis can be used to fit many curves at once, for example one for each species or site in a dataset. I show this functionality in the example application further below.
A C4 model of leaf photosynthesis [40] is also implemented (in AciC4), but at the moment it is only possible to fit the C3 model of leaf photosynthesis to A-C i curves.

Fitting stomatal conductance models
The straightforward fitBB function provides an interface to non-linear or linear regression to fit one of three stomatal conductance models [22][23][24]. This yields estimates of g 1 and (optionally) g 0 , which are necessary inputs to the coupled leaf gas exchange model. Note that the user must provide stomatal conductance to H 2 O (not CO 2 ) as input to the fitting process, which is the standard output of portable gas exchange instruments. This function is demonstrated in the example application further below.

Coupled leaf gas exchange model
The intersection of the supply and demand curves of photosynthesis (Fig 1B) gives the steadystate intercellular CO 2 concentration (C i ). This is solved by the Photosyn function. This flexible interface can be used to either 1) estimate A n when C i is known (Photosyn(Ci = . . .); equivalent to Aci(. . .)), 2) estimate A n when g s is known (Photosyn(GS = . . .)) (cf. Fig 1B) or c) solve for C i from the coupled leaf gas exchange model (Eqs 1,3 and 4).
To demonstrate the use of the coupled gas exchange model, I visualize the temperature response of A n when both T leaf and D are varying. In field conditions, D is always strongly positively related to T leaf . The consequence is that when studying D or T leaf responses in the field, both drivers have to be accounted for simultaneously [35,41]. Fig 3B shows simulated A-C i curves and the solutions of the coupled leaf gas exchange models at a range of T leaf and corresponding D (calculated following [35]). Both V cmax and J max have a peaked response to T leaf , so that at a given C i , A n first increases with T leaf and then decreases (lines, Fig 3A). As a result of increasing D, the modelled C i decreases (symbols, Fig 3A, as a consequence of Eq 4). The net result is a peaked response of A n as a function of D (Fig 3B).
The simplified code to produce Fig 3B, using the Photosyn function, is given below. Note that in this example the default values of many parameters (e.g. J max , g 1 ) are used in the call to Photosyn, but all of these can be set by the user.
# Set range of leaf temperature tleafs <-seq(5, 40, by = 5) # Define D as a function of Tleaf vpdfun <-function(tair)0.000605 Ã tair^2.39 # Simulate. run1 <-Photosyn(Tleaf = tleafs, VPD = vpdfun(tleafs)) # Plot (produces Fig 3b minus the special formatting) with(run1, plot(Tleaf, ALEAF)) The Photosyn function assumes that the boundary layer conductance (g bl ) is high compared to g s , so that T leaf is close to T air . As an alternative, the PhotosynEB function calculates T leaf from the leaf energy balance. Transpiration is calculated with the Penman-Monteith equation [42], which accounts for boundary layer effects. The details of PhotosynEB are not described here (see the built in help file for more information), because it is very similar to other implementations [28,43].

Numerical solution of optimal stomatal conductance
The FARAO function (FARquhar And Optimality) calculates optimal stomatal conductance based on the Cowan-Farquhar [25] hypothesis that stomata respond to environmental drivers in order to maximize photosynthesis while minimizing water loss. This implementation was used by Medlyn et al. [24] to compare a simplified model of optimal stomatal conductance to the full numerical solution.
To find optimal stomatal conductance, FARAO finds the C i for which the quantity A n −λE is maximal, where E is the leaf transpiration rate and λ is the marginal cost of water (an empirical parameter related to g 1 , see [24,25]). A n is calculated directly as a function of C i via the FvCB model (Eq 1), g s is calculated by rearranging Eq 3, and E is calculated assuming perfect coupling (thus E = g s D/P a , where P a is atmospheric pressure). This numerical routine does not need specification of an f(D) function as in Eq 4, instead, this function is an emergent property. In Fig 4A, I have calculated A n −λE across a range of C i values, and for different values of D. The FARAO function calculates the optima of these curves, which can for example be used to study the response of stomatal conductance to D (Fig 4B).
Optionally, the FARAO function accounts for the presence of a leaf boundary layer (when energybalance = TRUE). In that case it uses PhotosynEB (see description above) to calculate A n and E, and solves for T leaf . A very similar method was employed by Buckley et al. [43], who demonstrated that when a boundary layer is present, frequently an optimal g s cannot be found.

An example application
To demonstrate a practical application of the key functions in the package, I use field-collected data from Medlyn et al. [44,45] on Eucalyptus delegatensis. Both A-C i curves and 'spot gas exchange' data (i.e. leaf gas exchange measurements at prevailing environmental conditions) were collected. Using the fitacis function, it is straightforward to fit all 43 curves to the A-C i data, and make standard plots of the fitted curves (shown in Fig 5A). The fitted coefficients can   be extracted using the coef function, and used to plot a comparison of fitted V cmax and J max values, which show the typical correlation between the two (Fig 5B).
Next, I fit Eq 4 to the spot gas exchange data, yielding an estimate of g 1 (Fig 5C). In this example, I used the model of Medlyn et al. [24], which is given by Eq 6 (in this example, I assumed g 0 = 0) In Fig 5C, modelled g s is compared to measurements. To compare the model prediction of instantaneous transpiration efficiency (A n /E) to measurements along the variation in D ( Fig  5D), Eq 6 can be rearranged to give (cf. [34], where it is assumed that g 0 = 0) Because fitBB can fit a number of Ball-Berry type variants, the various models can be easily compared in terms of goodness of fit. This simple example application is available in the published repository (see Methods), and simplified code for this example (panels a-c only), omitting special formatting and a few minor settings, is given below.

Conclusions
We need an open source set of tools to analyse leaf gas exchange data, as these data form a cornerstone of plant physiological ecology. At the moment there are no publicly available tools to fit A-C i curves or perform simulations with the coupled leaf gas exchange model that can be used as part of a reproducable workflow. The plantecophys R package is implemented in widely used language for data analysis. The package includes a useful set of tools to perform standard, and more advanced, analyses of leaf gas exchange data. The open source framework combined with version control allows further development of the code.