Nothing
#' Fitting ACi curves
#'
#' @param data Dataframe for A-Ci curve fitting
#' @param varnames List of variable names. varnames = list(A_net = "A_net",
#' T_leaf = "T_leaf", C_i = "C_i", PPFD = "PPFD", g_mc = "g_mc"), where A_net
#' is net CO2 assimilation, T_leaf is leaf temperature in Celsius, C_i is
#' intercellular CO2 concentration in umol/mol, PPFD is incident irradiance
#' in umol m-2 s-1 (note that it is ASSUMED to be absorbed irradiance, so be
#' sure to adjust according to light absorbance and PSI/PSII partitioning
#' accordingly OR interpret the resultant values of J and J_max with caution),
#' g_mc is mesophyll conductance to CO2 in mol m-2 s-1 Pa-1.
#' @param P Atmospheric pressure in kPa
#' @param fitTPU Should triose phosphate utilization (V_TPU) be fit?
#' @param alpha_g Fraction of respiratory glycolate carbon that is not returned
#' to the chloroplast (von Caemmerer, 2000). If ACi curves show high-CO2
#' decline, then this value should be > 0.
#' @param R_d_meas Measured value of respiratory CO2 efflux in umol m-2 s-1.
#' Input value should be positive to work as expected with the equations.
#' @param useR_d Use a measured value of R_d? Set to TRUE if using R_d_meas.
#' @param useg_mc Use mesophyll conductance? Set to TRUE if specifying g_mc
#' in varnames above.
#' @param useg_mct Use mesophyll conductance temperature response? Set to TRUE
#' if using a temperature response of mesophyll conductance.
#' @param usegamma_star Specify gamma_star value? If FALSE, uses a temperature
#' response function with Nicotiana tabacum defaults from Bernacchi et al.
#' 2001.
#' @param useK_M Specify K_M? If FALSE, uses an Arrhenius temperature response
#' function with Nicotiana tabacum defaults from Bernacchi et al. 2001.
#' @param useK_C_K_O Use individual carboxylation/oxygenation constants for
#' rubisco? If TRUE, need to specify values at 25C and activation energy for
#' the Arrhenius temperature response function.
#' @param alpha Quantum yield of CO2 assimilation
#' @param theta_J Curvature of the photosynthetic light response curve
#' @param gamma_star25 gamma_star at 25C in umol mol-1
#' @param Ea_gamma_star Activation energy of gamma_star in J mol-1
#' @param K_M25 Michaelis-Menten constant for rubisco at 25C
#' @param Ea_K_M Activation energy for K_M in J mol-1
#' @param g_mc25 Mesophyll conductance at 25C in mol m-2 s-1
#' @param Ea_g_mc Activation energy of g_mc in J mol-1
#' @param K_C25 Michaelis-Menten constant for rubisco carboxylation at 25C
#' @param Ea_K_C Activation energy for K_C in J mol-1
#' @param K_O25 Michaelis-Menten constant for rubisco oxygenation at 25C
#' @param Ea_K_O Activation energy for K_O in J mol-2
#' @param Oconc O2 concentration in %. Used with P to calculate
#' intracellular O2 when using K_C_K_O
#' @param gamma_star_set Value of gamma_star to use (in ppm) if
#' usegamma_star = TRUE
#' @param K_M_set Value of K_M to use if useK_M = TRUE
#' @param ... Other arguments to pass on
#'
#' @return fit_aci_response fits ACi curves using an approach similar to
#' Gu et al. 2010. Iterates all possible C_i transition points and checks for
#' inadmissible curve fits. If no curves are admissible (either due to poor data
#' or poor assumed parameters), the output will include a dataframe of NA values.
#' Default parameters are all from Bernacchi et al. 2001, 2002.
#'
#' @references
#' Bernacchi CJ, Singsaas EL, Pimentel C, Portis AR, Long SP. 2001. Improved
#' temperature response functions for models of rubisco-limited photosynthesis.
#' Plant Cell Environment 24:253-259.
#'
#' Bernacchi CJ, Portis AR, Nakano H, von Caemmerer S, Long SP. 2002. 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.
#'
#' Gu L, Pallardy SG, Tu K, Law BE, Wullschleger SD. 2010. Reliable estimation
#' of biochemical parameters from C3 leaf photosynthesis-intercellular carbon
#' dioxide response curves. Plant Cell Environment 33:1852-1874.
#'
#' von Caemmerer S. 2000. Biochemical models of leaf photosynthesis. CSIRO
#' Publishing, Collingwood.
#'
#' @importFrom ggplot2 element_blank
#' @importFrom ggplot2 geom_line
#' @importFrom ggplot2 ggplot
#' @importFrom ggplot2 scale_color_manual
#' @importFrom ggplot2 scale_y_continuous
#' @importFrom stats coef
#' @importFrom stats lm
#' @importFrom stats sd
#' @export
#'
#' @examples
#' \donttest{
#' # Read in your data
#' # Note that this data is coming from data supplied by the package
#' # hence the complicated argument in read.csv()
#' # This dataset is a CO2 by light response curve for a single sunflower
#' data <- read.csv(system.file("extdata", "A_Ci_Q_data_1.csv",
#' package = "photosynthesis"
#' ))
#'
#' # Define a grouping factor based on light intensity to split the ACi
#' # curves
#' data$Q_2 <- as.factor((round(data$Qin, digits = 0)))
#'
#' # Convert leaf temperature to K
#' data$T_leaf <- data$Tleaf + 273.15
#'
#' # Fit ACi curve. Note that we are subsetting the dataframe
#' # here to fit for a single value of Q_2
#' fit <- fit_aci_response(data[data$Q_2 == 1500, ],
#' varnames = list(
#' A_net = "A",
#' T_leaf = "T_leaf",
#' C_i = "Ci",
#' PPFD = "Qin"
#' )
#' )
#'
#' # View fitted parameters
#' fit[[1]]
#'
#' # View graph
#' fit[[2]]
#'
#' # View data with modelled parameters attached
#' fit[[3]]
#'
#' # Fit many curves
#' fits <- fit_many(
#' data = data,
#' varnames = list(
#' A_net = "A",
#' T_leaf = "T_leaf",
#' C_i = "Ci",
#' PPFD = "Qin"
#' ),
#' funct = fit_aci_response,
#' group = "Q_2"
#' )
#'
#' # Print the parameters
#' # First set of double parentheses selects an individual group value
#' # Second set selects an element of the sublist
#' fits[[3]][[1]]
#'
#' # Print the graph
#' fits[[3]][[2]]
#'
#' # Compile graphs into a list for plotting
#' fits_graphs <- compile_data(fits,
#' list_element = 2
#' )
#'
#' # Compile parameters into dataframe for analysis
#' fits_pars <- compile_data(fits,
#' output_type = "dataframe",
#' list_element = 1
#' )
#' }
fit_aci_response <- function(data,
varnames = list(
A_net = "A_net",
T_leaf = "T_leaf",
C_i = "C_i",
PPFD = "PPFD",
g_mc = "g_mc"
),
P = 100,
fitTPU = TRUE,
alpha_g = 0,
R_d_meas = NULL,
useR_d = FALSE,
useg_mc = FALSE,
useg_mct = FALSE,
usegamma_star = FALSE,
useK_M = FALSE,
useK_C_K_O = FALSE,
alpha = 0.24,
theta_J = 0.85,
gamma_star25 = 42.75,
Ea_gamma_star = 37830,
K_M25 = 718.40,
Ea_K_M = 65508.28,
g_mc25 = 0.08701,
Ea_g_mc = 0,
K_C25 = NULL,
Ea_K_C = NULL,
K_O25 = NULL,
Ea_K_O = NULL,
Oconc = 21,
gamma_star_set = NULL,
K_M_set = NULL,
...) {
# Locally bind variables - avoids notes on check package
C_i <- NULL
A_model <- NULL
A_carbox <- NULL
A_regen <- NULL
A_tpu <- NULL
A_net <- NULL
PPFD <- NULL
# Set variable names
data$C_i <- data[, varnames$C_i]
data$A_net <- data[, varnames$A_net]
data$PPFD <- data[, varnames$PPFD]
data$T_leaf <- data[, varnames$T_leaf]
outputs <- vector("list", 3)
# Order data by increasing C_i, avoids calculation issues
data <- data[order(data$C_i), ]
# Convert O2 concentration to partial pressure
O <- Oconc * P / 100
# Create grid of possible C_i transition points
ci <- data[order(data$C_i), ]$C_i
nci <- length(ci)
citransitions <- diff(ci) / 2 + ci[-nci]
# Make sure there is a minimum of 3 points for V_cmax fitting
citransitions1 <- citransitions[3:length(citransitions)]
if (!fitTPU) {
citransitions2 <- max(ci) + 1
} else {
citransitions2 <- c(max(ci) + 1, rev(citransitions1))
}
# Create combinations of ci1 and ci2 to fit
citransdf <- expand.grid(ci1 = citransitions1, ci2 = citransitions2)
citransdf <- citransdf[citransdf$ci1 <= citransdf$ci2, ]
# Mesophyll conductance calculations
if (!useg_mc) {
# Assumes g_mc is infinite
data$C <- data$C_i * P / 100
} else {
# Uses measured values of g_mc
data$g_mc <- data[, varnames$g_mc]
data$C <- (data$C_i - data$A_net / data$g_mc) * P / 100
}
if (useg_mct) {
# Calculates g_mc based on a specified temperature response
data$g_mc <- g_mc25 * t_response_arrhenius(
T_leaf = data$T_leaf,
Ea = Ea_g_mc
)
data$C <- (data$C_i - data$A_net / data$g_mc) * P / 100
}
# gamma_star settings
if (!usegamma_star) {
# Calculates gamma_star based on temperature response function
gamma_star <- gamma_star25 * t_response_arrhenius(
T_leaf = mean(data$T_leaf),
Ea = Ea_gamma_star
) * P / 100
} else {
# Uses specified gamma_star, converts to partial pressure
gamma_star <- gamma_star_set * P / 100
}
# K_M settings
if (!useK_M) {
# Calculates K_M based on temperature response
K_M <- K_M25 * t_response_arrhenius(
T_leaf = mean(data$T_leaf),
Ea = Ea_K_M
)
} else {
# Uses specified K_M
K_M <- K_M_set
}
if (useK_C_K_O) {
# Calculates K_M based on temperature responses of K_C and K_O
K_C <- K_C25 * t_response_arrhenius(
T_leaf = mean(data$T_leaf),
Ea = Ea_K_C
)
K_O <- K_O25 * t_response_arrhenius(
T_leaf = mean(data$T_leaf),
Ea = Ea_K_O
)
K_M <- K_C * (1 + O / K_O)
}
# Generate x-variables for linearized prediction of V_cmax, J_max, V_TPU
# Note this is based on the Duursma (2015) approach eto Gu et al. 2010
data$V_cmax_pred <- (data$C - gamma_star) / (data$C + K_M)
data$J_max_pred <- (data$C - gamma_star) / (data$C + 2 * gamma_star)
data$V_TPU_part <- (data$C - gamma_star) / (data$C - (1 + 3 * alpha_g) *
gamma_star)
# Create dataframe for all possible curve fits
poss_fits <- data.frame(matrix(0,
nrow = nrow(citransdf),
ncol = 16
))
# Add column names
colnames(poss_fits) <- c(
"Num", "V_cmax", "V_cmax_se", "J_max",
"J", "J_se", "V_TPU", "V_TPU_se", "R_d", "R_d_se",
"cost", "citransition1", "citransition2",
"V_cmax_pts", "J_max_pts", "V_TPU_pts"
)
# Fit all possible citransition combinations
for (i in seq_len(nrow(citransdf))) {
# Locally bind variables
cost <- NULL
V_cmax_fit <- NULL
J_max_fit <- NULL
V_TPU <- NULL
datc <- NULL
datj <- NULL
datp <- NULL
datcomp <- NULL
# CO2-limited points
datc <- data[data$C_i < citransdf$ci1[i], ]
# RuBP regeneration-limited points
datj <- data[data$C_i > citransdf$ci1[i] &
data$C_i < citransdf$ci2[i], ]
# V_TPU-limited points
datp <- data[data$C_i > citransdf$ci2[i], ]
# Fits V_cmax
if (!useR_d) {
# Fit R_d and V_cmax
fitc <- lm(A_net ~ V_cmax_pred, data = datc)
R_d_fit <- coef(fitc)[[1]]
R_d_se <- summary(fitc)$coefficients[1, 2]
V_cmax_fit <- coef(fitc)[[2]]
V_cmax_se <- summary(fitc)$coefficients[2, 2]
datc$A_gross <- datc$A_net - R_d_fit
} else {
# Use R_d and fit V_cmax
R_d_fit <- -R_d_meas
datc$A_gross <- datc$A_net - R_d_fit
fitc <- lm(A_gross ~ V_cmax_pred - 1, data = datc)
V_cmax_fit <- coef(fitc)[[1]]
V_cmax_se <- summary(fitc)$coefficients[, 2]
}
# Fit J and J_max
if (nrow(datj) > 0) {
datj$A_gross <- datj$A_net - R_d_fit
if (nrow(datj) == 1) {
# Calculates J_max based on one point
J_fit <- 4 * datj$A_gross / datj$J_max_pred
J_se <- NA
J_max_fit <- suppressWarnings(calculate_jmax(mean(data$PPFD),
alpha,
J = J_fit, theta_J
))
} else {
# Calculates J_max based on a linear regression fit
fitj <- lm(A_gross ~ J_max_pred - 1, data = datj)
J_fit <- 4 * coef(fitj)[[1]]
J_se <- summary(fitj)$coefficients[2]
J_max_fit <- suppressWarnings(calculate_jmax(mean(data$PPFD),
alpha,
J = J_fit, theta_J
))
}
} else {
# Assign J_max a value of 10 ^ 6 if there's no RuBP limitation
J_max_fit <- 10^6
J_max_SS <- 0
J_se <- NA
}
# Calculating V_TPU limitations
if (nrow(datp) == 1 && nrow(datj) == 0) {
# Assign V_TPU a value of 1000 if there is only 1 assigned
# point and no RuBP-limited points
datp$A_gross <- datp$A_net - R_d_fit
V_TPU <- 1000
V_TPU_SS <- 0
V_TPU_se <- NA
} else {
# This section covers if there are no V_TPU-limited points
datp$A_gross <- datp$A_net - R_d_fit
V_TPU <- 1000 # same as default in Photosyn
V_TPU_SS <- 0
V_TPU_se <- NA
}
# Calculates V_TPU limitations if there are at least 3 points
# to ensure more reliable fit
if (nrow(datp) > 2) {
datp$A_gross <- datp$A_net - R_d_fit
V_TPU_vals <- (1 / 3) * datp$A_gross / datp$V_TPU_part
V_TPU <- mean(V_TPU_vals)
V_TPU_se <- sd(V_TPU_vals) / sqrt(length(V_TPU_vals))
}
# If V_TPU is fit to be < 0, assign value of 1000. Avoids
# strange issues.
if (V_TPU < 0) {
V_TPU <- 1000
V_TPU_SS <- 0
V_TPU_se <- NA
}
# Calculate rates of photosynthesis for each limitation state and sums of
# squares for model-wise cost function
# CO2 limitations
datc$A <- V_cmax_fit * (datc$C - gamma_star) / (datc$C + K_M) + R_d_fit
datc$SS <- (datc$A - datc$A_net)^2
V_cmax_SS <- sum(datc$SS)
# RuBP limitations
datj$A <- J_fit * (datj$C - gamma_star) /
(4 * datj$C + 8 * gamma_star) + R_d_fit
datj$SS <- (datj$A - datj$A_net)^2
J_max_SS <- sum(datj$SS)
# V_TPU limitations
datp$A <- 3 * V_TPU * (datp$C - gamma_star) /
(datp$C - (1 + 3 * alpha_g) * gamma_star / datp$C) + R_d_fit
datp$SS <- (datp$A - datp$A_net)^2
V_TPU_SS <- sum(datp$SS)
# Calculate cost functions
cost <- 1 / 2 * (V_cmax_SS + J_max_SS + V_TPU_SS)
# Bind calculated data
datcomp <- rbind(datc, datj, datp)
# Add outputs to possible fits
poss_fits$V_cmax[i] <- V_cmax_fit
poss_fits$V_cmax_se[i] <- V_cmax_se
poss_fits$J_max[i] <- J_max_fit
poss_fits$J[i] <- J_fit
poss_fits$J_se[i] <- J_se
poss_fits$V_TPU[i] <- V_TPU
poss_fits$V_TPU_se[i] <- V_TPU_se
poss_fits$R_d[i] <- R_d_fit
poss_fits$R_d_se[i] <- R_d_se
poss_fits$V_cmax_pts[i] <- nrow(datc)
poss_fits$J_max_pts[i] <- nrow(datj)
poss_fits$V_TPU_pts[i] <- nrow(datp)
poss_fits$cost[i] <- cost
poss_fits$citransition1[i] <- citransdf$ci1[i]
poss_fits$citransition2[i] <- citransdf$ci2[i]
} # End curve fitting of all possible C_i transitions
# Select fit with minimized cost function
best_fits <- poss_fits[poss_fits$cost == min(poss_fits$cost), ]
# New segment for sensitivity analysis
# Adds values for assumed constants to the output dataframe
best_fits$alpha <- alpha
best_fits$alpha_g <- alpha_g
best_fits$gamma_star25 <- gamma_star25
best_fits$Ea_gamma_star <- Ea_gamma_star
best_fits$K_M25 <- K_M25
best_fits$Ea_K_M <- Ea_K_M
best_fits$g_mc25 <- g_mc25
best_fits$Ea_g_mc <- Ea_g_mc
best_fits$K_C25 <- K_C25
best_fits$Ea_K_C <- Ea_K_C
best_fits$K_O25 <- K_O25
best_fits$Ea_K_O <- Ea_K_O
best_fits$Oconc <- Oconc
best_fits$theta_J <- theta_J
# Calculate net photosynthetic rates
data$A_carbox <- best_fits$V_cmax * data$C /
(data$C + K_M) * (1 - gamma_star / data$C) + best_fits$R_d
data$A_regen <- calculate_j(
PPFD = mean(data$PPFD), alpha = alpha,
J_max = best_fits$J_max, theta_J = theta_J
) *
(data$C - gamma_star) / (4 * data$C + 8 * gamma_star) + best_fits$R_d
data$A_tpu <- 3 * best_fits$V_TPU /
(1 - 0.5 * (1 + 3 * alpha_g) * (2 * gamma_star / data$C)) *
(1 - gamma_star / data$C) + best_fits$R_d
# Calculate gross photosynthetic rates
data$W_carbox <- data$A_carbox - best_fits$R_d
data$W_regen <- data$A_regen - best_fits$R_d
data$W_tpu <- data$A_tpu - best_fits$R_d
# Create empty variable for modelled CO2 assimilation
data$A_model <- rep(NA, nrow(data))
# To avoid issues with graphing and cases where W_j drops below W_c
# at very low CO2, for the first few points, A_model is calculated
# with W_c only - plantecophys took an approach where Aj was fixed
# on the graph until a certain C_i to avoid this same issue.
for (i in 1:(best_fits$V_cmax_pts - 2)) {
data$A_model[i] <- data$W_carbox[i] + best_fits$R_d
}
if (best_fits$V_TPU == 1000) {
for (i in (best_fits$V_cmax_pts - 1):nrow(data)) {
data$A_model[i] <- min(data$W_carbox[i], data$W_regen[i]) + best_fits$R_d
}
} else {
for (i in (best_fits$V_cmax_pts - 1):nrow(data)) {
data$A_model[i] <- min(
data$W_carbox[i], data$W_regen[i],
data$W_tpu[i]
) + best_fits$R_d
}
}
# Assign best fit to output element 1
outputs[[1]] <- best_fits
# Assign graph to output element 2
outputs[[2]] <- ggplot(data, aes(x = C_i, y = A_model)) +
scale_y_continuous(limits = c(
min(c(data$A_model, data$A_net)) - 3,
max(c(data$A_model, data$A_net)) + 3
)) +
geom_line(aes(color = "black"), linewidth = 4) +
geom_line(aes(y = A_carbox, color = "blue"), linewidth = 2) +
geom_line(aes(y = A_regen, color = "orange"), linewidth = 2) +
geom_line(aes(y = A_tpu, color = "red"), linewidth = 2) +
geom_point(aes(y = A_net),
color = "black", fill = "white",
size = 2, shape = 21
) +
scale_color_manual(
labels = c("Amod", "Ac", "Aj", "Ap", "Anet"),
values = c(
"black", "blue", "orange",
"red", "white"
)
) +
theme_bw() +
theme(legend.title = element_blank())
# Assign dataframe to output element 3
outputs[[3]] <- data
# Add names to list output
names(outputs) <- c("Fitted Parameters", "Plot", "Data")
# Return output list
return(outputs)
} # End function
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.