Introduction to the photosynthesis package

The photosynthesis package simulates photosynthetic rate given a set of leaf traits and environmental conditions by solving the Farquhar-von Caemmerer-Berry C3 biochemical model. There are two main steps to using photosynthesis:

define leaf parameters, environmental parameters, temperature response parameters, and physical constants; and
solve for the chloroplastic CO2 concentration that balances CO2 supply and demand (photo and photosynthesis for single and multiple parameter sets, respectively).

In this vignette, I’ll show you how to:

run a minimum worked example using default parameters
replace default parameters
simulate photosynthetic rate along a gradient of CO $_2$ concentrations ( $A-C_\mathrm{c}$ curve)

Minimum worked example

You can use the default parameter settings and simulate photosynthetic rate in a single leaf using the make_*() functions and photo().


library(dplyr)
library(magrittr)
library(photosynthesis)

# Leaving the make_* functions empty will automatically default to defaults
# parameters.
bake_par   = make_bakepar()                       # temperature response parameters
constants  = make_constants(use_tealeaves = FALSE) # physical constants
leaf_par   = make_leafpar(use_tealeaves = FALSE)   # leaf parameters
enviro_par = make_enviropar(use_tealeaves = FALSE) # environmental parameters

photo(leaf_par, enviro_par, bake_par, constants, quiet = TRUE,
      use_tealeaves = FALSE)
#>                 C_chl         value convergence                  g_tc
#> 1 258.1218 [umol/mol] -1.363227e-06           0 0.1726157 [mol/m^2/s]
#>                       A            g_mc          g_mc25            g_sc
#> 1 27.94273 [umol/m^2/s] 0.4 [mol/m^2/s] 0.4 [mol/m^2/s] 0.4 [mol/m^2/s]
#>               g_uc         gamma_star       gamma_star25            J_max
#> 1 0.01 [mol/m^2/s] 37.9258 [umol/mol] 37.9258 [umol/mol] 200 [umol/m^2/s]
#>            J_max25              K_C            K_C25  k_mc                 K_O
#> 1 200 [umol/m^2/s] 268.3 [umol/mol] 268.3 [umol/mol] 1 [1] 165084.2 [umol/mol]
#>                 K_O25  k_sc  k_uc leafsize     phi_J            R_d
#> 1 165084.2 [umol/mol] 1 [1] 1 [1]  0.1 [m] 0.331 [1] 2 [umol/m^2/s]
#>            R_d25     T_leaf   theta_J           V_cmax         V_cmax25
#> 1 2 [umol/m^2/s] 298.15 [K] 0.825 [1] 150 [umol/m^2/s] 150 [umol/m^2/s]
#>              V_tpu          V_tpu25          C_air              O
#> 1 200 [umol/m^2/s] 200 [umol/m^2/s] 420 [umol/mol] 0.21 [mol/mol]
#>                P              PPFD      RH    wind           Ds_gmc
#> 1 101.3246 [kPa] 1500 [umol/m^2/s] 0.5 [1] 2 [m/s] 487.29 [J/K/mol]
#>            Ds_Jmax     Ea_gammastar           Ea_gmc          Ea_Jmax
#> 1 388.04 [J/K/mol] 24459.97 [J/mol] 68901.56 [J/mol] 56095.18 [J/mol]
#>              Ea_KC            Ea_KO            Ea_Rd         Ea_Vcmax
#> 1 80989.78 [J/mol] 23719.97 [J/mol] 40446.75 [J/mol] 52245.78 [J/mol]
#>            Ea_Vtpu           Ed_gmc          Ed_Jmax             D_c0
#> 1 52245.78 [J/mol] 148788.6 [J/mol] 121244.8 [J/mol] 1.29e-05 [m^2/s]
#>              D_h0             D_m0             D_w0   epsilon       eT
#> 1 1.9e-05 [m^2/s] 1.33e-05 [m^2/s] 2.12e-05 [m^2/s] 0.622 [1] 1.75 [1]
#>             G                 R                sigma                 C_i
#> 1 9.8 [m/s^2] 8.31446 [J/K/mol] 5.67e-08 [W/K^4/m^2] 350.1432 [umol/mol]

Replace default parameters

You can look at default parameters settings in the manual (run ?make_parameters). These defaults are reasonable, but of course you will probably want to use different choices and allow some parameters to vary. Here, I’ll demonstrate how to replace a default. In the next section, I’ll show you how to set up a gradient of parameter values over which to solve for leaf temperature.


# Use the `replace` argument to replace defaults. This must be a named list, and
# each named element must have the proper units specified. See `?make_parameters`
# for all parameter names and proper units.

# Temperature response parameters can be updated (but we won't do that here)
bake_par = make_bakepar()

# Physical constants probably do not need to be replaced in most cases,
# that's why we call them 'constants'!
constants  = make_constants(use_tealeaves = FALSE)

# First, we'll change photosynthetic photon flux density to 1000 umol / (m^2 s)
enviro_par = make_enviropar(
  replace = list(
    PPFD = set_units(1000, "umol/m^2/s")
    ), use_tealeaves = FALSE
  )

# Next, we'll change stomatal conductance to 0.3 mol / m^2 / s.
leaf_par = make_leafpar(
  replace = list(
    g_sc = set_units(0.3, mol / m^2 / s)
    ), use_tealeaves = FALSE
  )

photo <- photo(leaf_par, enviro_par, bake_par, constants, quiet = TRUE,
               use_tealeaves = FALSE)

photo |>
  select(PPFD, C_chl, A) |>
  knitr::kable()

PPFD	C_chl	A
1000 [umol/m^2/s]	252.6879 [umol/mol]	25.63373 [umol/m^2/s]

Environmental gradients

In the previous two examples, I used the photo function to solve for a single parameter set. In most cases, you’ll want to solve for many parameter sets. The function photosynthesis generalizes photo and makes it easy to solve for multiple parameter sets using the same argument structure. All you need to do is specify multiple values for one or more leaf or environmental parameters and photosynthesis uses the purrr::cross function to fit all combinations¹.


# As before, use the `replace` argument to replace defaults, but this time we
# enter multiple values

bake_par = make_bakepar()
constants = make_constants(use_tealeaves = FALSE)

# First, we'll change the PPFD to 1000 and 1500 umol / (m^2 s)
enviro_par = make_enviropar(
  replace = list(
    PPFD = set_units(c(1000, 1500), umol / m^2 / s)
    ), use_tealeaves = FALSE
  )

# Next, we'll change stomatal conductance to to 0.2 and 0.4 mol / m^2 / s
leaf_par  = make_leafpar(
  replace = list(
    g_sc = set_units(c(0.2, 0.4), mol / m^2 / s)
    ), use_tealeaves = FALSE
  )

# Now there should be 4 combinations (high and low g_sc crossed with high and low PPFD)
ph = photosynthesis(leaf_par, enviro_par, bake_par, constants, 
                     use_tealeaves = FALSE, progress = FALSE, quiet = TRUE)

ph |> 
  select(g_sc, PPFD, A) |>
  knitr::kable()

g_sc	PPFD	A
0.2 [mol/m^2/s]	1000 [umol/m^2/s]	24.34817 [umol/m^2/s]
0.4 [mol/m^2/s]	1000 [umol/m^2/s]	26.27649 [umol/m^2/s]
0.2 [mol/m^2/s]	1500 [umol/m^2/s]	25.68222 [umol/m^2/s]
0.4 [mol/m^2/s]	1500 [umol/m^2/s]	27.94273 [umol/m^2/s]

Parallel processing

It can take a little while to simulate many different parameter sets. If you have multiple processors available, you can speed things up by running simulations in parallel. In the photosynthesis function, simply use the parallel = TRUE argument to simulate in parallel. You’ll need to set up a future plan(). See ?future::plan for more detail. Here I’ll provide an example simulating an $A-C_c$ curve.


# NOTE: parallel example is not evaluated because it was causing an issue with CRAN, but you can copy-and-paste the code to run on your own machine.

library(future)
plan("multisession") # Set up plan

# We'll use the `replace` argument to enter multiple atmospheric CO2 concentrations

bake_par = make_bakepar()
constants  = make_constants(use_tealeaves = FALSE)

enviro_par = make_enviropar(
  replace = list(
    C_air = set_units(seq(10, 2000, length.out = 20), umol / mol)
    ), use_tealeaves = FALSE
  )

leaf_par  = make_leafpar(use_tealeaves = FALSE)

ph = photosynthesis(leaf_par, enviro_par, bake_par, constants, 
                     use_tealeaves = FALSE, progress = FALSE, 
                     quiet = TRUE, parallel = TRUE)

# Plot C_c versus A
library(ggplot2)

## Drop units for plotting
ph %<>% mutate_if(~ is(.x, "units"), drop_units)
ggplot(ph, aes(C_chl, A)) +
  geom_line(size = 2) +
  xlab(expression(paste(C[chl], " [ppm]"))) +
  ylab(expression(paste("A [", mu, "mol ", m^-2~s^-1, "]"))) +
  theme_bw() +
  NULL

Incorporating leaf temperature using tealeaves

In experiments, leaf temperature can be kept close to air temperature, but in nature, leaf temperature can be quite a bit different than air temperature in the shade depending on environmental and leaf parameters. If use_tealeaves = TRUE, photo() and photosynthesis() will call on the tealeaves package to calculate leaf temperature using an energy balance model.


# NOTE: parallel example is not evaluated because it was causing an issue with CRAN, but you can copy-and-paste the code to run on your own machine.

# You will need to set use_tealeaves = TRUE when making parameters because additional parameters are needed for tealeaves.

bake_par = make_bakepar()
constants  = make_constants(use_tealeaves = TRUE)

enviro_par = make_enviropar(
  replace = list(
    T_air = set_units(seq(288.15, 313.15, 1), K)
    ), use_tealeaves = TRUE
  )

leaf_par = make_leafpar(replace = list(
    g_sc = set_units(c(0.2, 0.4), mol / m^2 / s)
    ), use_tealeaves = TRUE
  )

ph = photosynthesis(leaf_par, enviro_par, bake_par, constants, 
                     use_tealeaves = TRUE, progress = FALSE, 
                   quiet = TRUE, parallel = TRUE)

# Plot temperature and photosynthesis
library(ggplot2)

## Drop units for plotting
ph %<>% 
  mutate_if(~ is(.x, "units"), drop_units) %>%
  mutate(`g[s]` = ifelse(g_sc == 0.2, "low", "high"))

ggplot(ph, aes(T_air, T_leaf, color = `g[s]`)) +
  geom_line(size = 2, lineend = "round") +
  geom_abline(slope = 1, intercept = 0, linetype = "dotted") +
  scale_color_discrete(name = expression(g[s])) +
  xlab(expression(paste(T[air], " [K]"))) +
  ylab(expression(paste(T[leaf], " [K]"))) +
  theme_bw() +
  NULL

ggplot(ph, aes(T_air, A, color = `g[s]`)) +
  geom_line(size = 2, lineend = "round") +
  scale_color_discrete(name = expression(g[s])) +
  xlab(expression(paste(T[leaf], " [K]"))) +
  ylab(expression(paste("A [", mu, "mol ", m^-2~s^-1, "]"))) +
  theme_bw() +
  NULL

Chris Muir

2024-11-23

Minimum worked example

Replace default parameters

Environmental gradients

Parallel processing

Incorporating leaf temperature using tealeaves