authors: DS Falster, A O’Reilly-Nugent, I Towers, F Robinson date: 2022

Background

This page contains a brief introduction to get you started with plant.

At a high level, plant has two key use cases - modelling individuals or modelling patches. Modelling either is dependent on the species parameters, the environment and the intial conditions. By modifying these inputs you can grow different individuals or patches of species.

XXXX diagram here

More extensive information and tutorials are available, such as:

Details of the modelling approaches (concepts):

Modelling demography of individuals, patches and metapopulations
The FF16 physiological strategy model

Details of using plant in R:

Setup

First, load plant, tidyverse and purrr, which will be useful (for working with list objects):

library(plant)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)
library(purrr)

For a full list of available functions run library(help=plant), or see the online reference.

Strategy objects

Strategies are the corner stone of plant, describing the system of dynamical equations that determine an individuals vital rates – i.e. growth, reproduction and mortality. The plant package is setup to accept different stratgies. Think of these as sets of rules which determine how biological foundation of your model.

To get started, we’ll use the FF16 strategy.

You can feed these strategies into either individuals or patches to alter the solver outcomes.

For definitions of each of these parameters see Concepts: FF16 Strategy tables.

For further information on how to adjust the strategies see Modifying parameters of strategies

s <- FF16_Strategy()
str(s)

## List of 37
##  $ lma                   : num 0.198
##  $ rho                   : num 608
##  $ hmat                  : num 16.6
##  $ omega                 : num 3.8e-05
##  $ eta                   : num 12
##  $ theta                 : num 0.000214
##  $ a_l1                  : num 5.44
##  $ a_l2                  : num 0.306
##  $ a_r1                  : num 0.07
##  $ a_b1                  : num 0.17
##  $ r_s                   : num 6.6
##  $ r_b                   : num 13.2
##  $ r_r                   : num 217
##  $ r_l                   : num 198
##  $ a_y                   : num 0.7
##  $ a_bio                 : num 0.0245
##  $ k_l                   : num 0.457
##  $ k_b                   : num 0.2
##  $ k_s                   : num 0.2
##  $ k_r                   : num 1
##  $ a_p1                  : num 151
##  $ a_p2                  : num 0.205
##  $ a_f3                  : num 0.000114
##  $ a_f1                  : num 1
##  $ a_f2                  : num 50
##  $ S_D                   : num 0.25
##  $ a_d0                  : num 0.1
##  $ d_I                   : num 0.01
##  $ a_dG1                 : num 5.5
##  $ a_dG2                 : num 20
##  $ k_I                   : num 0.5
##  $ recruitment_decay     : num 0
##  $ control               :List of 18
##   ..$ function_integration_rule      : num 21
##   ..$ offspring_production_tol       : num 1e-08
##   ..$ offspring_production_iterations: num 1000
##   ..$ node_gradient_eps              : num 1e-06
##   ..$ node_gradient_direction        : int 1
##   ..$ node_gradient_richardson       : logi FALSE
##   ..$ node_gradient_richardson_depth : num 4
##   ..$ ode_step_size_initial          : num 1e-06
##   ..$ ode_step_size_min              : num 1e-06
##   ..$ ode_step_size_max              : num 0.1
##   ..$ ode_tol_rel                    : num 1e-06
##   ..$ ode_tol_abs                    : num 1e-06
##   ..$ ode_a_y                        : num 1
##   ..$ ode_a_dydt                     : num 0
##   ..$ schedule_nsteps                : num 20
##   ..$ schedule_eps                   : num 0.001
##   ..$ schedule_verbose               : logi FALSE
##   ..$ save_RK45_cache                : logi FALSE
##   ..- attr(*, "class")= chr "Control"
##  $ collect_all_auxiliary : logi FALSE
##  $ birth_rate_x          : num(0) 
##  $ birth_rate_y          : num 1
##  $ is_variable_birth_rate: logi FALSE
##  - attr(*, "class")= chr "FF16_Strategy"

Individuals

We can define one or more individuals for a given strategy.

These are accessed using the Individual class:

ind <- FF16_Individual()

Noting that FF16_Individual uses the FF16_Strategy by default (see ?FF16_Individual for more information).

Our individual also shares the FF16_Environment (more on that soon) and a number of rates and functions for us to explore. The Individuals vignette describes the nuts and bolts of all these functions, for now we’re only going to grow and plot our individual’s height.

str(ind)

## Classes 'Individual<FF16,FF16_Env>', 'Individual', 'R6' <Individual<FF16,FF16_Env>>
##   Inherits from: <Individual>
##   Public:
##     .ptr: externalptr
##     aux: function (name) 
##     aux_names: active binding
##     aux_size: active binding
##     clone: function (deep = FALSE) 
##     compute_competition: function (h) 
##     compute_rates: function (environment) 
##     establishment_probability: function (environment) 
##     initialize: function (ptr) 
##     internals: active binding
##     mortality_probability: active binding
##     net_mass_production_dt: function (environment) 
##     ode_names: active binding
##     ode_rates: active binding
##     ode_size: active binding
##     ode_state: active binding
##     rate: function (name) 
##     reset_mortality: function () 
##     resource_compensation_point: function () 
##     set_state: function (name, v) 
##     state: function (name) 
##     strategy: active binding
##     strategy_name: active binding

First we set a fixed environment (here 1.0 represents full light exposure in an open canopy) then use the grow_individual_to_time function to grow our individual for a range of time steps.

env <- FF16_fixed_environment(1.0)
times <- seq(0, 50, length.out = 101)
result <- grow_individual_to_time(ind, times, env) %>%
  tidy_individual()

The tidy_indivdual function puts the outputs from grow_individual_to_time into tidy format, as defined in the tidyverse. Examining our result, we see a matrix of our state variables at each timestep

result %>%
  head()

## # A tibble: 6 × 7
##    step  time height mortality fecundity area_heartwood mass_heartwood
##   <int> <dbl>  <dbl>     <dbl>     <dbl>          <dbl>          <dbl>
## 1     1   0    0.344   0        0         0                 0         
## 2     2   0.5  0.545   0.00500  7.69e-22  0.00000000611     0.00000153
## 3     3   1    0.819   0.0100   7.23e-21  0.0000000308      0.0000109 
## 4     4   1.5  1.16    0.0150   6.35e-20  0.000000115       0.0000569 
## 5     5   2    1.57    0.0200   5.66e-19  0.000000354       0.000237  
## 6     6   2.5  2.03    0.0250   5.10e-18  0.000000942       0.000817

Which we can plot against time

#need to create a function to convert individual results to tidy-format
#just do ggplot for now
result %>%
  ggplot() +
  geom_line(aes(x=time, y=height), linewidth =1.5) +
  theme_classic()

You can find out more about how to grow individuals under the Individual plants example OR if you’re wanting to learn more about the concept and equations underlying individuals you can find them under Concepts: demography of plants, patches and metacommunities

Patches

While we can torture our individual into all sorts of shapes and sizes, it’s far more interesting to see how many individuals interact.

The most important feature of plant is enabling easy modeling of size-distributions for multiple competing species within a patch.

A patch is simply an area in which a collection of individuals grow and compete for resources (see Concepts: patches). By size distribution we mean the distribution the relationship for density of individuals against their size (see Concepts-> XXX). By default, plant starts with a bare/empty patch, but you can also start with an existing size distribution (see XXXX).

(XXX figure: Example size distribution (stylised)

Modelling a patch involves

introducing individuals into the patch over time for each species, together which form a size-distribution for that species, and
stepping the size distributions for one or more species through time, integrating the rates of change on individuals to describe their state (such as height) and integrating the impact their state has on other individuals in a patch via the environment.

Plant uses a particular mathematical method to model size distributions, which assumes large populations and continuous size distribution. The dynamics are solved by integrating along so-called “characteristics” of a partial differential equation (see Concepts-> XXX).

The output from modelling a patch is size distribution at a series of times during patch development. At each time, the model reports states and rates for example individuals along the characteristics. From this, you can calculate a wide variety of metrics for the individual, population, community or environment within a patch at a given time.

Below we provide some examples to get you started. More worked examples are available at XXXXX.

Single species

First, let’s activate the logger (to keep us updated on what’s happening):

plant_log_console()

## [2025-03-13 07:23:07.074] Activating logging to console

Then start by specifying the model rules and load some basic controls. Here we’ll use the FF16 ruleset:

params <- scm_base_parameters("FF16")

We can then need to specify the particular species by specifying their traits and birth rates. Here we use lma but other traits are possible (see Examples -> XXX).

species_parameters <- expand_parameters(trait_matrix(0.0825, "lma"), params)

Then we run plant solver (SCM) to solve the size distribution over time. There are two steps

# 1. Resolve number of characteristics needed for specified accuracy
species_parameters <- build_schedule(species_parameters)

## [2025-03-13 07:23:07.846] schedule> 1: Splitting {20} times (141)
## [2025-03-13 07:23:08.594] schedule> 2: Splitting {9} times (161)

# 2. Rerun, gathering outputs at each time step
result <- 
  run_scm_collect(species_parameters) %>% 
  tidy_patch()

If you’re working through this yourself, the SCM solver should be blazingly fast. This is quite cool, given that this model has hundreds of coupled differential equations veing stepped through time, wtih complex compeititve interactions at each step.

The results is a list with a number of objects, which are explained in full at XXXX, and include

the times the size-distribution is reported (the exact times used are determined by the solver, via adaptive stepping of system. See XXXX),
details on the size distribution for each species in the patch at each timestep,
details on the environment in the patch at each at each timestep.

str(result, max.level = 1)

## List of 6
##  $ time                : num [1:171] 0e+00 1e-05 2e-05 3e-05 4e-05 ...
##  $ species             : tibble [29,070 × 13] (S3: tbl_df/tbl/data.frame)
##  $ env                 :List of 1
##  $ offspring_production: num 20
##  $ p                   :List of 9
##   ..- attr(*, "class")= chr [1:2] "Parameters<FF16,FF16_Env>" "Parameters"
##   ..- attr(*, "offspring_production")= num 20
##  $ n_spp               : int 1

For easy usage, the tidy_patch function puts the species and env outputs into tidy format, as defined in the tidyverse. The species object looks like this:

result %>% 
  pluck("species") %>%
  drop_na()

## # A tibble: 14,705 × 13
##     step    time patch_density species  node height mortality fecundity
##    <int>   <dbl>         <dbl> <chr>   <int>  <dbl>     <dbl>     <dbl>
##  1     1 0              0.0333 1           1  0.392    0.0107  0       
##  2     2 0.00001        0.0333 1           1  0.392    0.0107  7.36e-27
##  3     2 0.00001        0.0333 1           2  0.392    0.0107  0       
##  4     3 0.00002        0.0333 1           1  0.392    0.0107  1.47e-26
##  5     3 0.00002        0.0333 1           2  0.392    0.0107  7.36e-27
##  6     3 0.00002        0.0333 1           3  0.392    0.0107  0       
##  7     4 0.00003        0.0333 1           1  0.392    0.0107  2.21e-26
##  8     4 0.00003        0.0333 1           2  0.392    0.0107  1.47e-26
##  9     4 0.00003        0.0333 1           3  0.392    0.0107  7.36e-27
## 10     4 0.00003        0.0333 1           4  0.392    0.0107  0       
## # ℹ 14,695 more rows
## # ℹ 5 more variables: area_heartwood <dbl>, mass_heartwood <dbl>,
## #   offspring_produced_survival_weighted <dbl>, log_density <dbl>,
## #   density <dbl>

Let’s look at height. Each line represents the height of an example individual growing along a characteristic curve over time, beginning from the point at which the example individual recruited into the patch. Notably, the first example individuals follow much the same growth curve as our individual above, but subsequent individuals have a bumpier ride, with growth slowing as the canopy closes over. Another variable in the species output is the log density of the example individual, which we can use to weight the growth trajectories to show thinning. For example, by following the growth of an example individual we can see that the relative density of these individuals is initially high after recruitment into the patch and declines progressively through time.

result%>%
  pluck("species") %>%
  drop_na() %>%
  plot_size_distribution()

To see the change in canopy openness over time we can explore the patch Environment. Lets look at year 20 first, which corresponds to the 99th timestep in our model (not all timesteps are equal!) Our FF16 environment is described in terms of canopy openness, with 1.0 being completely open and 0.0 being completely shaded. We see that the shortest individuals experience intense shading while taller individuals enjoy full sunlight:

result %>%
  pluck("env", "light_availability") %>%
  filter(step == "99") %>%
  ggplot() +
  geom_line(aes(x=height, y = light_availability), linewidth =1.5) +
  theme_classic()

If we look at the light environment at the forest floor (height = 0.0) we can see that it varies through time as older individuals thin out and gaps form

result %>%
  pluck("env", "light_availability") %>%
  filter(height == 0) %>%
  ggplot() +
  geom_line(aes(x= time, y=light_availability), linewidth=1.5) +
  theme_classic()

We can now look at the size distributions of two species competing in a patch. We can use the trait matrix function to specify multiple values for a given trait, each representing a different species in the community, or, as below, we can build upon the one-patch example above by using the expand_parameters function, including species_parameters as an argument.

species_parameters_2sp <- expand_parameters(trait_matrix(0.2625, "lma"), species_parameters)

As above in the single-species patch, we then determine the appropriate number of characteristics and rerun the solver, collecting the output from the patch-level dynamics.

species_parameters_2sp <- build_schedule(species_parameters_2sp)

## [2025-03-13 07:23:12.976] schedule> 1: Splitting {1,7} times (170,170)

result <- 
  run_scm_collect(species_parameters_2sp) %>% 
  tidy_patch()

We can now plot the size distributions of these two species through time

result$species %>%
  drop_na() %>%
plot_size_distribution()

Meta-populations

🚧 Coming soon 🚧