Modelling Tradescantia fluminensis to assess long term survival

Field experiments

Two experimental sites were established on the North and South Islands of New Zealand. The South Island site, “Conway Flat” (43°S, 173°E, 50 m altitude), was established in late February 2007 and located in a 10.5 ha parcel of remnant podocarp-hardwood forest along a river flat, which had undergone a large flood in March 1988. To exclude large herbivores, mainly deer, temporary electric fences were erected at the beginning of the study. However, several months after commencement of the study, feral pigs entered the plots so rabbit fencing (1 m in height, 50 mm mesh size) was fixed around the bottom of the electric fence. The North Island site was established in June 2007 and was located in a 2 ha area of remnant podocarp-hardwood forest on a river floodplain southwest of Hamilton at Te Pahu (37.5°S, 175°E, 40 m altitude). Herbivory was not a problem at this site so fencing was not required.

At each site, five blocks (replicates) were selected in areas having similar cover of T. fluminensis with each block having regions of both high and low light levels (canopy edge and under canopy respectively). One plot per block was established in each of the high and low light areas giving 10 plots (5 high light and 5 low light) at each site. At Conway Flat, 12 T. fluminensis plants per plot were randomly selected by searching for the base ends of the plant. Only plants that had a minimum of five green leaves on five nodes plus a base end with roots were chosen; these were defined as “complete plants.” The total number of complete plants at this site was 120 (12 plants × 2 light levels × 5 replicates). At Hamilton, the mats of T. fluminensis were denser than those at Conway Flat, making it difficult to locate the base of a plant without damaging it and causing major disturbance to the surrounding sward. To avoid this, a combination of complete plants and partial plants was selected. Complete plants were selected as for the Conway site, except that only 5 plants were chosen from both high and low light plots in only two blocks, giving a total of 20 complete plants. In addition, 120 aerial tips, i.e., tips visible at the top of the plant mat, of plants were chosen (defined as “partial plants”): 5 replicates of 12 plants per plot at both light levels. To distinguish chosen plants from others in the sward, a jewellery tag was placed several nodes forward from the basal end (tag 1) on complete plants and on the internode immediately below the oldest leaf-bearing node (tag 2) on both complete and partial plants.

Data collection at both sites occurred once every 3 months for a 12 month period (late February 2007 to early March 2008 for Conway and June 2007 to June 2008 for Hamilton). For complete plants, during the first three censuses, a drawing was made for each plant noting which nodes had leaves present and the occurrence of any branching or flowering. The length between each successive node (i.e., the internode length) from the basal end to tag 1 was measured, as were the internode lengths nodes from tag 2 forward to the tip of the plant. Internode lengths at the tip were measured to the last measurable node (defined as ≥0.5 cm in length). Separate counts were made for mature (those expanded to ≥3 cm) and immature (<3 cm) leaves. When branching occurred, a wire tag was placed at the node where it joined the main stem and classified as either a green shoot (<1 cm long), a short branch (1–10 cm long), or a green branch (≥10 cm). On the fourth and fifth censuses, full measurements were taken only for two plants per plot and the main stem tips of the remaining 10 plants were tagged only. Where necessary, measurements were taken of any loss of stem length at bases and bases retagged. After the final (i.e., fifth) census, plants were harvested one block at a time. Individual plants were carefully pulled from the ground and placed in large numbered plastic bags for measurement at the lab. Before measurement, the roots of each plant were washed and surface water removed with a paper towel. All internode lengths were measured and the root position for the main stem, the existence of flower buds and the positions of branches were noted. The number of nodes was counted for individual branches. Leaf area, fresh weight and dry weight were recorded for the main stem section of each plant. The combined leaf area, fresh weight and dry weight of all branches were also recorded. A LiCor leaf area meter was used to measure leaf areas. Plants were split up into smaller lots when necessary and dried in brown paper bags at 65–70 °C for at least 4 days to measure dry weights.

Partial plants were measured identically to the complete plants except the measurements were only for the leaf-bearing stem forward to the tip. Because plants were growing quickly, at subsequent visits it became necessary to place extra tags on the internode adjacent to the last measurable tip. Notes were made on the plant drawing of the number of internodes between these tags to avoid problems should the tags slip (e.g. through the death of a leaf).

To summarise, data were collected from four treatments: the Conway Flat site under high and low light and the Hamilton site under high and low light. In each environment, 60 plants were measured at five time points over the course of one year period giving 300 plant measurements at each of the four treatments.

Model parameter estimates

From the data collected, we calculated the following quantities for each of the four environments.

1. Tip growth rate (λ_G): the mean rate of adding new nodes for a single stem (nodes stem⁻¹ year⁻¹)

   This is the number of nodes that was added to the front of the main stem between the first and final censuses, averaged across the 60 plants measured.

2. Death rate (λ_D): the mean death rate of nodes at the base of a stem (nodes base⁻¹ year⁻¹).

   This is the number of nodes which were lost from the rear of the main stem between the first and final censuses.

3. Branching age (A): the mean age of a node when it branches (years).

   For each node that branched during the study period, the time between the node first appearing and the subsequent side stem appearing was calculated. As censuses were taken at three-monthly intervals, this is an approximate value. Nodes that were already present at the initial census and branched subsequently were excluded as the ages of these parent nodes are unknown.

4. Branching probability (P): the probability of a node branching.

   The total number of main stem nodes across all plants that produced a side stem was divided by the total number of main stem nodes. Nodes and side stems that were already present at the initial census were excluded as the age of the side stem was unknown. The analysis of branching ages described above revealed that a negligible number of nodes branched at an age greater than 6 months. Nodes that appeared less than 6 months prior to the final census were, therefore, also excluded as they may have branched subsequently to the final census. Hence, only nodes that appeared between the initial and the third census were used in this calculation. All side stems that originated from one of these nodes at any census were included.

5. Total growth rate (λ_T): the net rate of increase in the number of nodes per plant (nodes year⁻¹).

   This is the net change in the total number of nodes between the initial and final census, averaged across all plants in each treatment. The net change consists of the sum of nodes that were added to the front of the main stem and side stems minus nodes lost from the rear of the plant.

6. Branching rate (B): the rate of increase in the number of stems per plant (stems year⁻¹).

   This was the average number of new stems (including sub-branches) that were added to the plant during the year.

Each of the continuous variables was tested for differences between sites using a two-tailed t-test. The probability variable was tested for differences using a two-sided binomial exact test. All variables were then tested for differences between high and low light environments. For variables with significant between-site differences, this was done separately for each site; for variables with no significant between-site differences, the data were aggregated across the two sites.

Field experiment

Table 1 gives the expected value, standard deviation and 95% confidence intervals (found by boot-strapping for non-normal data) for each variable for the four different treatments. Note that P, the probability of a node branching is calculated from the data as a whole, rather than as a mean of each plant and therefore does not have a standard deviation. The confidence interval in this case is found from the binomial distribution. Of the six variables examined, five were significantly different between sites (λ_G, p < 10⁻⁴; λ_D, p < 0.05; A, p < 0.01; λ_T, p < 10⁻⁴). The branching rate, B, showed no significant difference (p = 0.24) and the branching probability, P, showed a close to significant difference (p < 0.1) between the two sites.

Data from the Hamilton site showed no significant difference between the high and low light environments for four variables (λ_G, p = 0.67; λ_D, p = 0.38; A, p = 0.18; λ_T, p = 0.59), but there was a small difference in the branching probability (P, p < 0.1). At the Conway site, there was a significant difference between high and low light in the total growth rate (λ_T, p < 0.005) and the death rate (λ_D, p < 0.001) but no significant difference in the tip growth rate (λ_G, p < 0.1), the branching age (A, p > 0.1) or the branching probability (P, p > 0.1).

Figure 1(A, D, G, J) shows the time series for the total number of nodes on each plant (grey curves) and averaged across all plants in that treatment (black curve). Clearly, there was a wide range of growth rates for different plants in the same environment. Figure 1(B, E, H, K) shows the frequency distribution of the total growth rate, λ_T, for plants in each environment. In all four environments, some plants had a net decrease in number of nodes during the measurement period: this was due to the loss of nodes at the base of the plant (basal death) outstripping the increase of nodes at the tips. Figure 1(C, F, I, L) shows the distribution of the rate of increase in number of stems, B. In all environments, the majority of plants (>70%) did not branch during the measurement period, giving B = 0 for those plants.

Individual plant stochastic growth model

Rather than following each node and the relationships among nodes individually, the model describes the plant in terms of four key variables:

• The number of young nodes, N_Y. Young nodes have the potential to branch and create a new side stem.

• The number of old nodes, N_O. Old nodes no longer have the potential to branch, either because they are too old or because they have already branched previously. (Data indicate that it is extremely rare for a single node to be the parent for more than one side stem.)

• The number of tip nodes, T.

• The number of nodes (young or old) which have no live parent, i.e., basal nodes, N_B. Basal nodes are able to die. Note that basal nodes are already accounted for in the total number of nodes.

The total number of live nodes is N_Y + N_O + T. The basal nodes category and the total live nodes are not mutually exclusive: basal nodes are a subset of total live nodes. The model makes no distinction between the main stem and side stems, i.e., nodes on a side stem behave in the same fashion as nodes on the main stem.

The model is based on Poisson processes representing four types of event:

• Tip growth. Any stem on the plant can add an additional node at the tip. The tip growth rate is denoted λ_G per tip per unit time. The effect of this event is to increase the number of young nodes by 1 (there is no net change in the number of tips). Tip nodes are not able to produce side stems.

• Node branching. A young node can branch, creating a new side stem consisting of a single new tip. The node branching rate is denoted λ_B per young node per unit time. Since the parent node cannot branch again subsequently, it transitions from being a young node to an old node. The effect of this event is therefore to increase the number of tips by 1, to increase the number of old nodes by 1 and to decrease the number of young nodes by 1.

• Node ageing. A young node can transition into the old node category (and cannot subsequently create new side stems), representing a natural ageing process. The ageing rate is denoted λ_A per young node per unit time.

• Basal node death. The basal node death rate is denoted λ_D per basal node per unit time. Basal node death decreases the number of old nodes by 1 and exposes either 1 new basal node (if the dead node had not branched) or 2 new basal nodes (if the dead node had branched); hence the net effect is either to leave the number of basal nodes unchanged or to increase the number of basal nodes by 1. In the long term, the proportion of nodes that branch before transitioning into the old category is $\frac{{\lambda }_B}{{\lambda }_A+{\lambda }_B}$. Hence, this is taken as the probability that this event increases the number of basal nodes by 1. (Implicit in this is the assumption that the main stem is sufficiently long that all nodes have transitioned into the old node category before they become basal nodes.)

Figure 2 shows a schematic to illustrate the four types of event and Table 2 summarises their rates and effects. The model has no explicit spatial structure because individual nodes are not tracked; only the total number of each type of node is known. However, once these numbers are known, an example spatial structure can be inferred. The assumption that each event type is a Poisson process means that events occur non-simultaneously and the process is memoryless (i.e., depends only on the current state of the plant, not on its prior history). Each simulation of this stochastic model gives a different realisation of an individual plant. By doing repeated simulations, summary statistics, such as the mean rate of increase in the total number of nodes, can be calculated.

The four parameters of the stochastic model can be estimated from the data. The tip growth rate λ_G and the basal node death rate, λ_D, are taken directly from Table 1. By considering the expected time at which a node branches conditional on it not having transitioned from young to old before the branching event occurs, the mean age of a parent node when it branches is ${\displaystyle A={\int }_0^{\infty }{\lambda }_A{e}^{-{\lambda }_A\hspace{0.17em}{t}^{\prime }}{\int }_0^{t^{\prime }}t\hspace{0.17em}{\lambda }_B{e}^{-{\lambda }_B\hspace{0.17em}t}\hspace{0.17em}dt\hspace{0.17em}d{t}^{\prime }=\frac{{\lambda }_B}{{\left({\lambda }_B+{\lambda }_A\right)}^2}.}$ The probability of a node branching is ${\displaystyle P=\frac{{\lambda }_B}{{\lambda }_B+{\lambda }_A}.}$ Taking the values of A and P from Table 1, these two equations can be solved simultaneously to give values for λ_B and λ_A. The initial conditions for the model were chosen to match the state of the plants at the first field experiment census. The initial total number of nodes (N_Y + N_O + T) was the mean total number of nodes at the first census. Over a large number of simulations, the initial number of tips was given the same probability distribution as the number of tips on plants at the first census. As it is not possible to distinguish between young and old nodes in the data, the initial number of young nodes was estimated from the model prediction for the mean number of young nodes per stem, $\frac{{\lambda }_G}{{\lambda }_B+{\lambda }_A}$. Hence, each plant was initialised with $N_Y=T\frac{{\lambda }_G}{{\lambda }_B+{\lambda }_A}$, though changing this initial condition has very little effect on the model output. Since none of the event rates depend on N_O, the model is insensitive to the initial number of old nodes on the plant. Each plant was assumed to have one basal node initially.

Deterministic growth model

For a simple Poisson process of this type, the average behaviour of the stochastic simulations can be described by differential equations. For this, process the mean-field differential equations are given by ${\displaystyle \frac{dT}{dt}={\lambda }_B{N}_Y,}$ ${\displaystyle \frac{d{N}_Y}{dt}={\lambda }_{G}T-\left({\lambda }_A+{\lambda }_B\right)N_Y,}$ ${\displaystyle \frac{d{N}_o}{dt}=\left({\lambda }_A+{\lambda }_B\right)N_Y-{\lambda }_D{N}_B,}$ ${\displaystyle \frac{d{N}_B}{dt}=\frac{{\lambda }_D{\lambda }_B}{{\lambda }_A+{\lambda }_B}{N}_B.}$ Note that the number of basal nodes is independent of the total number of live nodes and can be solved separately to give exponential growth of basal nodes over time at rate ${\displaystyle \frac{{\lambda }_D{\lambda }_B}{{\lambda }_A+{\lambda }_B}.}$ Similarly, the equations for tips and young nodes are uncoupled from the rest of the system and can be solved to show that the number of tips and young nodes will increase exponentially over time at rate ${\displaystyle \frac{-\left({\lambda }_A+{\lambda }_B\right)+\sqrt{{\left({\lambda }_A+{\lambda }_B\right)}^2+4{\lambda }_G{\lambda }_B}}{2}.}$ Thus, the overall behaviour of system is given by the behaviour of the old nodes. If these are being lost by basal death faster than they are increasing by branching and aging of young nodes, then the plant will die in the long term (i.e., N_Y + N_O + T will eventually reach zero). The condition for the death of a plant is therefore ${\displaystyle \frac{{\lambda }_D{\lambda }_B}{{\lambda }_A+{\lambda }_B}>\frac{-{\left(\right.\lambda }_A+{\lambda }_B\left.\right)+\sqrt{{\left({\lambda }_A+{\lambda }_B\right)}^2+4{\lambda }_G{\lambda }_B}}{2}.}$

Model validation

The model, in either its stochastic or mean field format, only requires the first four statistics shown in Table 1 as inputs. This allows us to use the model to predict the final two statistics (total growth rate, λ_T, and branching rate, B) and compare these predictions with the data. The predictions are made using the behaviour of a large number of repeat simulations of the stochastic model or, in the case of the mean total growth rate, using the solution of the mean-field deterministic model. Figure 3 shows the same results as in Fig. 1, but using the model predictions rather than the observed data; the corresponding summary statistics are shown in Fig. 4 and Table 3. In all four environments, model predictions for the mean total growth rate and the mean branching rate lie within the 95% confidence intervals estimated from the data (Fig. 4).

Variance

Although the above simple Poisson process model predicts the mean plant characteristics with reasonable accuracy there are some differences in the distributions of these quantities (compare central column in Figs. 1 and 3). In particular, the variance in the net growth rate is higher in the data than in the stochastic model. Table 1 shows that, of the four model inputs, main stem growth rate λ_G and basal death rate λ_D vary widely between plants in the same treatment. The basal death rate has a coefficient of variation (CV) greater than one and the growth rate has CV ≈1 for all four cases. Conversely, branching age A has CV less than 1 for three of the four treatments.

Table 4 shows the variation seen in the model outputs and the data. The simple Poisson process model described above fails to capture the variance in the total growth rate, though it makes a reasonable estimate for the branching rate. As the underlying plant model is based on a stochastic process, additional sources of variance can be easily included. As the branching process appears to be adequately captured and the parameters that govern it have relatively low variance among plants, this part of the model was not changed. The high-variance model focuses on capturing the variance in the total growth rate by altering the main stem growth rate λ_G and basal death rate λ_D.

The main stem growth rate λ_G, is approximately normally distributed across plants at two of the four treatments (Lilliefors test: Con high p = 0.001, Con low p = 0.16, Ham high p = 0.21, Ham low p = 0.001). Interestingly, it has approximately the same coefficient of variation at each treatment, though the mean varies widely. The basal death rate λ_D is geometrically distributed at three of the four treatments (χ² goodness of fit test: Con high p = 0.004, Con low p = 0.82, Ham high p = 0.08, Ham low p = 0.23). To capture this additional variance we rerun the model but, in each realisation, generate new values for the growth and death rate parameters. We assume the growth parameter is normally distributed with mean λ_G for that treatment and coefficient of variation 0.8 ${\displaystyle {\lambda }_{G_i}\sim N\left({\lambda }_G,{\left(0.8{\lambda }_G\right)}^2\right).}$ In the data the death rate of an individual plant is discrete and well-approximated by a geometric distribution. In the model the death rate is a continuous variable so we assume it is exponentially distributed with mean λ_D for that treatment ${\displaystyle {\lambda }_D\sim Exp\left({\lambda }_D\right).}$ These additional sources of among-plant variability increase the variance of the model realisations (Fig. 5). The coefficient of variation of the model outputs for each treatment is shown in Table 4. This high-variance model is a better fit to the data than the original model, although it still fails in the case of the Conway low light treatment, which has an extremely high CV (66) in the total growth λ_T.

Probability of ultimate extinction

Using the deterministic model, we have calculated the conditions under which a plant is guaranteed to go extinct. Using a stochastic model allows this idea to be extended to predict the probability a plant will go extinct in cases where the criterion derived in the deterministic case does not hold. To assess the chance of survival for a plant at the different treatment levels, we calculate the probability of ultimate extinction (PUE). This is defined as the fraction of realisations (i.e., individual plants in the data or realisations of the stochastic process) that undergo a decrease in the total number of nodes after one year, i.e., λ_T < 0. Table 5 shows PUE for the data with the 95% confidence interval and the predictions from each model. The original Poisson model fails to capture PUE for any of the treatments whereas the high-variance model prediction falls within the confidence intervals in all four cases.