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(rsofun)
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
gamma <- 100 #  0.053       # 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 <- rsofun::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 
  )

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$ci )
## [1] 28.14209
print( out_analytical$ca - (out_analytical$gpp / c_molmass) / out_analytical$gs )
## [1] 28.14209
print( out_analytical$ca * out_analytical$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$gpp / c_molmass )
## [1] 10.68456
print( out_analytical$vcmax * (out_analytical$ci - out_analytical$gammastar) / (out_analytical$ci + out_analytical$kmm ))
## [1] 10.68456
print( out_analytical$gs * (out_analytical$ca - out_analytical$ci) )
## [1] 10.68456
print( kphio * ppfd * fapar * (out_analytical$ci - out_analytical$gammastar) / (out_analytical$ci + 2 * out_analytical$gammastar ))
## [1] 10.68456

Yes.

Numerical solution

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 <- 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.01,
    upper     = 0.99,
    fn        = maximise_this,
    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)
}

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$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, beta )
print(chi_opt)
## [1] 0.08229195

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

print(out_analytical$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$vcmax)
## [1] 31.98167
print( kphio * fapar * ppfd * (out_analytical$ci + out_analytical$kmm) / (out_analytical$ci + 2 * out_analytical$gammastar ) )
## [1] 31.98167

Performance benchmark

TODO: Check computation time for N calls of the numerical calc_optimal_chi_num() vs. the analytical calc_optimal_chi() functions.

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 ){

  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]
  
  ## Get assimilation and ci, given Vcmax and gs
  ## by solving the equation system
  ## assim = vcmax * (ci - gammastar)/(ci + kmm)
  ## assim = gs * (ca - ci)
  a_quad <- -1.0 * gs
  b_quad <- gs * ca - gs * kmm - vcmax
  c_quad <- gs * ca * kmm + vcmax * gammastar

  root_ci <- try( polyroot( c(c_quad, b_quad, a_quad) ) )

  if (class(root_ci)=="try-error"){

    return( list( vcmax=NA, gs=NA, ci=NA, chi=NA, a_c=NA, cost_transp=NA, cost_vcmax=NA, net_assim=NA  ) )

  } else {

    ## take only real part of the root 
    ci <- Re(root_ci)

    ## take only positive root
    ci <- ci[which(ci>0)]   # take positive root

    ## if both are positive, take the one that satisfies ci < ca (don't know if this is necessary)
    if (length(ci)>1) ci <- ci[which(ci<ca)]
    
    ## A_c
    a_c <- vcmax * (ci - gammastar) / (ci + kmm)

    ## 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
    
    ## Option B: This is equivalent to the P-model with its optimization of ci:ca. 
    net_assim <- -(cost_transp + cost_vcmax) / a_c
    
    if (maximize) net_assim <- -net_assim
    
    if (return_all){
      return( list( vcmax=vcmax, gs=gs, ci=ci, chi=ci/ca, a_c=a_c, cost_transp=cost_transp, cost_vcmax=cost_vcmax, net_assim=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 ){
  out_optim <- optimr::optimr(
    par        = c( out_analytical$vcmax, out_analytical$gs ), # starting values
    lower      = c( out_analytical$vcmax*0.0001, out_analytical$gs*0.001 ),
    upper      = c( out_analytical$vcmax*20,   out_analytical$gs*30 ),
    fn         = calc_optimal_gs_vcmax,
    args       = c( out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, beta ),
    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 )
  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.69435201913794"
## [1] "Optimal chi from P-model: 0.694352013202358"

Indeed, this is the same.

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 ) )
## [1] "optimal Vcmax: 31.9816667938782"
print( paste( "Vcmax^P:", kphio * fapar * ppfd * (varlist$ci + varlist$kmm) / (varlist$ci + 2 * varlist$gammastar )) )
## [1] "Vcmax^P: "
print( paste("Vcmax from P-model:", out_analytical$vcmax ) )
## [1] "Vcmax from P-model: 31.9816667938778"

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) \]

This is implemented as follows:

calc_optimal_gs_vcmax_ll <- function( par, args, iabs, kphio, maximize=FALSE, return_all=FALSE ){

  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]
  
  ## Get ci using Rubisco-limited assimilation
  ## Get assimilation and ci, given Vcmax and gs
  ## by solving the equation system
  ## assim = vcmax * (ci - gammastar)/(ci + kmm)
  ## assim = gs * (ca - ci)
  a_quad <- -1.0 * gs  # coefficient of quadratic term
  b_quad <- gs * ca - gs * kmm - vcmax
  c_quad <- gs * ca * kmm + vcmax * gammastar

  root_ci_c <- try( polyroot( c(c_quad, b_quad, a_quad) ) )

  ## Get ci using light-limited assimilation
  ## Get assimilation and ci, given Vcmax and gs
  ## by solving the equation system
  ## assim = kphio * iabs * (ci-gammastar)/(ci+2*gammastar)
  ## assim = gs * (ca - ci)
  a_quad <- gs  # coefficient of quadratic term
  b_quad <- kphio * iabs - gs * ca + 2 * gs * gammastar
  c_quad <- - kphio * iabs * gammastar - 2 * gs * ca * gammastar

  root_ci_j <- try( polyroot( c(c_quad, b_quad, a_quad) ) )

  if (class(root_ci_j)=="try-error" || class(root_ci_c)=="try-error"){

    return( list( vcmax=NA, gs=NA, ci=NA, chi=NA, a_c=NA, a_j=NA, ci_c=NA, ci_j=NA, cost_transp=NA, cost_vcmax=NA, net_assim=NA  ) )

  } else {

    ## take only real part of the root 
    ci_j <- Re(root_ci_j)
    ci_c <- Re(root_ci_c)

    ## take only positive root
    ci_j <- ci_j[which(ci_j>0)]   # take positive root
    ci_c <- ci_c[which(ci_c>0)]   # take positive root

    ## if both are positive, take the one that satisfies ci < ca (don't know if this is necessary)
    if (length(ci_j)>1) ci_j <- ci_j[which(ci_j<ca)]
    if (length(ci_c)>1) ci_c <- ci_c[which(ci_c<ca)]

    ## Rubisco-limited
    a_c <- vcmax * (ci_c - gammastar) / (ci_c + kmm)

    ## light-limited
    a_j <- kphio * iabs * (ci_j - gammastar) / (ci_j + 2 * gammastar)

    ## Take minimum of the two assimilation rates
    assim <- min( a_j, a_c )
    
    # if (return_all){
    #   if (a_j<a_c) print("Warning: A_J is lower than A_C")
    # }

    ## ... and consistently the maximum of the two ci
    ci <- max( ci_j, ci_c  )

    ## 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
    
    ## Option B: This is equivalent to the P-model with its optimization of ci:ca. 
    net_assim <- -(cost_transp + cost_vcmax) / assim
    
    if (maximize) net_assim <- -net_assim
    
    if (return_all){
      return( list( vcmax=vcmax, gs=gs, ci=ci, chi=ci/ca, a_c=a_c, a_j=a_j, ci_c=ci_c, ci_j=ci_j, cost_transp=cost_transp, cost_vcmax=cost_vcmax, net_assim=net_assim  ) )
    } else {
      return( net_assim )
    }

  }
}

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

wrap_calc_optimal_gs_vcmax_ll <- function( ppfd, fapar, kphio, kmm, gammastar, ns_star, ca, vpd, beta, vcmax_start, gs_start ){
  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*20, gs_start*30 ),
    fn         = calc_optimal_gs_vcmax_ll,
    args       = c(kmm, gammastar, ns_star, ca, vpd, beta),
    iabs       = (ppfd * fapar),
    kphio      = kphio,
    method     = "L-BFGS-B",
    maximize   = TRUE,
    control    = list( maxit = 1000 )
    )
  varlist <- calc_optimal_gs_vcmax_ll( par=out_optim$par, args=c(kmm, gammastar, ns_star, ca, vpd, beta), iabs=(fapar*ppfd), kphio=kphio, maximize=FALSE, return_all=TRUE )
  return(varlist)
}

varlist_ll <- wrap_calc_optimal_gs_vcmax_ll( ppfd, fapar, kphio, out_analytical$kmm, out_analytical$gammastar, out_analytical$ns_star, out_analytical$ca, vpd, beta, out_analytical$vcmax, out_analytical$gs )

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