Analyzing C3 A-Ci Curves

Overview

In this vignette, we will give an example showing how to analyze C₃ A-Ci data using the PhotoGEA package. The commands in this vignette can be used to initialize your own script, as described in Customizing Your Script.

Background

Understanding C3 A-Ci Curves

An A-Ci curve (or CO₂ response curve) is a particular type of gas exchange measurement where a living leaf is exposed to varying concentrations of CO₂. For each CO₂ concentration in the sequence, the net assimilation rate $A_n$, stomatal conductance to H₂O ($g_{sw}$), and other important quantities are measured and recorded. Typically, other environmental variables such as temperature, humidity, and incident photosynthetically-active photon flux density are held constant during the measurement sequence so that changes in photosynthesis can be attributed to CO₂ alone.

Because of their different cell structures and biochemical pathways, C₃ and C₄ plants have very different responses to CO₂. Here, we will only be discussing C₃ plants.

The full C₃ photosynthetic pathway is quite complicated, consisting of at least two hundred individual reactions, each of which may have an impact on a measured A-Ci curve. However, simplified models for photosynthesis are available and are much easier to understand and work with. These models tend to be based around Rubisco kinetics, describing how the net assimilation rate responds to the amount of CO₂ available in the chloroplast. The most widely-used model is typically referred to as the Farquhar-von-Caemmerer-Berry model (often shortened to the FvCB model or Farquhar model) after the three scientists who originally developed it. An excellent description of this model can be found in Biochemical Models of Leaf Photosynthesis (Caemmerer 2000), subject to a few additional considerations discussed in Lochocki and McGrath (in preparation).

This model provides a framework for understanding the changes in $A_n$ that occur as a C₃ plant is exposed to successively higher concentrations of CO₂. Overall, the photosynthetic response to CO₂ under high light conditions can be divided into three separate ranges:

At low levels of CO₂ in the chloroplast, CO₂ assimilation is primarily limited by Rubisco activity.
At intermediate levels of CO₂ in the chloroplast, CO₂ assimilation is primarily limited by RuBP regeneration.
At high levels of CO₂ in the chloroplast, CO₂ assimilation is primarily limited by triose phosphate utilization.

More specifically, the model provides equations for potential RuBP carboxylation rates that are each limited by one of three important processes:

The Rubisco-limited carboxylation rate ($W_c$).
The RuBP-regeneration-limited carboxylation rate ($W_j$).
The triose-phosphate-utilization-limited carboxylation rate ($W_p$).

Each of these rates depend on the CO₂ concentration in the chloroplast and other factors such as the incident light intensity. The actual RuBP carboxylation rate $V_c$ for a particular set of conditions is taken to be the smallest of the three potential rates. Then, the net CO₂ assimilation rate $A_n$ can be calculated from $V_c$ by including carbon losses due to photorespiration and day respiration. It is also possible to calculate the potential assimilation rates corresponding to each potential carboxylation rate: $A_c$, $A_j$, and $A_p$. Using the equations that compose this model, it is possible to calculate plots like the one below, which illustrates the three different assimilation ranges observed in a typical high-light A-Ci curve.

(Note: this figure was generated using the calculate_c3_assimilation function from the PhotoGEA package, and it represents the photosynethetic response of a C₃ leaf according to the FvCB model with Tp = 18 micromol m^(-2) s^(-1), J = 200 micromol m^(-2) s^(-1), Rd = 1.8 micromol m^(-2) s^(-1), Vcmax = 125 micromol m^(-2) s^(-1), and a leaf temperature of 30 degrees C. Arrhenius temperature response parameters were taken from Sharkey et al. (2007).)

Thus, one of the most common reasons to measure an A-Ci curve is to interpret it in the context of this model. In other words, by fitting the model’s equations to a measured curve, it is possible to estimate values for Tp, J, Rd, and Vcmax. See the documentation for calculate_c3_assimilation for more information about these important quantities.

One complication related to fitting is that the FvCB model calculates assimilation at the level of a single chloroplast; however, experimental gas exchange measurements are made by averaging the behavior of many different chloroplasts across a piece of leaf, which may each have different concentrations of CO₂. Thus, the sharp transitions between Rubisco-limited, RuBP-regeneration-limited, and triose-phosphate-utilization-limited assimilation predicted by the model are rarely observed in practice, and there are instead regions of co-limited assimilation. This can be represented mathematically via two empirical “curvature” parameters, as discussed in the documentation for calculate_c3_assimilation. The following figure shows a CO₂ response curve calculated with both curvatures set to 0.99 and other parameters identical to the ones used to generate the previous figure.

Note that An calculated with curvature < 1 is generally lower than An with curvature = 1 due to the presence of co-limitation; this can be seen by comparing this figure to the previous one (where curvature = 1).

Practicalities

There are a few important practicalities to keep in mind when thinking about CO₂ response curves.

One point is that photosynthesis models generally predict the response of assimilation to the CO₂ concentration in the chloroplast ($C_c$), but gas exchange measurements can only determine the CO₂ concentration in the leaf’s intercellular spaces ($C_i$). Thus, an extra step is required when interpreting A-Ci curves. If the mesophyll conductance $g_{mc}$ to CO₂ is known, then it is possible to calculate values of $C_c$ from $A_n$, $C_i$, and $g_{mc}$. Otherwise, it is typical to assume an infinite mesophyll conductance; in this case, $C_c = C_i$, and the estimated values of Vcmax and other parameters can be considered to be “effective values” describing the plant’s response to intercellular CO₂ rather than the true enzyme response to chloroplastic CO₂.

Another important point is that plants generally do not appreciate being starved of CO₂, so it is not usually possible to start a response curve at low CO₂ and proceed upwards. A more typical approach is to:

Begin at ambient atmospheric CO₂ levels.
Decrease towards a low value.
Return to ambient levels and wait for the plant to reacclimate; this waiting period is usually accomplished by logging several points at ambient CO₂ levels.
Increase to higher values.

When taking this approach, it therefore becomes necessary to remove the extra points measured at ambient CO₂ levels and to reorder the points according to their CO₂ values before plotting or analyzing them.

The Data

A-Ci curves are commonly measured using a Licor portable photosynthesis system such as the Li-6400 or Li-6800. These machines record values of $A_n$, $g_{sw}$, and many other important quantities. They produce two types of output files: plain-text and Microsoft Excel. It is often more convenient to work with the Excel files since the entries can be easily modified (for example, to remove an extraneous row or add a new column). On the other hand, it can be more difficult to access the file contents using other pieces of software such as R. However, the PhotoGEA package reduces this barrier by including tools for reading Licor Excel files in R, which will be demonstrated in the following section.

Loading Packages

As always, the first step is to load the packages we will be using. In addition to PhotoGEA, we will also use the lattice package for generating plots.

# Load required packages
library(PhotoGEA)
library(lattice)

If the lattice package is not installed on your R setup, you can install it by typing install.packages('lattice').

Loading Licor Data

The PhotoGEA package includes two files representing A-Ci curves measured using two Li-6800 instruments. The data is stored in Microsoft Excel files, and includes curves measured from two different crop species (tobacco and soybean) and several different plots of each. Each curve is a sixteen-point CO₂ response curve; in other words, the CO₂ concentration in the air surrounding the leaf was varied, and $A_n$ (among other variables) was measured at each CO₂ setpoint. Although these two files are based on real data, noise was added to it since it is unpublished, so these files should only be used as examples.

The files will be stored on your computer somewhere in your R package installation directory, and full paths to these files can be obtained with PhotoGEA_example_file_path:

# Define a vector of paths to the files we wish to load; in this case, we are
# loading example files included with the PhotoGEA package
file_paths <- c(
  PhotoGEA_example_file_path('c3_aci_1.xlsx'),
  PhotoGEA_example_file_path('c3_aci_2.xlsx')
)

(Note: When loading your own files for analysis, it is not advisable to use PhotoGEA_example_file_path as we have done here. Instead, file paths can be directly written, or files can be chosen using an interactive window. See Input Files below for more information.)

To actually read the data in the files and store them in R objects, we will use the read_gasex_file function from PhotoGEA. Since there are multiple files to read, we will call this function once for each file using lapply:

# Load each file, storing the result in a list
licor_exdf_list <- lapply(file_paths, function(fpath) {
  read_gasex_file(fpath, 'time')
})

The result from this command is an R list of “extended data frames” (abbreviated as exdf objects). The exdf class is a special data structure defined by the PhotoGEA package. In many ways, an exdf object is equivalent to a data frame, with the major difference being that an exdf object includes the units of each column. For more information, type ?exdf in the R terminal to access the built-in help menu entry, or check out the Working With Extended Data Frames vignette.

Generally, it is more convenient to work with a single exdf object rather than a list of them, so our next step will be to combine the objects in the list. This action can be accomplished using the rbind function, which combines table-like objects by their rows; in other words, it stacks two or more tables vertically. This action only makes sense if the tables have the same columns, so before we combine the exdf objects, we should make sure this is the case.

The PhotoGEA package includes a function called identify_common_columns that can be used to get the names of all columns that are present in all of the Licor files. Then, we can extract just those columns, and then combine the exdf objects into a single one.

# Get the names of all columns that are present in all of the Licor files
columns_to_keep <- do.call(identify_common_columns, licor_exdf_list)

# Extract just these columns
licor_exdf_list <- lapply(licor_exdf_list, function(x) {
  x[ , columns_to_keep, TRUE]
})

# Use `rbind` to combine all the data
licor_data <- do.call(rbind, licor_exdf_list)

Now we have a single R object called licor_data that includes all the data from several Licor Excel files. For more information about consolidating information from multiple files, see the Common Patterns section of the Working With Extended Data Frames vignette.

Validating the Data

Before attempting to fit the curves, it is a good idea to do some basic checks of the data to ensure it is organized properly and that it was measured properly.

Basic Checks

First, we should make sure there is a column in the data whose value uniquely identifies each curve. In this particular data set, several “user constants” were defined while making the measurements that help to identify each curve: instrument, species, and plot. However, neither of these columns alone are sufficient to uniquely identify each curve. We can solve this issue by creating a new column that combines the values from each of these:

# Create a new identifier column formatted like `instrument - species - plot`
licor_data[ , 'curve_identifier'] <-
  paste(licor_data[ , 'instrument'], '-', licor_data[ , 'species'], '-', licor_data[ , 'plot'])

The next step is to make sure that this column correctly identifies each response curve. To do this, we can use the check_response_curve_data function from PhotoGEA. Here we will supply the name of a column that should uniquely identify each response curve (curve_identifier), the expected number of points in each curve (16), the name of a “driving” column that should follow the same sequence in each curve (CO2_r_sp). If the data passes the checks, this function will have no output and will not produce any messages. (For more information, see the built-in help menu entry by typing ?check_response_curve_data.)

# Make sure the data meets basic requirements
check_response_curve_data(licor_data, 'curve_identifier', 16, 'CO2_r_sp')

However, if check_response_curve_data detects an issue, it will print a helpful message to the R terminal. For example, if we had specified the wrong number of points or the wrong identifier column, we would get error messages:

check_response_curve_data(licor_data, 'curve_identifier', 15)
#>       curve_identifier npts
#> 1 ripe1 - soybean - 5a   16
#> 2  ripe1 - tobacco - 1   16
#> 3  ripe1 - tobacco - 2   16
#> 4  ripe2 - soybean - 1   16
#> 5 ripe2 - soybean - 5b   16
#> 6  ripe2 - tobacco - 4   16
#> Error in check_response_curve_data(licor_data, "curve_identifier", 15): One or more curves does not have the expected number of points.


check_response_curve_data(licor_data, 'species', 16)
#>   species npts
#> 1 soybean   48
#> 2 tobacco   48
#> Error in check_response_curve_data(licor_data, "species", 16): One or more curves does not have the expected number of points.


check_response_curve_data(licor_data, 'curve_identifier', 16, 'Ci')
#> Error in check_response_curve_data(licor_data, "curve_identifier", 16, : Point 1 from curve `ripe1 - tobacco - 1` has value `Ci = 251.525996944688`, but the average value for this point across all curves is `Ci = 237.742020686325`
#>   Point 1 from curve `ripe1 - tobacco - 2` has value `Ci = 261.420818849899`, but the average value for this point across all curves is `Ci = 237.742020686325`
#>   Point 1 from curve `ripe2 - tobacco - 4` has value `Ci = 268.023119478639`, but the average value for this point across all curves is `Ci = 237.742020686325`

Plotting the A-Ci Curves

One qualitative way to check the data is to simply create a plot of the A-Ci curves. In this situation, the lattice library makes it simple to include each curve as its own separate subplot of a figure. For example:

# Plot all A-Ci curves in the data set
xyplot(
  A ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste('Net CO2 assimilation rate [', licor_data$units$A, ']')
)

Whoops! Why do these curves look so strange? Well, some of the issues are related to the sequence of CO₂ values that was used when measuring the curves. As discussed in Practicalities, there are several repeated points logged at the same CO₂ concentration, and the points are not logged in order of ascending or descending concentration. In fact, the sequence of CO₂ setpoints is as follows:

licor_data[licor_data[, 'curve_identifier'] == 'ripe2 - soybean - 1', 'CO2_r_sp']
#>  [1]  400  300  200  150  100   75   50   20  400  400  600  800 1000 1200 1500
#> [16] 1800

Ideally, we would like to remove the ninth and tenth points (where the setpoint has been reset to 400 to allow the leaf to reacclimate to ambient CO₂ levels), and reorder the data so it is arranged from low to high values of Ci. This can be done using the organize_response_curve function from PhotoGEA:

# Remove points with duplicated `CO2_r_sp` values and order by `Ci`
licor_data <- organize_response_curve_data(
    licor_data,
    'curve_identifier',
    c(9, 10),
    'Ci'
)

Now we can plot them again:

# Plot all A-Ci curves in the data set
xyplot(
  A ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste('Net CO2 assimilation rate [', licor_data$units$A, ']')
)

They still look a bit strange. Some of this is related to the noise that was intentionally added to the data. Nevertheless, there are a few points that we should probably exclude before attempting to fit the curves. One issue is that the model never predicts a decrease in $A$ when $Ci$ increases. So, it is usually a good idea to exclude any points at high $Ci$ where $A$ is observed to decrease. This dropoff in assimilation is due to one or more processes that are not captured by the model, and the fits will be unreliable. We will remove these points later, after we have made several other data quality checks.

Additional Plots for Qualitative Validation

Sometimes a Licor will override the temperature or humidity controls while making measurements; in this case, conditions inside the measurement chamber may not be stable, and we may wish to exclude some of these points. We can check for these types of issues by making more plots. In the following sections, we will generate several different plots to check each curve for quality.

Humidity Control

# Make a plot to check humidity control
xyplot(
  RHcham + `Humidifier_%` + `Desiccant_%` ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  ylim = c(0, 100),
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']')
)

Here, Humidifier_% and Desiccant_% represent the flow from the humidifier and desiccant columns, where a value of 0 indicates that the valve to the column is fully closed and a value of 100 indicates that the valve to the column is fully opened. RHcham represents the relative humidity inside the chamber as a percentage (in other words, as a value between 0 and 100).

When these curves were measured, a chamber humidity setpoint was specified. So, when looking at this plot, we should check that the relative humidity is fairly constant during each curve. Typically, this should be accompanied by relatively smooth changes in the valve percentages as they accomodate changes in ambient humidity and leaf photosynthesis. In this plot, all the data looks good.

Temperature Control

# Make a plot to check temperature control
xyplot(
  TleafCnd + Txchg ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  ylim = c(25, 40),
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste0('Temperature (', licor_data$units$TleafCnd, ')')
)

Here, TleafCnd is the leaf temperature measured using a thermocouple, and Txchg is the temperature of the heat exhanger that is used to control the air temperature in the measurement instrument. When these curves were measured, an exchanger setpoint was specified. So, when looking at this plot, we should check that Txchg is constant during each curve and that the leaf temperature does not vary in an erratic way. In this plot, all the data looks good.

CO₂ Control

# Make a plot to check CO2 control
xyplot(
  CO2_s + CO2_r + CO2_r_sp ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste0('CO2 concentration (', licor_data$units$CO2_r, ')')
)

Here, CO2_s is the CO₂ concentration in the sample cell, CO2_r is the CO₂ concentration in the reference cell, and CO2_r_sp is the setpoint for CO2_r. When these curves were measured, a sequence of CO2_r values was specified, so, when looking at this plot, we should check that CO2_r is close to CO2_r_sp. We also expect that CO2_s should be a bit lower than CO2_r because the leaf in the sample chamber is assimilating CO₂, which should reduce its concentration in the surrounding air. (An exception to this rule occurs at very low values of CO2_r_sp, since in this case there is not enough carbon available to assimilate, and the leaf actually releases CO₂ due to respiration.) In this plot, all the data looks good.

Stability

# Make a plot to check stability criteria
xyplot(
  `A:OK` + `gsw:OK` + Stable ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']')
)

When measuring response curves with a Licor, it is possible to specify stability criteria for each point in addition to minimum and maximum wait times. In other words, once the set point for the driving variable is changed, the machine waits until the stability criteria are met; there is a minimum waiting period, and also a maximum to prevent the machine from waiting for too long.

When these curves were measured, stability criteria were supplied for the net assimilation rate A and the stomatal conductance gsw. The stability status for each was stored in the log file because the appropriate logging option for stability was set. Now, for each point, it is possible to check whether stability was achieved or whether the point was logged because the maximum waiting period had been met. If the maximum waiting period is reached and the plant has still not stabilized, the data point may be unreliable, so it can be helpful to check this information.

In the plot, A:OK indicates whether A was stable (0 for no, 1 for yes), gsw:OK indicates whether gsw was stable (0 for no, 1 for yes), and Stable indicates the total number of stability conditions that were met. So, we are looking for points where Stable is 2. Otherwise, we can check the other traces to see whether A or gsw was unstable.

Comparing these plots with the ones in Plotting the A-Ci Curves, it seems that the unstable points correspond with some of the “odd-looking” points in the A-Ci curves, so it is probably a good idea to remove them before fitting the data.

Cleaning the Licor Data

While checking over the plots in the previous sections, two issues were noticed: (1) some points were logged before stability was achieved and (2) some of the curves exhibit a decrease in $A$ at high values of $C_i$. In this section, we will demonstrate how to remove the unstable and unusual points.

The following command will keep only the points where Stable is exactly 2; this condition means that all of the stability criteria were satisfied. Sometimes, following this procedure, a curve will have very few stable points remaining; it is usually a good idea to automatically exclude any curve with fewer than three stable points.

# Only keep points where stability was achieved
licor_data <- licor_data[licor_data[, 'Stable'] == 2, , TRUE]

# Remove any curves that have fewer than three remaining points
npts <- by(licor_data, licor_data[, 'curve_identifier'], nrow)
ids_to_keep <- names(npts[npts > 2])
licor_data <- licor_data[licor_data[, 'curve_identifier'] %in% ids_to_keep, , TRUE]

Next, we can use the remove_points function from PhotoGEA to exclude the points where there is a decrease in $A$ at high $C_i$. It just so happens that all of these points were measured by the ripe1 instrument and occur at the highest CO₂ setpoint value, so it is easy to specify them all at once:

# Remove points where `instrument` is `ripe1` and `CO2_r_sp` is 1800
licor_data <- remove_points(
  licor_data,
  list(instrument = 'ripe1', CO2_r_sp = 1800)
)

Fitting Licor Data

Now that we have checked the data quality, we are ready to perform the fitting. In order to fit the curves, there are several required pieces of information that are not included in the Licor data files as produced by the instrument: values of the total pressure, values of the chloroplast CO₂ concentration $C_c$, and temperature-dependent values of important photosynthetic parameters such as $\Gamma^*$. However, the PhotoGEA package includes three functions to help with these calculations: calculate_total_pressure, apply_gm, and calculate_arrhenius. Each of these requires an exdf object containing Licor data. The units for each required column will be checked in an attempt to avoid unit-related errors. More information about these functions can be obtained from the built-in help system by typing ?calculate_total_pressure, ?apply_gm, or ?calculate_arrhenius.

Another important consideration is that apply_gm requires values of the mesophyll conductance to CO₂. Values of mesophyll conductance can be specified using the set_variable function from PhotoGEA. Here we will set gmc = 1.0 mol / m^2 / s / bar for tobacco and gmc = Inf for soybean; in this case, the soybean parameters determined by the fitting process will be “effective” parameters. Then we will calculate values of $C_c$, and use Arrhenius parameters from Sharkey et al. (2007) to calculate the temperature-dependent inputs:

# Calculate total pressure in the Licor chamber
licor_data <- calculate_total_pressure(licor_data)

# Specify separate mesophyll conductance values for each species
licor_data <- set_variable(
  licor_data, 'gmc', 'mol m^(-2) s^(-1) bar^(-1)',
  id_column = 'species',
  value_table = list(tobacco = 1.0, soybean = Inf)
)

# Calculate Cc
licor_data <- apply_gm(licor_data)

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_data <- calculate_arrhenius(licor_data, c3_arrhenius_sharkey)

Together, these functions add several new columns to licor_data, including gmc, ‘Cc’, ‘Gamma_star’, and others. With this information, we are now ready to perform the fitting procedure. For this operation, we can use the fit_c3_aci function from the PhotoGEA package, which fits a single response curve to extract the values of key photosynthetic parameters. To apply this function to each curve in a larger data set and then consolidate the results, we can use it in conjunction with by and consolidate, which are also part of PhotoGEA. (For more information about these functions, see the built-in help menu entries by typing ?fit_c3_aci, ?by.exdf, or ?consolidate, or check out the Common Patterns section of the Working With Extended Data Frames vignette.) Together, these functions will split apart the main data using the curve identifier column we defined before (Basic Checks), fit each A-Ci curve using the FvCB model discussed in Understanding C3 A-Ci Curves, and return the resulting parameters and fits:

# Fit the C3 A-Ci curves
c3_aci_results <- consolidate(by(
  licor_data,                       # The `exdf` object containing the curves
  licor_data[, 'curve_identifier'], # A factor used to split `licor_data` into chunks
  fit_c3_aci,                       # The function to apply to each chunk of `licor_data`
  Ca_atmospheric = 420              # Additional argument passed to `fit_c3_aci`
))

Note that in this command, we specified fixed = c(NA, NA, NA, NA) when calling fit_c3_aci. By doing this, we are choosing to vary all four key photosynthetic parameters: Tp, J, Rd, and Vcmax. If fixed is not specified, it will take its default value: c(Inf, NA, NA, NA). In that case, Tp would be fixed to a high value that essentially disables triose-phosphate-utilization-limited assimilation, and only J, Rd, and Vcmax would be varied.

Viewing the Fitted Curves

Having made the fits, it is now a good idea to visually check them, making sure they look reasonable. As before, we can use xyplot from the lattice library to quickly make plots of all the fits, which are stored in c3_aci_results$fits, an exdf object:

# Plot the C3 A-Ci fits
xyplot(
  A + A_fit ~ Ci | curve_identifier,
  data = c3_aci_results$fits$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste0('Intercellular CO2 concentration (', c3_aci_results$fits$units$Ci, ')'),
  ylab = paste0('Net CO2 assimilation rate (', c3_aci_results$fits$units$A, ')')
)

Note that in addition to the fitted values of A, the output also includes values of the potential limiting rates Ac, Aj, and Ap, as well as the estimated operating point. We can take a look at these using a built-in plotting function:

# Plot the C3 A-Cc fits (including limiting rates)
plot_c3_aci_fit(c3_aci_results, 'curve_identifier', 'Ci', ylim = c(-10, 80))

Considering these two plots, we can see that although the model is able to reproduce the observed A values fairly well, some output parameters may be less reliable than others. For example, the triose-phosphate-utilization-limited rate Ap is determined by the Tp parameter; for some of the curves, the assimilation rate is never limited by Ap, so the resulting “best-fit” values of Tp may not be meaningful for these curves.

There is also one curve (ripe2 - soybean - 1) where the optimizer is possibly overestimating J, so Aj is never the limiting factor. To address this, it is possible to force the optimizer to only allow fits where Aj is less than Ac above a certain threshold value of Cc; essentially, this is a way to constrain the crossover from Rubisco-limited carboxylation to RuBP-regeneration-limited carboxylation to a limited range of Cc values. Adding this restriction can sometimes produce more reasonable fits.

# Fit the C3 A-Ci curves
c3_aci_results <- consolidate(by(
  licor_data,                       # The `exdf` object containing the curves
  licor_data[, 'curve_identifier'], # A factor used to split `licor_data` into chunks
  fit_c3_aci,                       # The function to apply to each chunk of `licor_data`
  Ca_atmospheric = 420,             # Additional argument passed to `fit_c3_aci`
  cj_crossover_min = 100,           # Wj must be > Wc when Cc < 100 ppm
  cj_crossover_max = 550            # Wj must be < Wc when Cc > 550 ppm
))

# Plot the C3 A-Cc fits (including limiting rates)
plot_c3_aci_fit(c3_aci_results, 'curve_identifier', 'Ci', ylim = c(-10, 80))

Here we can see that Aj is now significantly lower in the fit for ripe2 - soybean - 1; however, the fitted assimilation values are still never limited by Aj. There may be something pathological about this curve that will prevent us from ever getting a reasonable value for J. In other situations, this approach may be more successful. We can also notice that the fit for ripe1 - soybean - 5a is a little strange. Based on the values of Ci and Cc, it seems that this plant’s stomata did not fully open during the measurement, almost fully restricting the measured assimilation values to the Rubisco-limited range; there is nothing that can be done about this. It is also apparent that the Tp issue remains, and some curves are never triose-phosphate-utilization-limited; again, there is nothing that can really be done to address it besides measuring more A-Ci curves.

Checking the residuals is also a powerful way to gauge the quality of a fit. The output from fit_c3_aci includes the residuals (calculated as A - A_fit) in its output, so it is easy to plot them:

# Plot the residuals
xyplot(
  A_residuals ~ Ci | curve_identifier,
  data = c3_aci_results$fits$main_data,
  type = 'b',
  pch = 16,
  grid = TRUE,
  xlab = paste0('Intercellular CO2 concentration (', c3_aci_results$fits$units$Ci, ')'),
  ylab = paste0('Assimilation rate residual (measured - fitted)\n(', c3_aci_results$fits$units$A, ')')
)

A good fit should produce small and randomly-distributed residuals; here there is no clear pattern to the residuals, so this model is able to reproduce the measured values well.

Examining the Results

Visualizing Average Fit Parameters

Having performed the fits, we can now take a look at the fitted parameter values, which are stored in c3_aci_results$parameters, another exdf object. One way to do this is by using the barchart_with_errorbars function from PhotoGEA to create barcharts of the average values for each species, such as the following example showing the fitted values of Vcmax:

# Make a barchart showing average Vcmax values
barchart_with_errorbars(
  c3_aci_results$parameters[, 'Vcmax_at_25'],
  c3_aci_results$parameters[, 'species'],
  ylim = c(0, 180),
  xlab = 'Species',
  ylab = paste0('Vcmax at 25 degrees C (', c3_aci_results$parameters$units$Vcmax_at_25, ')')
)

Another option is to create box-whisper plots using the bwplot function from the lattice package. Here we use this function to show values of J:

# Make a boxplot showing the distribution of J values
bwplot(
  J_at_25 ~ species,
  data = c3_aci_results$parameters$main_data,
  ylim = c(0, 300),
  xlab = 'Species',
  ylab = paste0('J at 25 degrees C (', c3_aci_results$parameters$units$J_at_25, ')')
)

Accessing Raw and Average Values of Fit Parameters

We can also take a look at the raw numbers of the fitted parameters and their average values. The c3_aci_results$parameters object contains many columns but we only care about a few of them. We can view them as follows:

# View the C3 parameters
columns_for_viewing <-
  c('instrument', 'species', 'plot', 'Vcmax_at_25', 'J_at_25', 'Rd_at_25', 'Tp')

c3_aci_parameters <-
  c3_aci_results$parameters[ , columns_for_viewing, TRUE]

print(c3_aci_parameters)
#>   instrument [UserDefCon] (NA) species [UserDefCon] (NA) plot [UserDefCon] (NA)
#> 1                        ripe1                   soybean                     5a
#> 2                        ripe1                   tobacco                      1
#> 3                        ripe1                   tobacco                      2
#> 4                        ripe2                   soybean                      1
#> 5                        ripe2                   soybean                     5b
#> 6                        ripe2                   tobacco                      4
#>   Vcmax_at_25 [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                         103.76131
#> 2                                         152.22329
#> 3                                         126.80601
#> 4                                         108.91878
#> 5                                         122.01129
#> 6                                          93.87769
#>   J_at_25 [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                      202.7476
#> 2                                      240.2729
#> 3                                      211.6087
#> 4                                            NA
#> 5                                      160.9142
#> 6                                            NA
#>   Rd_at_25 [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                      0.8482917
#> 2                                      2.3279143
#> 3                                      1.9864791
#> 4                                      2.1702846
#> 5                                      2.7885614
#> 6                                      1.5785789
#>   Tp [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                       NA
#> 2                                       NA
#> 3                                       NA
#> 4                                 12.99728
#> 5                                       NA
#> 6                                 14.97514

Alternatively, we can extract average values of the parameters for each species using the basic_stats function from PhotoGEA:

# Compute the average and standard error of each parameter for each species
c3_aci_averages <- basic_stats(c3_aci_results$parameters, 'species')

# View the averages and errors
columns_to_view <- c(
  'species',
  'Vcmax_at_25_avg', 'Vcmax_at_25_stderr',
  'J_at_25_avg', 'J_at_25_stderr',
  'Rd_at_25_avg', 'Rd_at_25_stderr',
  'Tp_avg', 'Tp_stderr'
)
print(c3_aci_averages[ , columns_to_view, TRUE])
#>   species [UserDefCon] (NA)
#> 1                   soybean
#> 2                   tobacco
#>   Vcmax_at_25_avg [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                              111.5638
#> 2                                              124.3023
#>   Vcmax_at_25_stderr [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                                 5.431776
#> 2                                                16.889382
#>   J_at_25_avg [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                                NA
#> 2                                                NA
#>   J_at_25_stderr [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                                   NA
#> 2                                                   NA
#>   Rd_at_25_avg [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                           1.935713
#> 2                                           1.964324
#>   Rd_at_25_stderr [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                             0.5722557
#> 2                                             0.2165979
#>   Tp_avg [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                           NA
#> 2                                           NA
#>   Tp_stderr [fit_c3_aci] (micromol m^(-2) s^(-1))
#> 1                                              NA
#> 2                                              NA

Of course, more advanced statistical tests to check for differences between groups can also be performed within R using other packages such as onewaytests or DescTools. Alternatively, the parameter values can be exported to a comma-separated-value (CSV) file and analyzed in another software environment like jmp.

Customizing Your Script

Note that most of the commands in this vignette have been written in a general way so they can be used as the basis for your own analysis script (see Commands From This Document). In order to use them in your own script, some or all of the following changes may be required. There may also be others not specifically mentioned here.

Input Files

The file paths specified in file_paths will need to be modified so they point to your Licor files. One way to do this in your own script is to simply write out relative or absolute paths to the files you wish to load. For example, you could replace the previous definition of file_paths with this one:

# Define a vector of paths to the files we wish to load
file_paths <- c(
  'myfile1.xlsx',        # `myfile1.xlsx` must be in the current working directory
  'C:/documents/myfile2' # This is an absolute path to `myfile2`
)

You may also want to consider using the choose_input_licor_files function from PhotoGEA; this function will create a pop-up browser window where you can interactively select a set of files. Sometimes this is more convenient than writing out file paths or names. For example, you could replace the previous definition of file_paths with this one:

# Interactively define a vector of paths to the files we wish to load
file_paths <- choose_input_licor_files()

Unfortunately, choose_input_licor_files is only available in interactive R sessions running on Microsoft Windows, but there is also a platform-independent option: choose_input_files. See the Translation section of the Developing a Data Analysis Pipeline vignette for more details.

Curve Identifier

Depending on which user constants are defined in your Licor Excel files, you may need to modify the definition of the curve_identifier column.

Data Cleaning

Depending on the qualitative data checks, you may need to change the input arguments to remove_points. It might also not be necessary to remove the unstable points before performing the fits. Often, it is helpful to not perform any data cleaning at first, and then remove problematic points if they seem to cause problems with the fits.

Tp Limitations

In this example, we consider triose phosphate utilization as a factor that may possibly limit assimilation. However, as discussed in Fitting Licor Data, it is also possible to remove this potential limitation from the model by fixing Tp to a high value and only fitting J, Rd, and Vcmax. Often it is wise to first try fitting a set of curves without Tp limitations, and only consider allowing Tp limitations if the resulting fits are not good.

Averages and Standard Errors

Depending on how your data is organized, you may want to change the column used to divide the data when calculating averages and standard errors.

Plots

You may need to change the axis limits in some or all of the plots. Alternatively, you can remove them, allowing xyplot to automatically choose them for you.

Saving Results

You may want to use write.csv to save some or all of the fitting results as csv files. For example, the following commands will allow you to interactively choose output filenames for the resulting csv files:

write.csv(c3_aci_results$fits, file.choose(), row.names = FALSE)
write.csv(c3_aci_results$parameters, file.choose(), row.names = FALSE)
write.csv(c3_aci_averages, file.choose(), row.names = FALSE)

Commands From This Document

The following code chunk includes all the central commands used throughout this document. They are compiled here to make them easy to copy/paste into a text file to initialize your own script. Annotation has also been added to clearly indicate the four steps involved in data analysis, as described in the Developing a Data Analysis Pipeline vignette.

###
### PRELIMINARIES:
### Loading packages, defining constants, creating helping functions, etc.
###

# Load required packages
library(PhotoGEA)
library(lattice)

###
### TRANSLATION:
### Creating convenient R objects from raw data files
###

## IMPORTANT: When loading your own files, it is not advised to use
## `PhotoGEA_example_file_path` as in the code below. Instead, write out the
## names or use the `choose_input_licor_files` function.

# Define a vector of paths to the files we wish to load; in this case, we are
# loading example files included with the PhotoGEA package
file_paths <- c(
  PhotoGEA_example_file_path('c3_aci_1.xlsx'),
  PhotoGEA_example_file_path('c3_aci_2.xlsx')
)

# Load each file, storing the result in a list
licor_exdf_list <- lapply(file_paths, function(fpath) {
  read_gasex_file(fpath, 'time')
})

# Get the names of all columns that are present in all of the Licor files
columns_to_keep <- do.call(identify_common_columns, licor_exdf_list)

# Extract just these columns
licor_exdf_list <- lapply(licor_exdf_list, function(x) {
  x[ , columns_to_keep, TRUE]
})

# Use `rbind` to combine all the data
licor_data <- do.call(rbind, licor_exdf_list)

###
### VALIDATION:
### Organizing the data, checking its consistency and quality, cleaning it
###

# Create a new identifier column formatted like `instrument - species - plot`
licor_data[ , 'curve_identifier'] <-
  paste(licor_data[ , 'instrument'], '-', licor_data[ , 'species'], '-', licor_data[ , 'plot'])

# Make sure the data meets basic requirements
check_response_curve_data(licor_data, 'curve_identifier', 16, 'CO2_r_sp')

# Remove points with duplicated `CO2_r_sp` values and order by `Ci`
licor_data <- organize_response_curve_data(
    licor_data,
    'curve_identifier',
    c(9, 10),
    'Ci'
)

# Plot all A-Ci curves in the data set
xyplot(
  A ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste('Net CO2 assimilation rate [', licor_data$units$A, ']')
)

# Make a plot to check humidity control
xyplot(
  RHcham + `Humidifier_%` + `Desiccant_%` ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  ylim = c(0, 100),
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']')
)

# Make a plot to check temperature control
xyplot(
  TleafCnd + Txchg ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  ylim = c(25, 40),
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste0('Temperature (', licor_data$units$TleafCnd, ')')
)

# Make a plot to check CO2 control
xyplot(
  CO2_s + CO2_r + CO2_r_sp ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']'),
  ylab = paste0('CO2 concentration (', licor_data$units$CO2_r, ')')
)

# Make a plot to check stability criteria
xyplot(
  `A:OK` + `gsw:OK` + Stable ~ Ci | curve_identifier,
  data = licor_data$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste('Intercellular CO2 concentration [', licor_data$units$Ci, ']')
)

# Only keep points where stability was achieved
licor_data <- licor_data[licor_data[, 'Stable'] == 2, , TRUE]

# Remove any curves that have fewer than three remaining points
npts <- by(licor_data, licor_data[, 'curve_identifier'], nrow)
ids_to_keep <- names(npts[npts > 2])
licor_data <- licor_data[licor_data[, 'curve_identifier'] %in% ids_to_keep, , TRUE]

# Remove points where `instrument` is `ripe1` and `CO2_r_sp` is 1800
licor_data <- remove_points(
  licor_data,
  list(instrument = 'ripe1', CO2_r_sp = 1800)
)

###
### PROCESSING:
### Extracting new pieces of information from the data
###

# Calculate total pressure in the Licor chamber
licor_data <- calculate_total_pressure(licor_data)

# Specify separate mesophyll conductance values for each species
licor_data <- set_variable(
  licor_data, 'gmc', 'mol m^(-2) s^(-1) bar^(-1)',
  id_column = 'species',
  value_table = list(tobacco = 1.0, soybean = Inf)
)

# Calculate Cc
licor_data <- apply_gm(licor_data)

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_data <- calculate_arrhenius(licor_data, c3_arrhenius_sharkey)

# Fit the C3 A-Ci curves
c3_aci_results <- consolidate(by(
  licor_data,                       # The `exdf` object containing the curves
  licor_data[, 'curve_identifier'], # A factor used to split `licor_data` into chunks
  fit_c3_aci,                       # The function to apply to each chunk of `licor_data`
  Ca_atmospheric = 420              # Additional argument passed to `fit_c3_aci`
))

# Plot the C3 A-Ci fits
xyplot(
  A + A_fit ~ Ci | curve_identifier,
  data = c3_aci_results$fits$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  grid = TRUE,
  xlab = paste0('Intercellular CO2 concentration (', c3_aci_results$fits$units$Ci, ')'),
  ylab = paste0('Net CO2 assimilation rate (', c3_aci_results$fits$units$A, ')')
)

# Plot the residuals
xyplot(
  A_residuals ~ Ci | curve_identifier,
  data = c3_aci_results$fits$main_data,
  type = 'b',
  pch = 16,
  grid = TRUE,
  xlab = paste0('Intercellular CO2 concentration (', c3_aci_results$fits$units$Ci, ')'),
  ylab = paste0('Assimilation rate residual (measured - fitted)\n(', c3_aci_results$fits$units$A, ')')
)

###
### SYNTHESIS:
### Using plots and statistics to help draw conclusions from the data
###

# Make a barchart showing average Vcmax values
barchart_with_errorbars(
  c3_aci_results$parameters[, 'Vcmax_at_25'],
  c3_aci_results$parameters[, 'species'],
  ylim = c(0, 180),
  xlab = 'Species',
  ylab = paste0('Vcmax at 25 degrees C (', c3_aci_results$parameters$units$Vcmax_at_25, ')')
)

# Make a boxplot showing the distribution of J values
bwplot(
  J_at_25 ~ species,
  data = c3_aci_results$parameters$main_data,
  ylim = c(0, 300),
  xlab = 'Species',
  ylab = paste0('J at 25 degrees C (', c3_aci_results$parameters$units$J_at_25, ')')
)

# View the C3 parameters
columns_for_viewing <-
  c('instrument', 'species', 'plot', 'Vcmax_at_25', 'J_at_25', 'Rd_at_25', 'Tp')

c3_aci_parameters <-
  c3_aci_results$parameters[ , columns_for_viewing, TRUE]

print(c3_aci_parameters)

# Compute the average and standard error of each parameter for each species
c3_aci_averages <- basic_stats(c3_aci_results$parameters, 'species')

# View the averages and errors
columns_to_view <- c(
  'species',
  'Vcmax_at_25_avg', 'Vcmax_at_25_stderr',
  'J_at_25_avg', 'J_at_25_stderr',
  'Rd_at_25_avg', 'Rd_at_25_stderr',
  'Tp_avg', 'Tp_stderr'
)
print(c3_aci_averages[ , columns_to_view, TRUE])

References

Caemmerer, S. von. 2000. Biochemical Models of Leaf Photosynthesis. CSIRO Publishing. https://doi.org/10.1071/9780643103405.

Sharkey, Thomas D., Carl J. Bernacchi, Graham D. Farquhar, and Eric L. Singsaas. 2007. “Fitting Photosynthetic Carbon Dioxide Response Curves for C3 Leaves.” Plant, Cell & Environment 30 (9): 1035–40. https://doi.org/10.1111/j.1365-3040.2007.01710.x.