photofit: Estimate parameters from photosynthesis curves

Documented in fitACi

# Simple ACi model (no gm) ------------------------------------------------

ACi_limitations = function(pars, data) {
  Vcmax = pars[["Vcmax"]]
  KmC   = pars[["KmC"]]
  KmO   = pars[["KmO"]]
  CiStar = pars[["CiStar"]]
  J     = pars[["J"]]
  TPU   = pars[["TPU"]]
  Rd    = pars[["Rd"]]
  Ci = data$Ci
  Ac = Vcmax*(Ci - CiStar)/(Ci + KmC*(1 + 210/KmO))
  Aj = J*(Ci - CiStar)/(4*Ci + 8*CiStar)
  At = 3*TPU
  A = cbind(Ac, Aj, At) - Rd
}


ACi = function(pars, data) {
  lims = ACi_limitations(pars, data)
  pmin(pmin(lims[,1], lims[,2]), lims[,3])
}

gr_ACi = function(pars, data, parnames) {

  Vcmax = pars[["Vcmax"]]
  KmC   = pars[["KmC"]]
  KmO   = pars[["KmO"]]
  CiStar = pars[["CiStar"]]
  J     = pars[["J"]]
  TPU   = pars[["TPU"]]
  Rd    = pars[["Rd"]]

  # Need to know which process is limiting
  Ci = data$Ci
  n = length(Ci)
  O2 = 210
  n = length(Ci)
  Ac = Vcmax*(Ci - CiStar)/(Ci + KmC*(1 + 210/KmO))
  Aj = J*(Ci - CiStar)/(4*Ci + 8*CiStar)
  At = rep(3*TPU, n)

  # Figure out which are the limiting steps
  limiting = map_dbl(1:n, ~which.min(c(Ac[.x], Aj[.x], At[.x])))

  # Derivatives for Ac-limited photosynthesis
  .e1 = 1 + O2/KmO
  .e2 = Ci + KmC * .e1
  .e3 = Ci - CiStar
  .e4 = .e2^2
  dA_dVcmax = ifelse(limiting == 1, .e3/.e2, 0)
  dA_dKmC   = ifelse(limiting == 1, -(Vcmax*.e1 * .e3/.e4),  0)
  dA_dKmO   = ifelse(limiting == 1, O2*(KmC * Vcmax * .e3/(KmO^2*.e4)),  0)
  dA_dCiStar_Vcmax = -(Vcmax/.e2)

  # Derivatives for Aj-limited photosynthesis
  .e2 <- 4 * Ci + 8 * CiStar
  .e3 <- (Ci - CiStar)/.e2
  dA_dJ = ifelse(limiting == 2, .e3, 0)
  dA_dCiStar_J = -(J * (1 + 8 * .e3)/.e2)

  # Combine and others
  dA_dCiStar = ifelse(limiting == 1, dA_dCiStar_Vcmax, ifelse(limiting == 2, dA_dCiStar_J, 0))
  dA_dRd    = rep(-1, n)
  dA_dTPU   = ifelse(limiting == 3, 3, 0)

  cbind(Vcmax = dA_dVcmax, KmC = dA_dKmC, KmO = dA_dKmO, CiStar = dA_dCiStar,
        J = dA_dJ, TPU = dA_dTPU, Rd = dA_dRd)[,parnames]
}



# Generate list of parameters for prior distributions

#' @export
priors_ACi = function(Vcmax = prior(50, 50),
                      CiStar = prior(15, 15),
                      KmC = prior(230, 50),
                      KmO = prior(195, 40),
                      J = prior(100, 100),
                      TPU = prior(10, 10),
                      Rd = prior(1,1),
                      sigma = prior(1,1)) {
  list(Vcmax = Vcmax, CiStar = CiStar, KmC = KmC, KmO = KmO, J = J, TPU = TPU,
       Rd = Rd, sigma = sigma)
}


#' @export
initial_ACi = function(Vcmax = 50, CiStar = 15, KmC = 230, KmO = 195, J = 100,
                       TPU = 10, Rd = 0.75, sigma = 0.9) {
  c(Vcmax = Vcmax, CiStar = CiStar, KmC = KmC, KmO = KmO, J = J, TPU = TPU,
    Rd = Rd, sigma = sigma)
}


# Constraint on each parameter
parmodesACi = c(Vcmax = ">0", CiStar = ">0", KmC = ">0", KmO = ">0", J = ">0",
                TPU = ">0", Rd = ">0", sigma = ">0")

#' Fit A-Ci curve
#'
#' \code{fitACi} will fit a steady-state model of photosynthesis to a response curve
#' measuring CO₂ assimilation and intercellular CO₂ molar fraction.
#'
#'
#' A classic FvCB model that implicitly assumes infinite mesophyll conductance is used.
#' Net CO₂ assimilation for every measurement point is calculated as the minimum of three
#' potential rates limited by Rubisco kinetics (\code{Ac}), electron transport (\code{Aj})
#' and triose phosphate utilisation (\code{At}):
#'
#' \code{Ac = Vcmax*(Ci - CiStar)/(Ci + KmC*(1 + 210/KmO))}
#'
#' \code{Aj = J*(Ci - CiStar)/(4*Ci + 8*CiStar)}
#'
#' \code{At = 3*TPU}
#'
#' The parameters that can be fitted are:
#'
#' \itemize{
#'  \item \code{Vcmax} Maximum rate of carboxylation (μmol m⁻² s⁻¹).
#'  \item \code{CiStar} Intercellular CO₂ compensation point (μmol mol⁻¹).
#'  \item \code{Kmc} Michaelis-Menten constant of carboxylation with respect to CO₂ (μmol mol⁻¹).
#'  \item \code{Kmo} Michaelis-Menten constant of carboxylation with respect to O₂ (mmol mol⁻¹).
#'  \item \code{J} Maximum rate of electron transport under measurement conditions (μmol m⁻² s⁻¹).
#'  \item \code{TPU} Maximu rate of triose phosphate utilisation (μmol m⁻² s⁻¹).
#'  \item \code{Rd} Mitochondrial respiration in the light (μmol m⁻² s⁻¹).
#'  \item \code{sigma} Standard deviation of measurement error (μmol m⁻² s⁻¹).
#' }
#'
#' The default initial values are generated by \code{\link{initial_ACi}} whereas
#' the mean and standard deviations of Gaussian prior distributions are generated
#' by \code{\link{priors_ACi}}.
#'
#' Values allowed for the \code{method} argument are \code{"quad"} (quadratic
#' approximation) and \code{"sir"} (sampling importance resampling).
#'
#' All methods to approximate the posterior start with a calculated of the
#' maximum a posteriori. The argument \code{algorithm} corresponds to the \code{method}
#' argument of the function \code{\link{optimx}} (see help for this function for
#' a list of possible algorithms).
#'
#' @param data Measurements as a \code{data.frame} with two columns called \code{Ci} and \code{A}
#'   representing intercellular CO₂ molar fraction and net CO₂ assimilation, respectively.
#' @param priors List of parameters for the prior distributions as generated by
#'   \code{\link{priors_ACi}}.
#' @param pars Vector of initial values of parameters as generated by
#'   \code{\link{initial_ACi}}.
#' @param fixed Vector with names of parameters that will not be fitted to the
#'   data but fixed to their initial values. Same names as in \code{\link{priors_ACi}}
#'   and \code{\link{initial_ACi}} should be used.
#' @param method Method to approximate the posterior distribution of parameters. See Details.
#' @param algorithm Name of non-linear optimization algorithm used to calculate
#'   maximum a posteriori. See Details.
#'
#' @seealso initial_ACi, priors_ACi
#' @examples
#' data = data.frame(Ci = c(245.39, 164.70, 87.59, 50.22, 332.55, 333.15, 421.98,
#'                   515.53, 612.01, 707.88, 904.34, 1099.13, 1295.58, 1485.47),
#'            A = c(8.90, 5.72, 1.94, -0.18, 11.38, 11.49, 13.49, 14.66, 15.10,
#'                  15.61, 15.49, 15.43, 14.97, 15.51))
#'
#' # Quadratic approximation
#' fitQuad = fitACi(data = data)
#' summary(fitQuad)
#' plot(fitQuad)
#' plot(fitQuad, type = "limitations")
#' plot(fitQuad, type = "marginal")
#' plot(fitQuad, type = "pairs")
#'
#' # Quadratic approximation - Avoiding TPU
#' fitQuad2 = fitACi(data = data, pars = initial_ACi(TPU = 100),
#'                  fixed = c("KmC", "KmO", "TPU"))
#' summary(fitQuad2)
#' plot(fitQuad2)
#' plot(fitQuad2, type = "pairs")
#'
#' # SIR approximation
#' fitSIR = fitACi(data = data, pars = initial_ACi(TPU = 100),
#'                  fixed = c("KmC", "KmO", "TPU"))
#' summary(fitSIR)
#' plot(fitSIR)
#' plot(fitSIR, type = "marginal")
#' plot(fitSIR, type = "pairs")
#' @export
fitACi = function(data, priors = priors_ACi(), pars = initial_ACi(),
                  fixed = c("KmC", "KmO"), method = "quad", algorithm = "BFGS",
                  nStart = 25, ...) {
  # Fixed parameters not to be fitted to the data
  fixed_pars = list()
  parmodes = parmodesACi
  for(name in fixed) {
    fixed_pars[[name]] = pars[[name]]
    pars = pars[-which(names(pars) == name)]
    priors = priors[-which(names(priors) == name)]
    parmodes = parmodes[-which(names(parmodes) == name)]
  }
  # Create a structure with all the information required to fit the model
  model = structure(list(data = as.data.frame(data), pars = pars,
                         priors = priors, model  = ACi,
                         limitations = ACi_limitations, gradient = gr_ACi,
                         parmodes = parmodes, fixed = fixed_pars,
                         predictor = "Ci"), class = "ACi")
  fit = fitModel(model = model, pars = pars, method = method,
                 algorithm = algorithm, nStart = nStart, ...)
  class(fit) = c("ACi", "photofit")
  return(fit)
}