A Pair-Approximation Model for Spatial Patterns in Tree Populations with Asymmetrical Resource Competition

A pair-approximation model for the spatial dynamics of a height-structured tree population is defined on a regular lattice where each site can be in 1 of 3 states: empty (gap site), occupied by an immature tree, and occupied by a mature tree. The nonlinearities are associated with resource competition effects of mature trees on immature ones (asymmetric competition) affecting the mortality of the latter but not their growth. The survival–extinction transition of the forest is expressed; the early dynamics of colonization are described in terms of local densities. Predictions of the pair-approximation model are compared with results from numerical simulations of cellular automata.


INTRODUCTION
As in the study of other problems in biology and chemistry, in the 1970s partial differential equations (PDEs) were introduced to model the reaction and diffusion processes characterizing vegetation and forest dynamics. This classical framework offers a large-scale description of systems where population densities vary in a continuous space. However, reaction-diffusion PDE models neglect small-scale spatial correlations, and when differences among individuals are significant and local interactions matter, other frameworks are more appropriate (Cronhjort, 2000;Chen et al., 2002). Such are spatially explicit individual-based models, cellular automata, and moment-based models (Gratzer et al., 2004). When individual differences matter but the population spatial distribution, as in forest exploitations, does not, continuously size-structured population dynamics are described by hyperbolic PDEs (Metz and Diekmann, 1986). Goetz et al. (2011) considered this situation in the study of optimal control problems of forest management, an optimization issue which is usually formulated in terms of discrete age-class models (Tahvonen, 2004).
Garcia-Domingo and Saldañ na (2011) extended the pioneering works of Iwasa and coauthors on gap-occupancy models for forest spatial dynamics in order to include a simple vertical layering. This vertical structure of the forest consists of three layers: the canopy (the tallest layer, composed of mature trees), the understory (the layer composed of saplings of canopy trees), and the shrub layer (the lowest layer of woody vegetation). This vertical stratification was also considered in Adams et al. (2007) in the context of interspecific height-structured competition for light. Although simple, it allows for a new type of local interaction among individuals, namely the asymmetrical competition for nutrients and sunlight. Manrubia and Solé (1997) and Pagnutti et al. (2005) dealt with lattice models of forest dynamics with a richer vertical layering; Solé et al. (2005) addressed the interaction between dispersal strategies and vertical forest stratification. These authors used simulations of cellular automata. In contrast, Garcia-Domingo and Saldañ na (2011) theorize on a stochastic lattice forest with a vertical layering, using the so-called pair approximation (Harada et al., 1994(Harada et al., , 1995Rand, 1999;Keeling, 1999;Sato and Iwasa, 2000), which consists in deriving the differential equations describing the dynamics of the total number of ordered pairs ij of adjacent sites which are in state i and j.
The novelty of the model we present is the assumption of the existence of asymmetrical competition between mature and immature trees affecting the mortality of the latter. This sort of competition is based on the fact that mature trees have a larger root system than immature ones, are higher, and have larger crowns. In ecology, it is the so-called self-thinning process, whereby during a stand development, the density of trees decreases as the stand biomass increases. This process is most easily observed in even-aged stands. Li et al. (2000) claimed that it emerges from the ecological interactions among individuals (or local spatial field effects). In this sense, competition for light is one of the most important interactions in plant and tree populations, and the existence of an interspecific trade-off between high survivorship under low light availability (shade-tolerant species) against rapid growth under high light availability (shade-intolerant or light-demanding species) is established by Martin et al. (2010). However, in some tree species, mortality remains constant across different values of the (radial) growth rate in dense self-thinning stands, indicating that mortality can be driven by factors different from light competition in these species (Dekker et al., 2009). Martin et al. (2010) show the existence of exotic invading tree species combining very high growth rates with moderately high shade tolerance, diverging from the growth-survival trade-off pattern of the native species. We shall consider that competition primarily affects the survivorship of saplings but not their growth rate, which is assumed to remain constant for different levels of light. This can be the case of shade-tolerant species for which the effects of competition for light are small compared with those competing for resources (e.g., space, nutrients, water).
In terms of the model, the asymmetrical competition is modeled by an additional mortality of immature trees, which is assumed to be proportional to the total number of neighboring mature trees.

MODEL
The forest region is represented by a two-dimensional regular lattice of N stands or sites. Each site can be empty (O, a gap site), occupied by a mature tree (M) or an immature tree (I). Transitions among these states correspond to three vital processes: birth, growth, and death ( Figure 1). Mature trees produce seeds which germinate in nearby empty sites giving rise to immature trees (saplings). The latter can die, leaving an empty site, or grow and become a mature tree. Similarly, mature trees can die, creating new gaps.
Due to the local interactions among trees, the rates of these processes are not constant but rather depend on the states of the neighboring sites. For instance, the recruitment of immature trees depends on the total number of mature trees around empty sites because the more numerous mature trees surrounding a gap, the higher the probability of seeds germinating in this gap. In turn, the mortality of immature trees is increased by the presence of neighboring mature trees because of asymmetrical competition effects. In addition to natural mortality, mature trees can die due to wind disturbances which cause their fall. This additional mortality is proportional to the total number of gap sites surrounding a mature tree. Immature trees cannot produce seeds and are insensitive to wind stress thanks to a higher flexibility.
We denote by [O], [I], and [M] the total number of gaps, immature trees, and mature trees, respectively, in a two-dimensional regular lattice of N sites. The total number of nearest neighbors per site in the lattice is denoted by z. For a fixed site x, Q x (i) denotes the total number of the nearest neighbors of x which are in state i at time t (0 Q x (i) z). The transitions between states for every site x are symbolized by: where R x , G x , D I x , and D M x are the corresponding transition rates. The growth rate g of immature trees is assumed to be constant, whereas the rest of rates are defined as: where the recruitment is proportional to the total number of neighboring mature trees of a gap site, b=z is the per neighbor-pair fertility rate of a mature tree, d I is the natural mortality rate of immature trees, l=z is the competition effect caused by a neighboring mature tree on an immature tree. The natural mortality rate of mature trees is d M , and d=z is the induced mortality due to the presence of a gap in the neighborhood of a mature tree. Therefore, when all neighboring sites are gaps, the induced mortality is equal to d, the so-called wind disturbance in Kubo et al. (1996). We use the pair-approximation technique to analyze the dynamics of structured populations, which are spatially distributed in regular lattices (Rand, 1999;Keeling, 1999;Sato and Iwasa, 2000). The idea of this method is to extend the mean-field equations for the total number of sites in different states (O, M, and I) by incorporating pairwise interactions between neighboring sites which define the pairs. To approximate the total number of triples arising in the equation for pairs, we assume that the presence of a site in a given state at one end of a triple does not affect the probability of the state of the site at the other end of the triple. The system of equations for the total number of different types of pairs is closed under this triple closure.
Let [ij] be the total number of ordered pairs of adjacent sites in states i and j. Because only nearness in space defines pairs, [ij] ¼ [ji] and pairs of adjacent sites in the same state i are counted twice in [ii]. Let [ijk] be the total number of ordered triples whose sites are in states i, j, and k. Averaging over the lattice, the mean value of the transition rates R x , D I x , and D M x are given by: Using these averages, the differential equation for the total number of immature trees is: Similarly, one can compute the average total number of neighbors in a given state for the sites belonging to a given type of pair.  Proceeding along the same lines for the other site and pair variables, the system of differential equations governing the dynamics of sites and pairs is: This system of equations is closed by introducing a pair approximation for the total number of triples. Here we use the usual pair approximation (Kubo et al., 1996), which is given by: with k ¼ zÀ1 z , although other values of k are possible (Rand, 1999). Under this approximation, the expected total number of (i, j, l) triples in the lattice is equal to the total number of neighbors of those j-sites in (i, j) pairs, (z À 1)[ij], times the conditional probability that a j site has a neighbor in state l, [jl]=(z[j]). This pair approximation and the fact that, for amount to the approximations for the total number of triples of each type: The nonlinearities introduced by the pair approximation lead to the existence of singularities when the system is linearized around the trivial equilibrium To remove these singularities, we re-scale the state variables using the fact that the total number of immature trees approaches 0 at the same rate as that of mature trees (see Section 3.2).

EXTINCTION EQUILIBRIUM
One of the basic aspects of the dynamics of a population model is the stability of the trivial equilibrium. From an ecological point of view, this equilibrium is reached either when a population of trees experiences a change of (exogenous) environmental conditions and goes extinct, or when an initial colonization of an empty area fails to progress. In the first case, we think of a change in the value of one or more parameters and are interested in the extinction-survival transition curve. In the second case, we are interested in the early dynamics with initial conditions close to the trivial equilibrium. In either case, we refer to the trivial equilibrium as ''the extinction equilibrium.''

Existence of an Extinction Threshold
The main tuning parameter in the extinction-survival transition curve (or surface) is the wind disturbance d with d c denoting its critical value, corresponding to the extinction threshold of the system. Adding Eq. (4) and (5)  Moreover, when d is large enough, the competition-induced mortality l is also used as a tuning parameter. In this case, there exists a critical value l c such that for l > l c the tree population goes to extinction. Assuming the existence of a unique nonextinction equilibrium for each d < d c , we follow the branch of these equilibria toward the bifurcation point by taking d ! d À c . The value of the extinction threshold d c depends on other parameters. This dependence defines the survival-extinction transition surface in the parameter space: . To obtain this transition surface numerically, we proceed using the following steps (Garcia-Domingo and Saldañ na, 2011): and substitute their expressions into the other model equations. The condition bg > d(g þ d I ) is necessary for [IM] Ã > 0 and implies b > d. 2. Impose that [IM] Ã > 0 and [MM] Ã > 0 at any nontrivial equilibrium to obtain, from System (18), the lower and upper bounds of 3. Check that the inequality between the previous bounds amounts to which is a necessary (but not sufficient) condition on the values of the parameters for the existence of a nontrivial equilibrium. Here R MF 0 corresponds to the basic reproduction number obtained by decoupling Eq. (4) and (5) (7) and (8)

Stability of the Extinction Equilibrium
To study the linear stability of the extinction equilibrium we transform System {(4),(5),(6),(7),(8)} to get rid of the singularities around this equilibrium.

SPATIAL PATTERNS IN TREE POPULATIONS 183
the positivity of x 0 (t) follows because v(t) also tends to 0, and the first term in the right-hand side of Eq. (24) becomes the dominant one. This guarantees that x(t) > 0 for all t > 0, which implies that x Ã > 0 at any extinction equilibrium, and that Eq. (27) for v 0 (t) has no singularities close to this equilibrium. Whenever y Ã ¼ 0, any equilibrium of System {(24),(25),(26), (27) (27), v Ã > 0 because of the strict positivity of the first term in the right-hand side of this equation, which implies v 0 (t) > 0 for v % 0. At the survival-extinction transition, the equilibrium is given by ðx Ã ; y Ã ; u Ã ; v Ã ; w Ã Þ ¼ ðC; 0; 0; a c 2 ; a c 3 Þ with C given by the solution of Eq. (21) and (22), a c 2 :¼ lim Local stability of the extinction equilibrium is examined by linearizing System {(24), (25),(26),(27),(28)} around the equilibrium (x Ã , 0, 0, v Ã , w Ã ) and computing the eigenvalue of the corresponding Jacobian matrix J Ã with the largest real part, here denoted by k 1 . This eigenvalue turns out to be real and coincides with the stability modulus of J Ã . In particular, both factors on the right-hand side of Eq. (25) must be equal to 0 at the bifurcation equilibrium so that J Ã ðC; 0; 0; a c 2 ; a c 3 Þ has a row whose elements are all 0 and, subsequently an eigenvalue equals 0. Figure 4 shows k 1 parameterized by d and l. From Figure 4, as expected, k 1 ! 0 as d ! d c (left panel) or l ! l c (right panel). As k 1 crosses 0, the computation of the equilibria of System {(24),(25),(26),(27),(28)} shows that this system undergoes a transcritical bifurcation at the equilibrium ðx Ã ; y Ã ; u Ã ; v Ã ; w Ã Þ ¼ ðC; 0; 0; a c 2 ; a c 3 Þ with an exchange of stability between the only admissible (u Ã < x Ã ) nonextinction (y Ã > 0) equilibrium existing for d < d c and the extinction one. For d > d c , the second component of the non-trivial equilibrium becomes negative and this equilibrium has no biological meaning.

Early Dynamics of Colonization
To have an interpretation of why the rescaled extinction equilibrium has strictly positive components is useful to examine the early dynamics of some local densities at the beginning of the colonization of an empty area, when almost all the sites are gaps. As in Garcia-Domingo and Saldañ na (2011) From the phase portrait of the system given by Eq. (28) and (30) (with x ¼ n), this limit system has always a unique equilibrium ðv Ã O ; w Ã O Þ 2 fðv; wÞ 2 ð0; 1Þ Â ð0; 1Þ : v þ w < 1g, and v(t) and w(t) are always strictly positive. Even if we are very close to the extinction equilibrium, the densities [IM]=(z[M]) and [MM]=(z[M]) of immature and mature trees around a mature tree grow quickly. These local densities settle the environmental conditions that will determine the initial success or failure of the colonization of an empty area. Moreover, under a general condition on the parameters values (independent of n), the asymptotic stability of ðv Ã O ; w Ã O Þ is guaranteed (see Appendix). When the extinction equilibrium is unstable, the assumption [O]=N ! 1 is no longer true as the colonization progresses. This means that our previous analysis based on the limit equations is no longer valid and the full system must be considered.

SIMULATIONS
To check the validity of the model predictions, we performed simulations with cellular automata using the transition rates described in Section 2. The size of the lattice is 100 Â 100; periodic conditions are assumed at the boundaries and the size of the neighborhood of each site is z ¼ 4. The tuning parameters in the simulations are, alternatively, the competition-induced mortality rate l of immature trees and the additional mortality d of the mature trees due to wind disturbance. The other parameters are kept constant at: g ¼ b ¼ 0.2 and d I ¼ d M ¼ 0.01. All values of d and l used in the cellular automata satisfy the inequalities g þ d I þ l 1 and d M þ d 1. 1. These constraints and the additional one given by b 1 imply that, at each time step of the simulations, the probability that a given transition takes place in a given lattice site x is equal to the corresponding transition rate defined in Section 2. The initial configuration consists of approximately one third of the sites in each of the three possible states which are randomly distributed over the lattice.
The simulations show a good agreement among the local densities predicted by the pair approximation and those observed in the cellular automata when parameters are far from the survival-extinction transition (Figures 5 and 6). The disagreement between predictions and simulations of cellular automata for parameters close to   the critical values is known (Matsuda et al., 1992;Sato et al., 1994;Levin and Durrett, 1996;Tilman and Kareiva, 1997). It occurs because of long-range spatial correlations between the site states. When l is used as a tuning parameter, predictions about local densities close to extinction are slightly less accurate than those obtained when d is the tuning parameter (Figure 7).
The effect of the competition-induced mortality l on the spatial arrangement of trees is only noticeable for d > 0. In this case, there is a lower bound for l over which an empty site is more likely in the neighborhood of an immature tree than around a mature tree. Meanwhile, the probability of having another immature tree in the neighborhood tends to 0 smoothly. This effect of l is greater for large values of d (as the comparison between the left and right panels of Figure 5 shows).

CONCLUSION
Our pair-approximation model of the spatial distribution of a forest close to extinction shows that the competition-induced mortality l of immature trees caused by mature trees is less determinant for an initial success of colonization in an empty area than wind disturbance d.
Close to extinction, the agreement between model predictions and simulations from cellular automata is lower than in the case where mature trees affect the growth of immature trees by reducing sunlight availability but not survival (Garcia-Domingo and Saldañ na, 2011). The stronger interaction among mature and immature trees assumed in the present model and the fact that d must be large enough for having reasonable critical values of l might be the reason.
The effect of competition-induced mortality l is to reinforce the effect of wind disturbance d when the system approaches the survival-extinction transition. In particular, because we assumed very low natural mortalities for mature and immature trees in the simulations (d I ¼ d M ¼ 0.01), this effect, to be noticeable for realistic values of l c , requires that d ) 0 (top right panel in Figure 3). When l ! l c , as well as when d ! d c , immature trees show a trend to be located between mature trees and gaps ( Figure 8). In the first case, immature trees minimize the competition effects caused by mature trees while, in the second one, mature trees avoid the effect of wind disturbance. When l ! l c , the neighborhoods of mature and immature trees change less abruptly (Figures 5 and 6) than when d ! d c (Figure 7).