Theory

Theory

The P-model predicts an optimal ratio of \(c_i : c_a\), termed as \(\chi\), that balances the costs associated with maintaining the transpiration (\(E\)) stream and the carboxylation capacity \(V_{\text{cmax}}\). It can therefore be used to simulate the acclimation of the photosynthetic machinery to its environment - a mechanism that happens at a time scale of several days to months. At its core, it provides a solution for the optimality criterium \[ a \; \frac{\partial (E/A)}{\partial \chi} = -b \; \frac{\partial (V_{\mathrm{cmax}}/A)}{\partial \chi} \;\;\;\;\;\;\;\;\;\;\;\;(1) \] The optimal \(\chi\) solves the above equation and, with \(E = 1.6 g_s D\), \(A = g_s (1-\chi)\), and using the Rubisco-limited assimilation rate: \[ A = A_C = V_{\mathrm{cmax}} \; \frac{\chi\;c_a-\Gamma^{\ast}}{\chi\;c_a + K} \] is given by: \[ \chi = \frac{\Gamma^{\ast}}{c_a} + \left(1- \frac{\Gamma^{\ast}}{c_a}\right)\;\frac{\xi}{\xi + \sqrt{D}} \] with \[ \xi = \sqrt{\frac{b(K+\Gamma^{\ast})}{1.6\;a}} \] The unit cost ratio \(b/a\) is also referred to as \(\beta\).

P-model run

So much for the theory. Let’s run the P-model, without \(J_{\text{max}}\) limitation, for one set of inputs, being temperature, PPFD, VPD, CO\(_2\), elevation, and fAPAR.

To do so, run the rpmodel() function from the rsofun package:

library(rpmodel)
library(dplyr)
# modified seq() function to get a logarithmically spaced sequence
lseq <- function(from=1, to=100000, length.out=6) {
  exp(seq(log(from), log(to), length.out = length.out))
}

## Set parameters
beta <- 146          # unit cost ratio a/b
c_cost <- 0.41
gamma <- 0.105       # unit cost ratio c/b
kphio <- 0.05        # quantum yield efficiency
c_molmass <- 12.0107 # molar mass, g / mol

## Define environmental conditions
tc <- 20             # temperature, deg C
ppfd <- 300          # mol/m2/d
vpd  <- 1000         # Pa
co2  <- 400          # ppm
elv  <- 0            # m.a.s.l.
fapar <- 1           # fraction

out_analytical_none <- rpmodel::rpmodel(
  tc             = tc,
  vpd            = vpd,
  co2            = co2,
  elv            = elv,
  kphio          = kphio,
  beta           = beta,
  fapar          = fapar,
  ppfd           = ppfd,
  method_optci   = "prentice14",
  method_jmaxlim = "none",
  do_ftemp_kphio = FALSE
  )

out_analytical_jmaxlim <- rpmodel::rpmodel(
  tc             = tc,
  vpd            = vpd,
  co2            = co2,
  elv            = elv,
  kphio          = kphio,
  beta           = beta,
  fapar          = fapar,
  ppfd           = ppfd,
  method_optci   = "prentice14",
  method_jmaxlim = "wang17",
  do_ftemp_kphio = FALSE
  )

The function returns a list of variables (see also man page by ?rpmodel), including \(V_{\mathrm{cmax}}\), \(g_s\), and all the parameters of the photosynthesis model (\(K\), \(\Gamma^{\ast}\)), which are all internally consistent, as can be verified for… \[ c_i = c_a - A / g_s = \chi c_a \]

print( out_analytical_none$ci )

## [1] 28.14209

print( out_analytical_none$ca - (out_analytical_none$gpp / c_molmass) / out_analytical_none$gs )

## [1] 28.14209

print( out_analytical_none$ca * out_analytical_none$chi )

## [1] 28.14209

Yes.

And for… \[ A = V_{\text{cmax}} \frac{c_i-\Gamma^{\ast}}{c_i + K} = \phi_0 I_{\text{abs}} \frac{c_i-\Gamma^{\ast}}{c_i + 2 \Gamma^{\ast}} = g_s (c_a - c_i) \]

print( out_analytical_none$gpp / c_molmass )

## [1] 10.68456

print( out_analytical_none$vcmax * (out_analytical_none$ci - out_analytical_none$gammastar) / (out_analytical_none$ci + out_analytical_none$kmm ))

## [1] 10.68456

print( out_analytical_none$gs * (out_analytical_none$ca - out_analytical_none$ci) )

## [1] 10.68456

print( kphio * ppfd * fapar * (out_analytical_none$ci - out_analytical_none$gammastar) / (out_analytical_none$ci + 2 * out_analytical_none$gammastar ))

## [1] 10.68456

Yes.

Numerical solution without Jmax cost

Instead of formulating the optimality criterium (Eq. 1) with respect to equality in marginal costs (derivative w.r.t. \(\chi\)), one may also write this as a minimisation of the combined costs. With \(\beta=b/a\), the following criterium is equivalent to Eq. 1: \[ E/A + \beta V_{\mathrm{cmax}}/A = min. \;\;\;\;\;\;\;\;(2) \] This formulation has the advantage that we can easily apply a numerial search algorithm to find the solution in \(\chi\). Why bother? I think it has an advantage as it allows for a more modular treatment, e.g., of transpiration as a function of \(g_s\). But for now, let’s stick with \(E=1.6g_sD\). Eq. 2 can be written out as: \[ \frac{1.6\;D}{c_a(1-\chi)} + \beta \frac{\chi c_a + K}{\chi c_a - \Gamma^{\ast}} = min. \]

To find the minimum, let’s define the numerical optimization function using the L-BFGS-B algorithm implemented in the optimr package:

calc_optimal_chi_num <- function( kmm, gammastar, ns_star, ca, vpd, beta ){
  #-----------------------------------------------------------------------
  # Input:    - float, 'kmm' : Pa, Michaelis-Menten coeff.
  #           - float, 'ns_star'  : (unitless) viscosity correction factor for water
  #           - float, 'vpd' : Pa, vapor pressure deficit
  # Output:   float, ratio of ci/ca (chi)
  # Features: Returns an estimate of leaf internal to ambient CO2
  #           partial pressure following the "simple formulation".
  # Depends:  - kc
  #           - ns
  #           - vpd
  #-----------------------------------------------------------------------
  maximise_this_chi <- function( chi, kmm, gammastar, ns_star, ca, vpd, beta ){
    out <- 1.6 * ns_star * vpd / (ca * (1.0 - chi)) + beta * (chi * ca + kmm)/(chi * ca - gammastar)
    return(out)
  }

  out_optim <- optimr::optimr(
    par       = 0.7,
    lower     = 0.1,
    upper     = 0.99,
    fn        = maximise_this_chi,
    kmm       = kmm,
    gammastar = gammastar,
    ns_star   = ns_star,
    ca        = ca,
    vpd       = vpd,
    beta      = beta,
    method    = "L-BFGS-B",
    control   = list( maxit = 100, maximize = TRUE )
    )

  return(out_optim$par)
  return(list(jmax_opt = NA, chi_opt = out_optim$par, aj = NA))
}

Now, let’s find the optimal \(\chi\) numerically using above function with the same parameters as calculated inside and returned by the P-model.

chi_opt <- calc_optimal_chi_num( out_analytical_none$kmm, out_analytical_none$gammastar, out_analytical_none$ns_star, out_analytical_none$ca, vpd, beta)
print(chi_opt)

## [1] 0.6943512

As a check: is this the same as returned by the P-model?

print(out_analytical_none$chi)

## [1] 0.694352

Practically, yes. A slight difference is due to the limited number of iterations (maxit=100).

Let’s call this \(\chi^{\ast}\). For the P-model, it is then assumed that \(A=A_J\) with \[ A_J = \phi_0 \; I_{\mathrm{abs}}\;\frac{c_a \chi^{\ast} - \Gamma^{\ast}}{c_a \chi^{\ast} + 2\Gamma^{\ast}} \] It is further assumed that the light-limited and Rubisco-limited assimilation rates are equal for average conditions (which implies that the arguments to the P-model run should represent average conditions). This assumption is also called the “coordination hypothesis.” \[ A_J = A_C \] This allows for solving for \(V_{\text{cmax}}\): \[ V_{\text{cmax}} = \phi_0 \; I_{\mathrm{abs}}\;\frac{c_a \chi^{\ast} + K}{c_a \chi^{\ast} + 2 \Gamma^{\ast} } \] This is what the P-model returns as can be verified:

print( out_analytical_none$vcmax)

## [1] 31.98167

print( kphio * fapar * ppfd * (out_analytical_none$ci + out_analytical_none$kmm) / (out_analytical_none$ci + 2 * out_analytical_none$gammastar ) )

## [1] 31.98167

Response of \(\chi\)

Create data frames where one variable is varied at a time, holding all others constant, and apply rpmodel. Plot response of \(\chi\):

Response of \(A_J\)

Same as above, but showing \(A_J\) (no Jmax limitation):

## Warning: Removed 1 rows containing missing values (geom_path).

Numerical solution with Jmax cost

The method described by Wang et al. (2017) to account for \(J_{\mathrm{max}}\) limitation of the assimilation rate is to first calculate the optimal \(\chi\) as described above, and multiply the “unlimited” rate with a term \(L\) following Smith et al. (1937) to limit \(A_J\): \[ A_J = \varphi_0 \; I_{\mathrm{abs}} \; \frac{c_i - \Gamma^{\ast}}{c_i + 2\Gamma^{\ast}} \; \underbrace{ \frac{1}{\sqrt{1+ \left( \frac{4\;\varphi_0\;I_{\mathrm{abs}}}{J_{\mathrm{max}}} \right)^{2}}} }_{L} \] By using an additional optimality criterion \[ \frac{\partial A_J}{\partial J_{\mathrm{max}}} = c \;, \] we can derive optimal \(J_{\mathrm{max}}\) and \(A_J\) can be quantified. By taking the integral of the above optimality criterion, it can be “translated” into a equivalent maximisation criterion: \[ A_J - c \; J_{\mathrm{max}} = max. \] We can now implement a numerical version of the \(J_{\mathrm{max}}\)-limited P-model by…

• First (numerically) derive optimal \(\chi\) based on the original criterion (as describd in the previous section).
• Then calculate the \(J_{\mathrm{max}}\)-limited \(A_J\) using Smith (1937).
• Maximise \(A_J - c \; J_{\mathrm{max}}\).

calc_optimal_jmax_num <- function( kmm, gammastar, ns_star, ca, vpd, chi = NULL, ppfd, kphio, beta, c_cost){
  #-----------------------------------------------------------------------
  # Input:    - float, 'kmm' : Pa, Michaelis-Menten coeff.
  #           - float, 'ns_star'  : (unitless) viscosity correction factor for water
  #           - float, 'vpd' : Pa, vapor pressure deficit
  # Output:   float, ratio of ci/ca (chi)
  # Features: Returns an estimate of leaf internal to ambient CO2
  #           partial pressure following the "simple formulation".
  # Depends:  - kc
  #           - ns
  #           - vpd
  #-----------------------------------------------------------------------
  maximise_this_chi <- function( chi, kmm, gammastar, ns_star, ca, vpd, beta ){
    out <- 1.6 * ns_star * vpd / (ca * (1.0 - chi)) + beta * (chi * ca + kmm)/(chi * ca - gammastar)
    return(out)
  }
  
  maximise_this_jmax <- function( jmax, chi, kmm, gammastar, ns_star, ca, vpd, kphio, ppfd, beta, c_cost ){
    
    mj <- (chi * ca - gammastar) / (chi * ca + 2 * gammastar)
    assim_jmaxlim <- kphio * ppfd * mj * 1 / sqrt(1 + ((4 * kphio * ppfd)/jmax)^2)
    
    # ## extended optimality criterion
    # out <- 1.6 * ns_star * vpd / (ca * (1.0 - chi)) + beta * (chi * ca + kmm)/(chi * ca - gammastar) + c_cost * jmax / assim_jmaxlim
  
    ## second optimality criterion expressed as a maximisation
    out <- -(assim_jmaxlim - c_cost * jmax)
    
    return(out)
  }
  
  if (is.null(chi)){
    ## first get optimal chi independent of jmax
    out_optim_chi <- optimr::optimr(
      par       = 0.7,
      lower     = 0.1,
      upper     = 0.99,
      fn        = maximise_this_chi,
      kmm       = kmm,
      gammastar = gammastar,
      ns_star   = ns_star,
      ca        = ca,
      vpd       = vpd,
      beta      = beta,
      method    = "L-BFGS-B",
      control   = list( maxit = 100, maximize = TRUE )
      )
    chi <- out_optim_chi$par
  }
    
  ## Get optimal jmax given optimal chi
  out_optim_jmax <- optimr::optimr(
    par       = 10,
    lower     = 0.1,
    upper     = 500,
    fn        = maximise_this_jmax,
    chi       = chi,
    kmm       = kmm,
    gammastar = gammastar,
    ns_star   = ns_star,
    ca        = ca,
    vpd       = vpd,
    kphio     = kphio,
    ppfd      = ppfd,
    beta      = beta,
    c_cost    = c_cost,
    method    = "L-BFGS-B",
    control   = list( maxit = 100, maximize = TRUE )
    )

  jmax <- out_optim_jmax$par
  mj <- (chi * ca - gammastar) / (chi * ca + 2 * gammastar)
  assim_jmaxlim <- kphio * ppfd * mj * 1 / sqrt(1 + ((4 * kphio * ppfd)/jmax)^2)
  
  return(list(jmax_opt = jmax, chi_opt = chi, aj = assim_jmaxlim))
}

Now, let’s find the optimal \(\chi\) numerically using above function with the same parameters as calculated inside and returned by the P-model.

jmax_chi_opt <- calc_optimal_jmax_num( 
  kmm = out_analytical_none$kmm, 
  gammastar = out_analytical_none$gammastar, 
  ns_star = out_analytical_none$ns_star, 
  ca = out_analytical_none$ca, 
  vpd = vpd, 
  chi = chi_opt, 
  ppfd = ppfd, 
  kphio = kphio, 
  beta = beta, 
  c_cost = c_cost / 4.0
  )

print(jmax_chi_opt$aj)

## [1] 5.930098

print(out_analytical_jmaxlim$gpp / 12.0107)

## [1] 5.930102

Response to the environment

Create data frames where one variable is varied at a time, holding all others constant, and apply rpmodel.

## Warning: Removed 1 rows containing missing values (geom_path).

Optimising \(g_s\) and \(V_{\text{cmax}}\) instead of \(\chi\)

At the core of the assimilation-transpiration trade-off are the quantities \(g_s\) and \(V_{\text{cmax}}\). The combined effect of their magnitudes determines \(\chi\). Why not simutaneously optimize \(g_s\) and \(V_{\text{cmax}}\)? We can write Eq. 2 also as: \[ (1.6g_sD + \beta V_{\text{cmax}})/A_C = min. \;\;\;\;\;\;\;(3) \] and implement this in another numerical search function. Note that \(c_i\) can be calculated from solving the equation system \[ A = V_{\mathrm{cmax}} \; \frac{c_i-\Gamma^{\ast}}{c_i + K} \\ A = g_s(c_a - c_i) \] which leads to a quadratic equation for \(c_i\) where we will only take the positive real part (Re(root_ci)[which(Re(root_ci)>0)]). This can be implemented as follows.

calc_optimal_gs_vcmax <- function( par, args, maximize=FALSE, return_all=FALSE, method_opt = "LC", cost_scalar = 0.001717135){

  kmm       <- args[1]
  gammastar <- args[2]
  ns_star   <- args[3]
  ca        <- args[4]
  vpd       <- args[5]
  beta      <- args[6]

  vcmax <- par[1]
  gs    <- par[2]

  ## Rubisco is limiting
  ## Solve Eq. system
  ## A = gs (ca- ci)
  ## A = Vcmax * (ci - gammastar)/(ci + Kmm)
  
  ## This leads to a quadratic equation:
  ## A * ci^2 + B * ci + C  = 0
  ## 0 = a + b*x + c*x^2
  
  ## with
  A <- -1.0 * gs
  B <- gs * ca - gs * kmm - vcmax
  C <- gs * ca * kmm + vcmax * gammastar
  
  ci <- QUADM(A, B, C)
  a_c <- vcmax * (ci - gammastar) / (ci + kmm)

  # ## XXX from Wolfram Alpha:
  # a_c <- 1/2 * (-sqrt(4 * gs * vcmax * (gammastar - ca) + (ca * gs + gs * kmm + vcmax)^2) + ca * gs + gs * kmm + vcmax) 
  # ci  <- (sqrt(4 * gs * vcmax * (gammastar - ca) + (ca * gs + gs * kmm + vcmax)^2) + ca * gs - gs * kmm - vcmax)/(2 * gs)
  
  ## only cost ratio is defined. for this here we need absolute values. Set randomly
  cost_transp <- 1.6 * ns_star * gs * vpd
  cost_vcmax  <- beta * vcmax

  if (method_opt=="LC"){
    ## Option A: This is equivalent to the P-model with its optimization of ci:ca.
    net_assim <- -(cost_transp + cost_vcmax) / a_c
  } else if (method_opt=="PM"){
    net_assim <- a_c - cost_scalar * (cost_transp + cost_vcmax)
  }
    
  if (maximize) net_assim <- -net_assim

  if (return_all){
    return( tibble( vcmax_num = vcmax, gs_num = gs, ci_num = ci, chi_num = ci/ca, a_c_num = a_c, cost_transp_num = cost_transp, cost_vcmax_num = cost_vcmax, net_assim_num = net_assim  ) )
  } else {
    return( net_assim )
  }

}

Now that the target function is defined, we can actually numerically search for the optimal \(V_{\text{cmax}}\) and \(g_s\). We need to provide the algorithm with starting values, taken here as the \(V_{\text{cmax}}\) and \(g_s\) returned by our initial P-model call.

wrap_calc_optimal_gs_vcmax <- function( kmm, gammastar, ns_star, ca, vpd, beta, vcmax_start, gs_start, method_opt = "LC" ){
  
  out_optim <- optimr::optimr(
    par        = c( vcmax_start, gs_start ), # starting values
    lower      = c( vcmax_start*0.0001, gs_start*0.001 ),
    upper      = c( vcmax_start*1000,   gs_start*1000 ),
    fn         = calc_optimal_gs_vcmax,
    args       = c( kmm, gammastar, ns_star, ca, vpd, beta ),
    method_opt = method_opt,
    method     = "L-BFGS-B",
    maximize   = TRUE,
    control    = list( maxit = 1000 )
    )
  varlist <- calc_optimal_gs_vcmax( par=out_optim$par, args=c(kmm, gammastar, ns_star, ca, vpd, beta), maximize=FALSE, return_all=TRUE, method_opt = method_opt )
  return(varlist)
}

Note that, given that the optimality criterium implemented here (Eq. 3) embodies the same idea and is mathematically identical as Eq. 1, this should yield the same \(\chi\) (which can be back-calculated from the optimised \(V_{\text{cmax}}\) and \(g_s\)). Let’s check if this is the case.

## [1] "Optimal chi from Vcmax-gs optimization: 0.694352013366217"

## [1] "Optimal chi from P-model              : 0.694352013202358"

Indeed, this is the same.

Response to the environment

Ok. That’s identical.

Now, we have an optimised \(V_{\text{cmax}}\) (lets call it \(V_{\text{cmax}}^{\ast}\)), and we can have an optimised \(c_i\) which is consistent with the former. Given this, we can take the same step as is done in the P-model and assume that \(A_J=A_C\), which leads to \[ V_{\text{cmax}} = \phi_0 \; I_{\mathrm{abs}}\;\frac{c_a \chi^{\ast} + K}{c_a \chi^{\ast} + 2 \Gamma^{\ast} } \] To verify one again:

print( paste( "optimal Vcmax:", varlist$vcmax_num ) )

## [1] "optimal Vcmax: 31.9816667938778"

print( paste( "Vcmax^P:", kphio * fapar * ppfd * (varlist$ci_num + varlist$kmm) / (varlist$ci + 2 * varlist$gammastar )) )

## Warning: Unknown or uninitialised column: 'kmm'.

## Warning: Unknown or uninitialised column: 'ci'.

## Warning: Unknown or uninitialised column: 'gammastar'.

## [1] "Vcmax^P: "

print( paste( "Vcmax from P-model:", out_analytical_none$vcmax ) )

## [1] "Vcmax from P-model: 31.9816667938778"

It’s important to note that calculating \(V_{\text{cmax}}\) like this implies the assumption of \(A_J = A_C\), while solving for simultaneously optimal \(V_{\text{cmax}}\) and \(g_s\) doesn’t use this assumption. Nevertheless, the results are identical, which may not be immediately obvious.

Alternative Jmax limitation

Now, we can try to account for the cost of Jmax, gs, and Vcmax in one single optimality criterion: \[ (1.6g_sD + \beta V_{\text{cmax}} + c J_{\text{max}})/A_C = min. \;\;\;\;\;\;\;(3) \]

This is implemented as follows:

Let’s compare the results for \(\chi\) with this modified version to the original one:

source("~/rpmodel/vignettes_add/calc_optimal_gs_vcmax_jmax.R")
varlist_ll <- calc_optimal_gs_vcmax_jmax( out_analytical_none$kmm, out_analytical_none$gammastar, out_analytical_none$ns_star, out_analytical_none$ca, vpd, ppfd, fapar, kphio, beta, c_cost, out_analytical_none$vcmax, out_analytical_none$gs, 40 )

## Invoke function with optimised Vcmax and gs (in out_optim$par), now returning all variables
print( paste("Optimal chi from P-model:", out_analytical_none$chi ) )

## [1] "Optimal chi from P-model: 0.694352013202358"

print( paste("Optimal chi from Vcmax-gs optimization:", jmax_chi_opt$chi_opt ) )

## [1] "Optimal chi from Vcmax-gs optimization: 0.694351221313652"

print( paste("Optimal chi from Vcmax-gs opt. with light limit.:", varlist_ll$chi ) )

## Warning: Unknown or uninitialised column: 'chi'.

## [1] "Optimal chi from Vcmax-gs opt. with light limit.: "

Response to the environment

\(\chi\)

Create data frames where one variable is varied at a time, holding all others constant, and apply rpmodel.

\(A_J:A_C\)

This demonstrates that the coordination (for mean daily condition) emerges from the alternative optimality criterion (instead of being used as a principle).

\(J_\text{max}:V_\text{cmax}\)

This demonstrates that the coordination (for mean daily condition) emerges from the alternative optimality criterion (instead of being used as a principle).

Profit Maximisation vs. Least Cost

Maximise the following term numerically, again subject to \(V_{\mathrm{cmax}}\) and \(g_s\): \[ A_n = A - aE -bV_{\text{cmax}} = max. \;\;\;\;\;\;\;\;\;\;\;\;(4) \] We can divide this by \(a\) in order to define it with the cost ratio \(\beta\) which we have used above. However, an absolute cost remains in the first term: \[ A_n = \frac{1}{a} A - E -\beta V_{\text{cmax}} = max. \;\;\;\;\;\;\;\;\;\;\;\;(5) \]

Response to the environment

xxxxxxxxxx

Light limitation?

But what about the implicit assumption made in the optmisation that assimilation always follows the Rubisco-limited curve? Shouldn’t there be a limit given by light availability? We can spell this out by calculating \(A=min(A_C, A_J)\) or equivalently \(c_i=max(c_i^C, c_i^J)\) and modify calc_optimal_gs_vcmax() accordingly. We now need to additionally solve the equation system \[ A = \phi_0 I_{\text{abs}} \; \frac{c_i-\Gamma^{\ast}}{c_i + 2 \Gamma^{\ast}} \\ A = g_s(c_a - c_i) \]

–>

–>