# Fit the Farquhar-Berry-von Caemmerer model of leaf photosynthesis

### Description

Fits the Farquhar-Berry-von Caemmerer model of photosynthesis to measurements of photosynthesis and intercellular *CO2* concentration (Ci). Estimates Jmax, Vcmax, Rd and their standard errors. A simple plotting method is also included, as well as the function `fitacis`

which quickly fits multiple A-Ci curves (see its help page). Temperature dependencies of the parameters are taken into account following Medlyn et al. (2002), see `Photosyn`

for more details.

### Usage

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | ```
fitaci(data, varnames = list(ALEAF = "Photo", Tleaf = "Tleaf", Ci = "Ci", PPFD
= "PARi", Rd = "Rd"), Tcorrect = TRUE, Patm = 100, citransition = NULL,
quiet = FALSE, startValgrid = TRUE, fitmethod = c("default",
"bilinear"), algorithm = "default", fitTPU = FALSE, useRd = FALSE,
PPFD = NULL, Tleaf = NULL, alpha = 0.24, theta = 0.85, gmeso = NULL,
EaV = 82620.87, EdVC = 0, delsC = 645.1013, EaJ = 39676.89,
EdVJ = 2e+05, delsJ = 641.3615, GammaStar = NULL, Km = NULL,
id = NULL, ...)
## S3 method for class 'acifit'
plot(x, what = c("data", "model", "none"), xlim = NULL,
ylim = NULL, whichA = c("Ac", "Aj", "Amin", "Ap"), add = FALSE,
pch = 19, addzeroline = TRUE, addlegend = !add, legendbty = "o",
transitionpoint = TRUE, linecols = c("black", "blue", "red"), lwd = c(1,
2), ...)
``` |

### Arguments

`data` |
Dataframe with Ci, Photo, Tleaf, PPFD (the last two are optional). For |

`varnames` |
List of names of variables in the dataset (see Details). |

`Tcorrect` |
If TRUE, Vcmax and Jmax are corrected to 25C. Otherwise, Vcmax and Jmax are estimated at measurement temperature. |

`Patm` |
Atmospheric pressure (kPa) |

`citransition` |
If provided, fits the Vcmax and Jmax limited regions separately (see Details). |

`quiet` |
If TRUE, no messages are written to the screen. |

`startValgrid` |
If TRUE (the default), uses a fine grid of starting values to increase the chance of finding a solution. |

`fitmethod` |
Method to fit the A-Ci curve. Either 'default' (Duursma 2015), or 'bilinear'. See Details. |

`algorithm` |
Passed to |

`fitTPU` |
Logical (default FALSE). Attempt to fit TPU limitation (fitmethod set to 'bilinear' automatically if used). See Details. |

`useRd` |
If Rd provided in data, and useRd=TRUE (default is FALSE), uses measured Rd in fit. Otherwise it is estimatied from the fit to the A-Ci curve. |

`PPFD` |
Photosynthetic photon flux density ('PAR') (mu mol m-2 s-1) |

`Tleaf` |
Leaf temperature (degrees C) |

`alpha` |
Quantum yield of electron transport (mol mol-1) |

`theta` |
Shape of light response curve. |

`gmeso` |
Mesophyll conductance (mol m-2 s-1 bar-1). If not NULL (the default), Vcmax and Jmax are chloroplastic rates. |

`EaV, EdVC, delsC` |
Vcmax temperature response parameters |

`EaJ, EdVJ, delsJ` |
Jmax temperature response parameters |

`Km, GammaStar` |
Optionally, provide Michaelis-Menten coefficient for Farquhar model, and Gammastar. If not provided, they are calculated with a built-in function of leaf temperature. |

`id` |
Names of variables (quoted, can be a vector) in the original dataset to be stored in the result. Most useful when using |

`x` |
For plot.acifit, an object returned by |

`what` |
The default is to plot both the data and the model fit, or specify 'data' or 'model' to plot one of them, or 'none' for neither (only the plot region is set up) |

`xlim` |
Limits for the X axis, if left blank estimated from data |

`ylim` |
Limits for the Y axis, if left blank estimated from data |

`whichA` |
By default all three photosynthetic rates are plotted (Aj=Jmax-limited (blue), Ac=Vcmax-limited (red), Hyperbolic minimum (black)). Or, specify one or two of them. |

`add` |
If TRUE, adds to the current plot |

`pch` |
The plotting symbol for the data |

`addzeroline` |
If TRUE, the default, adds a dashed line at y=0 |

`addlegend` |
If TRUE, adds a legend (by default does not add a legend if add=TRUE) |

`legendbty` |
Box type for the legend, passed to argument bty in |

`transitionpoint` |
For plot.acifit, whether to plot a symbol at the transition point. |

`linecols` |
Vector of three colours for the lines (limiting rate, Ac, Aj), if one value provided it is used for all three. |

`lwd` |
Line widths, can be a vector of length 2 (first element for Ac and Aj, second one for the limiting rate). |

`...` |
Further arguments (ignored at the moment). |

### Details

**Fitting method - **
The default method to fit A-Ci curves (set by `fitmethod="default"`

) uses non-linear regression to fit the A-Ci curve. No assumptions are made on which part of the curve is Vcmax or Jmax limited. Normally, all three parameters are estimated: Jmax, Vcmax and Rd, unless Rd is provided as measured (when `useRd=TRUE`

, and Rd is contained in the data). This is the method as described by Duursma (2015, Plos One).

The 'bilinear' method to fit A-Ci curves (set by `fitmethod="bilinear"`

) linearizes the Vcmax and Jmax-limited regions, and applies linear regression twice to estimate first Vcmax and Rd, and then Jmax (using Rd estimated from the Vcmax-limited region). The transition point is found as the one which gives the best overall fit to the data (i.e. all possible transitions are tried out, similar to Gu et al. 2010). The advantage of this method is that it *always* returns parameter estimates, so it should be used in cases where the default method fails. Be aware, though, that the default method fails mostly when the curve contains bad data (so check your data before believing the fitted parameters).

When `citransition`

is set, it splits the data into a Vcmax-limited (where Ci < citransition), and Jmax-limited region (Ci > citransition). Both parameters are then estimated separately for each region (Rd is estimated only for the Vcmax-limited region). **Note** that the actual transition point as shown in the standard plot of the fitted A-Ci curve may be quite different from that provided, since the fitting method simply decides which part of the dataset to use for which limitation, it does not constrain the actual estimated transition point directly. See the example below. If `fitmethod="default"`

, it applies non-linear regression to both parts of the data, and when fitmethod="bilinear", it uses linear regression on the linearized photosynthesis rate. Results will differ between the two methods (slightly).

**TPU limitation - **
Optionally, the `fitaci`

function estimates the triose-phosphate utilization (TPU) rate. The TPU can act as another limitation on photosynthesis, and can be recognized by a 'flattening out' of the A-Ci curve at high Ci. When `fitTPU=TRUE`

, the fitting method used will always be 'bilinear'. The TPU is estimated by trying out whether the fit improves when the last n points of the curve are TPU-limited (where n=1,2,...). When TPU is estimated, it is possible (though rare) that no points are Jmax-limited (in which case estimated Jmax will be NA). A minimum of two points is always reserved for the estimate of Vcmax and Rd. See examples.

**Temperature correction - **
When `Tcorrect=TRUE`

(the default), Jmax and Vcmax are re-scaled to 25C, using the temperature response parameters provided (but Rd is always at measurement temperature). When `Tcorrect=FALSE`

, estimates of all parameters are at measurement temperature. If TPU is fit, it is never corrected for temperature. Important parameters to the fit are GammaStar and Km, both of which are calculated from leaf temperature using standard formulations. Alternatively, they can be provided as known inputs.

**Mesophyll conductance - **
It is possible to provide an estimate of the mesophyll conductance as input (`gmeso`

), in which case the fitted Vcmax and Jmax are to be interpreted as chloroplastic rates. When using gmeso, it is recommended to use the 'default' fitting method (which will use the Ethier&Livingston equations inside `Photosyn`

). It is also implemented with the 'bilinear' method but it requires more testing (and seems to give some strange results). When gmeso is set to a relatively low value, the resulting fit may be quite strange.

**Other parameters - **
The A-Ci curve parameters depend on the values of a number of other parameters. For Jmax, PPFD is needed in order to express it as the asymptote. If PPFD is not provided in the dataset, it is assumed to equal 1800 mu mol m-2 s-1 (in which case a warning is printed). It is possible to either provide PPFD as a variable in the dataset (with the default name 'PARi', which can be changed), or as an argument to the `fitaci`

directly.

**Plotting and summarizing - **
The default **plot** of the fit is constructed with `plot.acifit`

, see Examples below. When plotting the fit, the A-Ci curve is simulated using the `Aci`

function, with leaf temperature (Tleaf) and PPFD set to the mean value for the dataset. The **coefficients** estimated in the fit (Vcmax, Jmax, and usually Rd) are extracted with `coef`

. The summary of the fit is the same as the 'print' method, that is `myfit`

will give the same output as `summary(myfit)`

(where `myfit`

is an object returned by `fitaci`

).

Because fitaci returns the fitted `nls`

object, more details on statistics of the fit can be extracted with standard tools. The Examples below shows the use of the nlstools to extract many details of the fit at once. The fit also includes the **root mean squared error** (RMSE), which can be extracted as `myfit$RMSE`

. This is a useful metric to compare the different fitting methods.

**Atmospheric pressure correction - **
Note that atmospheric pressure (Patm) is taken into account, assuming the original data are in molar units (Ci in mu mol mol-1, or ppm). During the fit, Ci is converted to mu bar, and Km and Gammastar are recalculated accounting for the effects of Patm on the partial pressure of oxygen. When plotting the fit, though, molar units are shown on the X-axis. Thus, you should get (nearly) the same fitted curve when Patm was set to a value lower than 100kPa, but the fitted Vcmax and Jmax will be higher. This is because at low Patm, photosynthetic capacity has to be higher to achieve the same measured photosynthesis rate.

### References

Duursma, R.A., 2015. Plantecophys - An R Package for Analysing and Modelling Leaf Gas Exchange Data. PLoS ONE 10, e0143346. doi:10.1371/journal.pone.0143346

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 | ```
# Fit an A-Ci curve on a dataframe that contains Ci, Photo and optionally Tleaf and PPFD.
# Here, we use the built-in example dataset 'acidata1'.
f <- fitaci(acidata1)
# Note that the default behaviour is to correct Vcmax and Jmax for temperature,
# so the estimated values are at 25C. To turn this off:
f2 <- fitaci(acidata1, Tcorrect=FALSE)
# To use different T response parameters (see ?Photosyn),
f3 <- fitaci(acidata1, Tcorrect=TRUE, EaV=25000)
# Make a standard plot
plot(f)
# Look at a summary of the fit
summary(f)
# Extract coefficients only
coef(f)
# The object 'f' also contains the original data with predictions.
# Here, Amodel are the modelled (fitted) values, Ameas are the measured values.
with(f$df, plot(Amodel, Ameas))
abline(0,1)
# The fitted values can also be extracted with the fitted() function:
fitted(f)
# The non-linear regression (nls) fit is stored as well,
summary(f$nlsfit)
# Many more details can be extracted with the nlstools package
library(nlstools)
overview(f$nlsfit)
# The curve generator is stored as f$Photosyn:
# Calculate photosynthesis at some value for Ci, using estimated parameters and mean Tleaf,
# PPFD for the dataset.
f$Photosyn(Ci=820)
# Photosynthetic rate at the transition point:
f$Photosyn(Ci=f$Ci_transition)$ALEAF
# Set the transition point; this will fit Vcmax and Jmax separately. Note that the *actual*
# transition is quite different from that provided, this is perfectly fine :
# in this case Jmax is estimated from the latter 3 points only (Ci>800), but the actual
# transition point is at ca. 400ppm.
g <- fitaci(acidata1, citransition=800)
plot(g)
g$Ci_transition
# Use measured Rd instead of estimating it from the A-Ci curve.
# The Rd measurement must be added to the dataset used in fitting,
# and you must set useRd=TRUE.
acidata1$Rd <- 2
f2 <- fitaci(acidata1, useRd=TRUE)
f2
# Fit TPU limitation
ftpu <- fitaci(acidata1, fitTPU=TRUE, PPFD=1800, Tcorrect=TRUE)
plot(ftpu)
``` |