## 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)
b_quad <- gs * ca - gs * kmm - vcmax
c_quad <- gs * ca * kmm + vcmax * gammastar

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"