Version:	1.3.3
Date:	2025-06-11
Title:	Photosynthetic Gas Exchange Analysis
Description:	Read, process, fit, and analyze photosynthetic gas exchange measurements. Documentation is provided by several vignettes; also see Lochocki, Salesse-Smith, & McGrath (2025) <doi:10.1111/pce.15501>.
Depends:	R (≥ 3.6.0)
Imports:	openxlsx, lattice, dfoptim, DEoptim
Suggests:	knitr, rmarkdown, plantecophys, testthat (≥ 3.0.0)
VignetteBuilder:	knitr
License:	MIT + file LICENSE
Encoding:	UTF-8
LazyData:	true
URL:	https://github.com/eloch216/PhotoGEA, https://eloch216.github.io/PhotoGEA/
Config/testthat/edition:	3
NeedsCompilation:	no
Packaged:	2025-06-12 12:32:44 UTC; eloch
Author:	Edward B. Lochocki [cre, aut], Coralie E. Salesse-Smith [aut], Justin M. McGrath [aut], PhotoGEA authors [cph]
Maintainer:	Edward B. Lochocki <eloch@illinois.edu>
Repository:	CRAN
Date/Publication:	2025-06-12 14:30:09 UTC

The PhotoGEA R package

Description

PhotoGEA (short for photosynthetic gas exchange analysis) is an R package that provides a suite of tools for loading, processing, and analyzing photosynthetic gas exchange data. See Lochocki, Salesse-Smith, & McGrath (2025) [doi:10.1111/pce.15501] for more information.

The best way to learn about using PhotoGEA is to visit the PhotoGEA website and click the Get Started link in the top menu bar. The website includes documentation for all the functions and data sets included in the package, as well as articles that describe its general features and several important use cases.

Locate a PhotoGEA example file on your computer

Description

A convenience function that locates examples files included with the PhotoGEA package (see example_data_files). This function is intended for use in PhotoGEA examples and documentation, and users should not need to use it in their own analysis scripts.

Usage

  PhotoGEA_example_file_path(example_file_name)

Arguments

Details

The PhotoGEA package includes several instrument log files to use in examples and other documentation. A full list can be found in the article about example_data_files. When PhotoGEA is installed, these example files will be stored locally in the R package directory (in the PhotoGEA/extdata subdirectory), which will generally have a different path on every computer. The PhotoGEA_example_file_path function simply locates one of these files and returns its full file path.

When loading your own files for analysis, this function should not be used. Instead, either:

1. Directly write absolute file paths

2. Directly write relative file paths

3. Use one of the convenience functions from PhotoGEA to select files via a pop-up window, such as choose_input_licor_files

When directly writing relative file paths, consider using the file.path function from base R, which will ensure that the paths are properly formatted on any operating system. For example, instead of writing 'Documents\file.xlsx', write file.path('Documents', 'file.xlsx'). Doing this will make it easier to share your analysis scripts with other people who may be using different operating systems.

Value

A full path to a PhotoGEA example file.

Examples

PhotoGEA_example_file_path('c3_aci_1.xlsx')

Calculate CO2 concentration in the chloroplast or mesophyll

Description

Calculates CO2 concentration in the chloroplast or mesophyll, the CO2 drawdown across the stomata, and the CO2 drawdown across the mesophyll. This function can accomodate alternative column names for the variables taken from the Licor file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  apply_gm(
    exdf_obj,
    gmc_at_25 = '',
    photosynthesis_type = 'C3',
    calculate_drawdown = TRUE,
    a_column_name = 'A',
    ca_column_name = 'Ca',
    ci_column_name = 'Ci',
    gmc_norm_column_name = 'gmc_norm',
    total_pressure_column_name = 'total_pressure',
    perform_checks = TRUE,
    return_exdf = TRUE
  )

Arguments

Details

For a C3 plant, the mesophyll conductance to CO2 (gmc) is said to be the conductance satisfying the following one-dimensional flux-conductance equation:

(1) An = gmc * (PCi - PCc)

where An is the net CO2 assimilation rate, PCi is the partial pressure of CO2 in the intercellular spaces, and PCc is the partial pressure of CO2 in the chloroplast. A key underlying assumption for this equation is that the flow of CO2 has reached a steady state; in this case, the flow across the stomata is equal to the flow across the mesophyll.

This equation can be rearranged to calculate PCc:

(2) PCc = PCi - An / gmc

This version of the equation can be found in many places, for example, as Equation 4 in Sharkey et al. "Fitting photosynthetic carbon dioxide response curves for C3 leaves" Plant, Cell & Environment 30, 1035–1040 (2007) [doi:10.1111/j.1365-3040.2007.01710.x].

It is common to express the partial pressures in microbar and the assimilation rate in micromol m^(-2) s^(-1); in this case, the units of mesophyll conductance become mol m^(-2) s^(-1) bar^(-1).

Licor measurement systems provide CO2 levels as relative concentrations with units of parts per million (ppm), or equivalently, micromol mol^(-1). Concentrations and partial pressures are related by the total gas pressure according to:

(3) partial_pressure = total_pressure * relative_concentration

Thus, it is also possible to calculate the CO2 concentration in the choloroplast (Cc) using the following equation:

(4) Cc = Ci - An / (gmc * P)

where Ci is the intercellular CO2 concentration and P is the total pressure. In this function, Equation (4) is used to calculate Cc, where the total pressure is given by the sum of the atmospheric pressure and the chamber overpressure.

When a plant is photosynthesizing, it draws CO2 into its chloroplasts, and this flow is driven by a concentration gradient. In other words, as CO2 flows from the ambient air across the stomata to the intercellular spaces and then across the mesophyll into the chloroplast, there is a decrease in CO2 concentration at each step. Sometimes it is useful to calculate these changes, which are usually referred to as "CO2 drawdown" values. So, in addition to Ci, this function (optionally) calculates the drawdown of CO2 across the stomata (drawndown_cs = Ca - Ci) and the drawdown of CO2 across the mesophyll (drawdown_cm = Ci - Cc).

_Note_: mesophyll conductance is not specified in typical Licor files, so it usually must be added using set_variable before calling apply_gm.

For a C4 plant, mesophyll conductance instead refers to the conductance associated with the flow of CO2 from the intercellular spaces into the mesophyll (rather than into the chloroplast). In this case, the equations above just require a small modification where Pcc and Cc are replaced by PCm and Cm, the partial pressure and concentration of CO2 in the mesophyll.

Value

The return value depends on the value of return_exdf:

• If return_exdf is TRUE, the return value is an exdf object based on exdf_obj with the following columns, calculated as described above: Pci and Ci (for C3 plants) or PCm and Cm (for C4 plants), drawndown_s, and drawdown_cm. The category for each of these new columns is apply_gm to indicate that they were created using this function.

• If return_exdf is FALSE, the return value is a list with a single named element (internal_c), which contains values of Cc or Cm as a numeric vector.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

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

# Calculate temperature-dependent parameter values, including gmc_norm
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_sharkey)

# Calculate Cc and drawdowns assuming a mesophyll conductance of
# 1 mol / m^2 / s / bar at 25 degrees C
licor_file <- apply_gm(licor_file, 1)

licor_file$units$Cc      # View the units of the new `Cc` column
licor_file$categories$Cc # View the category of the new `Cc` column
licor_file[, 'Cc']       # View the values of the new `Cc` column

Convert an exdf object to a data frame

Description

Converts an exdf object to a data frame by appending the units and categories to the top of each column in the exdf object's main_data data frame. Typically this function is used for displaying the contents of an exdf object; in fact, it is used internally by View, write.csv, and other functions. The main_data of an exdf object x can be accessed directly (without including the units and categories in the first row) via x[['main_data']] as with any other list element.

Usage

  ## S3 method for class 'exdf'
as.data.frame(x, ...)

Arguments

Value

A data frame formed from x.

See Also

exdf

Examples

simple_exdf <- exdf(data.frame(A = 1), data.frame(A = 'u'), data.frame(A = 'c'))
as.data.frame(simple_exdf) # Includes units and categories in the first rows
simple_exdf[['main_data']] # Just returns the main data

Barcharts with error bars

Description

barchart_with_errorbars is a wrapper for lattice::barchart that includes error bars on the chart, while bwplot_wrapper is a simple wrapper for lattice::bwplot that gives it the same function signature as barchart_with_errorbars.

Usage

  barchart_with_errorbars(
    Y,
    X,
    eb_width = 0.2,
    eb_lwd = 1,
    eb_col = 'black',
    na.rm = TRUE,
    remove_outliers = FALSE,
    ...
  )

  bwplot_wrapper(Y, X, ...)

Arguments

Details

The barchart_with_errorbars function uses tapply to calculate the mean and standard error for each subset of Y as determined by the values of X. In other words, means <- tapply(Y, X, mean), and similar for the standard errors. The mean values are represented as bars in the final plot, while the standard error is used to create error bars located at mean +/- standard_error.

The bwplot_wrapper function is a simple wrapper for lattice::bwplot that gives it the same input arguments as barchart_with_errorbars. In other words, the same X and Y vectors can be used to create a barchart using barchart_with_errorbars or a box-whisker plot with bwplot_wrapper.

Value

A trellis object created by lattice::barchart or lattice::bwplot.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

# Plot the average assimilation value for each species. (Note: this is not a
# meaningful calculation since we are combining assimilation values measured
# at different PPFD.)
barchart_with_errorbars(
  licor_file[, 'A'],
  licor_file[, 'species'],
  ylim = c(0, 50),
  xlab = 'Species',
  ylab = paste0('Net assimilation (', licor_file$units$A, ')')
)

# Make a box-whisker plot using the same data. (Note: this is not a meaningful
# plot since we are combining assimilation values measured at different PPFD.)
bwplot_wrapper(
  licor_file[, 'A'],
  licor_file[, 'species'],
  ylim = c(0, 50),
  xlab = 'Species',
  ylab = paste0('Net assimilation (', licor_file$units$A, ')')
)

# Another way to create the plots. This method illustrates the utility of the
# bwplot_wrapper function.
plot_parameters <- list(
  Y = licor_file[, 'A'],
  X = licor_file[, 'species'],
  ylim = c(0, 50),
  xlab = 'Species',
  ylab = paste0('Net assimilation (', licor_file$units$A, ')')
)
do.call(barchart_with_errorbars, plot_parameters)
do.call(bwplot_wrapper, plot_parameters)

Calculate basic stats (mean and standard error)

Description

Calculates basic stats (mean and standard error) for each applicable column in an exdf object split up according to the values of one or more identifier columns.

Usage

  basic_stats(
    exdf_obj,
    identifier_columns,
    na.rm = TRUE
  )

Arguments

Details

This function first splits up exdf_obj into chunks according to the values of the identifier_columns. For each chunk, columns that have a single unique value are identified and excluded from the statistical calculations. For the remaining numeric columns, the mean and standard error are calculated.

Value

An exdf object including the mean and standard error for each applicable column, where each row represents one value of the identifier_columns. The column names are determined by appending '_avg' and '_stderr' to the original names.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

# Calculate the average assimilation and stomatal conductance values for each
# species. (Note: this is not a meaningful calculation!)
basic_stats(
  licor_file[ , c('species', 'K', 'A', 'gsw'), TRUE],
  'species'
)

Apply a function to an exdf object split by one or more factors

Description

Divides an exdf object into groups defined by one or more factors and applies a function to each group.

Usage

  ## S3 method for class 'exdf'
by(data, INDICES, FUN, ...)

Arguments

Value

Splits data into chunks x by the values of the INDICES and calls FUN(x, ...) for each chunk; returns a list where each element is the output from each call to FUN.

See Also

exdf

Examples

# Read a Licor file, split it into chunks according to the `species` column,
# and count the number of measurements for each species
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

by(licor_file, licor_file[, 'species'], nrow)

C3 temperature response parameters from Bernacchi et al.

Description

Parameters describing the temperature response of important C3 photosynthetic parameters, intended to be passed to the calculate_temperature_response function.

Usage

c3_temperature_param_bernacchi

Format

List with 12 named elements that each represent a variable whose temperature-dependent value can be calculated using an Arrhenius equation, Johnson-Eyring-Williams equation, or a polynomial equation:

• Gamma_star_at_25: The value of chloroplastic CO2 concentration at which CO2 gains from Rubisco carboxylation are exactly balanced by CO2 losses from Rubisco oxygenation (Gamma_star) at 25 degrees C.

• Gamma_star_norm: Gamma_star normalized to its value at 25 degrees C.

• gmc_norm: The mesophyll conductance to CO2 diffusion (gmc) normalized to its value at 25 degrees C.

• J_norm: The electron transport rate (J) normalized to its value at 25 degrees C.

• Kc_at_25: The Michaelis-Menten constant for rubisco carboxylation (Kc) at 25 degrees C.

• Kc_norm: Kc normalized to its value at 25 degrees C.

• Ko_at_25: The Michaelis-Menten constant for rubisco oxygenation (Ko) at 25 degrees C.

• Ko_norm: Ko normalized to its value at 25 degrees C.

• RL_norm: The rate of non-photorespiratory CO2 release in the light (RL) normalized to its value at 25 degrees C.

• Tp_norm: The maximum rate of triose phosphate utilization (Tp) normalized to its value at 25 degrees C.

• Vcmax_norm: The maximum rate of rubisco carboxylation (Vcmax) normalized to its value at 25 degrees C.

• Vomax_norm: The maximum rate of rubisco oxygenation (Vomax) normalized to Vcmax at 25 degrees C.

In turn, each of these elements is a list with at least 2 named elements:

• type: the type of temperature response

• units: the units of the corresponding variable.

Source

Many of these parameters are normalized to their values at 25 degrees C. Vomax is normalized to the value of Vcmax at 25 degrees C. These variables include _norm in their names to indicate this.

Arrhenius parameters for J were obtained from Bernacchi et al. (2003). Here, we use the values determined from chlorophyll fluorescence measured from plants grown at 25 degrees C (Table 1). Although Bernacchi et al. (2003) reports values of Jmax, here we assume that both Jmax and the light-dependent values of J follow the same temperature response function and refer to it as J for compatibility with c3_temperature_param_sharkey.

Johnson-Eyring-Williams parameters for gmc were obtained from Bernacchi et al. (2002).

The Bernacchi papers from the early 2000s do not specify a temperature response for Tp, so we instead use the Johnson-Eyring-Williams response from Sharkey et al. (2007). Another option would be to use a flat temperature response; in other words, to assume that Tp is constant with temperature. This could be done with the following code, which takes the flat response parameters from c3_temperature_param_flat: within(c3_temperature_param_bernacchi, {Tp_norm = c3_temperature_param_flat$Tp_norm})

The Arrhenius parameters for the other variables were obtained from Bernacchi et al. (2001).

References:

• Bernacchi, C. J., Singsaas, E. L., Pimentel, C., Jr, A. R. P. & Long, S. P. "Improved temperature response functions for models of Rubisco-limited photosynthesis" Plant, Cell & Environment 24, 253–259 (2001) [doi:10.1111/j.1365-3040.2001.00668.x].

• Bernacchi, C. J., Portis, A. R., Nakano, H., von Caemmerer, S. & Long, S. P. "Temperature Response of Mesophyll Conductance. Implications for the Determination of Rubisco Enzyme Kinetics and for Limitations to Photosynthesis in Vivo" Plant Physiology 130, 1992–1998 (2002) [doi:10.1104/pp.008250].

• Bernacchi, C. J., Pimentel, C. & Long, S. P. "In vivo temperature response functions of parameters required to model RuBP-limited photosynthesis" Plant, Cell & Environment 26, 1419–1430 (2003) [doi:10.1046/j.0016-8025.2003.01050.x].

• Sharkey, T. D., Bernacchi, C. J., Farquhar, G. D. & Singsaas, E. L. "Fitting photosynthetic carbon dioxide response curves for C3 leaves" Plant, Cell & Environment 30, 1035–1040 (2007) [doi:10.1111/j.1365-3040.2007.01710.x].

C3 temperature response parameters for a flat response

Description

Parameters that specify a flat temperature response (in other words, no dependence on temperature) for important C3 photosynthetic parameters, intended to be passed to the calculate_temperature_response function.

Usage

c3_temperature_param_flat

Format

List with 11 named elements that each represent a variable whose temperature-dependent value can be calculated using an Arrhenius equation or a polynomial equation:

• Gamma_star_at_25: The value of chloroplastic CO2 concentration at which CO2 gains from Rubisco carboxylation are exactly balanced by CO2 losses from Rubisco oxygenation (Gamma_star) at 25 degrees C.

• Gamma_star_norm: Gamma_star normalized to its value at 25 degrees C.

• gmc_norm: The mesophyll conductance to CO2 diffusion (gmc) normalized to its value at 25 degrees C.

• J_norm: The electron transport rate (J) normalized to its value at 25 degrees C.

• Kc_at_25: The Michaelis-Menten constant for rubisco carboxylation (Kc) at 25 degrees C.

• Kc_norm: Kc normalized to its value at 25 degrees C.

• Ko_at_25: The Michaelis-Menten constant for rubisco oxygenation (Ko) at 25 degrees C.

• Ko_norm: Ko normalized to its value at 25 degrees C.

• RL_norm: The rate of non-photorespiratory CO2 release in the light (RL) normalized to its value at 25 degrees C.

• Tp_norm: The maximum rate of triose phosphate utilization (Tp) normalized to its value at 25 degrees C.

• Vcmax_norm: The maximum rate of rubisco carboxylation (Vcmax) normalized to its value at 25 degrees C.

In turn, each of these elements is a list with at least 2 named elements:

• type: the type of temperature response

• units: the units of the corresponding variable.

Source

Many of these parameters are normalized to their values at 25 degrees C. These variables include _norm in their names to indicate this.

Here, the activation energy values (Ea) are all set to 0, which means that the values will not depend on temperature. Some parameters are specified at 25 degrees C; these values were obtained from Sharkey et al. (2007). (See c3_temperature_param_sharkey.)

References:

• Sharkey, T. D., Bernacchi, C. J., Farquhar, G. D. & Singsaas, E. L. "Fitting photosynthetic carbon dioxide response curves for C3 leaves" Plant, Cell & Environment 30, 1035–1040 (2007) [doi:10.1111/j.1365-3040.2007.01710.x].

C3 temperature response parameters from Sharkey et al.

Description

Parameters describing the temperature response of important C3 photosynthetic parameters, intended to be passed to the calculate_temperature_response function.

Usage

c3_temperature_param_sharkey

Format

List with 11 named elements that each represent a variable whose temperature-dependent value can be calculated using an Arrhenius equation or a polynomial equation:

• Gamma_star_at_25: The value of chloroplastic CO2 concentration at which CO2 gains from Rubisco carboxylation are exactly balanced by CO2 losses from Rubisco oxygenation (Gamma_star) at 25 degrees C.

• Gamma_star_norm: Gamma_star normalized to its value at 25 degrees C.

• gmc_norm: The mesophyll conductance to CO2 diffusion (gmc) normalized to its value at 25 degrees C.

• J_norm: The electron transport rate (J) normalized to its value at 25 degrees C.

• Kc_at_25: The Michaelis-Menten constant for rubisco carboxylation (Kc) at 25 degrees C.

• Kc_norm: Kc normalized to its value at 25 degrees C.

• Ko_at_25: The Michaelis-Menten constant for rubisco oxygenation (Ko) at 25 degrees C.

• Ko_norm: Ko normalized to its value at 25 degrees C.

• RL_norm: The rate of non-photorespiratory CO2 release in the light (RL) normalized to its value at 25 degrees C.

• Tp_norm: The maximum rate of triose phosphate utilization (Tp) normalized to its value at 25 degrees C.

• Vcmax_norm: The maximum rate of rubisco carboxylation (Vcmax) normalized to its value at 25 degrees C.

In turn, each of these elements is a list with at least 2 named elements:

• type: the type of temperature response

• units: the units of the corresponding variable.

Source

Many of these parameters are normalized to their values at 25 degrees C. These variables include _norm in their names to indicate this.

Response parameters were obtained from Sharkey et al. (2007). In this publication, gas concentrations are expressed as partial pressures (in Pa or kPa) rather than mole fractions (micromol / mol or mmol / mol). However, for consistency with c3_temperature_param_bernacchi, here we prefer to use mole fractions.

To convert a concentration expressed as a partial pressure (P; in Pa) to a concentration expressed as a mole fraction (C; in micromol / mol), we need a value for atmospheric pressure; we will use the typical value of 101325 Pa. Then C = P / 101325 * 1e6 or C = P * cf, where cf = 1e6 / 101325 is a conversion factor. The same correction can be used to convert kPa to mmol / mol. The value of cf can be accessed using PhotoGEA:::c_pa_to_ppm.

References:

• Sharkey, T. D., Bernacchi, C. J., Farquhar, G. D. & Singsaas, E. L. "Fitting photosynthetic carbon dioxide response curves for C3 leaves" Plant, Cell & Environment 30, 1035–1040 (2007) [doi:10.1111/j.1365-3040.2007.01710.x].

C4 temperature response parameters for a flat response

Description

Parameters that specify a flat temperature response (in other words, no dependence on temperature) for important C4 photosynthetic parameters, intended to be passed to the calculate_temperature_response function.

Usage

c4_temperature_param_flat

Format

List with 10 named elements that each represent a variable whose temperature-dependent value can be calculated using either an Arrhenius or Gaussian equation:

• Vcmax_norm: The maximum rate of rubisco carboxylation (Vcmax) normalized to its value at 25 degrees C.

• Vpmax_norm: The maximum rate of PEP carboxylase activity (Vpmax) normalized to its value at 25 degrees C.

• RL_norm: The respiration rate (RL) normalized to the value of Vcmax at 25 degrees C.

• Kc: The Michaelis-Menten constant for rubisco carboxylation.

• Ko: The Michaelis-Menten constant for rubisco oxygenation.

• Kp: The Michaelis-Menten constant of PEP carboxylase.

• gamma_star: Half the reciprocal of rubisco specificity.

• ao: The ratio of solubility and diffusivity of O2 to CO2.

• gmc_norm: The mesophyll conductance to CO2 diffusion normalized to its value at 25 degrees C.

• J_norm: The electron transport rate J normalized to its value at 25 degrees C.

Each of these is a list with 4 named elements:

• type: the type of temperature response ('Arrhenius')

• c: the (dimensionless) Arrhenius scaling factor.

• Ea: the activation energy in kJ / mol.

• units: the units of the corresponding variable.

Source

Some of these parameters (Vcmax, Vpmax, RL, gmc, and J) are normalized to their values at 25 degrees C. These variables include _norm in their names to indicate this.

The remaining parameters (Kc, Ko, Kp, gamma_star, ao, and gmc) are not normalized because they are assumed to not vary significantly between species.

Here, the activation energy values (Ea) are all set to 0, which means that the values will not depend on temperature. The Arrhenius scaling factors c are chosen to reproduce the parameter values at 25 degrees C as specified in von Caemmerer (2021). (See c4_temperature_param_vc.)

References:

• von Caemmerer, S. "Updating the steady-state model of C4 photosynthesis" Journal of Experimental Botany 72, 6003–6017 (2021) [doi:10.1093/jxb/erab266].

C4 temperature response parameters from von Caemmerer

Description

Temperature response parameters describing the temperature response of important C4 photosynthetic parameters, intended to be passed to the calculate_temperature_response function.

Usage

c4_temperature_param_vc

Format

List with 10 named elements that each represent a variable whose temperature-dependent value can be calculated using either an Arrhenius or Gaussian equation:

• Vcmax_norm: The maximum rate of rubisco carboxylation (Vcmax) normalized to its value at 25 degrees C.

• Vpmax_norm: The maximum rate of PEP carboxylase activity (Vpmax) normalized to its value at 25 degrees C.

• RL_norm: The respiration rate (RL) normalized to the value of Vcmax at 25 degrees C.

• Kc: The Michaelis-Menten constant for rubisco carboxylation.

• Ko: The Michaelis-Menten constant for rubisco oxygenation.

• Kp: The Michaelis-Menten constant of PEP carboxylase.

• gamma_star: Half the reciprocal of rubisco specificity.

• ao: The ratio of solubility and diffusivity of O2 to CO2.

• gmc_norm: The mesophyll conductance to CO2 diffusion normalized to its value at 25 degrees C.

• J_norm: The electron transport rate J normalized to its value at 25 degrees C.

The J_norm parameter is calculated using a Gaussian function and hence its corresponding list element is itself a list with 4 named elements:

• type: the type of temperature response ('Gaussian')

• optimum_rate: the largest value this parameter can take.

• t_opt: the temperature where the optimum occurs in degrees C.

• sigma: the width of the Gaussian in degrees C.

• units: the units of the corresponding variable.

Each of the remaining elements is a list with 4 named elements:

• type: the type of temperature response ('Arrhenius')

• c: the (dimensionless) Arrhenius scaling factor.

• Ea: the activation energy in kJ / mol.

• units: the units of the corresponding variable.

Source

Some of these parameters (Vcmax, Vpmax, RL, gmc, and J) are normalized to their values at 25 degrees C. These variables include _norm in their names to indicate this.

The remaining parameters (Kc, Ko, Kp, gamma_star, and ao) are not normalized because they are assumed to not vary significantly between species.

Here, the Arrhenius scaling factors (c; dimensionless) and activation energy values (Ea; kJ / mol) are obtained from von Caemmerer (2021). In that publication, the overall scaling for each parameter is specified by its value at 25 degrees C; the scaling factors are determined from this information as described in the documentation for calculate_temperature_response_arrhenius.

The Gaussian parameters (t_opt and sigma) for J_norm are also obtained from von Caemmerer (2021), assuming that J and Jmax follow the same temperature response. The value of optimum_rate is chosen such that J_norm is equal to 1 at a temperature of 25 degrees C.

References:

• von Caemmerer, S. "Updating the steady-state model of C4 photosynthesis" Journal of Experimental Botany 72, 6003–6017 (2021) [doi:10.1093/jxb/erab266].

Calculate the Ball-Berry index

Description

Calculates the Ball-Berry index. This function can accomodate alternative column names for the variables taken from the Licor file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_ball_berry_index(
    data_table,
    a_column_name = 'A',
    rhleaf_column_name = 'RHleaf',
    csurface_column_name = 'Csurface'
  )

Arguments

Details

The Ball-Berry index is defined as A * h_s / c_s, where A is the net assimilation rate, h_s is the relative humidity at the leaf surface, and c_s is the CO2 concentration at the leaf surface. This variable is a key part of the Ball-Berry model, which assumes that stomatal conductance is linearly related to the Ball-Berry index. For more information, please see the original publication describing the model: Ball, J. T., Woodrow, I. E. and Berry, J. A. "A Model Predicting Stomatal Conductance and its Contribution to the Control of Photosynthesis under Different Environmental Conditions." in "Progress in Photosynthesis Research: Volume 4" (1986) [doi:10.1007/978-94-017-0519-6_48].

Typically, the relative humidity and CO2 concentration at the leaf surface are not included in Licor output files. Instead, the output files only include the relative humidity and CO2 concentration in the sample chamber, and conditions at the leaf surface may be slightly different. These required inputs can be calculated using the calculate_gas_properties function.

Value

An object based on data_table that includes the Ball-Berry index as a new column called bb_index.

If data_table is an exdf object, the category of this new column will be calculate_ball_berry_index to indicate that it was created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package, calculate the
# total pressure, calculate additional gas properties, and finally calculate the
# Ball-Berry index.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_total_pressure(licor_file)

licor_file <- calculate_gas_properties(licor_file)

licor_file <- calculate_ball_berry_index(licor_file)

licor_file$units$bb_index      # View the units of the new `bb_index` column
licor_file$categories$bb_index # View the category of the new `bb_index` column
licor_file[,'bb_index']        # View the values of the new `bb_index` column

Calculate C3 assimilation rates

Description

Calculates C3 assimilation rates based on the Farquhar-von-Caemmerer-Berry model. This function can accomodate alternative colum names for the variables taken from Licor files in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c3_assimilation(
    data_table,
    alpha_g,
    alpha_old,
    alpha_s,
    alpha_t,
    Gamma_star_at_25,
    J_at_25,
    Kc_at_25,
    Ko_at_25,
    RL_at_25,
    Tp_at_25,
    Vcmax_at_25,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    cc_column_name = 'Cc',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    hard_constraints = 0,
    perform_checks = TRUE,
    return_table = TRUE,
    ...
  )

Arguments

Details

The Busch et al. (2018) and Busch (2020) model:

This function generally follows the Farquhar-von-Caemmerer-Berry model as described in Busch et al. (2018) and Busch (2020) with a few modifications described below. In this formulation, the steady-state net CO2 assimilation rate An is calculated according to

An = (1 - Gamma_star_agt / PCc) * Vc - RL,

where Gamma_star is the CO2 compensation point in the absence of non-photorespiratory CO2 release, Gamma_star_agt is the effective value of Gamma_star accounting for glycolate carbon remaining in the cytosol, PCc is the partial pressure of CO2 in the chloroplast, Vc is the RuBP carboxylation rate, and RL is the rate of non-photorespiratory CO2 release in the light. Gamma_star_agt is given by

Gamma_star_agt = (1 - alpha_g + 2 * alpha_t) * Gamma_star,

where alpha_g and alpha_t are the fractions of glycolate carbon leaving the photorespiratory pathway as glycine and CH2-THF, respectively.

The model considers three potential values of Vc that correspond to limitations set by three different processes: Rubisco activity, RuBP regeneration, and triose phopsphate utilization (TPU). The Rubisco-limited carboxylation rate Wc is given by

Wc = PCc * Vcmax / (PCc + Kc * (1.0 + POc / Ko)),

where Vcmax is the maximum rate of Rubisco carboxylation, Kc is the Michaelis-Menten constant for CO2, Ko is the Michaelis-Menten constant for O2, and POc is the partial pressure of O2 in the chloroplast.

The RuBP-regeneration-limited carboxylation rate Wj is given by

Wj = PCc * J / (4 * PCc + Gamma_star_agt * (8 + 16 * alpha_g - 8 * alpha_t + 8 * alpha_s)),

where J is the potential electron transport rate at a given light intensity and alpha_s is the fraction of glycolate carbon leaving the photorespiratory pathway as serine.

The TPU-limited carboxylation rate is given by

Wp = PCc * 3 * Tp / (PCc - Gamma_star_agt * (1 + 3 * alpha_g + 6 * alpha_t + 4 * alpha_s)),

where Tp is the maximum rate of triose phosphate utilization. Note that this equation only applies when PCc > Gamma_star_agt * (1 + 3 * alpha_g + 6 * alpha_t + 4 * alpha_s); for smaller values of PCc, TPU cannot limit the RuBP carboxylation rate and Wp = Inf. (Lochocki & McGrath, under review).

The actual carboxylation rate is typically chosen to be the smallest of the three potential rates:

Vc = min{Wc, Wj, Wp}.

In the equations above, several of the variables depend on the leaf temperature. In particular, the leaf-temperature-adjusted values of Gamma_star, J, Kc, Ko, RL, Tp, and Vcmax are determined from their base values at 25 degrees C and a temperature-dependent multiplicative factor.

Also note that PCc is calculated from the chloroplastic CO2 concentration Cc using the total pressure (ambient pressure + chamber overpressure).

In addition to the carboxylation and assimilation rates already mentioned, it is also possible to calculate the net CO2 assimilation rates determined by Rubisco activity, RuBP regeneration, and TPU as follows:

Ac = (1 - Gamma_star_agt / PCc) * Wc - RL

Aj = (1 - Gamma_star_agt / PCc) * Wj - RL

Ap = (1 - Gamma_star_agt / PCc) * Wp - RL

The Busch model with nitrogen restrictions:

Note that the implementation as described above does not currently facilitate the inclusion of nitrogen limitations (Equations 15-21 in Busch et al. (2018)).

The "old" model:

In an older version of the model, alpha_g, alpha_s, and alpha_t are replaced with a single parameter alpha_old. Most publications refer to this simply as alpha, but here we follow the notation of Busch et al. (2018) for clarity. In this version, there is no disctinction between Gamma_star_agt and Gamma_star. Other differences are described below.

The RuBP-regeneration-limited carboxylation rate Wj is given by

Wj = PCc * J / (Wj_coef_C * PCc + Wj_coef_Gamma_star * Gamma_star),

Here we have allowed Wj_coef_C and Wj_coef_Gamma_star to be variables rather than taking fixed values (as they do in many sources). This is necessary because not all descriptions of the FvCB model use the same values, where the different values are due to different assumptions about the NADPH and ATP requirements of RuBP regeneration.

The TPU-limited carboxylation rate is given by

Wp = PCc * 3 * Tp / (PCc - Gamma_star * (1 + 3 * alpha_old)),

Note that this equation only applies when PCc > Gamma_star * (1 + 3 * alpha_old); for smaller values of PCc, TPU cannot limit the RuBP carboxylation rate and Wp = Inf. (Lochocki & McGrath, under review).

Using either version of the model:

When using calculate_c3_assimilation, it is possible to use either version of the model. Setting alpha_g, alpha_s, and alpha_t to zero is equivalent to using the older version of the model, while setting alpha_old = 0 is equivalent to using the newer version of the model. If all alpha parameters are zero, there is effectively no difference between the two versions of the model. Attempting to set a nonzero alpha_old if either alpha_g, alpha_s, or alpha_t is nonzero is forbidden since it would represent a mix between the two models; if such values are passed as inputs, then an error will be thrown.

Hard constraints:

Most input parameters to the FvCB model have hard constraints on their values which are set by their biochemical or physical interpretation; for example, Vcmax cannot be negative and alpha_g must lie between 0 and 1. Yet, because of measurement noise, sometimes it is necessary to use values outside these ranges when fitting an A-Ci curve with fit_c3_aci or fit_c3_variable_j. To accomodate different potential use cases, it is possible to selectively apply these hard constraints by specifying different values of the hard_constraints input argument:

• hard_constraints = 0: Constraints are only placed on inputs that are user-supplied and cannot be fit, such as oxygen.

• hard_constraints = 1: Includes the same constraints as when hard_constraints is 0, with the additional constraint that all Cc values must be non-negative.

• hard_constraints = 2: Includes the same constraints as when hard_constraints is 1, which additional constraints on the parameters that can be fitted. For example, Vcmax_at_25 must be non-negative and alpha_g must lie between 0 and 1.

If any input values violate any of the specified constraints, an error message will be thrown.

Optional arguments:

• use_min_A: If an input argument called use_min_A is supplied and its value is TRUE, then the "minimum assimilation" variant of the FvCB model will be used. In this case, An will be calculated as An = min{Ac, Aj, Ap}. In general, using this variant is not recommended.It should only be used to investigate errors that may occur when using the minimal assimilation rate rather than the minimal carboxylation rate.

• TPU_threshold: If an input argument called TPU_threshold is supplied and its numeric value is not NULL, then TPU limitations will only be allowed for values of Cc above this threshold. This threshold will be used in place of the values discussed in the equations above. In general, using this option is not recommended. It should only be used to investigate errors that may occur when using a fixed TPU threshold.

• use_FRL: If an input argument called use_FRL is supplied and its value is TRUE, then An will always be set to Ac for Cc < Gamma_star_agt. This "forced Rubisco limitation" can only be used along with the "minimum assimilation" variant (use_min_A = TRUE).

• consider_depletion: If an input argument called consider_depletion is supplied and its value is TRUE, then RuBP depletion will be considered to be an additional potential limiting process. In this case, Vc will be calculated as Vc = min{Wc, Wj, Wp, Wd}, where Wd is zero when Cc < Gamma_star and Inf otherwise. Note that the value of Wd (and Ad = (1 - Gamma_star / PCc) * Wd - RL) will always be returned, regardless of whether RuBP depletion is considered when calculating An.

References:

• Busch, Sage, & Farquhar, G. D. "Plants increase CO2 uptake by assimilating nitrogen via the photorespiratory pathway." Nature Plants 4, 46–54 (2018) [doi:10.1038/s41477-017-0065-x].

• Busch "Photorespiration in the context of Rubisco biochemistry, CO2 diffusion and metabolism." The Plant Journal 101, 919–939 (2020) [doi:10.1111/tpj.14674].

• von Caemmerer, S. "Biochemical Models of Leaf Photosynthesis" (CSIRO Publishing, 2000) [doi:10.1071/9780643103405].

• Lochocki & McGrath "Widely Used Variants of the Farquhar-von-Caemmerer-Berry Model Can Cause Errors in Parameter Estimates and Simulations." submitted.

Value

The return value depends on the value of return_table:

• If return_table is TRUE, the return value is an exdf object with the following columns, calculated as described above: Tp_tl, Vcmax_tl, RL_tl, J_tl, Ac, Aj, Ap, An, Vc, and others. The category for each of these new columns is calculate_c3_assimilation to indicate that they were created using this function.

• If return_table is FALSE, the return value is a list with the following named elements: An, Ac, Aj, Ap, and J_tl. Each element is a numeric vector.

If data_table is not an exdf object, then the return value will be a data frame, and units and categories will not be reported.

Examples

# Simulate a C3 A-Cc curve with specified leaf temperature and photosynthetic
# parameters and plot the net assimilation rate along with the different
# enzyme-limited rates
inputs <- exdf(data.frame(
  Cc = seq(1, 601, by = 6),
  Tleaf = 30,
  total_pressure = 1,
  oxygen = 21
))

inputs <- document_variables(
  inputs,
  c('', 'Cc',             'micromol mol^(-1)'),
  c('', 'Tleaf',          'degrees C'),
  c('', 'total_pressure', 'bar'),
  c('', 'oxygen',         'percent')
)

inputs <- calculate_temperature_response(inputs, c3_temperature_param_sharkey, 'Tleaf')

assim <- calculate_c3_assimilation(inputs, 0, 0, 0, 0, '', 150, '', '', 1, 12, 120)

lattice::xyplot(
  Ac + Aj + Ap + An ~ inputs[, 'Cc'],
  data = assim$main_data,
  type = 'l',
  grid = TRUE,
  auto = TRUE,
  xlab = paste0('Chloroplast CO2 concentration (', inputs$units$Cc, ')'),
  ylab = paste0('Assimilation rate (', assim$units$An, ')')
)

Estimate the relative limiting factors to C3 photosynthesis

Description

Uses the method from Grassi & Magnani (2005) to estimate the relative limitations to C3 photosynthesis due to stomatal conductance, mesophyll conductance, and biochemistry. This function can accomodate alternative column names for the variables taken from the data file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c3_limitations_grassi(
    exdf_obj,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    cc_column_name = 'Cc',
    gamma_star_column_name = 'Gamma_star_tl',
    gmc_column_name = 'gmc_tl',
    gsc_column_name = 'gsc',
    kc_column_name = 'Kc_tl',
    ko_column_name = 'Ko_tl',
    oxygen_column_name = 'oxygen',
    total_pressure_column_name = 'total_pressure',
    vcmax_column_name = 'Vcmax_tl',
    j_column_name = NULL
  )

Arguments

Details

When analyzing or interpreting C3 gas exchange data, it is often useful to estimate the relative limitations to assimilation that are due to stomatal conductance, mesophyll conductance, and biochemistry. This can be done using a framework first introduced by Grassi & Magnani (2005). In this framework, the relative limitation due to stomatal conductance (ls) is

ls = [(g_t / g_sc) * (dAdC)] / [g_t + dAdC],

the relative limitation due to mesophyll conductance (lm) is

lm = [(g_t / g_mc) * (dAdC)] / [g_t + dAdC],

and the relative limitation due to biochemistry (lb) is

ln = [g_t] / [g_t + dAdC],

where g_sc is the stomatal conductance to CO2, g_mc is the mesophyll conductance to CO2, gt = 1 / (1 / g_mc + 1 / g_sc) is the total conductance to CO2, and dAdC is the partial derivative of the net CO2 assimilation rate (An) with respect to the chloroplast CO2 concentration (Cc). These can be found in Equation 7 from Grassi & Magnani (2005).

These equations were derived by assuming that CO2 assimilation is limited by Rubisco activity; in other words, that the net CO2 assimilation rate is given by

Ac = Vcmax * (Cc - Gamma_star) / (Cc + Km) - RL,

where Vcmax is the maximum Rubisco carboxylation rate, Gamma_star is the CO2 compensation point in the absence of day respiration, RL is the day respiration rate, Km is the effective Michaelis-Menten constant for Rubisco carboxylation. In turn, Km is given by Km = Kc * (1 + O / Ko), where Kc is the Michaelis-Menten constant for carboxylation, Ko is the Michaelis-Menten constant for oxygenation, and O is the oxygen concentration in the chloroplast.

Under this assumption, it is possible to analytically determine the partial derivative dAdC:

dAdC_rubisco = Vcmax * (Gamma_star + Km) / (Cc + Km)^2

In this case, the limitation due to "biochemistry" actually refers to limitation due to the value of Vcmax. Note that sometimes this derivative is estimated from the initial slope of a measured A-Ci curve rather than calculated analytically. (See, for example, Pathare et al. (2023).) However, we do not take that approach here. Also note that the value of Vcmax can be estimated using different approaches. For example, Xiong (2023) uses single-point gas exchange measurements. When possible, it would be better to use an estimate from fitting an entire A-Ci curve, as shown in the example below.

To understand the meaning of these limiting factors, note that simultaneously making small fractional increases to g_sc, g_mc, and Vcmax will generally cause an associated small fractional increase in An. The limiting factors describe the fraction of the increase in An that can be attributed to each of g_sc, g_mc, and Vcmax. For example, ls = 0.2, lm = 0.3, lb = 0.5 would mean that 20 percent of the increase in An would be due to an increase in stomatal conductance, 30 percent due to an increase in mesophyll conductance, and 50 percent due to an increase in Vcmax. Note that ls, lm, and lb always add up to 1.

Thus, when one of the factors is large, changes in the related parameter produce relatively larger changes in the assimilation rate. In that case, it can be said that that parameter is setting a large limit on the assimilation rate. On the other hand, if a factor is small, small changes in the related parameter produce relatively small changes in An, and therefore that parameter is not setting a large limit on the assimilation rate.

It is also possible to calculate dAdC when assimilation is limited by RuBP regeneration. In this case, we have

Aj = J * (Cc - Gamma_star) / (4 * Cc + 8 * Gamma_star) - RL,

where J is the RuBP regeneration rate, and the limitation due to "biochemistry" actually refers to limitation due to the value of J (rather than Vcmax. The same equations as before can be used to calculate the limiting factors (ls, lm, lb), but the partial derivative is now given by

dAdC_j = J * Gamma_star * 12 / (4 * Cc + 8 * Gamma_star)^2.

Most users will want the limitations assuming Rubisco-limited assimilation. However, if j_column_name is not NULL, values of J will be used to calculate the limiting factors assuming RuBP-regeneration-limited assimilation. For an example of how these additional factors can be used, see Sakoda et al. (2021).

References:

Grassi, G. & Magnani, F. "Stomatal, mesophyll conductance and biochemical limitations to photosynthesis as affected by drought and leaf ontogeny in ash and oak trees." Plant, Cell & Environment 28, 834–849 (2005) [doi:10.1111/j.1365-3040.2005.01333.x].

Pathare, V. S. et al. "Altered cell wall hydroxycinnamate composition impacts leaf- and canopy-level CO2 uptake and water use in rice." Plant Physiology kiad428 (2023) [doi:10.1093/plphys/kiad428].

Xiong, D. "Leaf anatomy does not explain the large variability of mesophyll conductance across C3 crop species." The Plant Journal 113, 1035–1048 (2023) [doi:10.1111/tpj.16098].

Sakoda, K., Yamori, W., Groszmann, M. & Evans, J. R. "Stomatal, mesophyll conductance, and biochemical limitations to photosynthesis during induction." Plant Physiology 185, 146–160 (2021) [doi:10.1093/plphys/kiaa011].

Value

This function returns an exdf object based on exdf_obj but with several new columns representing the partial derivatives and limiting factors discussed above: dAdC_rubisco, ls_rubisco_grassi, lm_rubisco_grassi, and lb_rubisco_grassi. If j_column_name is not NULL, the output will also include dAdC_j, ls_j_grassi, lm_j_grassi, and lb_j_grassi.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate additional gas properties
licor_file <- calculate_gas_properties(licor_file)

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Fit all curves in the data set. Here we use a faster optimizer than the
# default one to ensure the example runs quickly.
aci_results <- consolidate(by(
  licor_file,
  licor_file[, 'species_plot'],
  fit_c3_aci,
  Ca_atmospheric = 420,
  optim_fun = optimizer_nmkb(1e-7),
  fit_options = list(gmc_at_25 = 0.5)
))

# Get a subset of fitting results corresponding to the first measured point
# in each curve (where CO2_r_sp = 400 ppm)
aci_fit_subset <- aci_results$fits[aci_results$fits[, 'CO2_r_sp'] == 400, , TRUE]

# Calculate limiting factors
aci_fit_subset <- calculate_c3_limitations_grassi(aci_fit_subset)

# View the limiting factors for each species / plot
col_to_keep <- c(
  'species', 'plot',                                            # identifiers
  'ls_rubisco_grassi', 'lm_rubisco_grassi', 'lb_rubisco_grassi' # limitation info
)

aci_fit_subset[ , col_to_keep, TRUE]

# One of these fits has NA for all the limiting factors, which causes problems
# when making bar charts with some versions of the `lattice` package, so we
# exclude that curve for plotting
data_for_barchart <-
  aci_fit_subset$main_data[aci_fit_subset$main_data$species_plot != 'tobacco - 2', ]

# Display as a bar chart
lattice::barchart(
  ls_rubisco_grassi + lm_rubisco_grassi + lb_rubisco_grassi ~ species_plot,
  data = data_for_barchart,
  stack = TRUE,
  auto = TRUE,
  ylab = 'Factors limiting assimilation'
)

Estimate the relative limiting factors to C3 photosynthesis

Description

Uses the method from Warren et al. (2003) to estimate the relative limitations to C3 photosynthesis due to stomatal conductance and mesophyll conductance. This function can accomodate alternative column names for the variables taken from the data file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c3_limitations_warren(
    exdf_obj,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    ca_column_name = 'Ca',
    cc_column_name = 'Cc',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    hard_constraints = 0,
    ...
  )

Arguments

Details

When analyzing or interpreting C3 gas exchange data, it is often useful to estimate the relative limitations to assimilation that are due to stomatal conductance or mesophyll conductance. This can be done using a framework first introduced by Warren et al. (2003). In this framework, the relative limitation due to stomatal conductance (ls) is

ls = (An_inf_gsc - A_modeled) / An_inf_gsc

and the relative limitation due to mesophyll conductance (lm) is

lm = (An_inf_gmc - A_modeled) / An_inf_gmc. These are equations 10 and 11 in Warren et al. (2003).

In these equations A_modeled is the net assimilation rate calculated using the Farquhar-von-Caemmerer-Berry (FvCB) model at the measured value of the chloroplast CO2 concentration (Cc). The other two assimilation rates (An_inf_gsc and An_inf_gmc) are also calculated using the FvCB model, but under different assumptions: An_inf_gsc assumes that stomatal conductance is infinite while mesophyll conductance is as measured, while An_inf_gmc assumes that mesophyll conductance is infinite while stomatal conductance is as measured.

In other words, ls expresses the observed assimilation rate as a fractional decrease relative to a hypothetical plant with infinite stomatal conductance, while lm expresses the observed assimilation rate as a fractional decrease relative to a hypothetical plant with infinite mesophyll conductance.

For example, if lm = 0.4, this means that the observed assimilation rate is 40 percet lower than a hypothetical plant with infinite mesophyll conductance. If mesophyll conductance were to increase (all else remaining the same), then lm would decrease. This is not the case with other estimations of limiting factors, such as the one used in calculate_c3_limitations_grassi. (See Leverett & Kromdijk for more details.)

To actually calculate An_inf_gsc and An_inf_gmc, it is first necessary to estimate the corresponding values of Cc that would occur with infinite stomatal or mesophyll conductance. This can be done with a 1D diffusion equation expressed using drawdown values:

Cc = Ca - drawdown_cs - drawdown_cm,

where drawdown_cs = Ca - Ci is the drawdown of CO2 across the stomata (assuming infinite boundary layer conductance) and drawdown_cm = Ci - Cc is the drawdown of CO2 across the mesophyll. If one conductance is infinite, the corresponding drawdown becomes zero. Thus, we have:

Cc_inf_gsc = Ca - 0 - (Ci - Cc) = Ca - Ci + Cc

and

Cc_inf_gmc = Ca - (Ca - Ci) - 0 = Ci,

where Cc_inf_gsc is the value of Cc that would occur with infinite stomatal conductance and the measured mesophyll conductance, and Cc_inf_gmc is the value of Cc that would occur with infinite mesophyll conductance and the measured stomatal conductance.

Once values of Cc, Cc_inf_gsc, and Cc_inf_gmc, the corresponding assimilation rates are calculated using calculate_c3_assimilation, and then the limitation factors are calculated as described above.

References:

Warren, C. R. et al. "Transfer conductance in second growth Douglas-fir (Pseudotsuga menziesii (Mirb.)Franco) canopies." Plant, Cell & Environment 26, 1215–1227 (2003) [doi:10.1046/j.1365-3040.2003.01044.x].

Leverett, A. & Kromdijk, J. "The long and tortuous path towards improving photosynthesis by engineering elevated mesophyll conductance." [doi:10.22541/au.170016201.13513761/v1].

Value

This function returns an exdf object based on exdf_obj but with several new columns representing the quantities discussed above: Cc_inf_gsc, Cc_inf_gmc, An_inf_gsc, An_inf_gmc, ls_warren, and lm_warren.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate additional gas properties
licor_file <- calculate_gas_properties(licor_file)

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Fit all curves in the data set. Here we use a faster optimizer than the
# default one to ensure the example runs quickly.
aci_results <- consolidate(by(
  licor_file,
  licor_file[, 'species_plot'],
  fit_c3_aci,
  Ca_atmospheric = 420,
  optim_fun = optimizer_nmkb(1e-7)
))

# Get a subset of fitting results corresponding to the first measured point
# in each curve (where CO2_r_sp = 400 ppm)
aci_fit_subset <- aci_results$fits[aci_results$fits[, 'CO2_r_sp'] == 400, , TRUE]

# Calculate limiting factors
aci_fit_subset <- calculate_c3_limitations_warren(aci_fit_subset)

# View the limiting factors for each species / plot
col_to_keep <- c(
  'species', 'plot',       # identifiers
  'ls_warren', 'lm_warren' # limitation info
)

aci_fit_subset[ , col_to_keep, TRUE]

Calculate C3 variable J

Description

Calculates values of mesophyll conductance and chloroplast CO2 concentration using the "variable J" equation, as originally described in Harley et al. (1992) and modified in Moualeu-Ngangue, Chen, & Stutzel (2016). This function can accomodate alternative colum names for the variables taken from Licor files in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c3_variable_j(
    exdf_obj,
    alpha_g,
    alpha_s,
    alpha_t,
    Gamma_star_at_25,
    RL_at_25,
    tau,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    phips2_column_name = 'PhiPS2',
    qin_column_name = 'Qin',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    hard_constraints = 0,
    perform_checks = TRUE,
    return_exdf = TRUE
  )

Arguments

Details

The "Variable J" method is a way to estimate the chloroplast CO2 concentration Cc and the mesophyll conductance to CO2 gmc from combined gas exchange and chlorophyll fluorescence measurements, and was originally described in Harley et al. (1992). The main idea is that along with Cc, the net CO2 assimilation rate (An), day respiration rate (RL), and CO2 compensation point in the absence of day respiration (Gamma_star) determine the actual RuBP regeneration rate (J_actual) required to support the Calvin-Benson cycle:

J_actual = (A + RL) * (4 * Cc + 8 * Gamma_star) / (Cc - Gamma_star)

This is Equation 6 in Harley et al. (1992). (Note: this equation can be derived by solving the equation for Aj from the FvCB model for J. However, this relationship holds true even when CO2 assimilation is not limited by RuBP regeneration. Hence, we distinguish between the actual regeneration rate J_actual and the maximum regeneration rate for a given incident light level J.)

This equation can be rewritten by using a 1D diffusion equation to replace Cc with Cc = Ci - An / gmc and then solving for the mesophyll conductance. The result is Equation 7 in Harley et al. (1992), which we do not reproduce here. The importance of Equation 7 is that it calculates gmc from several quantities that can be measured using gas exchange (Ci, An, and RL), a quantity whose values can be known beforehand (Gamma_star), and J_actual (which can be estimated from chlorophyll fluorescence measurements). Here we update Equation 7 to include alpha_g and alpha_s following Busch et al. (2018) (also see calculate_c3_assimilation.)

The actual RuBP regeneration rate is related to the incident photosynthetically active flux density Qin and the operating efficiency of photosystem II PhiPSII according to:

J_actual = alpha_g * beta * Qin * PhiPSII,

where alpha_g is the leaf absorptance and beta is the fraction of absorbed light energy directed to photosystem II. Qin is set by the measurement conditions, while PhiPSII can be estimated from chlorophyll fluorescence. However, the values of alpha_g and beta are generally unknown; beta in particular is difficult or impossible to measure and is often assumed to be 0.5. Thus, while Equation 7 from Harley et al. (1992) can be used to estimate gmc, there is a practical uncertainty associated with determining a value of J_actual to be used in Equation 7.

Moualeu-Ngangue, Chen, & Stutzel (2016) developed a way to address this issue. The method from that paper replaces the product of alpha_g and beta by a single new parameter tau, and uses it to estimate the actual RuBP regeneration from fluoresence (J_F):

J_F = tau * Qin * PhiPSII.

This new parameter tau is assumed to be constant across an A-Ci curve, and is treated as an unknown whose value will be determined during a fitting procedure.

In this function, the supplied values of Qin, PhiPSII, and tau are used to calculate values of J_F. Then, the values of J_F are used along with Equation 7 from Harley et al. (1992) to calculate gmc. Finally, a 1D diffusion equation is used to calculate Cc.

Hard constraints:

Most input parameters to the Variable J equations have hard constraints on their values which are set by their biochemical or physical interpretation; for example, RL cannot be negative and tau must lie between 0 and 1. Yet, because of measurement noise, sometimes it is necessary to use values outside these ranges when fitting an A-Ci curve with fit_c3_variable_j. To accomodate different potential use cases, it is possible to selectively apply these hard constraints by specifying different values of the hard_constraints input argument:

• hard_constraints = 0: Constraints are only placed on inputs that are user-supplied and cannot be fit, such as Qin.

• hard_constraints = 1: Includes the same constraints as when hard_constraints is 0, with the additional constraint that all Ci values must be non-negative.

• hard_constraints = 2: Includes the same constraints as when hard_constraints is 1, which additional constraints on the parameters that can be fitted. For example, RL_at_25 must be non-negative and tau must lie between 0 and 1.

If any input values violate any of the specified constraints, an error message will be thrown.

References:

• Harley, P. C., Loreto, F., Di Marco, G. & Sharkey, T. D. "Theoretical Considerations when Estimating the Mesophyll Conductance to CO2 Flux by Analysis of the Response of Photosynthesis to CO2" Plant Physiology 98, 1429–1436 (1992) [doi:10.1104/pp.98.4.1429].

• Moualeu-Ngangue, D. P., Chen, T.-W. & Stutzel, H. "A new method to estimate photosynthetic parameters through net assimilation rate-intercellular space CO2 concentration (A-Ci) curve and chlorophyll fluorescence measurements" New Phytologist 213, 1543–1554 (2017) [doi:10.1111/nph.14260].

• Busch, Sage, & Farquhar, G. D. "Plants increase CO2 uptake by assimilating nitrogen via the photorespiratory pathway." Nature Plants 4, 46–54 (2018) [doi:10.1038/s41477-017-0065-x].

Value

The return value depends on the value of return_exdf:

• If return_exdf is TRUE, the return value is an exdf object with the following columns, calculated as described above: J_F, gmc, Cc, tau, and RL_tl. The category for each of these new columns is calculate_c3_variable_j to indicate that they were created using this function.

• If return_exdf is FALSE, the return value is a list with the following named elements: gmc, Cc, and J_F. Each element is a numeric vector.

Examples

# Read an example Licor file included in the PhotoGEA package. This file
# includes gas exchange and chlorophyll fluorescence data.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'])

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Calculate values of J_F, gmc, and Cc assuming alpha_g = alpha_s = alpha_t = 0,
# RL_at_25 = 1.5, and tau = 0.55.
vj_res <- calculate_c3_variable_j(licor_file, 0, 0, 0, '', 1.5, 0.55)

# Plot mesophyll conductance against Cc. Note: this information is not very
# meaningful since the values of Gamma_star, tau and RL used above are
# arbitrary.
lattice::xyplot(
  gmc ~ Cc | licor_file[, 'species_plot'],
  data = vj_res$main_data,
  type = 'b',
  pch = 16,
  auto = TRUE,
  xlab = paste0('Chloroplast CO2 concentration (', vj_res$units$Cc, ')'),
  ylab = paste0('Mesophyll conductance to CO2 (', vj_res$units$gmc, ')')
)

Calculate C4 assimilation rates

Description

Calculates C4 assimilation rates based on the von Caemmerer (2000) model. This function can accomodate alternative colum names for the variables taken from Licor files in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c4_assimilation(
    exdf_obj,
    alpha_psii,
    gbs,
    J_at_25,
    RL_at_25,
    Rm_frac,
    Vcmax_at_25,
    Vpmax_at_25,
    Vpr,
    x_etr = 0.4,
    ao_column_name = 'ao',
    gamma_star_column_name = 'gamma_star',
    j_norm_column_name = 'J_norm',
    kc_column_name = 'Kc',
    ko_column_name = 'Ko',
    kp_column_name = 'Kp',
    oxygen_column_name = 'oxygen',
    pcm_column_name = 'PCm',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    vcmax_norm_column_name = 'Vcmax_norm',
    vpmax_norm_column_name = 'Vpmax_norm',
    hard_constraints = 0,
    perform_checks = TRUE,
    return_exdf = TRUE
  )

Arguments

Details

General Description of the Model

This function generally follows Sections 4.2.1 and 4.2.2 from S. von Caemmerer (2000), which provides equations for calculating the enzyme-limited net assimilation rate Ac, the light- and electron-transport limited rate Aj, and the overall net assimilation rate An in a C4 leaf. (These equations are also reproduced in S. von Caemmerer (2021), although we use the equation numbers from the 2000 textbook here. Also note there is a typo in Equation 22 from the 2021 paper.) The enzyme-limited assimilation rate in this model is calculated according to Equation 4.21:

Ac = (-b - sqrt(b^2 - 4 * a * c)) / (2 * a)

where the parameters a, b, and c are determined by Equations 4.22, 4.23, and 4.24, respectively. These equations are fairly long, so we do not reproduce them here. Similarly, the light-limited rate Aj is also calculated according to a quadratic equation. Finally, the overall rate is calculated as the smaller of Ac and Aj:

An = min(Ac, Aj)

An Approximation to the Full Equations

The complicated equations above can be approximiated by simpler ones. For Ac, we can use Equation 4.25:

Ac = min(Vp + gbs * PCm - RLm, Vcmax - RL)

where Vp is the rate of PEP carboxylation, gbs is the bundle sheath conductance to CO2, PCm is the partial pressure of CO2 in the mesophyll, RLm is the rate of mitochondrial respiration occuring in the mesophyll, Vcmax is the maximum rate of Rubisco carboxylation, and RL is the rate of mitochondrial respiration occurring in the bundle sheath and mesophyll. Essentially, the first term in the equation above (Vp + gbs * PCm - RLm) can be thought of as a PEP-carboxylase-limited assimilation rate Ap, while the second term (Vcmax - RL) is a Rubisco-limited rate Ar.

The PEP carboxylation rate Vp is calculated according to Equation 4.19:

Vp = min(Pcm * Vpmax / (PCm + Kp), Vpr)

where Vpmax is the maximum rate of PEP carboxylation, Kp is a Michaelis-Menten constant for PEP carboxylation, and Vpr is the carboxylation rate when PEP carboxylase activity is limited by regeneration rather than carbon availability. Thus, we can see that the approximation above actually calculates the enzyme-limited rate as the smaller of three separate assimilation rates:

Ac = min(Apc, Apr, Ar)

where Apc = Pcm * Vpmax / (PCm + Kp) + gbs * PCm - RLm is the rate due to carbon-limited PEP carboxylation, Apr = Vpr + gbs * PCm - RLm is the rate due to regeneration-limited PEP carboxylation, and Ar = Vcmax - RL is the rate due to Rubisco-limited assimilation.

In the example at the end of this documentation page, we compare Apc, Apr, and Ar to Ac as calculated by Equation 4.21. From this example, it is clear that the approximation Ac = min(Apc, Apr, Ar) is quite accurate for low values of PCm, but introduces significant errors as PCm increases. Thus, while the approximation can be helpful for gaining an intuitive understanding of C4 photosynthesis, it should not be used for realistic calculations.

To be more precise, the approximation is only reliable when Vcmax is much larger than gbs * Kc * (1 + POm / Ko), which is rarely the case; otherwise, the limiting value of An at high PCm will be smaller than Ar = Vcmax - RL. Conversely, if gbs and alpha_psii are both set to zero, then the approximation is exact.

For Aj, the simplified version is Equation 4.45:

Aj = min(x_etr * J / 2 - RLm + gbs * PCm, (1 - x_etr) * J / 3 - RL)

where x_etr is the fraction of whole-chain electron transport occurring in the mesophyll and J is the electron transport rate. We can therefore think of this equation as

Aj = min(Ajm, Ajbs)

where Ajm is the mesophyll light-limited rate and Ajbs is the bundle sheath light-limited rate. These are given by Ajm = x_etr * J / 2 - RLm + gbs * PCm and (1 - x_etr) * J / 3 - RL As in the case with Ac, this approximation is not exact.

Combining these two simplifications, we can see that the overall net assimilation rate can be approximated as the smallest of five potential rates:

An = min(Apc, Apr, Ar, Ajm, Ajbs).

Here it is very important to note that some of these potential rates have identical or similar dependence on PCm. More specifically, Apr and Ajm have identical dependence, as do Ar and Ajbs. If gbs is zero, all four of these rates have no dependence on PCm. Thus, from a fitting point of view, it is not usually possible to distinguish between these potential limiting states. For this reason, it is not advisable to fit more than one of Vcmax, Vpr, and Jmax when estimating parameters from an experimentally measured curve.

Limiting Cases of the Approximate Equation

The bundle sheath conductance gbs is generally very small and can be ignored in a simple analysis of the above equations. In that case, when Pcm is very high, the approximate equation for Ac simplifies further to:

Ac = min(Vpmax - RLm, Vpr - RLm, Vcmax - RL)

Since respiration costs are also generally much smaller than the maximum enzyme activity and regeneration rates, the enzyme-limited assimilation rate at high levels of CO2 is therefore determined by the smaller of Vpmax, Vpr, and Vcmax. As shown in Table 4.1 of the textbook, Vpmax is typically much larger than the other two rates, so light- and CO2-saturated assimilation in C4 leaves is usually limited by either Vpr or Vcmax. The exact limiting factor can depend on many possible variables, such as the temperature. For example, see Wang (2008).

At lower values of PCm, enzyme-limited net assimilation is determined by CO2-limited PEP carboxylation according to:

An = PCm * Vpmax / Kp - RLm

where we have approximated gbs * PCm = 0 and PCm + Kp = Kp, as appropriate for small values of Pcm. Thus, we can see that for low CO2 levels, assimilation is linearly related to PCm with a slope of Vpmax / Kp and intercept of -RLm.

Respiration

Table 4.1 from von Caemmerer (2000) suggests that RL = 0.01 * Vcmax and RLm = 0.5 * RL. To allow more flexibility, we allow RL to be specified independently of Vcmax, and we also consider the ratio of RLm / RL = Rm_frac to be a variable (so that RLm is calculated from RL according to RLm = Rm_frac * RL). If Rm_frac is set to 1, then there is no distinction between RL and RLm.

Hard constraints:

Most input parameters to the C4 assimilation model have hard constraints on their values which are set by their biochemical or physical interpretation; for example, Vcmax cannot be negative and alpha_psii must lie between 0 and 1. Yet, because of measurement noise, sometimes it is necessary to use values outside these ranges when fitting an A-Ci curve with fit_c4_aci. To accomodate different potential use cases, it is possible to selectively apply these hard constraints by specifying different values of the hard_constraints input argument:

• hard_constraints = 0: Constraints are only placed on inputs that are user-supplied and cannot be fit, such as Kc.

• hard_constraints = 1: Includes the same constraints as when hard_constraints is 0, with the additional constraint that all PCm values must be non-negative.

• hard_constraints = 2: Includes the same constraints as when hard_constraints is 1, which additional constraints on the parameters that can be fitted. For example, Vcmax_at_25 must be non-negative and alpha_psii must lie between 0 and 1.

If any input values violate any of the specified constraints, an error message will be thrown.

References

• von Caemmerer, S. "Biochemical Models of Leaf Photosynthesis" (CSIRO Publishing, 2000) [doi:10.1071/9780643103405].

• von Caemmerer, S. "Updating the steady-state model of C4 photosynthesis." Journal of Experimental Botany 72, 6003–6017 (2021) [doi:10.1093/jxb/erab266].

• Wang, D., Portis, A. R., Jr., Moose, S. P. & Long, S. P. "Cool C4 Photosynthesis: Pyruvate Pi Dikinase Expression and Activity Corresponds to the Exceptional Cold Tolerance of Carbon Assimilation in Miscanthus × giganteus." Plant Physiology 148, 557–567 (2008) [doi:10.1104/pp.108.120709].

Value

The return value depends on the value of return_exdf:

• If return_exdf is TRUE, the return value is an exdf object with the following columns: alpha_psii, gbs, J_at_25, Jmax_tl, J_tl, Rm_frac, Vcmax_tl, Vpmax_tl, RL_tl, RLm_tl, Vpc, Vpr, Vp, Apc, Apr, Ap, Ar, Ajm, Ajbs, Ac, Aj, An, and c4_assimilation_msg. Most of these are calculated as described above, while several are copies of the input arguments with the same name. The c4_assimilation_msg is usually blank but may contain information about any issues with the inputs. The category for each of these new columns is calculate_c4_assimilation to indicate that they were created using this function.

• If return_exdf is FALSE, the return value is a numeric vector containing the calculated values of An.

Examples

# Simulate a C4 A-Cm curve with specified leaf temperature and photosynthetic
# parameters and plot the net assimilation rate.
npts <- 101

inputs <- exdf(data.frame(
  PCm = seq(0, 500, length.out = npts),
  Tleaf = 25,
  Qin = 1800,
  total_pressure = 1,
  oxygen = 21
))

inputs <- document_variables(
  inputs,
  c('', 'PCm',            'microbar'),
  c('', 'Tleaf',          'degrees C'),
  c('', 'Qin',            'micromol m^(-2) s^(-1)'),
  c('', 'total_pressure', 'bar'),
  c('', 'oxygen',         'percent')
)

inputs <- calculate_temperature_response(inputs, c4_temperature_param_vc, 'Tleaf')

assim <- calculate_c4_assimilation(inputs, 0, 0.003, 250, 1, 0.5, 40, 200, 80)

# Now we can plot Ac, Apr, Apc, and Ar. From this plot, we can see that
# replacing the complicated quadratic equation with a simple minimum yields
# very different results. Although this approximation is helpful for
# understanding C4 photosythesis, it should not be used for calculations.
lattice::xyplot(
  Apr + Apc + Ar + Ac ~ inputs[, 'PCm'],
  data = assim$main_data,
  type = 'l',
  grid = TRUE,
  auto = TRUE,
  ylim = c(-5, 100),
  xlab = paste0('Partial pressure of CO2 in the mesophyll (', inputs$units$PCm, ')'),
  ylab = paste0('Net CO2 assimilation rate (', assim$units$An, ')')
)

# Likewise, we can look at Ajm, Ajbs, and Aj
lattice::xyplot(
  Ajm + Ajbs + Aj ~ inputs[, 'PCm'],
  data = assim$main_data,
  type = 'l',
  grid = TRUE,
  auto = TRUE,
  ylim = c(-5, 45),
  xlab = paste0('Partial pressure of CO2 in the mesophyll (', inputs$units$PCm, ')'),
  ylab = paste0('Net CO2 assimilation rate (', assim$units$An, ')')
)

# Finally, we can see whether enzyme activity or light limits overall
# assimilation. In this case, assimilation is always enzyme-limited.
lattice::xyplot(
  Ac + Aj + An ~ inputs[, 'PCm'],
  data = assim$main_data,
  type = 'l',
  grid = TRUE,
  auto = TRUE,
  ylim = c(-5, 40),
  xlab = paste0('Partial pressure of CO2 in the mesophyll (', inputs$units$PCm, ')'),
  ylab = paste0('Net CO2 assimilation rate (', assim$units$An, ')')
)

Calculate C4 assimilation rates using a hyperbola

Description

Calculates C4 assimilation rates based on an empirical hyperbolic model. This function can accomodate alternative colum names for the variables taken from Licor files in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_c4_assimilation_hyperbola(
    exdf_obj,
    c4_curvature,
    c4_slope,
    rL,
    Vmax,
    ci_column_name = 'Ci',
    hard_constraints = 0,
    perform_checks = TRUE,
    return_exdf = TRUE
  )

Arguments

Details

General Description of the Model

In contrast to the mechanistic model implemented in calculate_c4_assimilation, this is a simple empirical model for C4 assimilation based on a four-parameter hyperbola. In this model, the net CO2 assimilation rate (An) is given by

An = Ag - rL,

where Ag is the gross assimilation rate and rL is the respiration rate. In turn, Ag is given by the smaller root of the following quadratic equation:

curvature * Ag^2 - (Vinitial + Vmax) * Ag + Vinitial * Vmax = 0,

where 0 <= curvature <= 1 is an empirical curvature factor, Vmax is the maximum gross assimilation rate, and Vinitial represents the initial response of Ag to increases in the intercellular CO2 concentration (Ci):

Vinitial = slope * Ci.

Here the slope is another empirical factor.

By including the respiration offset, it is also possible to define two other quantities: the maximum net CO2 assimilation rate (Amax) and the initial net CO2 assimilation rate (Ainitial). These are given by

Amax = Vmax - rL

and

Ainitial = Vinitial - rL.

Overall, this model exhibits a linear response of An to Ci at low Ci, a flat plateau of An at high Ci, and a smooth transition between these regions. The sharpess of the transition is set by the curvature. When curvature = 1, the model simplifies to

An = min{Vinitial, Vmax} - rL = min{Ainitial, Amax}.

As the curvature increases to 1, the transition becomes smoother. When the curvature is not zero, An approaches Amax asymptotically, and may not reach Amax at a reasonable value of Ci.

Code implementation

In this function, curvature and slope above are referred to as c4_curvature and c4_slope to avoid any potential ambiguity with other models that may also have curvature and slope parameters.

Temperature response

Because this model does not represent any photosynthetic mechanisms, temperature response functions are not applied.

Hard constraints

Most input parameters to the this model have hard constraints on their values which are set by their interpretation; for example, Vmax cannot be negative and c4_curvature must lie between 0 and 1. Yet, because of measurement noise, sometimes it is necessary to use values outside these ranges when fitting an A-Ci curve with fit_c4_aci_hyperbola. To accomodate different potential use cases, it is possible to selectively apply these hard constraints by specifying different values of the hard_constraints input argument:

• hard_constraints = 0: No constraints are applied.

• hard_constraints = 1: Checks whether all Ci values are non-negative.

• hard_constraints = 2: Includes the same constraints as when hard_constraints is 1, which additional constraints on the parameters that can be fitted. For example, Vmax must be non-negative and c4_curvature must lie between 0 and 1.

If any input values violate any of the specified constraints, an error message will be thrown.

Value

The return value depends on the value of return_exdf:

• If return_exdf is TRUE, the return value is an exdf object with the following columns: Ag, Ainitial, Amax, An, c4_curvature, c4_slope, rL, Vinitial, Vmax, and c4_assimilation_hyperbola_msg. Most of these are calculated as described above, while several are copies of the input arguments with the same name. The c4_assimilation_hyperbola_msg is usually blank but may contain information about any issues with the inputs. The category for each of these new columns is calculate_c4_assimilation_hyperbola to indicate that they were created using this function.

• If return_exdf is FALSE, the return value is a numeric vector containing the calculated values of An.

Examples

# Simulate a C4 A-Ci curve and plot the net assimilation rate.
npts <- 101

inputs <- exdf(data.frame(
  Ci = seq(0, 1000, length.out = npts),
  total_pressure = 1
))

inputs <- document_variables(
  inputs,
  c('', 'Ci',             'micromol mol^(-1)'),
  c('', 'total_pressure', 'bar')
)

assim <- calculate_c4_assimilation_hyperbola(inputs, 0.8, 0.5, 1.0, 55)

lattice::xyplot(
  Ainitial + Amax + An ~ inputs[, 'Ci'],
  data = assim$main_data,
  type = 'l',
  grid = TRUE,
  auto = TRUE,
  ylim = c(-5, 65),
  xlab = paste0('Intercellular CO2 concentration (', inputs$units$Ci, ')'),
  ylab = paste0('Net CO2 assimilation rate (', assim$units$An, ')')
)

Calculate Gamma_star from Rubisco specificity

Description

Calculates the CO2 compensation point in the absence of non-photorespiratory CO2 release (Gamma_star) from the Rubisco specificity (on a molarity basis), the oxygen concentration (as a percentage), and the temperature-dependent solubilities of CO2 and O2 in H2O.

Usage

  calculate_gamma_star(
    exdf_obj,
    alpha_pr = 0.5,
    oxygen_column_name = 'oxygen',
    rubisco_specificity_column_name = 'rubisco_specificity_tl',
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

The CO2 compensation point in the absence of non-photorespiratory CO2 release (Gamma_star) is the partial pressure of CO2 in the chloroplast at which CO2 gains from Rubisco carboxylation are exactly balanced by CO2 losses from Rubisco oxygenation; this quantity plays a key role in many photosynthesis calculations. One way to calculate its value is to use its definition, which can be found in many places, such as Equation 2.17 from von Caemmerer (2000):

Gamma_star = alpha_pr * O / S,

where O is the partial pressure (or mole fraction) of oxygen in the chloroplast, S is the Rubisco specificity on a gas basis, and alpha_pr is the number of CO2 molecules released by the photorespiratory cycle following each RuBP oxygenation (usually assumed to be 0.5).

The Rubisco specificity is often measured from an aqueous solution where the concentrations of O2 and CO2 are specified as molarities (moles of dissolved CO2 or O2 per mole of H2O). In this context, the equation above becomes

Gamma_star_aq = alpha_pr * O_aq / S_aq,

where Gamma_star_aq and O_aq are the molarities of CO2 and O2 corresponding to Gamma_star and O under the measurement conditions and S_aq is the specificity on a molarity basis.

Henry's law can be used to relate these two versions of the equation; Henry's law states that the concentration of dissolved gas is proportional to the partial pressure of that gas outside the solution. The proportionality factor H is called Henry's constant (or sometimes the solubility), and its value depends on the temperature, gas species, and other factors. Using Henry's law, we can write Gamma_star_aq = Gamma_star_aq * H_CO2 and O = O_aq * H_O2, where H_CO2 is Henry's constant for CO2 dissolved in H2O and H_O2 is Henry's constant for O2 dissolved in H2O. With these replacements, we can re-express the equation above as:

Gamma_star / H_CO2 = alpha_pr * (O / H_O2) / S_aq

Solving for Gamma_star, we see that:

Gamma_star = (alpha_pr * O / S_aq) * (H_CO2 / H_O2).

In other words, both the Rubisco specificity (as measured on a molarity basis) and the ratio of the two Henry's constants (H_CO2 / H_O2) play a role in determining Gamma_star. This equation also shows that it is possible to relate S (the specificity on a gas concentration basis) and S_aq as S = S_aq * H_O2 / H_CO2.

The values of H_O2 and H_CO2 can be calculated from the temperature using Equation 18 from Tromans (1998) and Equation 4 from Carroll et al. (1991), respectively.

In calculate_gamma_star, it is assumed that the value of specificity S_aq was was measured or otherwise determined at the leaf temperature; the leaf temperature is only used to determine the values of the two Henry's constants. Sometimes it is necessary to calculate the temperature-dependent value of the specificity using an Arrhenius equation; this can be accomplished via the calculate_temperature_response_arrhenius function from PhotoGEA.

Finally, it is important to note that Gamma_star can also be directly calculated using an Arrhenius equation, rather than using the oxygen concentration and the specificity. The best approach for determining a value of Gamma_star in any particular situation will generally depend on the available information and the measurement conditions.

References:

von Caemmerer, S. "Biochemical Models of Leaf Photosynthesis." (CSIRO Publishing, 2000) [doi:10.1071/9780643103405].

Carroll, J. J., Slupsky, J. D. and Mather, A. E. "The Solubility of Carbon Dioxide in Water at Low Pressure." Journal of Physical and Chemical Reference Data 20, 1201–1209 (1991) [doi:10.1063/1.555900].

Tromans, D. "Temperature and pressure dependent solubility of oxygen in water: a thermodynamic analysis." Hydrometallurgy 48, 327–342 (1998) [doi:10.1016/S0304-386X(98)00007-3].

Value

An exdf object based on exdf_obj that includes the following additional columns, calculated as described above: Gamma_star_tl (the value of Gamma_star at the leaf temperature), H_CO2, H_O2, and specificity_gas_basis. There are many choices for expressing Henry's constant values; here we express them as molalities per unit of pressure: (mol solute / kg H2O) / Pa. The category for each of these new columns is calculate_gamma_star to indicate that they were created using this function.

Examples

# Example 1: Calculate Gamma_star for each point in a gas exchange log file
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
)

licor_data <- get_oxygen_from_preamble(licor_data)

licor_data <- set_variable(
    licor_data,
    'rubisco_specificity_tl',
    'M / M',
    value = 90
)

licor_data <- calculate_gamma_star(licor_data)

licor_data[, c('specificity_gas_basis', 'oxygen', 'Gamma_star_tl'), TRUE]

# Example 2: Calculate Gamma_star at 21% and 2% oxygen for a Rubisco whose
# specificity was measured to be 100 M / M at 25 degrees C.

exdf_obj <- calculate_gamma_star(
  exdf(
    data.frame(
      oxygen = c(2, 21),
      rubisco_specificity_tl = c(100, 100),
      TleafCnd = c(25, 25)
    ),
    data.frame(
      oxygen = 'percent',
      rubisco_specificity_tl = 'M / M',
      TleafCnd = 'degrees C',
      stringsAsFactors = FALSE
    )
  )
)

exdf_obj[, c('specificity_gas_basis', 'oxygen', 'Gamma_star_tl'), TRUE]

# Example 3: Here we recreate Figure 1 from Long, S. P. "Modification of the
# response of photosynthetic productivity to rising temperature by atmospheric
# CO2 concentrations: Has its importance been underestimated?" Plant, Cell and
# Environment 14, 729–739 (1991). This is a fairly complicated example where
# Arrhenius constants for Rubisco parameters are determined by fitting
# published data and then used to determine the Rubisco specificity across a
# range of temperatures.

# Specify leaf temperature and oxygen concentration
leaf_temp <- seq(0, 50, by = 0.1)

exdf_obj <- exdf(
  data.frame(
    oxygen = rep_len(21, length(leaf_temp)),
    TleafCnd = leaf_temp
  ),
  data.frame(
    oxygen = 'percent',
    TleafCnd = 'degrees C',
    stringsAsFactors = FALSE
  )
)

# Get Arrhenius constants for Rubisco parameters using data from Table 2 of
# Jordan, D. B. and Ogren, W. L. "The CO2/O2 specificity of ribulose
# 1,5-bisphosphate carboxylase/oxygenase" Planta 161, 308–313 (1984).
rubisco_info <- data.frame(
  temperature = c(7,    12,   15,   25,   30,   35),
  Vc          = c(0.13, 0.36, 0.63, 1.50, 1.90, 2.90),
  Kc          = c(2,    3,    4,    11,   14,   19),
  Ko          = c(550,  510,  510,  500,  600,  540),
  Vo          = c(0.24, 0.48, 0.69, 0.77, 1.1,  1.3)
)

rubisco_info$x <- 1 / (8.314e-3 * (rubisco_info$temperature + 273.15))

lm_Vc <- stats::lm(log(Vc) ~ x, data = rubisco_info)
lm_Kc <- stats::lm(log(Kc) ~ x, data = rubisco_info)
lm_Ko <- stats::lm(log(Ko) ~ x, data = rubisco_info)
lm_Vo <- stats::lm(log(Vo) ~ x, data = rubisco_info)

arrhenius_info <- list(
  Vc = list(
    c = as.numeric(lm_Vc$coefficients[1]),
    Ea = -as.numeric(lm_Vc$coefficients[2]),
    units = 'micromol / mg / min'
  ),
  Kc = list(
    c = as.numeric(lm_Kc$coefficients[1]),
    Ea = -as.numeric(lm_Kc$coefficients[2]),
    units = 'microM'
  ),
  Ko = list(
    c = as.numeric(lm_Ko$coefficients[1]),
    Ea = -as.numeric(lm_Ko$coefficients[2]),
    units = 'microM'
  ),
  Vo = list(
    c = as.numeric(lm_Vo$coefficients[1]),
    Ea = -as.numeric(lm_Vo$coefficients[2]),
    units = 'micromol / mg / min'
  )
)

# Get temperature-dependent values of Rubisco parameters using Arrhenius
# equations
exdf_obj <- calculate_temperature_response_arrhenius(
  exdf_obj,
  arrhenius_info
)

# Calculate temperature-dependent specificity values
exdf_obj <- set_variable(
  exdf_obj,
  'rubisco_specificity_tl',
  units = 'M / M',
  value = exdf_obj[, 'Vc'] * exdf_obj[, 'Ko'] /
    (exdf_obj[, 'Vo'] * exdf_obj[, 'Kc'])
)

# Calculate Gamma_star and Henry constants
exdf_obj <- calculate_gamma_star(exdf_obj)

# Make a plot similar to Figure 1 from Long (1991)
lattice::xyplot(
  rubisco_specificity_tl + H_CO2 / H_O2 ~ TleafCnd,
  data = exdf_obj$main_data,
  auto = TRUE,
  grid = TRUE,
  type = 'l',
  xlim = c(0, 50),
  ylim = c(0, 250),
  xlab = "Temperature [ degrees C ]",
  ylab = "Rubisco specificity or ratio of Henry's constants (H_CO2 / H_O2)\n[ dimensionless ]"
)

# We can also make a plot of Gamma_star across this range
lattice::xyplot(
  Gamma_star_tl ~ TleafCnd,
  data = exdf_obj$main_data,
  auto = TRUE,
  grid = TRUE,
  type = 'l',
  xlim = c(0, 50),
  ylim = c(0, 120),
  xlab = "Temperature [ degrees C ]",
  ylab = paste('Gamma_star at leaf temperature [', exdf_obj$units$Gamma_star_tl, ']')
)

Calculate gas properties that are typically not included in Licor files

Description

Calculates gas properties that are typically not included in Licor files. This function can accomodate alternative column names for the variables taken from the Licor file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_gas_properties(
    licor_exdf,
    a_column_name = 'A',
    ca_column_name = 'Ca',
    total_pressure_column_name = 'total_pressure',
    e_column_name = 'E',
    gbw_column_name = 'gbw',
    gsw_column_name = 'gsw',
    h2o_s_column_name = 'H2O_s',
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

By default, a Licor file provides the following gas concentrations and conductances:

• Water vapor conductance to diffusion through the stomata (gsw).

• Water vapor conductance to diffusion through the boundary layer (gbw).

• Water vapor conductance to diffusion from the leaf's intercellular spaces to the ambient air; in other words, the total conductance to water vapor (gtw).

• Water vapor concentration in the sample cell (H2O_s).

• CO2 conductance to diffusion from the leaf's intercellular spaces to the ambient air; in other words, the total conductance to CO2 (gtc).

• CO2 concentration in the sample cell, corrected for any chamber leaks (Ca).

• CO2 concentration in the leaf's intercellular spaces (Ci).

However, it is sometimes helpful to know the "missing" conductances and concentrations, for example, when calculating mesophyll conductances or Ball-Berry parameters. This function adds these missing values, along with a few related water vapor properties:

• Water vapor concentration at the sample surface (H2O_surf).

• Water vapor concentration in the leaf's intercellular spaces (H2O_i).

• Saturation water vapor pressure at the leaf temperature (SVPleaf).

• Relative humidity at the leaf surface (RHleaf).

• CO2 conductance to diffusion through the stomata (gsc).

• CO2 conductance to diffusion through the boundary layer (gbc).

• CO2 concentration at the leaf surface (Cs).

Equations used for these calculations

The equations used to calculate these quantities can be found in the Licor Li-6800 manual (Appendix C), which relies heavily on Appendix 2 of the following paper: von Caemmerer, S. & Farquhar, G. D. "Some relationships between the biochemistry of photosynthesis and the gas exchange of leaves" Planta 153, 376–387 (1981) [doi:10.1007/BF00384257]

Equation C-79 in the Licor manual describes the total flow of water vapor from the leaf interior to the ambient air using gtw, H2O_i, H2O_s, and the transpiration rate E:

(1) gtw = E * (1000 - (H2O_i + H2O_s) / 2) / (H2O_i - H2O_s)

In steady-state conditions, the flux of H2O molecules across any portion of the gas flow is identical to E, so we can also apply this equation to the flow of water vapor from the leaf surface to the ambient air:

(2) gbw = E * (1000 - (H2O_surf + H2O_s) / 2) / (H2O_surf - H2O_s)

Equation (2) can be solved for H2O_surf:

(3) H2O_surf = (E * (1000 - H2O_s / 2) + gbw * H2O_s) / (gbw + E / 2)

Equation C-70 in the Licor manual describes how to calculate saturation water vapor pressure from air temperature. At the leaf surface, the air temperature should be the same as the leaf temperature (Tleaf; in degrees C), so we can determine SVPleaf using Equation C-70 as follows:

(4) SVPleaf = 0.6135 * e^((17.502 * Tleaf) / (240.97 + Tleaf))

For gas exchange measurements, we assume that water vapor is saturated in the leaf's intecellular spaces, so we can determine H2O_i from SVPleaf and the relationship between partial pressure and molar gas concentration:

(5) H2O_i = SVPleaf / Pcham = SVPleaf / (Pa + deltaPcham)

where Pcham is th total pressure in the sample chamber, Pa is the atmospheric pressure, and deltaPcham is the chamber overpressure. These are related by Pcham = Pa + deltaPcham.

The relative humidity at the leaf surface RHleaf can be determined from H2O_surf and SVPleaf using the definitions of relative humidity and partial pressure:

(6) RHleaf = Pwl / SVPleaf = H2O_surf * (Pa + deltaPcham) / SVPleaf

where Pwl, the partial pressure of H2O at the leaf surface, is given by H2O_surf * Pcham.

The CO2 conductances through the stomata and boundary layer can be determined from the corresponding H2O conductances using the ratios of molecular diffusivities for the two molecules, as explained in the vicinty of Equation C-106 in the Licor manual:

(7) gsc = gsw / 1.6

(8) gbc = gbw / 1.37

Equation C-105 in the Licor manual describes the flow of CO2 from the ambient air to the intercellular spaces:

(9) C_i = ((gtc - E / 2) * Ca - A) / (gtc + E / 2)

where we have replaced C_s (the CO2 concentration in the sample chamber) with Ca for clarity. In steady state conditions, the flows of H2O and CO2 are identical to E and A, respectively, so we can also apply this equation to the flow of CO2 from the ambient air to the leaf surface:

(10) Csurface = ((gbc - E / 2) * Ca - A) / (gbc + E / 2)

This function uses Equations (3)-(8) and (10) to calculate the desired values.

Value

An exdf object based on licor_exdf that includes the following additional columns, calculated as described above: H2O_surf, SVPleaf, H2O_i, RHleaf, gsc, gbc, and Csurface. The category for each of these new columns is calculate_gas_properties to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package, calculate the
# total pressure, and calculate additional gas properties.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_total_pressure(licor_file)

licor_file <- calculate_gas_properties(licor_file)

licor_file$units$RHleaf      # View the units of the new `RHleaf` column
licor_file$categories$RHleaf # View the category of the new `RHleaf` column
licor_file[,'RHleaf']        # View the values of the new `RHleaf` column

Calculate mesophyll conductance to CO2 diffusion

Description

Calculates mesophyll conductance to CO2 diffusion (gmc) from combined gas exchange and isotope discrimination measurements as described in Busch et al. (2020). This function can accomodate alternative colum names for the variables taken from exdf_obj; it also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_gm_busch(
    exdf_obj,
    e = -3,
    f = 11,
    e_star_equation = 20,
    gm_type = 'dis',
    a_bar_column_name = 'a_bar',
    a_column_name = 'A',
    ci_column_name = 'Ci',
    co2_s_column_name = 'CO2_s',
    csurface_column_name = 'Csurface',
    delta_c13_r_column_name = 'delta_C13_r',
    delta_obs_growth_column_name = 'Delta_obs_growth',
    delta_obs_tdl_column_name = 'Delta_obs_tdl',
    gamma_star_column_name = 'Gamma_star_tl',
    rl_column_name = 'RL',
    total_pressure_column_name = 'total_pressure',
    t_column_name = 't'
  )

Arguments

Details

This function uses a model for photosynthetic discrimination against 13C in C3 plants to determine mesophyll conductance values as described in Busch et al. (2020). That paper provides two alternate ways to calculate e_star, and two alternate ways to calculate mesophyll conductance gmc; this function allows the user to choose between them. In more detail:

• Isotopic fractionation due to day respiration (e_prime = e + e_star) is calculated with e_star given by either Equation 19 or 20 depending on the value of e_star_equation.

• Isotopic discrimination assuming infinite mesophyll conductance (Delta_i) is calculated by setting Cc = Ci in either Equation 2 or 13, depending on the value of gm_type.

• Mesophyll conductance to CO2 (gmc) is calculated using either Equation 21 or 22, depending on the value of gm_type.

Note 1: Setting e_star_equation = 19 and gm_type = 'con' should produce identical or similar results to calculate_gm_ubierna.

Note 2: Using e_star_equation = 20 and gm_type = 'dis' is expected to be more accurate, as discussed in Busch et al. (2020); however, be aware that this method requires a value for Delta_obs_growth, which may not always be available unless it is intentionally measured.

References:

Busch, F. A., Holloway-Phillips, M., Stuart-Williams, H. and Farquhar, G. D. "Revisiting carbon isotope discrimination in C3 plants shows respiration rules when photosynthesis is low." Nat. Plants 6, 245–258 (2020) [doi:10.1038/s41477-020-0606-6].

Value

An exdf object based on exdf_obj that includes the following additional columns, calculated as described above: e_prime, e_star, Delta_i, and gmc, as well as the values of a few intermediate calculations such as Delta_i_term_1 and Delta_i_term_2. The category for each of these new columns is calculate_gm_busch to indicate that they were created using this function.

Examples

## In this example we load gas exchange and TDL data files, calibrate the TDL
## data, pair the data tables together, and then calculate mesophyll conductance

# Read the TDL data file, making sure to interpret the time zone as US Central
# time
tdl_data <- read_gasex_file(
  PhotoGEA_example_file_path('tdl_for_gm.dat'),
  'TIMESTAMP',
  list(tz = 'America/Chicago')
)

# Identify cycles within the TDL data
tdl_data <- identify_tdl_cycles(
  tdl_data,
  valve_column_name = 'valve_number',
  cycle_start_valve = 20,
  expected_cycle_length_minutes = 2.7,
  expected_cycle_num_valves = 9,
  timestamp_colname = 'TIMESTAMP'
)

# Use reference tanks to calibrate the TDL data
processed_tdl <- consolidate(by(
  tdl_data,
  tdl_data[, 'cycle_num'],
  process_tdl_cycle_erml,
  noaa_valve = 2,
  calibration_0_valve = 20,
  calibration_1_valve = 21,
  calibration_2_valve = 23,
  calibration_3_valve = 26,
  noaa_cylinder_co2_concentration = 294.996,
  noaa_cylinder_isotope_ratio = -8.40,
  calibration_isotope_ratio = -11.505
))

# Read the gas exchange data, making sure to interpret the time stamp in the US
# Central time zone
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
  'time',
  list(tz = 'America/Chicago')
)

# Get TDL valve information from Licor file name; for this TDL system, the
# reference valve is 12 when the sample valve is 11
licor_data <- get_sample_valve_from_filename(licor_data, list('11' = 12))

# Get oxygen info from the Licor file preamble (needed for calculate_gamma_star)
licor_data <- get_oxygen_from_preamble(licor_data)

# Pair the Licor and TDL data by locating the TDL cycle corresponding to each
# Licor measurement
licor_data <- pair_gasex_and_tdl(licor_data, processed_tdl$tdl_data)

# Calculate total pressure (needed for calculate_gas_properties)
licor_data <- calculate_total_pressure(licor_data)

# Calculate Csurface (needed for calculate_ternary_correction)
licor_data <- calculate_gas_properties(licor_data)

# Calculate ternary correction
licor_data <- calculate_ternary_correction(licor_data)

# Set Rubisco specificity (needed for calculate_gamma_star)
licor_data <- set_variable(
    licor_data,
    'rubisco_specificity_tl',
    'M / M',
    value = 90
)

# Calculate Gamma_star (needed for calculate_gm_busch)
licor_data <- calculate_gamma_star(licor_data)

# Calculate isotope discrimination (needed for calculate_gm_busch)
licor_data <- calculate_isotope_discrimination(licor_data)

# Set Delta_obs_growth to the average of Delta_obs_tdl over the first 6 points,
# where the ambient CO2 concentration was set to the atmospheric value (420 ppm)
# (needed for calculate_gm_busch).
licor_data <- set_variable(
  licor_data,
  'Delta_obs_growth',
  'ppt',
  value = mean(licor_data[1:6, 'Delta_obs_tdl'])
)

# Set respiration (needed for calculate_gm_busch)
licor_data <- set_variable(
  licor_data,
  'RL',
  'micromol m^(-2) s^(-1)',
  value = 1.2
)

# Calculate mesophyll conductance
licor_data <- calculate_gm_busch(licor_data)

# Calculate Cc using the new values of mesophyll conductance
licor_data <- calculate_temperature_response(
  licor_data,
  c3_temperature_param_flat['gmc_norm']
)

licor_data <- set_variable(
  licor_data,
  'gmc_at_25',
  units = licor_data$units$gmc,
  value = licor_data[, 'gmc']
)

licor_data <- apply_gm(licor_data)

# View some of the results
licor_data[, c('replicate', 'CO2_s', 'Delta_obs_tdl', 'e_prime', 'gmc', 'Ci', 'Cc')]

Calculate mesophyll conductance to CO2 diffusion

Description

Calculates mesophyll conductance to CO2 diffusion (gmc) from combined gas exchange and isotope discrimination measurements as described in Ubierna et al. (2018). This function can accomodate alternative colum names for the variables taken from exdf_obj; it also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_gm_ubierna(
    exdf_obj,
    e = -3,
    f = 11,
    a_bar_column_name = 'a_bar',
    a_column_name = 'A',
    ci_column_name = 'Ci',
    co2_s_column_name = 'CO2_s',
    csurface_column_name = 'Csurface',
    delta_c13_r_column_name = 'delta_C13_r',
    delta_obs_tdl_column_name = 'Delta_obs_tdl',
    gamma_star_column_name = 'Gamma_star_tl',
    rl_column_name = 'RL',
    total_pressure_column_name = 'total_pressure',
    t_column_name = 't'
  )

Arguments

Details

This function uses the comprehensive model for photosynthetic discrimination against 13C in C3 plants to calculate mesophyll conductance, as described in Ubierna et al. (2018). In particular, the following equations from that source are implemented in the code:

• Isotopic fractionation due to day respiration (e_prime) is calculated using Equations 28 and 30.

• Isotopic discrimination due to photorespiration (Delta_f), due to day respiration (Delta_e), and that would occur if Ci = Cc in the absence of any respiratory fractionation (Delta_i) are calculated using Equations 34, 33, and 31, respectively.

• Mesophyll conductance to CO2 diffusion (gmc) is calculated using Equation 44. This equation is broken up into two factors called Delta_difference and equation_top which are separately returned in the output from calculate_gm_ubierna.

For an alternative method for calculating gmc, see calculate_gm_busch.

References:

Ubierna, N., Holloway-Phillips, M.-M. and Farquhar, G. D. "Using Stable Carbon Isotopes to Study C3 and C4 Photosynthesis: Models and Calculations." in Photosynthesis: Methods and Protocols (ed. Covshoff, S.) 155–196 (Springer, 2018) [doi:10.1007/978-1-4939-7786-4_10].

Value

An exdf object based on exdf_obj that includes the following additional columns, calculated as described above: e_prime, Delta_i, Delta_e, Delta_f, Delta_difference, equation_top, and gmc. The category for each of these new columns is calculate_gm_ubierna to indicate that they were created using this function.

Examples

## In this example we load gas exchange and TDL data files, calibrate the TDL
## data, pair the data tables together, and then calculate mesophyll conductance

# Read the TDL data file, making sure to interpret the time zone as US Central
# time
tdl_data <- read_gasex_file(
  PhotoGEA_example_file_path('tdl_for_gm.dat'),
  'TIMESTAMP',
  list(tz = 'America/Chicago')
)

# Identify cycles within the TDL data
tdl_data <- identify_tdl_cycles(
  tdl_data,
  valve_column_name = 'valve_number',
  cycle_start_valve = 20,
  expected_cycle_length_minutes = 2.7,
  expected_cycle_num_valves = 9,
  timestamp_colname = 'TIMESTAMP'
)

# Use reference tanks to calibrate the TDL data
processed_tdl <- consolidate(by(
  tdl_data,
  tdl_data[, 'cycle_num'],
  process_tdl_cycle_erml,
  noaa_valve = 2,
  calibration_0_valve = 20,
  calibration_1_valve = 21,
  calibration_2_valve = 23,
  calibration_3_valve = 26,
  noaa_cylinder_co2_concentration = 294.996,
  noaa_cylinder_isotope_ratio = -8.40,
  calibration_isotope_ratio = -11.505
))

# Read the gas exchange data, making sure to interpret the time stamp in the US
# Central time zone
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
  'time',
  list(tz = 'America/Chicago')
)

# Get TDL valve information from Licor file name; for this TDL system, the
# reference valve is 12 when the sample valve is 11
licor_data <- get_sample_valve_from_filename(licor_data, list('11' = 12))

# Get oxygen info from the Licor file preamble (needed for calculate_gamma_star)
licor_data <- get_oxygen_from_preamble(licor_data)

# Pair the Licor and TDL data by locating the TDL cycle corresponding to each
# Licor measurement
licor_data <- pair_gasex_and_tdl(licor_data, processed_tdl$tdl_data)

# Calculate total pressure (needed for calculate_gas_properties)
licor_data <- calculate_total_pressure(licor_data)

# Calculate Csurface (needed for calculate_ternary_correction)
licor_data <- calculate_gas_properties(licor_data)

# Calculate ternary correction
licor_data <- calculate_ternary_correction(licor_data)

# Set Rubisco specificity (needed for calculate_gamma_star)
licor_data <- set_variable(
    licor_data,
    'rubisco_specificity_tl',
    'M / M',
    value = 90
)

# Calculate Gamma_star (needed for calculate_gm_ubierna)
licor_data <- calculate_gamma_star(licor_data)

# Calculate isotope discrimination (needed for calculate_gm_ubierna)
licor_data <- calculate_isotope_discrimination(licor_data)

# Set respiration (needed for calculate_gm_ubierna)
licor_data <- set_variable(
  licor_data,
  'RL',
  'micromol m^(-2) s^(-1)',
  value = 1.2
)

# Calculate mesophyll conductance
licor_data <- calculate_gm_ubierna(licor_data)

# Calculate Cc using the new values of mesophyll conductance
licor_data <- calculate_temperature_response(
  licor_data,
  c3_temperature_param_flat['gmc_norm']
)

licor_data <- set_variable(
  licor_data,
  'gmc_at_25',
  units = licor_data$units$gmc,
  value = licor_data[, 'gmc']
)

licor_data <- apply_gm(licor_data)

# View some of the results
licor_data[, c('replicate', 'CO2_s', 'Delta_obs_tdl', 'gmc', 'Ci', 'Cc')]

Calculate photosynthetic isotope discrimination

Description

Calculates photosynthetic carbon isotope discrimination from combined gas exchange and tunable diode laser absorption spectroscopy measurements.

Usage

  calculate_isotope_discrimination(
    exdf_obj,
    co2_r_column_name = 'CO2_r',
    co2_s_column_name = 'CO2_s',
    delta_C13_r_column_name = 'delta_C13_r',
    delta_C13_s_column_name = 'delta_C13_s',
    h2o_r_column_name = 'H2O_r',
    h2o_s_column_name = 'H2O_s',
    tdl_12C_r_column_name = 'calibrated_12c_r',
    tdl_12C_s_column_name = 'calibrated_12c_s'
  )

Arguments

Details

As described in Ubierna et al. (2018), photosynthetic 13C discrimination can be determined from combined gas exchange and tunable diode laser (TDL) absorption spectroscopy measurements according to:

Delta_obs = xsi * (delta_out - delta_in) / (1 + delta_out - xsi * (delta_out - delta_in)),

where Delta_obs is the observed discrimination, delta_in and delta_out are the carbon isotope ratios in dry air flowing in and out of the leaf chamber. xsi is given by

xsi = C_in / (C_in - C_out),

where C_in and C_out are the mole fractions of 12CO2 in dry air flowing in and out of the leaf chamber. (See equations 5 and 6 in Ubierna et al. (2018)).

In practice, there are multiple options for calculating Delta_obs and xsi because CO2 concentrations are measured by both the gas exchange system and the TDL. For example, we can alternately calculate xsi as xsi_tdl = C_in_tdl / (C_in_tdl - C_out_tdl) or xsi_gasex = C_in_gasex / (C_in_gasex - C_out_gasex). Likewise, we can also calculate Delta_obs_tdl using xsi_tdl or Delta_obs_gasex using xsi_gasex. The TDL values are typically preferred in subsequent calculations, but it can be useful to compare the two different versions as a consistency check; the TDL and gas exchange values should be similar to each other.

There are two subtelties associated with xsi_gasex. One is that the gas exchange system generally measures the total CO2 concentration, not just the 12CO2 concentration. Typically there is much less 13CO2 than 12CO2 so this is usually not a large source of error.

The other issue is that the gas exchange system generally measures CO2 concentrations in wet air. Thus, it is important to use "corrected" values of CO2 concentrations that account for the "dilution effect" due to water vapor in the air. This effect is described in the Licor LI-6400 manual: "This is a correction we don’t do, at least when computing CO2 concentration in the LI-6400. The dilution effect is simply this: as you add molecules of a gas (water vapor, for example) to a mixture, the fraction of that mixture that is made up of something else (mole fraction of CO2, for instance) has to decrease, since the total number of molecules in the mixture has increased. Now for an airsteam flowing though a chamber containing a transpiring leaf (or in a chamber sitting on moist soil), there very definitely is dilution. However, we ignore that effect when computing CO2 concentration, but account for it when computing photosynthetic rate (or soil CO2 efflux). Thus, the LI-6400 IRGA is always indicating the actual CO2 concentration, not what the CO2 concentration would be if there were no water vapor in it."

To account for the dilution effect, we define a "corrected" CO2 concentration as CO2_corrected = CO2 / (1 - H2O), where H2O is the water vapor concentration in the air. Note: the TDL always measures concentrations in dry air, so no correction is required.

References:

Ubierna, N., Holloway-Phillips, M.-M. and Farquhar, G. D. "Using Stable Carbon Isotopes to Study C3 and C4 Photosynthesis: Models and Calculations." in Photosynthesis: Methods and Protocols (ed. Covshoff, S.) 155–196 (Springer, 2018) [doi:10.1007/978-1-4939-7786-4_10].

Value

An exdf object based on exdf_obj that includes several new columns: CO2_r_corrected, CO2_s_corrected, Delta_obs_gasex, Delta_obs_tdl, xsi_gasex, and xsi_tdl.

Examples

## In this example we load gas exchange and TDL data files, calibrate the TDL
## data, pair the data tables together, and then calculate isotope
## discrimination

# Read the TDL data file, making sure to interpret the time zone as US Central
# time
tdl_data <- read_gasex_file(
  PhotoGEA_example_file_path('tdl_for_gm.dat'),
  'TIMESTAMP',
  list(tz = 'America/Chicago')
)

# Identify cycles within the TDL data
tdl_data <- identify_tdl_cycles(
  tdl_data,
  valve_column_name = 'valve_number',
  cycle_start_valve = 20,
  expected_cycle_length_minutes = 2.7,
  expected_cycle_num_valves = 9,
  timestamp_colname = 'TIMESTAMP'
)

# Use reference tanks to calibrate the TDL data
processed_tdl <- consolidate(by(
  tdl_data,
  tdl_data[, 'cycle_num'],
  process_tdl_cycle_erml,
  noaa_valve = 2,
  calibration_0_valve = 20,
  calibration_1_valve = 21,
  calibration_2_valve = 23,
  calibration_3_valve = 26,
  noaa_cylinder_co2_concentration = 294.996,
  noaa_cylinder_isotope_ratio = -8.40,
  calibration_isotope_ratio = -11.505
))

# Read the gas exchange data, making sure to interpret the time stamp in the US
# Central time zone
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
  'time',
  list(tz = 'America/Chicago')
)

# Get TDL valve information from Licor file name; for this TDL system, the
# reference valve is 12 when the sample valve is 11
licor_data <- get_sample_valve_from_filename(licor_data, list('11' = 12))

# Pair the Licor and TDL data by locating the TDL cycle corresponding to each
# Licor measurement
licor_data <- pair_gasex_and_tdl(licor_data, processed_tdl$tdl_data)

# Calculate isotope discrimination
licor_data <- calculate_isotope_discrimination(licor_data)

# View some of the results
licor_data[, c('A', 'xsi_gasex', 'xsi_tdl', 'Delta_obs_gasex', 'Delta_obs_tdl')]

Calculate maximum electron transport rate

Description

Calculates maximum electron transport rates (Jmax) from estimates of the electron transport rate (J) at particular values of incident light (Qin).

This function is typically used after fit_c3_aci, fit_c3_variable_j, or fit_c4_aci is used to estimate values of J.

Usage

  calculate_jmax(
    data_table,
    alpha_j_at_25 = 'column',
    theta_j_at_25 = 'column',
    alpha_j_norm_column_name = 'alpha_j_norm',
    qin_column_name = 'Qin_avg',
    theta_j_norm_column_name = 'theta_j_norm',
    tleaf_column_name = 'TleafCnd_avg',
    ...
  )

Arguments

Details

Basic Requirements:

This function requires that data_table contains columns called J_at_25 and J_tl_avg, as would be included in the output from one of the PhotoGEA fitting functions (fit_c3_aci, fit_c3_variable_j, and fit_c4_aci). These will be used to calculate values of Jmax at 25 degrees C and at leaf temperature.

If any columns for the J confidence intervals are included in data_table (J_at_25_upper, J_at_25_lower, J_tl_avg_upper, or J_tl_avg_lower), the corresponding confidence intervals for Jmax will also be calculated.

By default, this function will take values of alpha_j and theta_j from columns of data_table with the same names.

If data_table is an exdf object, units will be checked for any columns used in the calculations.

Overview of Jmax Calculations:

The potential electron transport rate going to support RuBP regeneration (J) depends on the available light energy. J quickly increases with the incident photosynthetically active photon flux density (Qin) at low light levels, gradually reaching a plateau at high values of Qin. Although other mathematical representations have been used (Walker et al. 2021), this dependence is typically represented as a non-rectangular hyperbola:

J = (I2 + Jmax - sqrt[(I2 + Jmax)^2 - 4 * theta_j * I2 * Jmax]) / (2 * theta_j), (Eq. 1)

where Jmax is the maximum value of J that would be achieved at infinitely large Qin, 0 < theta_j <= 1 is an empirical curvature parameter, and I2 is the useful energy absorbed by photosystem II. In turn, I2 is calculated by

I2 = alpha_j * Qin,

where alpha_j is the apparent quantum efficiency of electron transport. alpha_j is often defined as

alpha_j = absorptance * phi_psii,max * beta_psii,

where absorptance is the leaf absorptance, phi_psii,max is the maximum quantum yield of photosytem II, and beta_psii is the fraction of light energy partitioned to photosystem II.

Equation 1 can be understood as a "smooth minimum" of two potential rates of electron transport: I2 (which increases linearly with Qin) and Jmax (which is independent of Qin). For lower light levels, I2 is the smaller rate, and J is approximately equal to I2; for very high light levels, Jmax is the smaller rate, and J is approximately equal to Jmax. For intermediate values of Qin, J smoothly transitions from I2 to Jmax.

This equation is often solved for Jmax, and thus it is necessary to consider the conditions for which the solution is appropriate. One key property of Equation 1 is that the largest possible value of J at a given Qin is I2, which only occurs when Jmax is much larger than I2. In other words, when considered as a function of Jmax, the range of the function in Equation 1 is 0 <= J <= I2.

Equation 1 can be solved for Jmax, enabling calculations of Jmax from estimates of J:

Jmax = J * (I2 - theta_j * J) / (I2 - J) (Eq. 2)

Because the range of the function in Equation 1 is 0 <= J <= I2, the domain of its inverse function (defined in Equation 2) is also 0 <= J <= I2. In other words, Jmax can only be calculated using Equation 2 when J < I2. Otherwise, there is no value of Jmax that can reproduce the value of J for the given value of alpha_j. This restriction can also be derived more rigorously; see the Detailed algebra section below for more information.

If J >= I2, the calculate_jmax function will return NA for the value of Jmax. This behavior can be bypassed by setting the optional input argument ignore_restriction to TRUE, but this is not recommended outside of pedagogical purposes. See Example 2 below for a demonstration of what goes wrong when Equation 2 is used for J >= I2.

Note that this issue is more significant at lower light levels. For example, assuming a typical value of alpha_j (0.293), I2 for Qin = 1800 micromol / m^2 / s would be 527.4 micromol / m^2 / s. Values of J are typically smaller than this, so an estimate of Jmax can almost always be made. But if a curve were measured at Qin = 300, I2 would only be 87.9 micromol / m^2 / s, placing a stronger restriction on the values of J where Jmax can be estimated. Say the best-fit value of J was 88.9 micromol / m^2 / s for a curve measured with Qin = 300 micromol / m^2 / s; in this case, it would not be possible to estimate Jmax, potentially indicating that the assumed value of alpha_j was not correct.

Typical values:

According to von Caemmerer (2000), typical values of absorptance, phi_psii,max, and beta_psii are 0.85, 1 - 0.15, and 0.5, respectively, leading to alpha_j = 0.36125, and the curvature parameter theta_j is typically 0.7.

Bernacchi et al. (2003) reports that phi_psii,max is 0.6895 for light-adapted leaves at 25 degrees C, while theta_j at 25 degrees C is 0.97875. Using this value of of phi_psii,max with typical values of absorptance and beta_psii results in an alpha_j estimate of 0.2930375.

It is not clear whether the temperture response defined in Bernacchi et al. (2003) is applicable to C4 leaves. For C4 leaves, it may be better to use the temperature-independent estimates from von Caemmerer (2000).

PhotoGEA provides two Jmax parameter lists that can be passed to calculate_temperature_response: jmax_temperature_param_bernacchi (implements the Bernacchi et al. 2003 values) and jmax_temperature_param_flat (implements the von Caemmerer 2000 values). Each of these parameter lists will calculate values of alpha_j_at_25, alpha_j_norm, theta_j_at_25, and theta_j_norm.

Absorbed light basis:

Values of Jmax can also be estimated from the absorbed photosynthetically active photon flux density (Qabs). In that case, we can regroup the terms in the definition of I2 as follows:

I2 = (Qin * absorptance) * (phi_psii,max * beta_psii) = Qabs * alpha_j_abs,

where alpha_j_abs is given by phi_psii,max * beta_psii. When working in this basis, the default value of alpha_j at 25 degrees C should be divided by the assumed absorptance (0.85). For example, the default value of alpha_j_at_25 used with the Bernacchi et al. (2003) parameters is 0.2930375, so dividing this by 0.95 would yielding an alpha_j_abs value of about 0.345. This value could be passed directly to calculate_jmax via the alpha_j_at_25 input argument, overriding the default value. Along with this change, it would also be necessary to change the name of the light column, likely to Qabs_avg.

Why PhotoGEA Uses a Separate Function for Jmax:

In principle, values of Jmax could be estimated by the fitting functions that estimate J: fit_c3_aci, fit_c3_variable_j, and fit_c4_aci. Instead, PhotoGEA requires users to use a separate function (calculate_jmax) to estimate Jmax. This serves several purposes:

• It highlights that estimates of Jmax are made using the same equations for C3 and C4 leaves.

• It leaves open the possibility of other estimates of Jmax, such as those based on a rectangular hyperbola instead of the non-rectangular hyperbola used here.

• It emphasizes that sometimes it is not possible to provide an estimate for Jmax, depending on the values of Qin, alpha_j, and J, because of the requirement that J < I2 = alpha_j * Qin.

The last point is especially important. If Jmax were varied during the fitting process, and J was estimated from Jmax using Equation 1, there would be a restriction on the possible values of J that could be obtained: J < alpha_j * Qin. This could potentially bias the fitting results, since it may be the case that the best fit would be found for J outside this range.

In other words, keeping estimates of Jmax separate from the fitting process ensures that the values of alpha_j and theta_j have no influence on the fits or best-fit values of J. This is important since the true values of these parameters for a particular leaf are difficult or impossible to determine.

Detailed algebra:

Here we will solve Equation 1 for Jmax, arriving at Equation 2. This algebra is reproduced here to highlight the important restriction that J < I2.

First, multiply both sides of Equation 1 by 2 * theta_j:

2 * theta_j * J = I2 + Jmax - sqrt[(I2 + Jmax)^2 - 4 * theta_j * I2 * Jmax]. (Eq. 3)

Next, isolate the square root term on one side:

I2 + Jmax - 2 * theta_j * J = sqrt[(I2 + Jmax)^2 - 4 * theta_j * I2 * Jmax]. (Eq. 4)

A key point here is that the right hand side cannot be negative, since the square root of a real number is never negative. Thus, the left hand side also cannot be negative. In other words,

I2 + Jmax - 2 * theta_j * J >= 0. (Eq. 5)

We will return to this restriction later. For now, we square both sides of Equation 4:

(I2 + Jmax)^2 - 4 * theta_j * J * (I2 + Jmax) + 4 * theta_j^2 * J^2 = (I2 + Jmax)^2 - 4 * theta_j * I2 * Jmax. (Eq. 6)

The term (I2 + Jmax)^2 appears on both sides of Equation 6 and can therefore be cancelled out. Grouping the remaining terms that contain Jmax on one side, we have:

4 * theta_j * Jmax * (I2 - J) = 4 * theta_j * J * (I2 - theta_j * J) (Eq. 7)

Finally, provided that I2 - J is not zero (in other words, that I2 is not equal to J), we can divide both sides of Equation 7 by 4 * theta_j * (I2 - J) to obtain Equation 2 above.

Now, we can use this expression (Equation 2) to replace Jmax in Equation 5:

I2 + J * (I2 - theta_j * J) / (I2 - J) - 2 * theta_j * J >= 0. (Eq. 8)

This can be converted to a single ratio as follows:

[(I2 - 2 * theta_j * J) * (I2 - J) + J * (I2 - theta_j * J)] / (I2 - J) >= 0. (Eq. 9)

Multiplying out the factors in the numerator and collecting like terms, Equation 9 becomes

[I2^2 - 2 * theta_j * I2 * J + theta_j * J^2] / (I2 - J) >= 0. (Eq. 10)

Because theta_j must lie between 0 and 1, theta_j^2 is always less than or equal to theta_j. This allows us to place a lower bound on the value of the numerator of the left hand side of Equation 10:

I2^2 - 2 * theta_j * I2 * J + theta_j * J^2 >= I2^2 - 2 * theta_j * I2 * J + theta_j^2 * J^2. (Eq. 11)

The right hand side of Equation 11 can be refactored:

I2^2 - 2 * theta_j * I2 * J + theta_j * J^2 >= (I2 - theta_j J)^2. (Eq. 12)

The right hand side of Equation 12 can never be negative, so from this we can also conclude that the numerator of the left hand side of Equation 10 can also never be negative. Thus, the inequality in Equation 10 is satisfied whenever its denominator is positive. In other words, whenever I2 - J > 0, or, equivalently, J < I2.

Thus, we have shown that Equation 2 holds whenever J < I2, since, when this inequality is satisfied, Equation 5 is also satisfied.

Although we do not do so here, it can be shown that when I2 < J, the value of Jmax that would be calculated by Equation 2 is the inverse of

J = (I2 + Jmax + sqrt[(I2 + Jmax)^2 - 4 * theta_j * I2 * Jmax]) / (2 * theta_j) (Eq. 13)

rather than the inverse of Equation 1. Note the difference: in Equation 13, the square root term is added to I2 + Jmax rather than subtracted. This is a "smooth maximum" function, rather than a smooth minimum. In fact, whenever I2 > Jmax, Equation 13 would predict J > Jmax, clearly a nonsensical result. Likewise, the inverse of the function in Equation 13 would predict some values of Jmax that are smaller than J. Example 2 below shows that it can even return negative values of Jmax, which is clearly not reasonable from a biological perspective.

References:

• von Caemmerer, S. "Biochemical Models of Leaf Photosynthesis" (CSIRO Publishing, 2000) [doi:10.1071/9780643103405].

• Walker, A. P. et al. "Multi-hypothesis comparison of Farquhar and Collatz photosynthesis models reveals the unexpected influence of empirical assumptions at leaf and global scales." Global Change Biology 27, 804–822 (2021) [doi:10.1111/gcb.15366].

• Bernacchi, C. J., Pimentel, C. & Long, S. P. "In vivo temperature response functions of parameters required to model RuBP-limited photosynthesis" Plant, Cell & Environment 26, 1419–1430 (2003) [doi:10.1046/j.0016-8025.2003.01050.x].

Value

The return value is a table based on data_table that includes several new columns: I2_at_25, Jmax_at_25, Jmax_at_25_msg, I2_tl, Jmax_tl, and Jmax_tl_msg. The _msg columns indicate when the error condition J >= I2 has occurred.

If J confidence intervals were provided in the inputs, then there will be correspoding columns for the related Jmax, and msg values; for example, Jmax_at_25_lower and Jmax_at_25_lower_msg.

Examples

## Example 1: Estimating Jmax after fitting several C3 A-Ci curves

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data; we will need average values of leaf temperature and
# incident PPFD in order to calculate Jmax later
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp',
    columns_to_average = c('TleafCnd', 'Qin')
)

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# For these examples, we will use a faster (but sometimes less reliable)
# optimizer so they run faster
optimizer <- optimizer_nmkb(1e-7)

# Fit all curves in the data set (it is more common to do this)
aci_results <- consolidate(by(
  licor_file,
  licor_file[, 'species_plot'],
  fit_c3_aci,
  Ca_atmospheric = 420,
  optim_fun = optimizer
))

# Calculate temperature-dependent values of Jmax-related parameters
aci_results$parameters <- calculate_temperature_response(
    aci_results$parameters,
    jmax_temperature_param_bernacchi,
    'TleafCnd_avg'
)

# Calculate Jmax
aci_results$parameters <- calculate_jmax(aci_results$parameters)

# Print a few columns
col_to_view <- c('species_plot', 'J_at_25', 'J_tl_avg', 'Jmax_at_25', 'Jmax_tl')

print(aci_results$parameters[, col_to_view, TRUE])

## Example 2: Illustrating the importance of requiring I2 > J

# Define a data frame with input values
npts <- 200
J_seq <- seq_len(npts)

jmax_df <- data.frame(
  J_at_25 = J_seq,
  J_tl_avg = J_seq,
  alpha_j_norm = 1,
  Qin_avg = 300,
  theta_j_norm = 1,
  TleafCnd_avg = 25
)

# Calculate Jmax values, overriding the default behavior so that values of Jmax
# are returned even when I2 < J.
jmax_df <- calculate_jmax(
  jmax_df, alpha_j_at_25 = 0.293, theta_j_at_25 = 0.979,
  ignore_restriction = TRUE
)

# Plot the Jmax values, distinguishing between cases where J < I2 and where
# J > I2. Here we can see that when J > I2, values of Jmax are smaller than J,
# and can even be negative, which is clearly unreasonable from a biological
# perspective. To highlight these considerations, J = I2 is plotted as a dashed
# black line, Jmax = J is plotted as a black long-dashed line, and Jmax = 0 is
# plotted as a solid black line.
ymin <- -50
ymax <- 250

xmin <- min(J_seq)
xmax <- max(J_seq)

I2 <- jmax_df$I2_at_25[1]

jmax_df$Jmax_at_25_msg[jmax_df$Jmax_at_25_msg == ''] <- 'J < I2'

lattice::xyplot(
  Jmax_at_25 ~ J_at_25,
  group = Jmax_at_25_msg,
  data = jmax_df,
  auto = TRUE,
  type = 'l',
  xlim = c(xmin, xmax),
  ylim = c(ymin, ymax),
  xlab = 'J (micromol / m^2 / s)',
  ylab = 'Jmax (micromol / m^2 / s)',
  panel = function(x, y, ...) {
    lattice::panel.lines(c(0, 0) ~ c(xmin, xmax), lty = 1, col = 'black')
    lattice::panel.lines(c(ymin, ymax) ~ c(I2, I2), lty = 2, col = 'black')
    lattice::panel.lines(J_seq ~ J_seq, lty = 5, col = 'black')
    lattice::panel.xyplot(x, y, ...)
  }
)

Calculate leakiness

Description

Calculates leakiness (phi) from combined gas exchange and isotope discrimination measurements as described in Ubierna et al. (2013). This function can accomodate alternative colum names for the variables taken from exdf_obj; it also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_leakiness_ubierna(
    exdf_obj,
    e = -3,
    a_bar_column_name = 'a_bar',
    a_column_name = 'A',
    ci_column_name = 'Ci',
    co2_s_column_name = 'CO2_s',
    csurface_column_name = 'Csurface',
    delta_c13_r_column_name = 'delta_C13_r',
    delta_obs_tdl_column_name = 'Delta_obs_tdl',
    rl_column_name = 'RL',
    t_column_name = 't'
  )

Arguments

Details

This function uses the model for photosynthetic discrimination against 13C in C4 plants to determine leakiness values, as described in Ubierna et al. (2013). In particular, the following equations from that source are implemented in the code:

• Isotopic fractionation due to day respiration (e_prime) is calculated using Equation 21.

• Leakiness including respiratory and photorespiratory fractionations under high light (phi_i) is calculated using Equation 16.

• Leakiness including respiratory and photorespiratory fractionations and Cs under high light (phi_is) is calculated using Equation 15.

• Leakiness ignoring respiratory and photorespiratory fractionations and Cs (phi_sim) is calculated using Equation 17.

References:

Ubierna, N., Sun, W., Kramer, D. M. and Cousins, A. B. "The efficiency of C4 photosynthesis under low light conditions in Zea mays, Miscanthus x giganteus and Flaveria bidentis." Plant, Cell & Environment 36, 365–381 (2013) [doi:10.1111/j.1365-3040.2012.02579.x].

Value

An exdf object based on exdf_obj that includes the following additional columns, calculated as described above: e_prime, phi_i, phi_is, and phi_sim. The category for each of these new columns is calculate_leakiness_ubierna to indicate that they were created using this function.

Examples

## In this example we load gas exchange and TDL data files, calibrate the TDL
## data, pair the data tables together, and then calculate leakiness. The
## results from this example are not meaningful because these measurements
## were not collected from C4 plants.

# Read the TDL data file, making sure to interpret the time zone as US Central
# time
tdl_data <- read_gasex_file(
  PhotoGEA_example_file_path('tdl_for_gm.dat'),
  'TIMESTAMP',
  list(tz = 'America/Chicago')
)

# Identify cycles within the TDL data
tdl_data <- identify_tdl_cycles(
  tdl_data,
  valve_column_name = 'valve_number',
  cycle_start_valve = 20,
  expected_cycle_length_minutes = 2.7,
  expected_cycle_num_valves = 9,
  timestamp_colname = 'TIMESTAMP'
)

# Use reference tanks to calibrate the TDL data
processed_tdl <- consolidate(by(
  tdl_data,
  tdl_data[, 'cycle_num'],
  process_tdl_cycle_erml,
  noaa_valve = 2,
  calibration_0_valve = 20,
  calibration_1_valve = 21,
  calibration_2_valve = 23,
  calibration_3_valve = 26,
  noaa_cylinder_co2_concentration = 294.996,
  noaa_cylinder_isotope_ratio = -8.40,
  calibration_isotope_ratio = -11.505
))

# Read the gas exchange data, making sure to interpret the time stamp in the US
# Central time zone
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
  'time',
  list(tz = 'America/Chicago')
)

# Get TDL valve information from Licor file name; for this TDL system, the
# reference valve is 12 when the sample valve is 11
licor_data <- get_sample_valve_from_filename(licor_data, list('11' = 12))

# Get oxygen info from the Licor file preamble (needed for calculate_gamma_star)
licor_data <- get_oxygen_from_preamble(licor_data)

# Pair the Licor and TDL data by locating the TDL cycle corresponding to each
# Licor measurement
licor_data <- pair_gasex_and_tdl(licor_data, processed_tdl$tdl_data)

# Calculate total pressure (needed for calculate_gas_properties)
licor_data <- calculate_total_pressure(licor_data)

# Calculate Csurface (needed for calculate_ternary_correction)
licor_data <- calculate_gas_properties(licor_data)

# Calculate ternary correction
licor_data <- calculate_ternary_correction(licor_data)

# Calculate isotope discrimination (needed for calculate_leakiness_ubierna)
licor_data <- calculate_isotope_discrimination(licor_data)

# Set respiration (needed for calculate_leakiness_ubierna)
licor_data <- set_variable(
  licor_data,
  'RL',
  'micromol m^(-2) s^(-1)',
  value = 1.2
)

# Calculate leakiness
licor_data <- calculate_leakiness_ubierna(licor_data)

# View some of the results
licor_data[, c('replicate', 'CO2_s', 'Delta_obs_tdl', 'phi_i', 'phi_sim')]

Calculate temperature-dependent parameter values

Description

Calculate leaf-temperature-dependent values of various parameters using various temperature response functions.

Usage

  calculate_temperature_response(
    exdf_obj,
    temperature_response_parameters,
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

Some key photosynthetic parameters are known to vary with temperature according to well-established temperature response functions such as the Arrhenius equation. The calculate_temperature_response function can be used to calculate such temperature-dependent parameter values at leaf temperature.

Depending on the type value supplied in each element of temperature_response_parameters, one of several possible functions will be used to calculate the temperature response:

• When type is 'Arrhenius', the calculate_temperature_response_arrhenius function will be used.

• When type is 'Gaussian', the calculate_temperature_response_gaussian function will be used.

• When type is 'Johnson', the calculate_temperature_response_johnson function will be used.

• When type is 'Polynomial', the calculate_temperature_response_polynomial function will be used.

Values of type are not case-sensitive.

It is rare to directly specify these parameters; instead, it is more typical to use one of the pre-set values such as those included in c3_temperature_param_sharkey.

Value

An exdf object based on exdf_obj that includes one new column for each element of temperature_response_parameters, where the temperature-dependent values of these new columns are determined using the temperature values specified by the tleaf_column_name column. The category of each of these new columns is calculate_temperature_response to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

# In this example we will calculate temperature-dependent values of two
# parameters:
#
# - The `Kc` parameter (in units of `micromol mol^(-1)`) will be calculated
#   using an Arrhenius function with scaling constant `c` = 38.05 and activation
#   energy `Ea` = 79.43 kJ / mol.
#
# - The `Jmax` parameter (in units of `micromol m^(-2) s^(-1)) will be
#   using a Gaussian function with optimal temperature `t_opt` = 43 degrees C
#   and width `sigma` = 16 degrees C.
#
# So the `temperature_response_parameters` list will contain two elements,
# defined as follows:

trp <- list(
  Kc = list(
    type = 'Arrhenius',
    c = 38.05,
    Ea = 79.43,
    units = 'micromol mol^(-1)'
  ),
  Jmax = list(
    type = 'Gaussian',
    optimum_rate = 4,
    t_opt = 43,
    sigma = 16,
    units = 'micromol m^(-2) s^(-1)'
  )
)

# Now we can calculate the values of Kc and Jmax at the measured leaf
# temperatures recorded in the log file
licor_file <- calculate_temperature_response(licor_file, trp)

licor_file$units$Kc      # View the units of the new `Kc` column
licor_file$categories$Kc # View the category of the new `Kc` column
licor_file[,'Kc']        # View the values of the new `Kc` column

licor_file$units$Jmax      # View the units of the new `Jmax` column
licor_file$categories$Jmax # View the category of the new `Jmax` column
licor_file[,'Jmax']        # View the values of the new `Jmax` column

Calculate temperature-dependent values using Arrhenius equations

Description

Calculate leaf-temperature-dependent values of various parameters using Arrhenius equations. It is rare for users to call this function directly; instead, it is used internally by calculate_temperature_response.

Usage

  calculate_temperature_response_arrhenius(
    exdf_obj,
    arrhenius_parameters,
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

The Arrhenius equation is often used to calculate the temperature dependence of the rate of a chemical reaction. It is often stated as follows:

(1) rate = A * exp(-Ea / (R * T))

where A is the "pre-exponential factor" that sets the overall scaling, Ea is the activation energy, R is the ideal gas constant, and T is the temperature in Kelvin. See, for example, the Wikipedia page for the equation.

In photosynthesis research, it is common to use an alternative form of the equation, where the pre-exponential factor A is rewritten as an exponent A = exp(c), where c is a "scaling factor" whose value can be calculated from A according to c = ln(A)). In this formulation, the equation becomes:

(2) rate = exp(c) * exp(-Ea / (R * T)) = exp(c - Ea / (R * T))

The advantage of this version is that the natural logarithm of the rate is equal to c - Ea / (R * T). This means that the Arrhenius paramerer values can be easily determined from a linear fit of log(rate) against 1 / (R * T); c is the y-intercept and -Ea is the slope.

In calculate_temperature_response_arrhenius, the scaling factor (c), activation energy (Ea), and units (units) for a variable must be specified as elements of a list, which itself is a named element of arrhenius_parameters. For example, if a variable called Kc has c = 38.05, Ea = 79.43, and units of micromol mol^(-1), the arrhenius_parameters argument could be specified as follows: list(Kc = list(c = 38.05, Ea = 79.43, units = 'micromol mol^(-1)')).

It is rare to directly specify the Arrhenius parameters; instead, it is more typical to use one of the pre-set values such as those included in c3_temperature_param_sharkey.

Sometimes a publication will specify the value of a variable at 25 degrees C instead of the Arrhenius scaling factor c. In this case, there is a "trick" for determining the value of c. For example, if the Arrhenius exponent should be X at 25 degrees C, then we have the following: X = exp(c - Ea / (R * (25 + 273.15))), which we can solve algebraically for c as follows: c = ln(X) + Ea / f, where f = R * (25 + 273.15). As a special case, for parameters normalized to 1 at 25 degrees C, we have c = Ea / f. The value of f can be accessed as PhotoGEA:::f.

Another common scenario is that we may wish to convert the units of a variable defined by Arrhenius exponents. For example, let's say Y is determined by an Arrhenius exponent, i.e., that Y = exp(c - Ea / (R * T)), and we want to convert Y to different units via a multiplicative conversion factor cf. Then, in the new units, Y becomes Y_new = cf * Y = cf * exp(c - (R * T)). Through algebra, it is possible to combine cf with the original value of c as c_new = c + ln(cf). Then we can continue calculating Y_new using an Arrhenius factor as Y_new = exp(c_new - Ea / (R * T)).

Value

An exdf object based on exdf_obj that includes one new column for each element of arrhenius_parameters, where the temperature-dependent values of these new columns are determined using the temperature values specified by the tleaf_column_name column. The category of each of these new columns is calculate_temperature_response_arrhenius to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_temperature_response_arrhenius(
  licor_file,
  list(Kc_norm = c3_temperature_param_sharkey$Kc_norm)
)

licor_file$units$Kc_norm      # View the units of the new `Kc_norm` column
licor_file$categories$Kc_norm # View the category of the new `Kc_norm` column
licor_file[,'Kc_norm']        # View the values of the new `Kc_norm` column

Calculate temperature-dependent values using Gaussian equations

Description

Calculate leaf-temperature-dependent values of various parameters using Gaussian equations. It is rare for users to call this function directly; instead, it is used internally by calculate_temperature_response.

Usage

  calculate_temperature_response_gaussian(
    exdf_obj,
    gaussian_parameters,
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

A Gaussian equation is sometimes used to model the temperature dependence of a biochemical rate parameter. Typically this is expressed by

rate = optimal_rate * exp(-(T - T_opt)^2 / sigma^2)

where optimal_rate is the highest rate which occurs at the optimal temperature T_opt, T is the current temperature, and sigma represents the "width" of the peak. More technically, it can be described as the difference in temperature away from the optimal value at which the rate falls to 37 percent (1/e) of its maximum.

In calculate_temperature_response_gaussian, the optimal rate (optimal_rate), optimal temperature (t_opt), width (sigma), and units (units) for a variable must be specified as elements of a list, which itself is a named element of gaussian_parameters. For example, if a variable called Jmax has optimal_rate = 1, t_opt = 43, sigma = 26, and units of micromol mol^(-1), the gaussian_parameters argument could be specified as follows: list(Jmax = list(optimal_rate = 1, t_opt = 43, sigma = 26, units = 'micromol mol^(-1)')).

It is rare to specify these parameters directly; instead, it is more typical to use one of the pre-set values such as those included in c4_temperature_param_vc.

Value

An exdf object based on exdf_obj that includes one new column for each element of gaussian_parameters, where the temperature-dependent values of these new columns are determined using the temperature values specified by the tleaf_column_name column. The category of each of these new columns is calculate_temperature_response_gaussian to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_temperature_response_gaussian(
  licor_file,
  list(J_norm = c4_temperature_param_vc$J_norm)
)

licor_file$units$J_norm      # View the units of the new `J_norm` column
licor_file$categories$J_norm # View the category of the new `J_norm` column
licor_file[,'J_norm']        # View the values of the new `J_norm` column

Calculate temperature-dependent values using Johnson-Eyring-Williams equations

Description

Calculate leaf-temperature-dependent values of various parameters using Johnson-Eyring-Williams equations. It is rare for users to call this function directly; instead, it is used internally by calculate_temperature_response.

Usage

  calculate_temperature_response_johnson(
    exdf_obj,
    johnson_parameters,
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

The Johnson-Eyring-Williams equation is often used to calculate the temperature dependence of the rate of a chemical reaction. It can be stated as follows:

rate = exp(c - Ha / (R * T)) / (1 + exp(S / R - Hd / (R * T)))

where c is the scaling factor that sets the overall magnitude of the rate, Ha is the enthalpy of activation, Hd is the enthalpy of deactivation, S is the entropy, R is the ideal gas constant, and T is the temperature in Kelvin.

This equation exhibits a peak; in other words, there is a particular temperature (the optimal temperature) where the rate is maximized. Thus, it is often used in place of an Arrhenius equation (see calculate_temperature_response_arrhenius) for photosynthetic parameters that exhibit a decrease at high temperatures.

This equation was originally published in Johnson, Eyring, & Williams (1942) and has been used to model the temperature dependence of key photosynthetic parameters, as in Harley et al. (1992), Bernacchi et al. (2003), Sharkey et al. (2007), and others.

In calculate_temperature_response_johnson, the scaling factor (c), enthalpy of activation (Ha), enthalpy of deactivation (Hd), entopy (S), and units (units) for a variable must be specified as elements of a list, which itself is a named element of johnson_parameters. For example, if a variable called Tp has c = 21.46, Ha = 53.1, Hd = 201.8, S = 0.65, and units of micromol mol^(-1), the johnson_parameters argument could be specified as follows: list(Tp = list(c = 21.46, Ha = 53.1, Hd = 201.8, S = 0.65, units = 'micromol mol^(-1)')).

It is rare to directly specify these parameters; instead, it is more typical to use one of the pre-set values such as those included in c3_temperature_param_sharkey.

References:

• Johnson, F. H., Eyring, H. & Williams, R. W. "The nature of enzyme inhibitions in bacterial luminescence: Sulfanilamide, urethane, temperature and pressure." Journal of Cellular and Comparative Physiology 20, 247–268 (1942) [doi:10.1002/jcp.1030200302].

• Harley, P. C., Thomas, R. B., Reynolds, J. F. & Strain, B. R. "Modelling photosynthesis of cotton grown in elevated CO2." Plant, Cell & Environment 15, 271–282 (1992) [doi:10.1111/j.1365-3040.1992.tb00974.x].

• Bernacchi, C. J., Pimentel, C. & Long, S. P. "In vivo temperature response functions of parameters required to model RuBP-limited photosynthesis." Plant, Cell & Environment 26, 1419–1430 (2003) [doi:10.1046/j.0016-8025.2003.01050.x].

• Sharkey, T. D., Bernacchi, C. J., Farquhar, G. D. & Singsaas, E. L. "Fitting photosynthetic carbon dioxide response curves for C3 leaves." Plant, Cell & Environment 30, 1035–1040 (2007) [doi:10.1111/j.1365-3040.2007.01710.x].

Value

An exdf object based on exdf_obj that includes one new column for each element of johnson_parameters, where the temperature-dependent values of these new columns are determined using the temperature values specified by the tleaf_column_name column. The category of each of these new columns is calculate_temperature_response_johnson to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_temperature_response_johnson(
  licor_file,
  list(Tp_norm = c3_temperature_param_sharkey$Tp_norm)
)

licor_file$units$Tp_norm      # View the units of the new `Tp_norm` column
licor_file$categories$Tp_norm # View the category of the new `Tp_norm` column
licor_file[,'Tp_norm']        # View the values of the new `Tp_norm` column

Calculate temperature-dependent values using polynomial equations

Description

Calculate leaf-temperature-dependent values of various parameters using polynomial equations. It is rare for users to call this function directly; instead, it is used internally by calculate_temperature_response.

Usage

  calculate_temperature_response_polynomial(
    exdf_obj,
    polynomial_parameters,
    tleaf_column_name = 'TleafCnd'
  )

Arguments

Details

Polynomial equations are often used to calculate the temperature dependence of the rates of chemical reactions. For example, a second-order polynomial could be given as follows:

(1) rate = R_0 + R_1 * T + R_2 * T^2

where R_0, R_1, and R_2 are the zeroth, first, and second order coefficients and T is the temperature. Higher order polynomials can also be defined, where an order-N polynomial is given by

(2) rate = R_0 + R_1 * T + R_2 * T^2 + ... + R_N * T^N

In general, an order-N polynomial has N coefficients, although some of them may be zero.

In calculate_temperature_response_polynomial, the coefficients (coef) and units (units) for a variable must be specified as elements of a list, which itself is a named element of polynomial_parameters. The coefficients must be specified as a numeric vector, where the ith element represents the ith coefficient. For example, if a dimensionless variable called theta is calculated according to theta = 0.352 + 0.022 * T - 3.4e-4 * T^2, the polynomial_parameters argument could be supplied as follows: list(theta = list(coef = c(0.352, 0.022, -3.4e-4), units = 'dimensionless')).

It is rare to directly specify the polynomial parameters; instead, it is more typical to use one of the pre-set values such as those included in jmax_temperature_param_bernacchi.

Value

An exdf object based on exdf_obj that includes one new column for each element of polynomial_parameters, where the temperature-dependent values of these new columns are determined using the temperature values specified by the tleaf_column_name column. The category of each of these new columns is calculate_temperature_response_polynomial to indicate that they were created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_temperature_response_polynomial(
  licor_file,
  list(theta = list(coef = c(0.352, 0.022, -3.4e-4), units = 'dimensionless'))
)

licor_file$units$theta      # View the units of the new `theta` column
licor_file$categories$theta # View the category of the new `theta` column
licor_file[,'theta']        # View the values of the new `theta` column

Calculate ternary correction factor

Description

Calculates the ternary correction factor t that is used in many carbon isotope discrimination calculations.

Usage

  calculate_ternary_correction(
    exdf_obj,
    ci_column_name = 'Ci',
    co2_s_column_name = 'CO2_s',
    csurface_column_name = 'Csurface',
    e_column_name = 'E',
    gtc_column_name = 'gtc'
  )

Arguments

Details

During photosynthetic gas exchange, there are separate fluxes of CO2 and H2O flowing in and out of the leaf. These gases interact with each other and with air, forming a ternary mixture. These interactions must be taken into account when modeling carbon isotope discrimination. Typically this is done via t, a ternary correction factor first introduced by Farquhar and Cernusak (2012). Here we calculate t as described in Equations 9 and 10 from Ubierna et al. (2018):

t = alpha_ac * E / (2 * g_ac)

and

a_bar = (a_b * (C_a - C_s) + a_s * (C_s - C_i)) / (C_a - C_i),

where E is the transpiration rate, g_ac is the total conductance to CO2 diffusion across the boundary layer and stomata in series, a_bar is the weighted fractionation across the boundary layer and stomata in series, a_b is the fractionation during diffusion through the boundary layer, a_s is the fractionation during diffusion through the stomata, C_a is the ambient CO2 concentration (in wet air), C_s is the CO2 concentration (in wet air) at the leaf surface, and C_i is the CO2 concentration (in wet air) in the intercellular spaces.

alpha_ac is the overall fractionation during diffusion through air; alpha_ac and a_bar are related according to an un-numbered equation in Ubierna et al. (2018) that appears just after Equation 9:

alpha_ac = 1 + a_bar

References:

Farquhar, G. D. and Cernusak, L. A. "Ternary effects on the gas exchange of isotopologues of carbon dioxide." Plant, Cell & Environment 35, 1221–1231 (2012) [doi:10.1111/j.1365-3040.2012.02484.x].

Ubierna, N., Holloway-Phillips, M.-M. and Farquhar, G. D. "Using Stable Carbon Isotopes to Study C3 and C4 Photosynthesis: Models and Calculations." in Photosynthesis: Methods and Protocols (ed. Covshoff, S.) 155–196 (Springer, 2018) [doi:10.1007/978-1-4939-7786-4_10].

Value

An exdf object based on exdf_obj that includes values of t, a_bar, and alpha_ac calculated as described above. The category of each new column is calculate_ternary_correction to indicate that it was created using this function.

Examples

## In this example we load a gas exchange data file and then calculate the
## ternary correction factor

# Read the gas exchange data
licor_data <- read_gasex_file(
  PhotoGEA_example_file_path('licor_for_gm_site11.xlsx'),
  'time'
)

# Calculate total pressure (needed for calculate_gas_properties)
licor_data <- calculate_total_pressure(licor_data)

# Calculate Csurface (needed for calculate_ternary_correction)
licor_data <- calculate_gas_properties(licor_data)

# Calculate ternary correction
licor_data <- calculate_ternary_correction(licor_data)

# View some of the results
licor_data[, c('replicate', 'A', 'E', 'Csurface', 't', 'a_bar', 'alpha_ac')]

Calculate the total pressure in bar

Description

Calculates the total pressure in bar. Licor gas exchange measurement systems report both the abient air pressure (Pa) and the chamber overpressure (DeltaPcham) in kPa; the total pressure in the chamber is therefore given by the sum of these two columns. This function can accomodate alternative column names for the variables taken from Licor log files in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_total_pressure(
    exdf_obj,
    pa_column_name = 'Pa',
    deltapcham_column_name = 'DeltaPcham'
  )

Arguments

Details

If deltapcham_column_name is NA, this function will simply convert the values in the pa_column_name to units of bar. Otherwise, the values from the pa_column_name and deltapcham_column_name columns will be added together and converted to bar.

Value

An exdf object based on exdf_obj that includes the total pressure values in a new column called total_pressure. The category of this new column is calculate_total_pressure to indicate that it was created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package and calculate the
# total pressure.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_total_pressure(licor_file)

licor_file$units$total_pressure      # View the units of the new `total_pressure` column
licor_file$categories$total_pressure # View the category of the new `total_pressure` column
licor_file[, 'total_pressure']       # View the values of the new `total_pressure` column

Calculate intrinsic water use efficiency

Description

Calculates the intrinsic water use efficiency (iWUE). This function can accomodate alternative column names for the variables taken from the data file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  calculate_wue(
    exdf_obj,
    calculate_c3 = FALSE,
    a_column_name = 'A',
    ca_column_name = 'Ca',
    cc_column_name = 'Cc',
    ci_column_name = 'Ci',
    e_column_name = 'E',
    gmc_column_name = 'gmc_tl',
    gsw_column_name = 'gsw',
    h2o_a_column_name = 'H2O_s',
    h2o_i_column_name = 'H2O_i',
    total_pressure_column_name = 'total_pressure'
  )

Arguments

Details

Leaf-level water use efficiency (lWUE) is defined as the ratio of net CO2 assimilation (An) to transpiration (E):

lWUE = An / E.

This quantity can also be expressed in terms of water and CO2 concentrations:

lWUE = 0.6 * Ca * (1 - Ci / Ca) / (H2Oi - H2Oa).

Here, Ca and Ci are the atmospheric and intercellular CO2 concentrations, and H2Oa and H2Oi are the atmospheric and intercellular water vapor concentrations. If differences in lWUE are measured between different groups of plants, it can be helpful to separately investigate Ci / Ca and H2Oi - H2Oa to see which factor is driving the differences.

The intrinsic water use efficiency iWUE is a measure of leaf-level water use efficiency, and it is defined to be the ratio An and the stomatal conductance to H2O diffusion (gsw):

iWUE = An / gsw.

For C3 plants, iWUE can be reexpressed as

iWUE = (gmc / gsw) / (1 + (gmc / gsw)) * (Ca - Cc),

where gmc is the mesophyll conductance to CO2 diffusion and Cc is the chloroplast CO2 concentration. If differences in iWUE are measured between different groups of plants, it can be helpful to separately investigate gmc / gsw and Ca - Cc to see which factor is driving the differences.

Note: both measures of water use efficiency depend directly or indirectly on stomatal conductance. Stomata are notoriously slow to reach steady-state, but water use efficiency is only reliable at steady-state. For this reason, it is recommended to only analyze water use efficiency for gas exchange measurements where stomatal conductance has stabilized. For an A-Ci or A-Q curve, only the first measured point has typically reached steady-state stomatal conductance. On the other hand, for a Ball-Berry curve, all measured points should have reached steady-state stomatal conductance.

For more details about these quantities, see Leakey et al. "Water Use Efficiency as a Constraint and Target for Improving the Resilience and Productivity of C3 and C4 Crops." Annual Review of Plant Biology 70 (1): 781–808 (2019) [doi:10.1146/annurev-arplant-042817-040305].

In this function, the following variables are calculated:

• lWUE, given by iWUE = An / E

• Cia_ratio, given by Cia_ratio = Ci / Ca

• drawdown_sw, given by drawdown_sw = H2Oi - H2Oa (this is the drawdown of water vapor across the stomata)

• iWUE, given by iWUE = An / gsw

• g_ratio, given by g_ratio = gmc / gsw

• drawdown_ct, given by drawdown_ct = Ca - Cc (this is the total drawdown of CO2 from the ambient air to the chloroplast)

Note: g_ratio and drawdown_ct are only calculated if calculate_c3 is TRUE.

Value

An exdf object based on exdf_obj that includes the quantities listed above, along with their units. The category of each of these new columns is calculate_wue to indicate that it was created using this function.

Examples

# Read an example Licor file included in the PhotoGEA package and calculate the
# water use efficiency.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

licor_file <- calculate_total_pressure(licor_file)

licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_sharkey)

licor_file <- calculate_gas_properties(licor_file)

licor_file <- apply_gm(licor_file, gmc_at_25 = 0.5)

licor_file <- calculate_wue(licor_file, calculate_c3 = TRUE)

licor_file$units$iWUE      # View the units of the new `iWUE` column
licor_file$categories$iWUE # View the category of the new `iWUE` column
licor_file[, 'iWUE']       # View the values of the new `iWUE` column

Combine exdf objects by columns or rows

Description

Combines one or more exdf objects by the columns or rows of their main_data. For rbind, errors will occur if column names are not the same in all of the exdf objects, and if all units and categories are not identical.

Usage

  ## S3 method for class 'exdf'
cbind(..., deparse.level = 1)

  ## S3 method for class 'exdf'
rbind(
    ...,
    deparse.level = 1,
    make.row.names = TRUE,
    stringsAsFactors = FALSE
  )

Arguments

Value

Returns a new exdf object.

See Also

exdf

Examples

# Make some simple exdf objects. 1 and 2 have the same number of rows but
# different columns, while 1 and 3 have the same columns but different rows.
simple_exdf_1 <- exdf(data.frame(A = 1), data.frame(A = 'au'), data.frame(A = 'ac'))
simple_exdf_2 <- exdf(data.frame(B = 2), data.frame(B = 'bu'), data.frame(B = 'bc'))
simple_exdf_3 <- exdf(data.frame(A = 2), data.frame(A = 'au'), data.frame(A = 'ac'))

cbind(simple_exdf_1) # will just return simple_exdf_1
cbind(simple_exdf_1, simple_exdf_2)

rbind(simple_exdf_1) # will just return simple_exdf_1
rbind(simple_exdf_1, simple_exdf_3)

Make sure required variables exist

Description

Checks whether the input table has the required variables.

Usage

  check_required_variables(x, required_variables, check_NA = TRUE)

  ## S3 method for class 'data.frame'
check_required_variables(x, required_variables, check_NA = TRUE)

  ## S3 method for class 'exdf'
check_required_variables(x, required_variables, check_NA = TRUE)

Arguments

Details

check_required_variables is generic, with methods defined for data frames and exdf objects.

When x is an exdf, the required_variables input argument must be a list of named strings, where the name of each element specifies the name of a column that must be included in x, while the value of each column specifies the corresponding units for that column. If the value is NA, no unit checking will be performed.

When x is a data.frame, the required_variables input argument can be specified as a list (as if x were an exdf object) or as a character vector specifying the names of columns that should be included in x.

The required variables will be checked as follows:

• If any required variable columns are missing from the table, an informative error message will be thrown.

• If check_NA is TRUE and any required variable columns are entirely NA, an informative error message will be thrown.

• If any required variable colums have incorrect units, an informative error message will be thrown. (Only applies to exdf objects.)

Otherwise, check_required_variables will have no output and produce no messages.

This function is used internally by many other functions from the PhotoGEA package to check for important columns and make sure they have the correct units. For example, see the code for apply_gm by typing PhotoGEA::apply_gm in the R terminal.

Value

The check_required_variables function does not return anything.

See Also

exdf

Examples

# Create a simple exdf object
simple_exdf <- exdf(
  data.frame(A = c(3, 2, 7, 9), B = c(4, 5, 1, 8)),
  data.frame(A = 'm', B = 's', stringsAsFactors = FALSE),
  data.frame(A = 'Cat1', B = 'Cat2', stringsAsFactors = FALSE)
)

# Confirm that columns named `A` and `B` are in the object, and that they have
# units of `m` and `s`, respectively.
check_required_variables(simple_exdf, list(A = 'm', B = 's'))

# Confirm that columns named `A` and `B` are in the object, but only check units
# for the `A` column.
check_required_variables(simple_exdf, list(A = 'm', B = NA))

# Use the data frame method on `simple_exdf$main_data` to confirm that columns
# named `A` and `B` are present
check_required_variables(simple_exdf$main_data, c('A', 'B'))

Check response curve data for common issues

Description

Checks to make sure an exdf object representing multiple response curves meets basic expectations.

Usage

  check_response_curve_data(
    exdf_obj,
    identifier_columns,
    expected_npts = 0,
    driving_column = NULL,
    driving_column_tolerance = 1.0,
    col_to_ignore_for_inf = 'gmc',
    constant_col = list(),
    error_on_failure = TRUE,
    print_information = TRUE
  )

Arguments

Details

Basic Behavior:

This function makes a few basic checks to ensure that the response curve data includes the expected information and does not include any mistakes. If no problems are detected, this function will be silent with no return value. If a problem is detected, then the user will be notified in one or more ways:

• If error_on_failure is TRUE, then this function will throw an error with a short message. If print_information is also TRUE, then additional information will be printed to the R terminal.

• If error_on_failure is FALSE and print_information is also FALSE, then this function will throw a warning with a short message.

• If error_on_failure is FALSE and print_information is true, information about the problem will be printed to the R terminal.

This function will (optionally) perform several checks:

• Checking for infinite values: If col_to_ignore_for_inf is not NULL, no numeric columns in exdf_obj should have infinite values, with the exception of columns designated in col_to_ignore_for_inf.

• Checking required columns: All elements of identifier_columns should be present as columns in exdf_obj. If driving_column is not NULL, it should also be present as a column in exdf_obj. If constant_col is not empty, then these columns must also be present in exdf_obj.

• Checking the number of points in each curve: The general idea is to ensure that each curve has the expected number of points. Several options can be specified via the value of expected_npts, as discussed below.

• Checking the driving column: If driving_column is not NULL, then each curve should have the same sequence of values in this column. To allow for small variations, a nonzero driving_column_tolerance can be specified.

• Checking the constant columns: If constant_col is not empty, then each specified column should either be constant, or only vary by a specified amount. See details below.

By default, most of these are not performed (except the simplest ones like checking for infinite values or checking that key columns are present). This enables an "opt-in" use style, where users can specify just the checks they wish to make.

More Details:

There are several options for checking the number of points in each curve:

• If expected_npts is a single negative number, no check will be performed.

• If expected_npts is 0, then each curve is expected to have the same number of points.

• If expected_npts is a single positive number, then each curve is expected to have that many points. For example, if expected_npts is 7, then each curve must have 7 points.

• If expected_npts is a pair of positive numbers, then each curve is expected to have a number of points lying within the range defined by expected_npts. For example, if expected_npts is c(6, 8), then each curve must have no fewer than 6 points and no more than 8 points.

• If expected_npts is a pair of numbers, one of which is zero and one of which is positive, then the positive number specifies a range; each curve must differ from the average number of points by less than the range. For example, if expected_npts is c(0, 3), then every curve must have a number of points within 3 of the average number of points.

There are two options for checking columns that should be constant:

• A value of NA indicates that all values of that column must be exactly identical; this check applies for numeric and character columns.

• A numeric value indicates that the range of values of that column must be smaller than the specified range; this range applies for numeric columns only.

For example, setting constant_col = list(species = NA, Qin = 10) means that each curve must have only a single value of the species column, and that the value of the Qin column cannot vary by more than 10 across each curve.

Use Cases:

Using check_response_curve_data is not strictly necessary, but it can be helpful both to you and to anyone else reading your analysis code. Here are a few typical use cases:

• Average response curves: It is common to calculate and plot average response curves, either manually or by using xyplot_avg_rc. But, it only makes sense to do this if each curve followed the same sequence of the driving variable. In this case, check_response_curve_data can be used to confirm that each curve used the same values of CO2_r_sp (for an A-Ci curve) or Qin (for an A-Q curve).

• Removing recovery points: It is common to wish to remove one or more recovery points from a set of curves. The safest way to do this is to confirm that all the curves use the same sequence of setpoints; then you can be sure that, for example, points 9 and 10 are the recovery points in every curve.

• Making a statement of expectations: If you measured a set of A-Ci curves where each curve has 16 points and used the same sequence of CO2_r setpoints, you could record this somewhere in your notes. But it would be even more meaningful to use check_response_curve_data in your script with expected_npts set to 16. If this check passes, then it means not only that your claim is correct, but also that the identifier columns are being interpreted properly.

• Checking identifiers: If the data set includes some identifying metadata, such as a species or location, it may be helpful to confirm that each curve has a single value of these "identifier" columns. Otherwise, the data set may be difficult to interpret.

• Checking measurement conditions: If the response curves are expected to be measured under constant temperature, humidity, light, or other environmental variables, it may be helpful to confirm that these variables do not vary too much across each individual curve. Otherwise, parameter values estimated from the curves may not be meaningful.

Sometimes the response curves in a large data set were not all measured with the same sequence of setpoints. If only a few different sequences were used, it is possible to split them into groups and separately run check_response_curve_data on each group. This scenario is discussed in the Frequently Asked Questions vignette.

Even if none of the above situations are relevant to you, it may still be helpful to run run check_response_curve_data but with expected_npts set to 0 and error_on_failure set to FALSE. With these settings, if there are curves with different numbers of points, the function will print the number of points in each curve to the R terminal, but won't stop the rest of the script from running. This can be useful for detecting problems with the curve_identifier column. For example, if the longest curves in the set are known to have 17 points, but check_response_curve_data identifies a curve with 34 points, it is clear that the same identifier was accidentally used for two different curves.

Value

The check_response_curve_data function does not return anything.

Examples

# Read an example Licor file included in the PhotoGEA package and check it.
# This file includes several 7-point light-response curves that can be uniquely
# identified by the values of its 'species' and 'plot' columns. Since these are
# light-response curves, each one follows a pre-set sequence of `Qin` values.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

# Make sure there are no infinite values and that all curves have the same
# number of points
check_response_curve_data(licor_file, c('species', 'plot'))

# Make sure there are no inifinite values and that all curves have 7 points
check_response_curve_data(licor_file, c('species', 'plot'), expected_npts = 7)

# Make sure there are no infinite values, that all curves have 7 points, and
# that the values of the `Qin` column follow the same sequence in all curves
# (to within 1.0 micromol / m^2 / s)
check_response_curve_data(
  licor_file,
  c('species', 'plot'),
  expected_npts = 7,
  driving_column = 'Qin',
  driving_column_tolerance = 1.0
)

# Make sure that there are no infinite values and that all curves have between
# 8 and 10 points; this will intentionally fail
check_response_curve_data(
  licor_file,
  c('species', 'plot'),
  expected_npts = c(8, 10),
  error_on_failure = FALSE
)

# Split the data set by `species` and make sure all curves have similar numbers
# of points (within 3 of the mean value); this will intentionally fail
check_response_curve_data(
  licor_file,
  'species',
  expected_npts = c(0, 3),
  error_on_failure = FALSE
)

# Split the data set by `species` and make sure all curves have a constant value
# of `plot` and a limited range of `TLeafCnd`; this will intentionally fail
check_response_curve_data(
  licor_file,
  'species',
  constant_col = list(plot = NA, TleafCnd = 0.001),
  error_on_failure = FALSE
)

Choosing input files

Description

Tools for choosing input files via dialog windows.

Usage

  choose_input_files()

  choose_input_licor_files()

  choose_input_tdl_files()

Details

These functions are only available in interactive sessions; moreover, choose_input_licor_files and choose_input_tdl_files are only available in Microsoft Windows.

• choose_input_files will prompt the user to select a single file, and will return full file paths for all files in the same directory that have the same extension.

• choose_input_licor_files can be used to select one or more Microsoft Excel files (with extension *.xlsx) or plaintext files (with no extension).

• choose_input_tdl_files can be used to select one or more TDL data files (with extension *.dat).

The outputs from these functions are typically passed to read_gasex_file via lapply.

Value

A character vector of full file paths.

Examples

# Interactively select a single file and get full file paths to all
# other files in the same directory that have the same extension
if (interactive()) {
  file_paths <- choose_input_files()
}


# Interactively select one or more Licor Excel files and read each one to create
# a list of exdf objects
if (interactive() && .Platform$OS.type == "windows") {
  lapply(choose_input_licor_files(), function(fname) {
    read_gasex_file(fname, 'time')
  })
}

# Interactively select one or more TDL data files and read each one to create a
# list of exdf objects
if (interactive() && .Platform$OS.type == "windows") {
  lapply(choose_input_tdl_files(), function(fname) {
    read_gasex_file(fname, 'TIMESTAMP')
  })
}

Calculate confidence intervals for C3 A-Ci fitting parameters

Description

Calculates confidence intervals for parameters estimated by a C3 A-Ci curve fit. It is rare for users to call this function directly, because it can be automatically applied to each curve when calling fit_c3_aci.

Usage

  confidence_intervals_c3_aci(
    replicate_exdf,
    best_fit_parameters,
    lower = list(),
    upper = list(),
    fit_options = list(),
    sd_A = 1,
    relative_likelihood_threshold = 0.147,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    gmc_norm_column_name = 'gmc_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    cj_crossover_min = NA,
    cj_crossover_max = NA,
    hard_constraints = 0,
    ...
  )

Arguments

Details

In maximum likelihood fitting, each set of parameter values has an associated likelihood value. If the maximum likelihood is known, then it is also possible to define a relative likelihood p according to p = L / L_max. The set of all parameter values where p exceeds a threshold value p_0 defines a region in parameter space called like a "relative likelihood region." When taking one-dimensional cuts through parameter space, the boundaries of the relative likelihood region define a relative likelihood interval.

Here we calculate the upper and lower limits of the relative likelihood intervals for each parameter. This is done by fixing the other parameters to their best-fit values, and varying a single parameter to find the interval where the relative likelihood is above the threshold value (set by the relative_likelihood_threshold input argument). If the threshold is set to 0.147, then these intervals are equivalent to 95% confidence intervals in most situations. See the Wikipedia page about relative likelihood for more information.

Internally, this function uses error_function_c3_aci to calculate the negative logarithm of the likelihood (-ln(L)). It varies each fitting parameter independendently to find values where ln(L) - ln(p_0) - ln(L_max) = 0.

If the upper limit of a confidence interval is found to exceed ten times the upper limit specified when fitting that parameter, then the upper limit of the condfidence interval is taken to be infinity.

Value

An exdf object based on best_fit_parameters that contains lower and upper bounds for each parameter; for example, if Vcmax_at_25 was fit, best_fit_parameters will contain new columns called Vcmax_at_25_lower and Vcmax_at_25_upper.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Fit just one curve from the data set
one_result <- fit_c3_aci(
  licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE],
  calculate_confidence_intervals = FALSE
)

# Calculate confidence limits for the fit parameters
parameters_with_limits <- confidence_intervals_c3_aci(
    licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE],
    one_result$parameters
)

# View confidence limits and best estimate for Vcmax_at_25
parameters_with_limits[, c('Vcmax_at_25_lower', 'Vcmax_at_25', 'Vcmax_at_25_upper')]

Calculate confidence intervals for C3 Variable J fitting parameters

Description

Calculates confidence intervals for parameters estimated by a C3 A-Ci curve fit. It is rare for users to call this function directly, because it can be automatically applied to each curve when calling fit_c3_variable_j.

Usage

  confidence_intervals_c3_variable_j(
    replicate_exdf,
    best_fit_parameters,
    lower = list(),
    upper = list(),
    fit_options = list(),
    sd_A = 1,
    relative_likelihood_threshold = 0.147,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    phips2_column_name = 'PhiPS2',
    qin_column_name = 'Qin',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    cj_crossover_min = NA,
    cj_crossover_max = NA,
    hard_constraints = 0,
    require_positive_gmc = 'positive_a',
    gmc_max = Inf,
    check_j = TRUE,
    ...
  )

Arguments

Details

In maximum likelihood fitting, each set of parameter values has an associated likelihood value. If the maximum likelihood is known, then it is also possible to define a relative likelihood p according to p = L / L_max. The set of all parameter values where p exceeds a threshold value p_0 defines a region in parameter space called like a "relative likelihood region." When taking one-dimensional cuts through parameter space, the boundaries of the relative likelihood region define a relative likelihood interval.

Here we calculate the upper and lower limits of the relative likelihood intervals for each parameter. This is done by fixing the other parameters to their best-fit values, and varying a single parameter to find the interval where the relative likelihood is above the threshold value (set by the relative_likelihood_threshold input argument). If the threshold is set to 0.147, then these intervals are equivalent to 95% confidence intervals in most situations. See the Wikipedia page about relative likelihood for more information.

Internally, this function uses error_function_c3_variable_j to calculate the negative logarithm of the likelihood (-ln(L)). It varies each fitting parameter independendently to find values where ln(L) - ln(p_0) - ln(L_max) = 0.

If the upper limit of a confidence interval is found to exceed ten times the upper limit specified when fitting that parameter, then the upper limit of the condfidence interval is taken to be infinity.

Value

An exdf object based on best_fit_parameters that contains lower and upper bounds for each parameter; for example, if Vcmax_at_25 was fit, best_fit_parameters will contain new columns called Vcmax_at_25_lower and Vcmax_at_25_upper.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Fit just one curve from the data set, using a less reliable optimizer so the
# example runs faster
one_result <- fit_c3_variable_j(
  licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE],
  optim_fun = optimizer_nmkb(1e-7),
  calculate_confidence_intervals = FALSE
)

# Calculate confidence limits for the fit parameters
parameters_with_limits <- confidence_intervals_c3_variable_j(
    licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE],
    one_result$parameters
)

# View confidence limits and best estimate for Vcmax_at_25
parameters_with_limits[, c('Vcmax_at_25_lower', 'Vcmax_at_25', 'Vcmax_at_25_upper')]

Calculate confidence intervals for C4 A-Ci fitting parameters

Description

Calculates confidence intervals for parameters estimated by a C4 A-Ci curve fit. It is rare for users to call this function directly, because it can be automatically applied to each curve when calling fit_c4_aci.

Usage

  confidence_intervals_c4_aci(
    replicate_exdf,
    best_fit_parameters,
    lower = list(),
    upper = list(),
    fit_options = list(),
    sd_A = 1,
    relative_likelihood_threshold = 0.147,
    x_etr = 0.4,
    a_column_name = 'A',
    ao_column_name = 'ao',
    ci_column_name = 'Ci',
    gamma_star_column_name = 'gamma_star',
    gmc_norm_column_name = 'gmc_norm',
    j_norm_column_name = 'J_norm',
    kc_column_name = 'Kc',
    ko_column_name = 'Ko',
    kp_column_name = 'Kp',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    vcmax_norm_column_name = 'Vcmax_norm',
    vpmax_norm_column_name = 'Vpmax_norm',
    hard_constraints = 0
  )

Arguments

Details

In maximum likelihood fitting, each set of parameter values has an associated likelihood value. If the maximum likelihood is known, then it is also possible to define a relative likelihood p according to p = L / L_max. The set of all parameter values where p exceeds a threshold value p_0 defines a region in parameter space called like a "relative likelihood region." When taking one-dimensional cuts through parameter space, the boundaries of the relative likelihood region define a relative likelihood interval.

Here we calculate the upper and lower limits of the relative likelihood intervals for each parameter. This is done by fixing the other parameters to their best-fit values, and varying a single parameter to find the interval where the relative likelihood is above the threshold value (set by the relative_likelihood_threshold input argument). If the threshold is set to 0.147, then these intervals are equivalent to 95% confidence intervals in most situations. See the Wikipedia page about relative likelihood for more information.

Internally, this function uses error_function_c4_aci to calculate the negative logarithm of the likelihood (-ln(L)). It varies each fitting parameter independendently to find values where ln(L) - ln(p_0) - ln(L_max) = 0.

If the upper limit of a confidence interval is found to exceed ten times the upper limit specified when fitting that parameter, then the upper limit of the condfidence interval is taken to be infinity.

Value

An exdf object based on best_fit_parameters that contains lower and upper bounds for each parameter; for example, if Vcmax_at_25 was fit, best_fit_parameters will contain new columns called Vcmax_at_25_lower and Vcmax_at_25_upper.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c4_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Calculate temperature-dependent values of C4 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c4_temperature_param_vc)

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

# Fit just one curve from the data set
one_result <- fit_c4_aci(
  licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE],
  calculate_confidence_intervals = FALSE
)

# Calculate confidence limits for the fit parameters
parameters_with_limits <- confidence_intervals_c4_aci(
    licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE],
    one_result$parameters
)

# View confidence limits and best estimate for Vcmax_at_25
parameters_with_limits[, c('Vcmax_at_25_lower', 'Vcmax_at_25', 'Vcmax_at_25_upper')]

Calculate confidence intervals for C4 A-Ci hyperbola fitting parameters

Description

Calculates confidence intervals for parameters estimated by a C4 A-Ci curve fit. It is rare for users to call this function directly, because it can be automatically applied to each curve when calling fit_c4_aci_hyperbola.

Usage

  confidence_intervals_c4_aci_hyperbola(
    replicate_exdf,
    best_fit_parameters,
    lower = list(),
    upper = list(),
    fit_options = list(),
    sd_A = 1,
    relative_likelihood_threshold = 0.147,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    hard_constraints = 0
  )

Arguments

Details

In maximum likelihood fitting, each set of parameter values has an associated likelihood value. If the maximum likelihood is known, then it is also possible to define a relative likelihood p according to p = L / L_max. The set of all parameter values where p exceeds a threshold value p_0 defines a region in parameter space called like a "relative likelihood region." When taking one-dimensional cuts through parameter space, the boundaries of the relative likelihood region define a relative likelihood interval.

Here we calculate the upper and lower limits of the relative likelihood intervals for each parameter. This is done by fixing the other parameters to their best-fit values, and varying a single parameter to find the interval where the relative likelihood is above the threshold value (set by the relative_likelihood_threshold input argument). If the threshold is set to 0.147, then these intervals are equivalent to 95% confidence intervals in most situations. See the Wikipedia page about relative likelihood for more information.

Internally, this function uses error_function_c4_aci_hyperbola to calculate the negative logarithm of the likelihood (-ln(L)). It varies each fitting parameter independendently to find values where ln(L) - ln(p_0) - ln(L_max) = 0.

If the upper limit of a confidence interval is found to exceed ten times the upper limit specified when fitting that parameter, then the upper limit of the condfidence interval is taken to be infinity.

Value

An exdf object based on best_fit_parameters that contains lower and upper bounds for each parameter; for example, if Vmax was fit, best_fit_parameters will contain new columns called Vmax_lower and Vmax_upper.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c4_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Fit just one curve from the data set
one_result <- fit_c4_aci_hyperbola(
  licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE],
  calculate_confidence_intervals = FALSE
)

# Calculate confidence limits for the fit parameters
parameters_with_limits <- confidence_intervals_c4_aci_hyperbola(
    licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE],
    one_result$parameters
)

# View confidence limits and best estimate for Vmax
parameters_with_limits[, c('Vmax_lower', 'Vmax', 'Vmax_upper')]

Consolidate a list of lists

Description

Consolidates a list of lists into a regular list by combining like-named elements.

Usage

  consolidate(x)

  ## S3 method for class 'data.frame'
consolidate(x)

  ## S3 method for class 'exdf'
consolidate(x)

Arguments

Details

consolidate is generic, with methods defined for nested lists of data frames and exdf objects.

Value

A list with elements named name_1, name_2, ..., name_M, where each element was created by combining all elements of x with the same name using rbind; for example, the element with name name_1 will be created by calling rbind(list_1$name_1, list_2$name_1, ..., list_N$name_1). Before calling rbind, each element will be limited to the columns that are common to all elements with the same name.

See Also

exdf

Examples

# Example 1: Create a nested list of data frames and then consolidate them into
# a regular list by combining the like-named elements
nested_df_list <- list(
  list_1 = list(
    name_1 = data.frame(A = c(1, 2), B = c(0, 0)),
    name_2 = data.frame(A = c(3, 4), B = c(0, 0)),
    name_3 = data.frame(A = c(5, 6), B = c(0, 0))
  ),
  list_2 = list(
    name_1 = data.frame(A = c(7, 8), B = c(0, 0)),
    name_2 = data.frame(A = c(9, 10), B = c(0, 0)),
    name_3 = data.frame(A = c(11, 12), B = c(0, 0))
  ),
  list_3 = list(
    name_1 = data.frame(A = c(13, 14), B = c(0, 0)),
    name_2 = data.frame(A = c(15, 16), B = c(0, 0)),
    name_3 = data.frame(A = c(17, 18), B = c(0, 0))
  )
)

str(nested_df_list)

consolidated_df_list <- consolidate(nested_df_list)

str(consolidated_df_list)

# Example 2: Create a nested list of `exdf` objects and then consolidate them
# into a regular list by combining the like-named elements. Here, some of the
# elements have columns not present in the others (for example,
# `nested_exdf_list$list_3$name_1`). However, these "extra" columns are removed
# before calling `rbind` and they do not appear in `consolidated_exdf_list`.
nested_exdf_list <- list(
  list_1 = list(
    name_1 = exdf(data.frame(A = c(1, 2), B = c(0, 0))),
    name_2 = exdf(data.frame(A = c(3, 4), B = c(0, 0))),
    name_3 = exdf(data.frame(A = c(5, 6), B = c(0, 0)))
  ),
  list_2 = list(
    name_1 = exdf(data.frame(A = c(7, 8), B = c(0, 0))),
    name_2 = exdf(data.frame(A = c(9, 10), B = c(0, 0))),
    name_3 = exdf(data.frame(A = c(11, 12), B = c(0, 0)))
  ),
  list_3 = list(
    name_1 = exdf(data.frame(A = c(13, 14), B = c(0, 0), C = c(-1, -2))),
    name_2 = exdf(data.frame(A = c(15, 16), B = c(0, 0), C = c(-1, -2))),
    name_3 = exdf(data.frame(A = c(17, 18), B = c(0, 0), C = c(-1, -2)))
  )
)

str(nested_exdf_list)

consolidated_exdf_list <- consolidate(nested_exdf_list)

str(consolidated_exdf_list)

Read and write CSV files representing an exdf object

Description

Functions for reading and writing CSV files that represent an exdf object.

Usage

  read.csv.exdf(file, ...)

  write.csv.exdf(x, file, ...)

Arguments

Details

An exdf object can be written to a CSV file by directly calling write.csv, but this approach causes some column names to be unintentionally modified. For example, any spaces will be replaced by periods. This can potentially cause problems when reloading the data from the CSV file.

Instead, it is preferred to use write.csv.exdf, which will not modify any column names. When writing the CSV file, it is saved with the column names in the first row, the categories in the second row, the units in the third row, and the data in the remaining rows.

The resulting file can be read using read.csv.exdf. Here, the names, categories, and units are read from the first three rows of the specified file, and the data values from the remaining rows. An exdf object is then created from this information.

Value

The write.csv.exdf function does not return anything. The read.csv.exdf function returns an exdf object representing the contents of file.

Examples

# Read a CSV file included with the PhotoGEA package; this file was created
# using `write.csv.exdf`.
licor_file <- read.csv.exdf(
  PhotoGEA_example_file_path('ball_berry_1.csv')
)

# Now rewrite this to a temporary CSV file
tf <- tempfile(fileext = ".csv")
tf

write.csv.exdf(licor_file, tf)

Deprecated functions

Description

Deprecated functions that will be fully removed in future releases. Each of these functions will produce an error when called that will redirect the user to a suitable replacement.

Usage

  read_tdl_file(...)

  read_licor_file(...)

  check_licor_data(...)

  calculate_arrhenius(...)

  calculate_peaked_gaussian(...)

Arguments

Value

None of the deprecated functions return anything.

Examples

# These functions all throw errors, so we will wrap them in `tryCatch` here

tryCatch(
  read_tdl_file(),
  error = function(e) {print(e)}
)

tryCatch(
  read_licor_file(),
  error = function(e) {print(e)}
)

tryCatch(
  check_licor_data(),
  error = function(e) {print(e)}
)

tryCatch(
  calculate_arrhenius(),
  error = function(e) {print(e)}
)

tryCatch(
  calculate_peaked_gaussian(),
  error = function(e) {print(e)}
)

Retrieve the dimension of an exdf object

Description

Returns the dimensions of an exdf object's main_data. Also enables nrow and ncol for exdf objects.

Usage

  ## S3 method for class 'exdf'
dim(x)

Arguments

Value

Returns dim(x[['main_data']]).

See Also

exdf

Examples

simple_exdf <- exdf(data.frame(A = 1), data.frame(A = 'u'), data.frame(A = 'c'))

dim(simple_exdf)
dim(simple_exdf[['main_data']]) # An equivalent command

nrow(simple_exdf)
ncol(simple_exdf)

Retrieve or set the dimension names of an exdf object

Description

Returns or sets the dimension names of an exdf object's main_data. When setting names, the column names of the exdf object's units and categories are also set. Also enables colnames and rownames for exdf objects.

Usage

  ## S3 method for class 'exdf'
dimnames(x)

  ## S3 replacement method for class 'exdf'
dimnames(x) <- value

Arguments

Value

Returns dimnames(x[['main_data']]).

See Also

exdf

Examples

simple_exdf <- exdf(data.frame(A = 1), data.frame(A = 'u'), data.frame(A = 'c'))

dimnames(simple_exdf)
dimnames(simple_exdf[['main_data']]) # An equivalent command

colnames(simple_exdf) <- "B"
rownames(simple_exdf) <- 2

colnames(simple_exdf)
rownames(simple_exdf)

Document exdf columns by specifying units and categories

Description

Adds new columns to a table-like object, and sets/modifies the units or categories of columns in an exdf object.

Usage

  document_variables(x, ...)

  ## S3 method for class 'data.frame'
document_variables(x, ...)

  ## S3 method for class 'exdf'
document_variables(x, ...)

Arguments

Value

An object based on x with new and/or modified columns.

See Also

exdf

Examples

# Create a simple exdf object with two columns (`A` and `B`) and default values
# for its units and categories.
simple_exdf <- exdf(data.frame(A = c(3, 2, 7, 9), B = c(4, 5, 1, 8)))

print(simple_exdf)

# Specify units and categories for the `A` and `B` columns, and add a new `C`
# column.
document_variables(
  simple_exdf,
  c('cat1', 'A', 'm'), # The category of `A` is `cat1` and its units are `m`
  c('cat2', 'B', 's'), # The category of `B` is `cat2` and its units are `s`
  c('cat3', 'C', 'g')  # The category of `C` is `cat3` and its units are `g`
)

# Do the same but for a data frame; in this case columns A and B will not be
# altered, but a new column C will be added (and initialized to NA)
document_variables(
  simple_exdf$main_data,
  c('cat1', 'A', 'm'), # The category of `A` is `cat1` and its units are `m`
  c('cat2', 'B', 's'), # The category of `B` is `cat2` and its units are `s`
  c('cat3', 'C', 'g')  # The category of `C` is `cat3` and its units are `g`
)

Generate an error function for C3 A-Ci curve fitting

Description

Creates a function that returns an error value (the negative of the natural logarithm of the likelihood) representing the amount of agreement between modeled and measured An values. When this function is minimized, the likelihood is maximized.

Internally, this function uses apply_gm to calculate Cc, and then uses link{calculate_c3_assimilation} to calculate assimilation rate values that are compared to the measured ones.

Usage

  error_function_c3_aci(
    replicate_exdf,
    fit_options = list(),
    sd_A = 1,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    gmc_norm_column_name = 'gmc_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    cj_crossover_min = NA,
    cj_crossover_max = NA,
    hard_constraints = 0,
    debug_mode = FALSE,
    ...
  )

Arguments

Details

When fitting A-Ci curves using a maximum likelihood approach, it is necessary to define a function that calculates the likelihood of a given set of alpha_g, alpha_old, alpha_s, alpha_t, Gamma_star_at_25, gmc_at_25, J_at_25, Kc_at_25, Ko_at_25, RL_at_25, Tp_at_25, and Vcmax_at_25 values by comparing a model prediction to a measured curve. This function will be passed to an optimization algorithm which will determine the values that produce the largest likelihood.

The error_function_c3_aci returns such a function, which is based on a particular A-Ci curve and a set of fitting options. It is possible to just fit a subset of the available fitting parameters; by default, the fitting parameters are alpha_old, J_at_25, RL_at_25, Tp_at_25, and Vcmax_at_25. This behavior can be changed via the fit_options argument.

For practical reasons, the function actually returns values of -ln(L), where L is the likelihood. The logarithm of L is simpler to calculate than L itself, and the minus sign converts the problem from a maximization to a minimization, which is important because most optimizers are designed to minimize a value.

Sometimes an optimizer will choose biologically unreasonable parameter values that nevertheless produce good fits to the supplied assimilation values. A common problem is that the fit result may not indicate Ac-limited assimilation at low CO2 values, which should be the case for any A-Ci curves measured at saturating light. In this case, the optional cj_crossover_min and cj_crossover_max can be used to constrain the range of Cc values (in ppm) where Aj is allowed to be the overall rate limiting factor. If the crossover from Rubisco-limited to RuBP-regeneration limited assimilation occurs outside these bounds (when they are supplied), a heavy penalty will be added to the error function, preventing the optimizer from choosing those parameter values.

A penalty is also added for any parameter combination where An is not a number, or where calculate_c3_assimilation produces an error.

Value

A function with one input argument guess, which should be a numeric vector representing values of the parameters to be fitted (which are specified by the fit_options input argument.) Each element of guess is the value of one parameter (arranged in alphabetical order.) For example, with the default settings, guess should contain values of alpha_old, J_at_25, RL_at_25, Tp_at_25, and Vcmax_at_25 (in that order).

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Define an error function for one curve from the set
error_fcn <- error_function_c3_aci(
  licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE]
)

# Evaluate the error for:
#  alpha_old = 0
#  J_at_25 = 236
#  RL_at_25 = 4e-8
#  Tp_at_25 = 22.7
#  Vcmax_at_25 = 147
error_fcn(c(0, 236, 4e-8, 22.7, 147))

# Make a plot of likelihood vs. Vcmax when other parameters are fixed to the
# values above.
vcmax_error_fcn <- function(Vcmax) {error_fcn(c(0, 236, 4e-8, 22.7, Vcmax))}
vcmax_seq <- seq(135, 152, length.out = 41)

lattice::xyplot(
  exp(-sapply(vcmax_seq, vcmax_error_fcn)) ~ vcmax_seq,
  type = 'b',
  xlab = 'Vcmax_at_25 (micromol / m^2 / s)',
  ylab = 'Negative log likelihood (dimensionless)'
)

Generate an error function for C3 Variable J curve fitting

Description

Creates a function that returns an error value (the negative of the natural logarithm of the likelihood) representing the amount of agreement between modeled and measured An values. When this function is minimized, the likelihood is maximized.

Internally, this function uses link{calculate_c3_variable_j} and link{calculate_c3_assimilation} to calculate assimilation rate values that are compared to the measured ones.

Usage

  error_function_c3_variable_j(
    replicate_exdf,
    fit_options = list(),
    sd_A = 1,
    Wj_coef_C = 4.0,
    Wj_coef_Gamma_star = 8.0,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    gamma_star_norm_column_name = 'Gamma_star_norm',
    j_norm_column_name = 'J_norm',
    kc_norm_column_name = 'Kc_norm',
    ko_norm_column_name = 'Ko_norm',
    oxygen_column_name = 'oxygen',
    phips2_column_name = 'PhiPS2',
    qin_column_name = 'Qin',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    tp_norm_column_name = 'Tp_norm',
    vcmax_norm_column_name = 'Vcmax_norm',
    cj_crossover_min = NA,
    cj_crossover_max = NA,
    hard_constraints = 0,
    require_positive_gmc = 'positive_a',
    gmc_max = Inf,
    check_j = TRUE,
    debug_mode = FALSE,
    ...
  )

Arguments

Details

When fitting A-Ci + chlorophyll fluorescence curves using the Variable J method, it is necessary to define a function that calculates the likelihood of a given set of alpha_g, alpha_old, alpha_s, alpha_t, Gamma_star, J_at_25, RL_at_25, tau, Tp_at_25, and Vcmax_at_25 values by comparing a model prediction to a measured curve. This function will be passed to an optimization algorithm which will determine the values that produce the smallest error.

The error_function_c3_variable_j returns such a function, which is based on a particular replicate and a set of fitting options. It is possible to just fit a subset of the available fitting parameters; by default, the fitting parameters are alpha_old, J_at_25, RL_at_25, Tp_at_25, tau, and Vcmax_at_25. This behavior can be changed via the fit_options argument.

For practical reasons, the function actually returns values of -ln(L), where L is the likelihood. The logarithm of L is simpler to calculate than L itself, and the minus sign converts the problem from a maximization to a minimization, which is important because most optimizers are designed to minimize a value.

Sometimes an optimizer will choose biologically unreasonable parameter values that nevertheless produce good fits to the supplied assimilation values. There are several options for preventing an optimizer from choosing such parameter values:

• Enforcing Rubisco limitations: A common problem is that the fit result may not indicate Rubisc-limited assimilation at low CO2 values, which should be the case for any A-Ci curves measured at saturating light. In this case, the optional cj_crossover_min and cj_crossover_max can be used to constrain the range of Cc values (in ppm) where Wj is allowed to be the overall rate limiting factor. If the crossover from Rubisco-limited to RuBP-regeneration limited carboxylation occurs outside these bounds (when they are supplied), a heavy penalty will be added to the error function, preventing the optimizer from choosing those parameter values.

• Requiring positive gmc: The Variable J method sometimes predicts negative values of the mesophyll conductance (gmc). Such values are probably not realistic, especially when Cc is above the CO2 compensation point. The require_positive_gmc input argument can be used to penalize negative values of gmc. When require_positive_gmc is 'all', a heavy penalty will be added to the error function if there are any negative gmc values. When require_positive_gmc is 'positive_a', a heavy penalty will be added to the error function if there are any negative gmc values when A is positive; negative gmc for negative A will be allowed. When require_positive_gmc is 'none', these restrictions are disabled and no penalties are added for negative gmc.

• Preventing large values of gmc: The Variable J method sometimes produces unreasonably high values of gmc. The gmc_max argument can be used to address this. If any predicted gmc values exceed gmc_max when A is positive, a heavy penalty will be added to the error function.

• Enforcing consistent RuBP regeneration rates: In principle, the actual RuBP regeneration rate (J_F) should always be less than or equal to its maximum value for a given Qin and leaf temperature (J_tl), with equality only occuring when assimilation is RuBP-regeneration-limited. When check_j is TRUE, a heavy penalty will be added to the error function for any parameter values where J_F is greater than J_tl at any point in the curve.

A penalty is also added for any parameter combination where An is not a number, or where calculate_c3_variable_j or calculate_c3_assimilation produce an error.

Value

A function with one input argument guess, which should be a numeric vector representing values of the parameters to be fitted (which are specified by the fit_options input argument.) Each element of guess is the value of one parameter (arranged in alphabetical order.) For example, with the default settings, guess should contain values of alpha_old, J_at_25, RL_at_25, tau, Tp_at_25, and Vcmax_at_25 (in that order).

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

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

# Calculate temperature-dependent values of C3 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_bernacchi)

# Define an error function for one curve from the set
error_fcn <- error_function_c3_variable_j(
  licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE]
)

# Evaluate the error for:
#  alpha_old = 1.9
#  J_at_25 = 270
#  RL_at_25 = 1.9
#  tau = 0.42
#  Tp_at_25 = 8.7
#  Vcmax_at_25 = 258
error_fcn(c(1.9, 270, 1.9, 0.42, 8.7, 258))

# Make a plot of error vs. Tp_at_25 when the other parameters are fixed to the
# values above. As Tp_at_25 increases, it eventually stops limiting the
# assimilation rate and its value stops influencing the error.
tpu_error_fcn <- function(Tp_at_25) {error_fcn(c(1.9, 270, 1.9, 0.42, Tp_at_25, 258))}
tpu_seq <- seq(5, 12, by = 0.25)

lattice::xyplot(
  sapply(tpu_seq, tpu_error_fcn) ~ tpu_seq,
  type = 'b',
  xlab = 'Tp at 25 degrees C (micromol / m^2 / s)',
  ylab = 'Negative log likelihood (dimensionless)'
)

Generate an error function for C4 A-Ci curve fitting

Description

Creates a function that returns an error value (the negative of the natural logarithm of the likelihood) representing the amount of agreement between modeled and measured An values. When this function is minimized, the likelihood is maximized.

Internally, this function uses apply_gm to calculate Cc, and then uses link{calculate_c4_assimilation} to calculate assimilation rate values that are compared to the measured ones.

Usage

  error_function_c4_aci(
    replicate_exdf,
    fit_options = list(),
    sd_A = 1,
    x_etr = 0.4,
    a_column_name = 'A',
    ao_column_name = 'ao',
    ci_column_name = 'Ci',
    gamma_star_column_name = 'gamma_star',
    gmc_norm_column_name = 'gmc_norm',
    j_norm_column_name = 'J_norm',
    kc_column_name = 'Kc',
    ko_column_name = 'Ko',
    kp_column_name = 'Kp',
    oxygen_column_name = 'oxygen',
    rl_norm_column_name = 'RL_norm',
    total_pressure_column_name = 'total_pressure',
    vcmax_norm_column_name = 'Vcmax_norm',
    vpmax_norm_column_name = 'Vpmax_norm',
    hard_constraints = 0,
    debug_mode = FALSE
  )

Arguments

Details

When fitting A-Ci curves, it is necessary to define a function that calculates the likelihood of a given set of alpha_psii, gbs, gmc_at_25, J_at_25, RL_at_25, Rm_frac, Vcmax_at_25, Vpmax_at_25, and Vpr values by comparing a model prediction to a measured curve. This function will be passed to an optimization algorithm which will determine the values that produce the smallest error.

The error_function_c4_aci returns such a function, which is based on a particular A-Ci curve and a set of fitting options. It is possible to just fit a subset of the available fitting parameters; by default, the fitting parameters are RL_at_25, Vcmax_at_25, and Vpmax_at_25. This behavior can be changed via the fit_options argument.

For practical reasons, the function actually returns values of -ln(L), where L is the likelihood. The logarithm of L is simpler to calculate than L itself, and the minus sign converts the problem from a maximization to a minimization, which is important because most optimizers are designed to minimize a value.

A penalty is added to the error value for any parameter combination where An is not a number, or where calculate_c4_assimilation produces an error.

Value

A function with one input argument guess, which should be a numeric vector representing values of the parameters to be fitted (which are specified by the fit_options input argument.) Each element of guess is the value of one parameter (arranged in alphabetical order.) For example, with the default settings, guess should contain values of RL_at_25, Vcmax_at_25, and Vpmax_at_25 (in that order).

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c4_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Calculate temperature-dependent values of C4 photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c4_temperature_param_vc)

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

# Define an error function for one curve from the set
error_fcn <- error_function_c4_aci(
  licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE]
)

# Evaluate the error for RL_at_25 = 0, Vcmax_at_25 = 35, Vpmax_at_25 = 180
error_fcn(c(0, 35, 180))

# Make a plot of error vs. Vcmax_at_25 when the other parameters are fixed to
# the values above.
vcmax_error_fcn <- function(Vcmax_at_25) {error_fcn(c(0, Vcmax_at_25, 180))}
vcmax_seq <- seq(20, 50)

lattice::xyplot(
  sapply(vcmax_seq, vcmax_error_fcn) ~ vcmax_seq,
  type = 'b',
  xlab = 'Vcmax at 25 degrees C (micromol / m^2 / s)',
  ylab = 'Negative log likelihood (dimensionless)'
)

Generate an error function for C4 A-Ci curve fitting with a hyperbola

Description

Creates a function that returns an error value (the negative of the natural logarithm of the likelihood) representing the amount of agreement between modeled and measured An values. When this function is minimized, the likelihood is maximized.

Internally, this function uses link{calculate_c4_assimilation_hyperbola} to calculate assimilation rate values that are compared to the measured ones.

Usage

  error_function_c4_aci_hyperbola(
    replicate_exdf,
    fit_options = list(),
    sd_A = 1,
    a_column_name = 'A',
    ci_column_name = 'Ci',
    hard_constraints = 0,
    debug_mode = FALSE
  )

Arguments

Details

When fitting A-Ci curves, it is necessary to define a function that calculates the likelihood of a given set of c4_curvature, c4_slope, rL, and Vmax values by comparing a model prediction to a measured curve. This function will be passed to an optimization algorithm which will determine the values that produce the smallest error.

The error_function_c4_aci_hyperbola returns such a function, which is based on a particular A-Ci curve and a set of fitting options. It is possible to just fit a subset of the available fitting parameters; by default, all are fit. This behavior can be changed via the fit_options argument.

For practical reasons, the function actually returns values of -ln(L), where L is the likelihood. The logarithm of L is simpler to calculate than L itself, and the minus sign converts the problem from a maximization to a minimization, which is important because most optimizers are designed to minimize a value.

A penalty is added to the error value for any parameter combination where An is not a number, or where calculate_c4_assimilation_hyperbola produces an error.

Value

A function with one input argument guess, which should be a numeric vector representing values of the parameters to be fitted (which are specified by the fit_options input argument.) Each element of guess is the value of one parameter (arranged in alphabetical order.) For example, with the default settings, guess should contain values of c4_curvature, c4_slope, rL, and Vmax (in that order).

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c4_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Define an error function for one curve from the set
error_fcn <- error_function_c4_aci_hyperbola(
  licor_file[licor_file[, 'species_plot'] == 'maize - 5', , TRUE]
)

# Evaluate the error for c4_curvature = 0.8, c4_slope = 0.5, rL = 1.0, Vmax = 65
error_fcn(c(0.8, 0.5, 1.0, 65))

# Make a plot of error vs. Vmax when the other parameters are fixed to
# the values above.
vmax_error_fcn <- function(Vmax) {error_fcn(c(0.8, 0.5, 1.0, Vmax))}
vmax_seq <- seq(55, 75)

lattice::xyplot(
  sapply(vmax_seq, vmax_error_fcn) ~ vmax_seq,
  type = 'b',
  xlab = 'Vmax (micromol / m^2 / s)',
  ylab = 'Negative log likelihood (dimensionless)'
)

Estimate variance of measured Licor values

Description

Estimates variance and standard deviation of the net CO2 assimilation rate as measured by a Licor Li-6800 or similar portable photosynthesis system.

Usage

  estimate_licor_variance(
    exdf_obj,
    sd_CO2_r,
    sd_CO2_s,
    sd_flow,
    sd_H2O_r,
    sd_H2O_s,
    a_column_name = 'A',
    co2_r_column_name = 'CO2_r',
    co2_s_column_name = 'CO2_s',
    corrfact_column_name = 'CorrFact',
    flow_column_name = 'Flow',
    h2o_r_column_name = 'H2O_r',
    h2o_s_column_name = 'H2O_s',
    s_column_name = 'S'
  )

Arguments

Details

Uses the error propogation formula to calculate the influence of the variance in CO2_r, CO2_s, etc on the variance of A, as calculated by a Licor LI-6800.

Value

An exdf object based on exdf_obj that includes additional columns representing the standard deviation of A measurements (sd_A), and the individual terms comprising the total variance of A, such as var_CO2_r, var_CO2_s, etc.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Estimate variance in measured A values
licor_file <- estimate_licor_variance(
  licor_file,
  sd_CO2_r = 1,
  sd_CO2_s = 0.1,
  sd_flow = 0.2,
  sd_H2O_r = 0.5,
  sd_H2O_s = 0.1
)

# Plot each component of the total variance of A
lattice::xyplot(
    var_CO2_r + var_CO2_s + var_flow + var_H2O_r + var_H2O_s + var_A ~ Ci | species_plot,
    data = licor_file$main_data,
    type = 'b',
    pch = 16,
    auto = TRUE
)

# Plot the standard deviation of A
lattice::xyplot(
    sd_A ~ Ci,
    group = species_plot,
    data = licor_file$main_data,
    type = 'b',
    pch = 16,
    auto = TRUE
)

Estimate the operating point from an A-Ci curve

Description

Uses linear interpolation to estimate Cc, Ci, and An at atmospheric CO2 concentration from the data in the exdf object, which should represent a single A-Ci curve. This function can accomodate alternative column names for the variables taken from the data file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  estimate_operating_point(
    aci_exdf,
    Ca_atmospheric,
    type = 'c3',
    a_column_name = 'A',
    ca_column_name = 'Ca',
    cc_column_name = 'Cc',
    ci_column_name = 'Ci',
    pcm_column_name = 'PCm',
    return_list = FALSE
  )

Arguments

Details

When analyzing or interpreting A-Ci curves, it is often useful to determine the values of Ci and An that correspond to typical growth conditions (where Ca is set to the atmospheric value). Together, these special values of Ci and An specify the "operating point" of the leaf.

However, for a variety of practical reasons, most A-Ci curves do not actually contain a measurement point where Ca is at the atmospheric value. Nevertheless, it is possible to apply linear interpolation to the observed Ci - Ca and An - Ca relations to estimate the operating point. This function automates that procedure. It also calculates the operating values of Cc (for c3 A-Ci curves) and PCm (for c4 A-Ci curves).

This function assumes that aci_exdf represents a single A-Ci curve. Typically, this function is not directly called by users because the fitting functions fit_c3_aci and fit_c4_aci automatically use this function to determine the operating point.

Value

The return value depends on return_list and type.

When return_list is FALSE, this function returns an exdf object based on aci_exdf that includes its identifier columns as well as values of Ca_atmospheric, operating_Ci, operating_An, and operating_Cc (or operating_PCm) in columns with those names.

When return_list is TRUE, this function returns a list with the following named elements: Ca_atmospheric, operating_Ci, operating_An, operating_Cc (or operating_PCm), and operating_exdf. The first four are numeric values as described above, while operating_exdf is an exdf object with one row that can be passed to calculate_c3_assimilation or calculate_c4_assimilation in order to estimate the operating An from a photosynthesis model.

If Ca_atmospheric is outside the range of Ca values in aci_exdf, or if all provided values of Ca are NA, then the operating point cannot be reasonably estimated; in this case, an explanation is returned as the operating_point_msg column or list element, and all other calculated return values are set to NA. Otherwise, the operating_point_msg is an empty string.

If Ca_atmospheric is NA, all calculated return values are set to NA without any additional explanation.

Examples

# Read an example Licor file included in the PhotoGEA package
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('c3_aci_1.xlsx')
)

# Define a new column that uniquely identifies each curve
licor_file[, 'species_plot'] <-
  paste(licor_file[, 'species'], '-', licor_file[, 'plot'] )

# Organize the data
licor_file <- organize_response_curve_data(
    licor_file,
    'species_plot',
    c(9, 10, 16),
    'CO2_r_sp'
)

# Calculate temperature-dependent values of photosynthetic parameters
licor_file <- calculate_temperature_response(licor_file, c3_temperature_param_sharkey)

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

# Calculate Cc, assuming an infinite mesophyll conductance (so `Cc` = `Ci`)
licor_file <- apply_gm(licor_file, Inf)

# Determine the operating point for just one curve from the data set
one_result <- estimate_operating_point(
  licor_file[licor_file[, 'species_plot'] == 'tobacco - 1', , TRUE],
  Ca_atmospheric = 420
)

one_result[, 'operating_Cc']
one_result[, 'operating_Ci']
one_result[, 'operating_An']
one_result[, 'operating_point_msg']

Example data files

Description

The PhotoGEA package includes several data files that can be used to demonstrate different functions and analysis techniques.

Details

The following files are included with the package:

• ball_berry_1.xlsx and ball_berry_2.xlsx: Two log files created by Licor Li-6800 portable gas exchange measurement systems. These log files each contain several Ball-Berry curves. Several user constants were defined in these logs that can be used to identify individual curves or subsets of curves: species, plot, and instrument. These files are used in the "Analyzing Ball-Berry Data" vignette and in other examples.

• ball_berry_1.csv: A CSV version of ball_berry_1.xlsx, which was created by reading the Excel file with read_gasex_file and then saving it using write.csv.exdf. This can be done as follows: tmp <- read_gasex_file(PhotoGEA_example_file_path('ball_berry_1.xlsx')); write.csv.exdf(tmp, 'ball_berry_1.csv')

• c3_aci_1.xlsx and c3_aci_2.xlsx: Two log files created by Licor Li-6800 portable gas exchange measurement systems. These log files each contain several C3 CO2 response (or A-Ci) curves. Several user constants were defined in these logs that can be used to identify individual curves or subsets of curves: species, plot, and instrument. These files are used in the "Analyzing C3 A-Ci Curves" vignette and in other examples. The Remarks sheet of c3_aci_2.xlsx was deleted from the original version as a test for read_licor_6800_Excel.

• c4_aci_1.xlsx and c4_aci_2.xlsx: Two log files created by Licor Li-6800 portable gas exchange measurement systems. These log files each contain several C4 CO2 response (or A-Ci) curves. Several user constants were defined in these logs that can be used to identify individual curves or subsets of curves: species, plot, and instrument. These files are used in the "Analyzing C4 A-Ci Curves" vignette and in other examples.

• tdl_sampling_1.dat and tdl_sampling_2.dat: Two log files created by a Campbell Scientific CR3000 data logger, representing data from a tunable diode laser (TDL) system. These files are used in the "Analyzing TDL Data" vignette and in other examples.

• plaintext_licor_file: A log file created by a Licor Li-6800 portable gas exchange measurement system. This file contains several CO2 response (or A-Ci) curves. Several user constants were defined in this log that can be used to identify individual curves or subsets of curves: species, plot, and instrument.

• plaintext_licor_file_v2: A log file based on plaintext_licor_file that has two separate [Data] and [Header] sections, as if the log file had been closed and reopened halfway through the measurement sequence. It also has an extra blank line at the end.

• licor_for_gm_site11.xlsx, licor_for_gm_site13.xslsx, and tdl_for_gm: Two Licor Li-6800 log files and a CR3000 TDL log file, respectively. These files are used as an example of loading and processing combined gas exchange and isotope discrimination measurements. Each Licor log file includes 6 points measured with the CO2_r setpoint set to 715 ppm and 6 points with the setpoint set to 450 ppm.

Since none of these data files have been published, noise has been added to the original data. Thus, they are similar to real measurements, but no useful conclusions can be drawn from them.

After installing 'PhotoGEA', copies of these files will be stored in the R package directory (in the PhotoGEA/extdata subdirectory). This location will be unique to your computer, but full paths to these files can be obtained using the PhotoGEA_example_file_path function.

Examples

# Print full paths to the example files
PhotoGEA_example_file_path('ball_berry_1.xlsx')
PhotoGEA_example_file_path('ball_berry_2.xlsx')
PhotoGEA_example_file_path('c3_aci_1.xlsx')
PhotoGEA_example_file_path('c3_aci_2.xlsx')
PhotoGEA_example_file_path('c4_aci_1.xlsx')
PhotoGEA_example_file_path('c4_aci_2.xlsx')
PhotoGEA_example_file_path('licor_for_gm_site11.xlsx')
PhotoGEA_example_file_path('licor_for_gm_site13.xlsx')
PhotoGEA_example_file_path('plaintext_licor_file')
PhotoGEA_example_file_path('plaintext_licor_file_v2')
PhotoGEA_example_file_path('tdl_for_gm.dat')
PhotoGEA_example_file_path('tdl_sampling_1.dat')
PhotoGEA_example_file_path('tdl_sampling_2.dat')

Exclude outliers from a data set

Description

Excludes outliers from a data set using the "1.5 interquartile range" rule.

Usage

  exclude_outliers(x, col_for_analysis, INDICES, method = 'exclude')

  ## S3 method for class 'data.frame'
exclude_outliers(x, col_for_analysis, INDICES, method = 'exclude')

  ## S3 method for class 'exdf'
exclude_outliers(x, col_for_analysis, INDICES, method = 'exclude')

Arguments

Details

exclude_outliers is generic, with methods defined for data frames and exdf objects. This function uses a simple rule to detect outliers, where any point that deviates from the mean by more than 1.5 * IQR, where IQR is the interquartile range, is said to be an outlier. This method is also sometimes referred to as "Tukey's Fences," as seen in the Wikipedia page about outliers.

For data sets with extreme outliers, it may be necessary to exclude outliers more than once to actually remove them all.

Value

This function returns an object formed from x, where the results depend on on the value of method.

When method is 'remove', the returned object is a modified copy of x where all rows in which the value of col_for_analysis is an outlier have been removed.

When method is 'exclude', the returned object is a modified copy of x where all outlier values of col_for_analysis have been replaced with NA.

See Also

exdf

Examples

# Read a Licor file included with the PhotoGEA package; this file includes
# several light response curves that can be identified by the 'species' and
# 'plot' columns.
licor_file <- read_gasex_file(
  PhotoGEA_example_file_path('ball_berry_1.xlsx')
)

# Remove points from each response curve in the data where the leaf temperature
# is determined to be an outlier
licor_file_clean <- exclude_outliers(
  licor_file,
  'TleafCnd',
  list(licor_file[, 'species'], licor_file[, 'plot']),
  method = 'remove'
)

# Check to see how many points remain after removing outliers
str(list('original' = nrow(licor_file), 'clean' = nrow(licor_file_clean)))

Extended data frame

Description

An "extended data frame" (exdf) is an object similar to a data frame, but which also contains information about the units and categories of each column.

Usage

  exdf(
    main_data = data.frame(),
    units = NULL,
    categories = NULL,
    ...
  )

Arguments

Details

The exdf class was originally created as a way to represent the contents of a Licor Excel file in an R structure. In Licor Excel files, each column has a name, units, and a category; for example, the column for values of net assimilation rate is called A, has units of micromol / m^2 / s, and is categorized as a GasEx variable.

From a technical point of view, an exdf object is simply a list with three required elements: main_data, units, and categories. Each of these should be a data frame with the same column names, as described above. It is also possible for an exdf object to have additional entries such as a filename that stores the name of the file that was used to create the exdf.

Several S3 methods have been defined for exdf objects, following the general guidance from Advanced R on S3 classes:

• is.exdf

• as.data.frame.exdf

• print.exdf

• str.exdf

• length.exdf

• dim.exdf

• dimnames.exdf

• [.exdf

• [<-.exdf

• rbind.exdf

• cbind.exdf

• split.exdf

• by.exdf

Note that the column names of main_data, units, and categories must be unique; the make.unique function can be useful for ensuring this.

Value

An exdf object as described above.

Examples

# Example 1: Creating a simple exdf object with two columns (`A` and `B`) and
# default values for its units and categories. There are four values of each
# variable.
exdf(data.frame(A = c(3, 2, 7, 9), B = c(4, 5, 1, 8)))

# Example 2: Creating a simple exdf object with two columns (`A` and `B`) that
# have units of `m` and `s`, respectively, and categories of `Cat1` and `Cat2`,
# respectively. There are four values of each variable.
exdf(
  data.frame(A = c(3, 2, 7, 9), B = c(4, 5, 1, 8)),
  data.frame(A = 'm', B = 's'),
  data.frame(A = 'Cat1', B = 'Cat2')
)

Access or modify exdf elements

Description

Returns or sets the values of elements in an exdf object.

Usage

  ## S3 method for class 'exdf'
x[i, j, return_exdf = FALSE]

  ## S3 replacement method for class 'exdf'
x[i, j] <- value

Arguments

Details

Since an exdf object is actually a list of named elements, those elements can be accessed using the [[ or $ operators, and a list of all named elements can be obtained by calling names.

Elements of the main_data data frame of an exdf object can be accessed and set using the [ and [<- operators. When applied to an exdf object, these operators are essentially shortcuts to calling the same operators on the object's main_data data frame.

To create a new exdf object with a subset of the data contained in another exdf object, the [ operator with return_exdf = TRUE can be used.

Value

When return_exdf is FALSE, the access operator will return either a vector or a data frame, depending on the dimension of j. When return_exdf is TRUE, the access operator will return an exdf object.

See Also

exdf

Examples

# Create a small exdf object that includes an extra element in addition to the
# required ones (`main_data`, `units`, and `categories`).
small_exdf <- exdf(
  data.frame(A = c(3, 2, 7, 9), B = c(4, 5, 1, 8)),
  data.frame(A = 'm', B = 's'),
  data.frame(A = 'Cat1', B = 'Cat2'),
  extra_exdf_element = "This is an example of an extra exdf element"
)

# Accessing elements of `small_exdf`
names(small_exdf)     # Get the names of all elements of small_exdf
small_exdf[['units']] # View the units using the `[[` operator
small_exdf$categories # View the categories using the `$` operator

# Accessing elements of `small_exdf$main_data`
small_exdf[,1]   # Access the first column
small_exdf[1,]   # Access the first row
small_exdf[,'B'] # Access the column named 'B'
small_exdf[1,2]  # Access element 1 of column 2

# Equivalent (but longer) commands for accessing elements of `small_exdf$main_data`
small_exdf$main_data[,1]   # Access the first column
small_exdf$main_data[1,]   # Access the first row
small_exdf$main_data[,'B'] # Access the column named 'B'
small_exdf$main_data[1,2]  # Access element 1 of column 2

# Replacing elements of `small_exdf$main_data`
small_exdf[,'A'] <- seq_len(4)             # Replace column A with new values
small_exdf[small_exdf[,'A'] > 2, 'B'] <- 0 # Replace some rows of column B with new values

# Creating a new exdf object with a subset of the data from small_exdf. Here we
# specify `return_exdf = TRUE` so that the `[` operator returns an exdf object
# instead of a data frame
new_exdf <- small_exdf[small_exdf[,'A'] > 2, , TRUE]
names(new_exdf) # Check that the `extra_exdf_element` is still present
print(new_exdf) # Check that only the rows with A > 2 are included

Convert ID column to a factor with a suitable ordering

Description

Converts an ID column to a factor with a suitable ordering. In particular, this function will ensure that any IDs beginning with WT (or any other control group name, case-insensitive) will be ordered before other values. This is helpful when plotting results according to genotype.

Usage

  factorize_id_column(x, ...)

  ## S3 method for class 'character'
factorize_id_column(x, control_group_name = 'WT', ...)

  ## S3 method for class 'data.frame'
factorize_id_column(x, id_column_name, control_group_name = 'WT', ...)

  ## S3 method for class 'exdf'
factorize_id_column(x, id_column_name, control_group_name = 'WT', ...)

Arguments

Details

To choose an ordering, each unique identifier is split into three components: an initial control_group_name (if present), a final numeric value, and any other content in between these two. Then, the identifiers are sorted according to these three values, in order of control_group_name -> other content -> numeric value. Note that capitalization of any initial control_group_name values will be standardized to match the user-specified version.

This system works well with identifiers that represent genotypes/events, or that combine genotype/event with a replicate number.

Value

factorize_id_column.character returns the character vector as a factor with an appropriate ordering.

factorize_id_column.data.frame and factorize_id_column.exdf return a copy of the original table, where one column (specified by id_column_name) has been converted to a factor with an appropriate ordering.

See Also

exdf

Examples

# Identifiers that represent genotypes
genotype_ids <- c('4', 'control', '2', 'CONTROL', '8')

factorize_id_column(genotype_ids, control_group_name = 'control')

# Identifiers that represent `genotype - replicate` values
replicate_ids <- c('4 - 4', 'wT - 2', 'a - 2', 'WT - 1', '4 - 8', 'wt - 9')

factorize_id_column(replicate_ids)

# Data frame
dat <- data.frame(replicate_id = replicate_ids, val = seq_along(replicate_ids))

# Display data in bar chart - note the order of the replicates
lattice::barchart(val ~ replicate_id, data = dat)

# Display factorized data in bar chart - note the order of the replicates
lattice::barchart(val ~ replicate_id, data = factorize_id_column(dat, 'replicate_id'))

# Extended data frame
exdf_obj <- exdf(dat, units = data.frame(replicate_id = '', val = 'm / s'))

exdf_obj <- factorize_id_column(exdf_obj, 'replicate_id')

exdf_obj[, 'replicate_id']

Fits the Ball-Berry model to an experimental curve

Description

Calculates a linear fit of stomatal conductance vs. the Ball-Berry index using the data in the exdf object. This function can accomodate alternative column names for the variables taken from the Licor file in case they change at some point in the future. This function also checks the units of each required column and will produce an error if any units are incorrect.

Usage

  fit_ball_berry(
    replicate_exdf,
    bb_index_column_name = 'bb_index',
    gsw_column_name = 'gsw'
  )

Arguments

Details

The Ball-Berry model is a simple way to describe the response of a leaf's stomata to its assimilation rate and local environmental conditions. Specifically, it predicts stomatal conductance to water vapor using the following equation:

gsw = bb_0 + bb_1 * A * h_s / C_s

where gsw is the stomatal conductance, A is the net assimilation rate, h_s is the relative humidity at the leaf surface, and C_s is the CO2 concentration at the leaf surface. The term A * h_s / C_s is commonly referred to as the Ball-Berry index, while the intercept (bb_0) and slope (bb_1) of the linear relationship are the Ball-Berry parameters which describe the stomatal response.