Background

The aim here is to use the plant package to investigate dynamics within a patch of competing individuals, focusing on emergent patch-level properties, rather than properties of individuals within the patch. We’ll work with a single species at demographic equilibrium:

library(plant)

params <- scm_base_parameters("FF16")
patch <- expand_parameters(trait_matrix(0.0825, "lma"), params, mutant = FALSE)
patch_eq <- equilibrium_seed_rain(patch)

## [2021-02-24 21:43:13.133] equilibrium> Solving seed rain using iteration
## [2021-02-24 21:43:13.133] equilibrium> Solving seed rain using iteration
## [2021-02-24 21:43:15.316] schedule> 1: Splitting {20} times (141)
## [2021-02-24 21:43:15.316] schedule> 1: Splitting {20} times (141)
## [2021-02-24 21:43:17.673] schedule> 2: Splitting {9} times (161)
## [2021-02-24 21:43:17.673] schedule> 2: Splitting {9} times (161)
## [2021-02-24 21:43:20.185] equilibrium> eq> 1: {1} -> {19.97185} (delta = {18.97185})
## [2021-02-24 21:43:20.185] equilibrium> eq> 1: {1} -> {19.97185} (delta = {18.97185})
## [2021-02-24 21:43:23.918] schedule> 1: Splitting {34} times (141)
## [2021-02-24 21:43:23.918] schedule> 1: Splitting {34} times (141)
## [2021-02-24 21:43:28.684] schedule> 2: Splitting {27} times (175)
## [2021-02-24 21:43:28.684] schedule> 2: Splitting {27} times (175)
## [2021-02-24 21:43:34.365] schedule> 3: Splitting {9} times (202)
## [2021-02-24 21:43:34.365] schedule> 3: Splitting {9} times (202)
## [2021-02-24 21:43:39.920] schedule> 4: Splitting {1} times (211)
## [2021-02-24 21:43:39.920] schedule> 4: Splitting {1} times (211)
## [2021-02-24 21:43:45.452] equilibrium> eq> 2: {19.97185} -> {17.15125} (delta = {-2.820598})
## [2021-02-24 21:43:45.452] equilibrium> eq> 2: {19.97185} -> {17.15125} (delta = {-2.820598})
## [2021-02-24 21:43:51.681] schedule> 1: Splitting {1} times (212)
## [2021-02-24 21:43:51.681] schedule> 1: Splitting {1} times (212)
## [2021-02-24 21:43:57.905] equilibrium> eq> 3: {17.15125} -> {17.31508} (delta = {0.1638286})
## [2021-02-24 21:43:57.905] equilibrium> eq> 3: {17.15125} -> {17.31508} (delta = {0.1638286})
## [2021-02-24 21:44:04.226] equilibrium> eq> 4: {17.31508} -> {17.30373} (delta = {-0.01134973})
## [2021-02-24 21:44:04.226] equilibrium> eq> 4: {17.31508} -> {17.30373} (delta = {-0.01134973})
## [2021-02-24 21:44:04.227] equilibrium> Reached target accuracy (delta 1.13497e-02, 6.55482e-04 < 1.00000e-03 eps)
## [2021-02-24 21:44:04.227] equilibrium> Reached target accuracy (delta 1.13497e-02, 6.55482e-04 < 1.00000e-03 eps)

result <- run_scm_collect(patch_eq)

Cohort densities

There was a bit of hassle in Patches vignette in plotting all trajectories along with a couple of lines corresponding to focal years. We’re going to do the same thing here:

closest <- function(t, time) {
  which.min(abs(time - t))
}
last <- function(x) {
  x[[length(x)]]
}

times <- c(5, 10, 20, 40, last(result$time))
i <- vapply(times, closest, integer(1), result$time)
blues <- c("#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6",
           "#4292C6", "#2171B5", "#08519C", "#08306B")
cols <- colorRampPalette(blues[-(1:2)])(length(i))

height      <- t(result$species[[1]]["height", , ])
log_density <- t(result$species[[1]]["log_density", , ])

As with height in the Patches vignette, cohort density is a matrix of time and cohort identity. However, to show the profiles for a given time slice using matplot it is convenient to transpose this matrix (matplot works column-by-column so we want each column to to be the state of the system at a given time):

density <- exp(log_density)
matplot(height, density, type="l", lty=1,
        col=make_transparent("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y")

Note that the densities here can be extremely low, and yet individuals within the cohorts continue to grow in size; this is the “atto-fox problem”, though here we drop down as low as 8 yocto individuals / m / metre squared (a yocto-x being 1 millionth of an atto-x). Anyway….

The trajectories are easier to understand if a few are highlighted.

xlim <- c(0, max(height, na.rm=TRUE) * 1.05)
matplot(height, density, type="l", lty=1,
        col=make_transparent("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y", xlim=xlim)
matlines(height[, i], density[, i], col=cols, lty=1, type="l")
points(height[1, i], density[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), density[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

Early on there is a low density “dip” caused by self thinning (5 years). As the stand develops that dip broadens and deepens and moves forward in height (these very low density cohorts are still traveling the characteristic equations of the SCM). Later on (20 years) additional an second wave of recruitment gives a second pulse of high density (around 4 m tall), which can be seen traveling along in the 40 year profile to about 13 m. When the patch is very old the stand approaches a stable age distribution, though a very narrow window of “missing” heights exists just below the top of the canopy.

Leaf area

It’s also possible see where the leaf area in the patch is coming from; a profile of leaf area with respect to height. Again, this requires reconstructing the patches, and using an unexported function from plant to put this into a matrix:

patches <- lapply(seq_along(result$time), scm_patch, result)
leaf_area <- lapply(patches, function(x) x$species[[1]]$competition_effects)
leaf_area <- plant:::pad_list_to_array(leaf_area)

matplot(height, leaf_area, type="l", lty=1, col="lightgrey",
        xlim=xlim, xlab="Height (m)",
        ylab="Leaf area density (m2 / m2 / m)", las=1)
matlines(height[, i], leaf_area[, i], col=cols, lty=1, type="l")
points(height[1, i], leaf_area[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), leaf_area[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

# Growth rate

Finally, we can see where height growth rate is concentrated in the population. This differs from the profile in the Individuals vignette because it is reduced depending on the light environment, but that environment is the one constructed by the individual itself.

ode_size <- patches[[1]]$ode_size

growth_rate <- lapply(patches, function(x)
                      matrix(x$species[[1]]$ode_rates, ode_size)[1, ])
growth_rate <- plant:::pad_list_to_array(growth_rate)

matplot(height, growth_rate, type="l", lty=1, col="lightgrey",
        xlim=xlim, xlab="Height (m)",
        ylab="Height growth rate (m / year)", las=1)
matlines(height[, i], growth_rate[, i], col=cols, lty=1, type="l")
points(height[1, i], growth_rate[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), growth_rate[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))

# Integrating over a metapopulation

The above plots show relationships with patches of a given age What about the average relationship across the entire metapopulation? To get that, we average (integrate) over the distribution over patch (for formula see demography vignette). To achieve this we first need the patch-level relationships to a series of fixed heights

hh <- seq_log_range(range(height, na.rm=TRUE), 500)

We’ll use a spline interpolation, on log-log-scaled data, clamped so that for x values outside the observed range are set to zero.

f <- function(height, density, hout) {
  r <- range(height, na.rm=TRUE)
  clamp_domain(splinefun_loglog(height, density), r, 0)(hout)
}

Now interpolate the height-density relationship in each patch to the series of specified heights

xx <- lapply(seq_along(result$time),
             function(i) f(height[, i], density[, i], hh))
n_hh <- plant:::pad_list_to_array(xx)

For each of these heights, we can now integrate across patches (i.e. across rows), weighting by patch abundance

trapezium <- plant:::trapezium
n_av <- apply(n_hh, 1,
              function(x) trapezium(result$time, x * result$patch_density))

Add this average to the plot (red line):

xlim <- c(0, max(height, na.rm=TRUE) * 1.05)
matplot(height, density, type="l", lty=1,
        col=make_transparent("black", 0.15),
        xlab="Height (m)", ylab="Density (1 / m / m2)", las=1,
        log="y", xlim=xlim)
matlines(height[, i], density[, i], col=cols, lty=1, type="l")
points(height[1, i], density[1, i], pch=19, col=cols)
text(height[1, i] + strwidth("x"), density[1, i],
     paste0(round(times), c(" years", rep("", length(times) - 1))),
     adj=c(0, 0))
points(hh, n_av, col="red", type='l')

Patch-level emergent properties

Rich FitzJohn

Daniel Falster

2016

Background

Cohort densities

Leaf area