Strategies & traits

guide

How to inspect and modify the parameters that define a plant strategy, and how trait values propagate through the hyper-parameterisation. For what the strategies mean physiologically, see the model docs.

library(plant)

There are a large number of parameters to the physiological model, but all are changeable. The default strategy is detailed in the model docs:

s <- FF16_Strategy()
names(s)

[1] "pars"                   "control"                "collect_all_auxiliary" 
[4] "birth_rate_x"           "birth_rate_y"           "is_variable_birth_rate"

The strategy object here is a special object of class FF16_Strategy. In contrast with most of the “reference” objects used by plant, this is simply a list with a class attribute. Some validation will be done on the parameters every time it is passed into the C++ bits of the model.

The physiological parameters now live in the nested pars element (the top-level strategy also carries control, birth-rate settings, and a few flags). All of the pars values are floating point (decimal) numbers:

unlist(s$pars)

              lma               rho              hmat             omega 
     1.978791e-01      6.080000e+02      1.659587e+01      3.800000e-05 
              eta             theta              a_l1              a_l2 
     1.200000e+01      2.141786e-04      5.440000e+00      3.060000e-01 
             a_r1              a_b1               r_s               r_b 
     7.000000e-02      1.700000e-01      6.598684e+00      1.319737e+01 
              r_r               r_l               a_y             a_bio 
     2.170000e+02      1.984545e+02      7.000000e-01      2.450000e-02 
              k_l               k_b               k_s               k_r 
     4.565855e-01      2.000000e-01      2.000000e-01      1.000000e+00 
             a_p1              a_p2              a_f3              a_f1 
     1.511778e+02      2.047162e-01      1.140000e-04      1.000000e+00 
             a_f2               S_D              a_d0               d_I 
     5.000000e+01      2.500000e-01      1.000000e-01      1.000000e-02 
            a_dG1             a_dG2               k_I recruitment_decay 
     5.500000e+00      2.000000e+01      5.000000e-01      0.000000e+00 

All of these parameters can be directly changed. For example, we can double the leaf mass per unit leaf area (lma) value:

s$pars$lma <- s$pars$lma * 2

and then from this construct a plant:

pl <- FF16_Individual(s)
pl$strategy$pars$lma

[1] 0.3957582

lma affects a few places in the model; see the source code, but only really for converting from leaf area to leaf mass (leaf mass being leaf area multiplied by leaf mass per unit leaf area). However, as a component of the leaf economic spectrum we imagine LMA as affecting a number of other components of the model — see FF16 for the trait trade-offs this functional-balance model builds in.

We capture this through what we call a “hyper-parameterisation”; additional parameters and functions that mean that changing one parameter might affect a number of other lower-level parameters. Our default hyper-parameterisation is via the FF16_hyperpar function:

FF16_hyperpar

function(m, s, filter=TRUE) {
    with_default <- function(name, default_value=s$pars[[name]]) {
      rep_len(if (name %in% colnames(m)) m[, name] else default_value,
              nrow(m))
    }
    lma       <- with_default("lma")
    rho       <- with_default("rho")
    omega     <- with_default("omega")
    narea     <- with_default("narea", narea)

    ## Light extinction coefficient: source from the strategy so the
    ## assimilation integral below stays consistent with the canopy
    ## model, rather than carrying a separate hard-coded default here.
    k_I       <- s$pars$k_I

    ## lma / leaf turnover relationship:
    k_l   <- B_kl1 * (lma / lma_0) ^ (-B_kl2)

    ## rho / mortality relationship:
    d_I  <- B_dI1 * (rho / rho_0) ^ (-B_dI2)

    ## rho / wood turnover relationship:
    k_s  <- B_ks1 *  (rho / rho_0) ^ (-B_ks2)

    ## rho / sapwood respiration relationship:

    ## Respiration rates are per unit mass, so this next line has the
    ## effect of holding constant the respiration rate per unit volume.
    ## So respiration rates per unit mass vary with rho, respiration
    ## rates per unit volume don't.
    r_s <- B_rs1 / rho
    # bark respiration follows from sapwood
    r_b <- B_rb1 / rho

    ## omega / accessory cost relationship
    a_f3 <- B_f1 * omega

    ## Narea, photosynthesis, respiration

    assimilation_rectangular_hyperbolae <- function(I, Amax, theta, QY) {
      x <- QY * I + Amax
      (x - sqrt(x^2 - 4 * theta * QY * I * Amax)) / (2 * theta)
    }

    ## Photosynthesis  [mol CO2 / m2 / yr]
    approximate_annual_assimilation <- function(narea, latitude) {
      E <- seq(0, 1, by=0.02)
      ## Only integrate over half year, as solar path is symmetrical
      D <- seq(0, 365/2, length.out = 10000)
      I <- PAR_given_solar_angle(solar_angle(D, latitude = abs(latitude)))

      Amax <- B_lf1 * (narea/narea_0) ^  B_lf5
      theta <- B_lf2
      QY <- B_lf3

      AA <- NA * E

      for (i in seq_len(length(E))) {
        AA[i] <- 2 * trapezium(D, assimilation_rectangular_hyperbolae(
                                    k_I * I * E[i], Amax, theta, QY))
      }
      if(all(diff(AA) < 1E-8)) {
        # line fitting will fail if all have are zero, or potentially same value
        ret <- c(dplyr::last(AA), 0)
        names(ret) <- c("p1","p2")
      } else {
        fit <- nls(AA ~ p1 * E/(p2 + E), data.frame(E = E, AA = AA), start = list(p1 = 100, p2 = 0.2))
        ret <- coef(fit)
      }
      ret
    }

    # Guard against narea having length zero (e.g. an empty trait matrix),
    # in which case trapezium() fails.
    a_p1 <- a_p2 <- 0 * narea
    if (length(narea) > 0) {
      i <- match(narea, unique(narea))
      y <- vapply(unique(narea), approximate_annual_assimilation,
                  numeric(2), latitude)
      a_p1  <- y["p1", i]
      a_p2  <- y["p2", i]
    }

    ## Respiration rates are per unit mass, so convert to mass-based
    ## rate by dividing with lma
    ## So respiration rates per unit mass vary with lma, while
    ## respiration rates per unit area don't.
    r_l  <- B_lf4 * narea / lma

    extra <- cbind(k_l,                # lma
                   d_I, k_s, r_s, r_b, # rho
                   a_f3,               # omega
                   a_p1, a_p2,         # narea
                   r_l)                # lma, narea

    overlap <- intersect(colnames(m), colnames(extra))
    if (length(overlap) > 0L) {
      stop("Attempt to overwrite generated parameters: ",
           paste(overlap, collapse=", "))
    }

    ## Check for infitinte values - these cause issues
    if(any(is.infinite(extra))) {
      stop("Attempt to use infinite value in derived parameters: ",
           paste(colnames(extra)[is.infinite(extra)], collapse=", "))
    }

    ## Filter extra so that any column where all numbers are with eps
    ## of the default strategy are not replaced:
    if (filter) {
      if (nrow(extra) == 0L) {
        extra <- NULL
      } else {
        pos <- diff(apply(extra, 2, range)) == 0
        if (any(pos)) {
          eps <- sqrt(.Machine$double.eps)
          x1 <- extra[1, pos]
          x2 <- unlist(s$pars[names(x1)])
          drop <- abs(x1 - x2) < eps & abs(1 - x1/x2) < eps
          if (any(drop)) {
            keep <- setdiff(colnames(extra), names(drop)[drop])
            extra <- extra[, keep, drop=FALSE]
          }
        }
      }
    }

    if (!is.null(extra)) {
      m <- cbind(m, extra)
    }
    m
  }
<bytecode: 0x1185a99a8>
<environment: 0x11859ffa8>

You will rarely need to call this function directly (see below) but note how setting of lma affects parameters k_l, a_p1 and r_l

FF16_hyperpar(trait_matrix(0.1, "lma"), s)

   lma      k_l     a_p1   r_l
p1 0.1 1.466783 151.1778 392.7

These are:

• k_l: Turnover rate for leaves
• a_p1: Leaf photosynthesis per area
• r_l: Leaf respiration per mass

This means that it is possible to implement different trade-offs between parameters relatively easily by modifying the hyper-parameterisation function in R, rather than having to modify the underlying physiological model in C++.

The trait_matrix function is a simple wrapper that just makes sure the trait matrix has the right format:

trait_matrix(0.1, "lma")

     lma
[1,] 0.1

trait_matrix(c(0.1, 500), "rho")

       rho
[1,]   0.1
[2,] 500.0

To make use of the hyper-parameterisation, the preferred way of setting parameters is through the utility functions generate_strategy and add_strategies. These take a Parameters object:

p <- FF16_Parameters()

The Parameters object mostly contains information about the patch:

names(p)

 [1] "patch_area"                  "n_patches"                  
 [3] "patch_type"                  "max_patch_lifetime"         
 [5] "strategies"                  "strategy_default"           
 [7] "node_schedule_times_default" "node_schedule_times"        
 [9] "ode_times"                   "initial_state"              
[11] "n_initial_cohorts"           "initial_node_times"         
[13] "initial_patch_density"       "initial_pr_patch_survival"  
[15] "initial_time"               

The default hyper-parameterisation function associated with a Parameters object is retrieved with param_hyperpar():

identical(param_hyperpar(p), FF16_hyperpar)

[1] TRUE

max_patch_lifetime sets the patch disturbance interval (in years). The physiological parameters (including k_I, the light extinction coefficient) now live on the strategy under pars; the strategy_default member is the Strategy that all others are built from:

class(p$strategy_default)

[1] "FF16_Strategy"

This is the Strategy object that all others will be built from by difference. Running

s <- generate_strategy(p, trait_matrix(0.1, "lma"), birth_rate = 1)[[1]]

will create a strategy s where lma is set but also all the parameters that depend on lma.

The function generate_strategy can be used to create a list of strategies:

lma <- trait_matrix(seq(0.1, 0.5, length.out=5), "lma")
FF16_hyperpar(trait_matrix(lma, "lma"), s)

     lma        k_l     r_l
[1,] 0.1 1.46678341 392.700
[2,] 0.2 0.44833712 196.350
[3,] 0.3 0.22412414 130.900
[4,] 0.4 0.13703875  98.175
[5,] 0.5 0.09356799  78.540

ss <- generate_strategy(p, lma, birth_rate = rep(1, 5))
length(ss)

[1] 5

We can then use standard R commands to extract variable from this list

sapply(ss, function(x) x$pars$lma)

[1] 0.1 0.2 0.3 0.4 0.5

sapply(ss, function(x) x$pars$k_l)

[1] 1.46678341 0.44833712 0.22412414 0.13703875 0.09356799

sapply(ss, function(x) x$pars$a_p1)

[1] 151.1778 151.1778 151.1778 151.1778 151.1778

sapply(ss, function(x) x$pars$r_l)

[1] 392.700 196.350 130.900  98.175  78.540

In addition to the physiological parameters there are a large number of “control” parameters that affect the behaviour of the various numerical algorithms used (note that in contrast to the physiological parameters these have a variety of types). These live in a separate Control object rather than on the Parameters:

Control()

$function_integration_rule
[1] 21

$shading_model
[1] ""

$ppa_layer_optical_depth
[1] 0.5

$ppa_layer_smoothing
[1] 0.3

$offspring_production_tol
[1] 1e-08

$offspring_production_iterations
[1] 1000

$node_gradient_eps
[1] 1e-06

$node_gradient_direction
[1] -1

$node_gradient_richardson
[1] FALSE

$node_gradient_richardson_depth
[1] 4

$ode_step_size_initial
[1] 1e-06

$ode_step_size_min
[1] 1e-06

$ode_step_size_max
[1] 5

$ode_tol_rel
[1] 1e-04

$ode_tol_abs
[1] 1e-04

$ode_a_y
[1] 1

$ode_a_dydt
[1] 0

$fixed_time_step
[1] 0

$schedule_nsteps
[1] 20

$schedule_eps
[1] 0.02

$schedule_verbose
[1] FALSE

$save_RK45_cache
[1] FALSE

$GSS_tol_abs
[1] 0.001

$vulnerability_curve_ncontrol
[1] 100

$ci_abs_tol
[1] 0.001

$ci_niter
[1] 1000

attr(,"class")
[1] "Control"

The Control() defaults are the pragmatic, fast-ish settings used for essentially all of plant’s runs (control() is a lowercase alias for the same constructor). When you need a high-accuracy run, control_accurate() tightens the ODE and schedule tolerances (see ODE stepping & control and The node-spacing algorithm). We can see which fields it changes:

ctrl_accurate <- control_accurate()
ctrl_accurate[unlist(control()) != unlist(ctrl_accurate)]

$ode_step_size_max
[1] 0.1

$ode_tol_rel
[1] 1e-06

$ode_tol_abs
[1] 1e-06

$schedule_eps
[1] 0.001