Approximating spatial interactions in a model of forest dynamics as a means of understanding spatial patterns

ecological complexity 3 (2006) 209–218

available at www.sciencedirect.com

journal homepage: http://www.elsevier.com/locate/ecocom

Approximating spatial interactions in a model of forest dynamics as a means of understanding spatial patterns Nicolas Verzelen a, Nicolas Picard b,*, Sylvie Gourlet-Fleury b a b

E´cole Normale Supe´rieure de Paris, 45 rue d’Ulm, 75230 Paris Cedex 5, France Cirad, Forest Department, Campus International de Baillarguet, TA 10/D, 34398 Montpellier Cedex 5, France

article info

abstract

Article history:

Understanding the behaviour of individual-based models of forest dynamics becomes

Received 28 June 2005

difficult as their complexity increases. A useful strategy consists in simplifying parts of

Received in revised form

the model, then scrutinizing for possible differences between the original model’s predic-

16 February 2006

tions and those of the simplified version. This strategy is adopted to understand the role of

Accepted 1 March 2006

space in a complex individual-based space-dependent model of tropical rain-forest

Published on line 20 July 2006

dynamics in French Guiana. The model is made of three components: growth, recruitment and mortality. Two aggregation techniques are used to simplify space-dependent interac-

Keywords:

tions into space-independent interactions: the mean-field approximation and an ad hoc

Aggregation

approximation that considers pair interactions. Each technique is applied to one of the

Complex models

components of the model, the other two being left untouched, or to all components. Nine

Individual-based models

versions of the model (including the original one) are thus compared on the basis of the

Forest dynamics

predicted diameter distribution at stationary state. The fully distance-independent models

Spatial interactions

are consistent with the original one, whereas growth and recruitment significantly modify the model’s predictions when simplified by the aggregation techniques. The interplay between the spatial distribution of trees at the population level and the dynamic processes at the individual level is thus understood, and the excessive impact of the mortality component that favours the establishment of clusters of big trees is diagnosed. # 2006 Elsevier B.V. All rights reserved.

1.

Introduction

Developments in forest dynamics modelling have led to complex individual-based distance-dependent models (Deutschman et al., 1997). This tendency is not limited to forestry but more generally occurs in ecology (Huston et al., 1988; Grimm, 1999; Van Nes and Marten Scheffer, 2005). In this kind of model, the level of description is the tree (or, more generally, the individual) and the environment is spatially explicit so that a local neighbourhood is defined around each tree. A tree then interacts with its neighbours.

However, analysing and understanding such complex models is difficult (Van Nes and Marten Scheffer, 2005). As mathematical analysis often is intractable, the study of the model relies on tedious simulation runs. Moreover, we can lose sight of the important phenomena. If it is easy to see qualitatively the effects of spatial structure on the dynamics, it becomes harder to understand the influence of spatial interactions on the formation of patterns (Parrott, 2005). To better understand the behaviour of complex models, model simplification has arisen as a tool, together with sensitivity analysis, bifurcation analysis, etc. (Murray, 2001; Garcı´a, 2003;

* Corresponding author. Tel.: +33 467 593 939; fax: +33 467 593 733. E-mail address: [email protected] (N. Picard). 1476-945X/$ – see front matter # 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ecocom.2006.03.001

210

ecological complexity 3 (2006) 209–218

Van Nes and Marten Scheffer, 2005). Model simplification consists in replacing the model with a simplified version, then scrutinizing the possible differences between the original model’s predictions and those of the simplified version (Van Nes and Marten Scheffer, 2005). It thus permits the identification of processes, variables or parameters that play a significant role in the model. Moreover the simplified model may sometimes be amenable to mathematical analysis. Model simplification thus has been used to answer the question: How much detail is actually required in a model to maintain its predictive abilities? (Fulton, 1991; Deutschman, 1996; Dieckmann et al., 1997; Lischke et al., 1998; Berec, 2002). This question should not be confused with identifying the processes that maintain the forest ecosystem function (Pacala and Deutschman, 1995), as the complex model may depart from reality (Van Nes and Marten Scheffer, 2005). Two independent models, one simple and the other complex, may actually bring as reliable predictions. Picard et al. (2004) or Gourlet-Fleury et al. (2005), for instance, showed that a sizestructured matrix model (or distribution-based model) could predict the diameter distribution of a forest after disturbance as reliably as an individual-based distance-dependent model (the two models being fitted independently of each other). In that case, Occam’s razor would suggest to use the simpler one. Conversely, the detailed model could be preferred if it includes processes of interest that cannot be included in the simple model. For instance, the size of canopy gaps due to treefall can be predicted with an individual-based spatially explicit model, but it cannot with a Usher matrix model. This means that the level of description used in the model will be also chosen according to the objectives for which the model is used. In the context of tropical forest management for instance, the diameter distribution is the quantity of interest. In this study, model simplification will be treated in the context of aggregation theory (Iwasa et al., 1987; Ritchie and Hann, 1997). It means shifting between two levels of description, a detailed level and an aggregated one. Here, the detailed level is the individual level and the aggregated level is the distribution level. It is known that the aggregation between an individual-based distance-independent model and a distribution-based model is perfect as the number of trees tend to infinity (Garcı´a Vidal, 1974; Uchman´ski and Grimm, 1996; Picard and Franc, 2001). Aggregating from an individual-based distance-dependent model to a distributionbased model then is equivalent with transforming distancedependent tree interactions into distance-independent interactions. The goal of this study is to aggregate distance-dependent interactions into distance-independent interactions for an individual-based model of forest dynamics, to better understand the role of space in this model. This study thus can be seen as a contribution to the strategy of model simplification (Van Nes and Marten Scheffer, 2005). The model used, named Selva, predicts the dynamics of a tropical rain-forest in French Guiana (Gourlet-Fleury, 1997; Gourlet-Fleury and Houllier, 2000). The forest dynamics in Selva is split into three components that are: growth, mortality and recruitment. Three versions of each component will be defined: the original Selva version and two distance-independent versions. By erasing the spatial dependence of one of the components and

leaving the other two components untouched, we will infer the influence of the spatial dependence of this component. The main techniques used in ecological studies to simplify spatial dependence are the mean field-approximation, the moment method, and a technique that we shall call ‘‘remodelling’’ (we did not find any name in the literature for this technique). The mean-field approximation (Levin et al., 1997) is the simplest technique but it relies on strong hypotheses. The moment method is a fruitful technique (Bolker and Pacala, 1997; Dieckmann et al., 1997) but it is tractable only for some type of models (and not for Selva). The ‘‘re-modelling’’ technique consists in using the distance-dependent model as a data generator, then using simulated data to fit empirical distance-independent equations (Acevedo et al., 1995, 1996; Wilson, 1996, 1998; Garcı´a, 2003). In this study we used two aggregation techniques: the mean-field approximation and an ad hoc approximation that is a mix of the moment method and the ‘‘re-modelling’’ technique.

2.

Material and methods

2.1.

The Selva model

The material is the grey-species version of the Selva model, that has been described in detail in Gourlet-Fleury (1997), and Gourlet-Fleury and Houllier (2000). Selva models the dynamics of a tropical rain-forest at Paracou (5
150 N, 52
550 W) in French Guiana (Gourlet-Fleury et al., 2004). A tree in Selva is characterized by its diameter (denoted D) and its spatial coordinates (denoted x 2 R2 ). Time is discrete with a time step Dt ¼ 3 years. By default, Selva models the dynamics of a 250 m  250 m domain (its area, denoted A, thus equals 6.25 ha). Selva is made of three components: recruitment, growth, and mortality. All of them depend on distancedependent indices. In the subsequent sections the main equations and indices of each component of Selva are reviewed.

2.2.

Aggregation techniques

Each distance-dependent component of Selva is separately aggregated into a distance-independent component. Two aggregation techniques are used: the mean-field approximation and an ad hoc approximation. Both primarily rely on the same hypotheses. As in the moment equation method (Bolker and Pacala, 1997, 1999; Levin and Pacala, 1997; Law and Dieckmann, 1999; Dieckmann et al., 2000), a stationary state predicted by the model is seen as a realization of an ergodic marked point process (Cressie, 1991; Rathbun and Cressie, 1994; Stoyan and Penttinen, 2000). Under this hypothesis, the spatial average, denoted hi, of any quantity equals its expectation over all possible realizations of the model. We further assume that the point process is homogeneous and isotropic. We shall denote C the marked point process, and F the corresponding non-marked point process. Given an individual characteristic f (such as growth rate, or recruitment or mortality probability) that depends on distance-dependent indices X1 ; . . . ; X p , aggregation is performed by taking the spatial average h f ðX1 ; . . . ; X p Þi. Depending on the

ecological complexity 3 (2006) 209–218

linearity of f, this will be approximated by f ðhX1 i; . . . ; hX p iÞ or by: p

f ðhX1 i; . . . ; hX p iÞ þ

p

1 XX hðXi  hXi iÞðX j  hX j iÞi 2 i¼1 j¼1

211

As the indicator function I equals 0 or 1, and interpreting the spatial average as the expectation over all possible realizations of an ergodic process, this may be rewritten as: hGi ¼ N PrC ðD0 > D and kx0  xk  RÞ

(1)

j>i

@2 f ðhX1 i; . . . ; hX p iÞ @xi @x j The former approximation is exact if f is linear. The latter approximation ignores the third and higher derivatives of f. Moreover the spatial variance hðXi  hXi iÞ2 i of Xi will be supposed to be related to its spatial average by a power relationship (Fairfield Smith, 1938; Taylor, 1961). The subsequent steps consist in interpreting the spatial averages as expectations over all possible realizations of the model, and in re-expressing these expectations using conditional expectations. The ad hoc approximation consists in using data simulated by Selva in the stationary state to fit empirical expressions to well-chosen quantities, whereas the mean-field approximation additionally assumes that there is no spatial structure in the forest, which means: (i) trees are located at random, and (ii) the position of a tree is independent of its diameter (random labelling hypothesis).

where ðD; xÞ is the subject tree, ðD0 ; x0 Þ a tree taken at random, and the probability refers to the law of the marked point process C. It should be clear that ðD; xÞ is fixed whereas ðD0 ; x0 Þ is random. We then define a function g from the conditional probability: PrC ðkx  x0 k  R j D0 > DÞ ¼ gðD; xÞ PrF ðkx  x0 k  RÞ

(2)

The g function characterizes the competition stress that the stand imposes on the subject tree. As the point process is homogeneous and isotropic, g is actually independent of x (and so does PrF ðkx  x0 k  RÞ). If diameters are independent of tree locations (random labelling hypothesis), g ¼ 1. If the number of trees larger than the subject tree in its neighbourhood is greater than expected under the random labelling hypothesis, g > 1; the subject tree then undergoes higher competition. Introducing the conditional probability (2) in (1) yields:

2.3.

Growth

hGi ¼ NgðDÞ PrF ðkx  x0 k  RÞ PrðD0 > DÞ

2.3.1.

Growth in Selva

The quantity NPrF ðkx  x0 k  RÞ represents the average number of trees at a distance less than R from a subject tree. As the point process is homogeneous and isotropic, it is equal to ðN=AÞ KðRÞ, where A is the area of the domain and K is Ripley’s K-function (Cressie, 1991). We further introduce the distanceindependent competition index: G1 ¼ ðpR2 N=AÞ PrðD0 > DÞ. G1 simply is the stocking of trees larger than the subject tree scaled by the area of the competition neighbourhood. In practise, in simulations, G1 will be computed as its empirical estimate:

The growth equation gives the diameter increment of a tree during one time step. It is written as an average growth term, times a modification factor, times a random error. The average growth term is a Korf growth curve (Zeide, 1993) and depends on the diameter of the subject tree only. The modification factor, that may inflate (if > 1) or reduce (if < 1) the average growth, non-linearly depends on one competition index G (G for growth), and on its variation DG between two consecutive time steps. The random error follows a log-normal distribution such that two consecutive errors for a given tree are positively correlated (Gourlet-Fleury and Houllier, 2000). We shall write in short: DD ¼ aðG; DGÞDt, where DD is the diameter increment and a is the diameter growth rate, to underline the dependence of the growth equation on the competition indices. The competition index G is defined as the number of trees larger than the subject tree within 30 m from it: N X IðDk > Di Þ Iðkxk  xi k  RÞ Gi ¼ k¼1

where N is the number of trees, R ¼ 30 m is the diameter of the competition neighbourhood, i 2 f1; . . . ; Ng is the index of the subject tree, Ið pÞ the indicator function of proposition p (=1 is p is true and 0 otherwise), and k  k denotes the Euclidean norm in R2 .

2.3.2.

Aggregated growth components

The spatial average of the diameter growth rate, haðG; DGÞi, is approximated by aðhGi; DhGiÞ. The spatial average of the competition index is written: hGi i ¼

N X hIðDk > Di Þ Iðkxk  xi k  RÞi k¼1

G1 ¼

N pR2 X IðDk > DÞ A k¼1

Finally: hGi ¼ gðDÞ

KðRÞ G1 pR2

(3)

The computation of DhGi follows from that of hGi. The ad hoc approximation consists in using data simulated by Selva in the stationary state to fit an empirical expression for KðRÞ and gðDÞ. Replacing these empirical expressions in (3) gives a distanceindependent competition index. The mean-field approximation permits us to write gðDÞ ¼ 1 (random labelling hypothesis) and KðRÞ ¼ pR2 (tree are located at random; Cressie, 1991). The competition index then simplifies to G1 .

2.4.

Mortality

2.4.1.

Mortality in Selva

Four types of death are distinguished in Selva: standing death, primary treefall, secondary treefall, and complex treefall. At

212

ecological complexity 3 (2006) 209–218

each time step, the occurrence of dying by primary treefall is drawn for each tree. The probability of primary treefall is computed as the inverse logit of a linear combination of diameter and a distance-dependent competition index M. This index is defined as the cumulative basal area in the circle of radius 30 m centred on the subject tree, and is written:

Mi ¼

N X p k¼1

4

D2k Iðkxk  xi k  RÞ

The probability of primary treefall decreases as M increases. If the tree does not die by primary treefall, the occurrence of its standing death is then drawn. The probability of standing death for a tree is computed as the inverse logit of a linear function of its diameter increment. Secondary treefall correspond to the trees destroyed by the fall of a tree killed by primary treefall. If a tree dies by primary treefall, each tree in the area affected by its fall undergoes a death probability that depends only on the size of the fallen tree. The area affected by the fall of a tree is a 10 m  10 m plot whose location depends on the location of the fallen tree and on its falling direction (Gourlet-Fleury, 1997). Complex treefall corresponds to the trees at the border of a canopy gap opened by a primary treefall, and who fall by loss of balance. Complex treefall occurs if the gap is large enough, that is if the tree died by primary treefall has a diameter greater than 40 cm. The death probability by complex treefall is computed as the inverse logit of a linear function of M. It decreases as M increases. It applies to all trees in a square 30 m  30 m neighbourhood centred on the primary treefall. The process is recursive, which means that a tree that died by a complex treefall can be the origin of another complex treefall, and so on.

2.4.2.

Aggregated mortality components

Aggregating the mortality component into a distance-independent component requires: (1) computing the four death probabilities with distance-independent indices, and (2) choosing trees for secondary and complex treefalls in a neighbourhood-free manner. The latter is obtained by choosing trees in the entire stand rather than in a neighbourhood centred on the primary treefall. Beforehand, death probabilities are scaled to the entire stand: if p is a death probability in a neighbourhood with area S centred on a primary treefall (S ¼ 100 m2 for secondary treefall and S ¼ 800 m2 for complex treefall), the corresponding probability in the entire stand is pS=A. The former point is achieved as follows: death probability by secondary treefall depends on the diameter of the fallen tree only, and is thus already distance-independent. Standing death probability depends on the diameter increment only, and aggregating diameter increments into distance-independent increments was already presented in the growth section. Death probabilities by primary or complex treefall, denoted p, depend in a similar way on the competition index M:

diameter D of the subject tree. The spatial average of the death probability is computed as invlogitðaD þ bhMiÞ in the meanfield approximation, and it is computed using the second derivative of the inverse logit function in the ad hoc approximation: VarðMÞ 2 h pi ¼ invlogitðaD þ bhMiÞ þ b invlogit00 ðaD þ bhMiÞ 2 and VarðMÞ ¼ where invlogit00 ðxÞ ¼ ex ðex  1Þ=ð1 þ ex Þ3 hðM  hMiÞ2 i is the spatial variance of M. The spatial variance of M is assumed to be related to its spatial average by a power relationship: VarðMÞ ¼ k1 hMij1 , where k1 and j1 are coefficients that are estimated from data simulated by Selva at the stationary state. We are now left with computing the spatial average of the competition index M. As its spatial average equals its expectation over all possible realizations of the process under the hypothesis of ergodicity, it may be written as: o p p n 2 2 hMi ¼ N E½D0 Iðkx0  xk  RÞ ¼ N E D0 PrC ½kx0  xk  R j D0  4 4 where ðD; xÞ is the subject tree (fixed) and ðD0 ; x0 Þ is a tree taken at random. We now introduce z such that: 1 PrC ½kx0  xk  R jD0 ¼ z  dz ¼ zðx; D; zÞ PrF ½kx0  xk  R dz where dz is an infinitesimally small diameter interval. z does depend on D as the points in the left hand side of the equation are marked points and thus implicitly include a position and a diameter (this misuse was made to simplify notations). As the process is homogeneous and isotropic, z is independent of x. It characterizes the spatial interaction of trees with diameters D and z. For instance, under the hypothesis of random labelling, z ¼ 1. The spatial average of M may now be written as: p 2 hMi ¼ N E½D0 zðD; D0 Þ PrF ðkx0  xk  RÞ 4 N p KðRÞ X p 2 2 KðRÞE½D0 zðD; D0 Þ ¼ D zðD; Dk Þ ¼N 4A A k¼1 4 k

(4)

The mean-field approximation implies z ¼ 1 and KðRÞ ¼ pR2 , so that the competition index then simplifies to:

M1 ¼

N pR2 X p 2 D A k¼1 4 k

The sum represents the total basal area of the forest stand. The ad hoc approximation consists in using data simulated by Selva in the stationary state to fit an empirical expression for KðRÞ and zðD; D0 Þ. Replacing these empirical expressions in (4) gives a distance-independent competition index.

p ¼ invlogitðaD þ bMÞ

2.5.

Recruitment

where invlogit: x 7! 1=ð1 þ ex Þ is the inverse logit function, and aD and b are coefficients. The probability p decreases as M increases, which implies b > 0 (Gourlet-Fleury, 1997, p. 177). aD depends (primary treefall) or not (complex treefall) on the

2.5.1.

Recruitment in Selva

Recruited trees have diameter D ¼ 10 cm. Recruitment in Selva relies on a regular grid of 10 m  10 m square cells that is superimposed on the plot domain. Let P be the number of

213

ecological complexity 3 (2006) 209–218

recruited trees in a subject cell. It is computed at each time step according to the following rule: if the cell is empty, then P equals 5; otherwise it is randomly drawn according to the distribution:

neighbourhood used for recruitment is not centred on a tree, the spatial average of B is readily computed as the total basal area of the plot scaled down to the recruitment neighbourhood:

PrðP  pÞ ¼ invlogitðm p þ nBÞ for PrðP  5Þ ¼ 1

hBi ¼

p ¼ 0; . . . ; 4 (5)

where m p and n are empirical parameters and B is the cumulative basal area in the subject cell plus its eight closest neighbours (which gives a 30 m  30 m square centred on the subject cell). The probabilities PrðP  pÞ increase as B increases, which implies n < 0 (Gourlet-Fleury, 1997, p.186). Newly recruited trees are placed at random in the cell.

2.5.2.

Aggregated recruitment components

Aggregating the recruitment component into a distanceindependent component requires: (1) computing the law of P using a distance-independent index, (2) computing the number of empty cells, and (3) locating newly recruited trees independently of cells. This latter requirement is simply met by placing newly recruited trees at random in the entire domain. The number of empty cells for the aggregated component is computed as the expectation of the number of empty cells under the hypothesis that trees are located at random. Under this hypothesis, the number of trees in a cell follows a Poisson distribution with parameter NS=A, where S ¼ 100 m2 is the cell area and A is the domain area. Then the expected number of empty cells, denoted E, is: E ¼ ðA=SÞexp ðNS=AÞ. The first requirement is met by taking the spatial averages of the recruitment probabilities (5). These are computed as: (6)

hPrðP  pÞi ¼ invlogitðm p þ nhBiÞ

in the mean-field approximation, and they are computed using the second derivative of the inverse logit function in the ad hoc approximation: hPrðP  pÞi ¼ invlogitðm p þ nhBiÞ þ

VarB 2 n invlogit00 ðm p þ nhBiÞ 2

(7)

for p ¼ 0; . . . ; 4. The spatial variance of B is assumed to be a power function of its spatial average: VarB ¼ k2 hBij2 . As the

N 9S X p 2 D A k¼1 4 k

where S ¼ 100 m2 is the area of a cell (the recruitment neighbourhood being made of nine cells). Finally, the total number of recruited trees at each time step PA=SE can be computed as 5E þ k¼1 Xk , where the Xk are independent and identically distributed random variables that follow the distribution given by (6)(mean-field approximation) or (7)(ad hoc approximation). In this expression, A=S  E simply is the number of non-empty cells.

2.6.

Model comparison

Thus, three versions of each dynamics component (growth, mortality, recruitment) are obtained: the original Selva version, and two distance-independent versions. By combining these different versions of each component, 27 models can be obtained, among which 9 are compared in this study. These nine combinations are listed in Table 1. Models are compared on the basis of (1) the predicted number of trees, (2) the predicted cumulative basal area (that can be seen as a non-centred second-order moment of the diameter distribution), and (3) the predicted diameter distribution in the stationary state. The diameter distribution is here defined as an histogram based on 13 diameter classes: the first 12 classes are 5 cm wide and range from 10–15 cm to 65– 70 cm; the last class gathers all trees with a diameter greater than 70 cm. As the models are stochastic, the stationary diameter distribution (and any variable that is computed from it) is considered as a random vector. As the numbers of trees in the classes are not mutually independent, they cannot be compared separately from each other between models. We thus used a global test of comparison based on a dissimilarity index between diameter histograms. This index is defined as:

0

dðH; H Þ ¼

"
 #1=2 13 X Hi  H0 2 i

i¼1

Hi

Table 1 – Components of the nine models (Selva and eight aggregated models) that are compared Model name

Selva G-mf G-ah M-mf M-ah R-mf R-ah m.f. ad hoc

Component Growth

Mortality

Recruitment

As in Selva m.f. approximation ad hoc As in Selva – – – m.f. approximation ad hoc

As in Selva – – m.f. approximation ad hoc As in Selva – m.f. approximation ad hoc

As in Selva – – – – m.f. approximation ad hoc m.f. approximation ad hoc

The abbreviation ‘‘m.f.’’ stands for ‘‘mean-field’’.

214

ecological complexity 3 (2006) 209–218

Fig. 1 – Estimated g function (solid line) from 600 predictions of the Selva model at stationary state. Dashed lines are confidence intervals at 10% level.

where Hi (respectively H0i ) is the number of trees in the ith diameter class of the diameter histogram H  ðHi Þi¼1;...;13 (respectively H0 ). Let s be the number of repeated simulations, m one of the models listed in Table 1, and Hm ji the stationary number of trees in the ith diameter class predicted by the jth simulation of model m (where i ¼ 1; . . . ; 13 and j ¼ 1; . . . ; s). Let Hm j  ðHm ji Þi¼1;...;13 be the stationary diameter histogram preP ˆ m ¼ ð1 sj¼1 Hm ji Þ dicted by the jth simulation of model m, and H s i the empirical average of the diameter histogram predicted by model m in the stationary state. To compare a model m to ˆ m Þ with the 90% quantile of ˆ Selva ; H Selva, we compared dðH ˆ Selva ; HSelva j ÞÞ ðdðH . j¼1;...;s Each model was run three times starting from an arbitrary initial state. For each simulation, a burn-in period of 1000 time steps was thrown away, and the predicted state was recorded every 10 time steps from time step 1000–3000. We checked beforehand that the duration of the burn-in period was much longer than that required to reach the stationary state, and that 10 time steps was longer than the range of temporal autocorrelations. Thus, s ¼ 600 independent observations were obtained for each model.

3.

Results

3.1.

Empirical relationships

The ad hoc approximation relied on empirical relationships that were fitted to data simulated by Selva in the stationary state. We here specify these relationships. Ripley’s K-function for the stationary state simulated by Selva was estimated as the empirical mean over 600 replicates from Selva. The stationary spatial pattern simulated by Selva is significantly clustered, and the fitted expression for KðRÞ was: KðRÞ ¼ 1:021pR2 . The g function was estimated after dividing diameters in 42 classes, and again as the empirical mean over 600 replicates from Selva. The resulting estimate of g is shown in Fig. 1. g is always greater than one, which indicates that competition is stronger than it would be if trees were located independently of their diameter. Moreover, g increases as D increases: the larger the tree, the stronger the competition. The parameters k1 and j1 of the power relationship that relates the spatial variance of the competition index M to its spatial average were estimated using a linear regression of the logarithm of the variance of M with respect to the logarithm of hMi. Their values are given in Table 2. The same operation was repeated with the competition index B to estimate k2 and j2 (Table 2). Finally, exploratory analysis of simulated data showed that zðD; D0 Þ was approximately independent of D0 . This means that the conditional probability of finding a tree with diameter D0 within 30 m from a tree with diameter D approximately depends on D only. The spatial average of M then can be expressed as a function of its mean-field approximation M1 : KðRÞ hMi zðDÞ M1 pR2 As KðRÞ=ðpR2 Þ was estimated by a constant, z can be assessed by plotting hMi against M1 . A linear relationship was fitted: zðDÞ ¼ a þ bD. The coefficients a and b were estimated by linear regression (Table 2). z linearly decreases from 1.05 to 0.90 as the diameter increases from 10 to 92 cm. It equals one for a diameter of 37 cm. Small trees (D < 37 cm) thus tend to have more trees in their neighbourhood than expected under random labelling, whereas large trees (D > 37 cm) tend to have fewer trees.

Table 2 – Estimates of the coefficients of the empirical relationships used in the ad hoc approximation Coefficient

Estimate (Unit)

Standard error

Model: ln ðVarMÞ ¼ ln k1 þ j1 ln hMi (750 observations, R2 ¼ 0:164; F ¼ 146:3, p-value < 1016 ) k1 1:55  102 (m2(2j1)) 0.414a j1 2.39 0.194 Model: ln ðVarBÞ ¼ ln k2 þ j2 ln hBi (750 observations, R2 ¼ 0:345; F ¼ 394:8, p-value < 1016 ) 1:48  102 (m2(2j2)) 0.126a k2 j2 3.74 0.122 Model: zðDÞ ¼ a þ bD (1080 observations, R2 ¼ 0:75; F ¼ 4424, p-value < 1016 ) a 1.07 9:97  104 b 1:82  103 (cm1) 2:68  105 The equations are fitted to data simulated by Selva in its stationary state. The number relates to the log-transformed value of the parameter.

a

t-Value 9.881a 12.09

Prð > jtjÞ < 1016 < 1016

a

22.6a 19.9

< 1016 < 1016

a

920.5 66.5

< 1016 < 1016

215

ecological complexity 3 (2006) 209–218

Table 3 – Dissimilarity index d between Selva and the other models Model

d

Significance

G-mf G-ah M-mf M-ah R-mf R-ah m.f. ad hoc

0.499 0.324 0.121 0.235 1.94 2.65 0.194 0.235

* *

* *

ˆ Selva j Þ is 0.287. The last column ˆ Selva ; H The 90% quantile of dðH indicates whether the difference between Selva and the model is (asterisk) or not (no asterisk) significant on the basis of the comparison between d and 0.287.

Fig. 2 – Basal area predicted by each of the nine models listed in Table 1 in the stationary state. The solid line, the box and the whiskers represent respectively the median, the 50% and the 95% confidence intervals.

3.2.

Numerical results

We first compared the predicted cumulative basal area and the predicted number of trees in the stationary state for each of the nine models. Figs. 2 and 3 present the median and the confidence intervals for each of the distributions. The meanfield and the ad hoc models both present the nearest cumulative basal area to Selva’s one. On the contrary, the intermediary models present worse basal areas, especially for the models with distance-independent recruitment (R-mf and R-ah). Whereas the mean difference between the basal areas of the ad hoc model and Selva is 0.11 m2 ha1, the nearest intermediate model’s error (M-mf) is 0.8 m2 ha1 and it reaches 20 m2 ha1 for the worst model (R-ah). The results for the number of trees are less clear-cut. Admittedly, the ad hoc model is the closest model to Selva for this characteristic but the intermediate models (except models with distance-

independent recruitment) have results quite similar to Selva. Meanwhile, the mean-field model predicts fewer trees than Selva in the stationary state. Table 3 shows the dissimilarity index d between Selva and the other eight models. We compared them with the 90% ˆ Selva j Þ. If the dissimilarity index is larger ˆ Selva ; H quantile of dðH than the quantile, we say that the difference is significant. It follows that the models G-mf, G-ah, R-mf and R-ah present a diameter distribution significantly different from Selva. Both the mean-field and the ad hoc model are consistent with Selva, but the mean-field model seems a bit better than the ad hoc model. Finally, Fig. 4 shows the diameter distribution in the stationary state for Selva, the mean-field model, and the ad hoc model. The mean-field model tends to underestimate the number of trees for the small diameter classes, whereas the ad hoc model overestimates the number of trees in the largest classes. This is consistent with the larger cumulative basal area for the ad hoc model and with the smaller number of trees for the mean-field model.

4.

Discussion

The stationary state predicted by the mean-field model and the ad hoc model, where all three components (growth, mortality, recruitment) are distance-independent, is consistent with that predicted by Selva. One may thus conclude that distance-dependence can be removed from Selva without degrading its predictive ability in the stationary state. However when one of the component is aggregated using the mean-field or the ad hoc technique while the two other components are left untouched, the predictions of the model can be significantly modified. Let us discuss first the impact of each component.

4.1.

Fig. 3 – Density of trees predicted by each of the nine models listed in Table 1 in the stationary state. The solid line, the box and the whiskers represent respectively the median, the 50% and the 95% confidence intervals.

Aggregation of the growth component

When the growth component is made distance-independent while recruitment and mortality are left untouched (models G-mf and G-ah in Table 1), the predictions of the model are significantly modified (Table 3). The diameter distributions predicted by G-mf and G-ah (not shown here to save space) have an inverse-J shape (like those shown in Fig. 4) that is

216

ecological complexity 3 (2006) 209–218

of primary or complex treefall decreases as the basal area M in the neighbourhood increases. Hence, the probability for a tree to die decreases as the local stocking of trees in its neighbourhood increases. This favours the formation of clusters of big trees.

4.3.

Fig. 4 – Diameter distribution at stationary state predicted by Selva and two distance-independent models (the mean-field model and the ad hoc model, see Table 1).

flatter than that predicted by Selva: large trees are overpredicted and small trees are under-predicted with respect to Selva, which results in a greater basal area (Fig. 2) and a smaller number of trees (Fig. 3). Meanwhile, the spatial pattern of trees gets more clustered than in Selva (Ripley’s K-function was computed for the G-mf and the G-ah models but is not shown here to save space). The distance-dependent competition index G thus draws the spatial pattern of trees towards regularity. As G is replaced by hGi or G1 , there is less restraint to clustering, and the overstocking of big trees may be a secondary consequence of the stronger clustering bound to the mortality component (see next section). The relationship between an asymmetric distance-dependent competition index for growth and a tendency towards the regularity of the spatial pattern has already been stressed (Ford and Diggle, 1981; Hara and Wyszomirski, 1994; Picard et al., 2001). The pressure towards regularity is all the stronger as the size of the competition neighbourhood is small (Picard et al., 2001; Picard and Franc, 2004). The radius of the neighbourhood in Selva is large (R ¼ 30 m), so that the distance-dependent competition index is already a spatial average. Then the discrepancy of the distance-independent growth component is not that important, and the spatial pattern of trees in Selva remains clustered. It could be possible to get a regular pattern for big trees by artificially decreasing the size of the competition neighbourhood (Picard and Franc, 2004).

4.2.

Aggregation of the mortality component

When the mortality component is made distance-independent while growth and recruitment are left untouched (models M-mf and M-ah in Table 1), the predictions of the model are not significantly modified (Table 3). Nevertheless, the spatial pattern of trees gets regular. The distance-dependent mortality component thus draws the spatial pattern of trees towards clustering. This can be understood as the probability

Aggregation of the recruitment component

The recruitment component has a qualitative impact that is similar to that of the growth component (and that is opposed to that of the mortality component): as the recruitment component is made distance-independent while growth and mortality are left untouched (models R-mf and R-ah in Table 1), the predictions of the model are significantly modified (Table 3). Large trees are over-predicted and small trees are under-predicted with respect to Selva, which results in a greater basal area (Fig. 2) and a smaller number of trees (Fig. 3). Meanwhile, the spatial pattern of trees gets much clustered than in Selva. The impact of the recruitment component is the strongest within the three components, with large discrepancies from Selva of the predicted characteristics of the forest stand. Two approximations in the distance-independent recruitment component may explain this behaviour. First, the number of empty cells is under-estimated. Indeed, the number of empty cells was computed using a Poisson distribution for the count of trees. In fact, as the spatial pattern in Selva is clustered, the count of trees is over-dispersed with a greater occurrence of zeros than expected with the Poisson distribution. As recruitment is maximum in empty cells, the distanceindependent recruitment component under-estimates the number of recruited trees. This partly explains the small predicted number of trees. The second approximation relates to the location of newly recruited trees. In Selva, places with a high local stocking of trees have a smaller recruitment, and gaps are filled with recruits, which draws the spatial pattern towards randomness. This is no longer true with the distanceindependent recruitment component, as recruited trees are located at random: highly stocked places thus recruit as much as places with low basal areas. This favours the formation of clusters of big trees, that experience less mortality as seen before.

4.4.

Interactions between the three components

The three components (growth, mortality, recruitment) thus have different, and sometimes opposite, effects on the spatial pattern of the predicted forest stand. The interplay between this spatial pattern and the distance-dependent competition indices may result in significant changes of the number of trees and basal area. However, when all spatial dependencies are removed, the effects of the three components seem to balance and the resulting distance-independent model (meanfield or ad hoc model) is consistent with Selva. We may eventually wonder whether the consistency between Selva and the distance-independent models (mean-field or ad hoc model) is due to the efficiency of the aggregation techniques, or to a coincidence that balances the effects of the three components. Other results that were obtained with the ad hoc

ecological complexity 3 (2006) 209–218

technique (Picard and Franc, 2004) would advocate for the former alternative. Moreover, the discrepancy between Selva and the distance-independent models should be put into perspective with respect to the discrepancy between Selva and the forest at Paracou. On the other hand, we cannot decide which of the mean-field or the ad hoc is the best technique in the present case. In any case, the aggregation techniques allowed us to understand the role of each component (Van Nes and Marten Scheffer, 2005). In particular, the role of the clusters of big trees, that are created by the mortality component despite the growth component, was identified. This could help to improve the Selva model. Other advantages may then be gained by turning back to the observed field data at Paracou. For instance Ripley’s K-function and the g function can be estimated at Paracou and compared to their simulated values. Both K and g functions estimated at Paracou are strongly deviant from their values predicted by Selva. For instance, the g function at Paracou is decreasing and is less than one for large trees. The K function at Paracou indicates that the spatial pattern of big trees is regular till a distance of 10 m. These observations have a biological interpretation: as trees grow, they become regularly spaced and escape to competition (Ford and Diggle, 1981; Picard et al., 2001; Neeff et al., 2005). Thus the clusters of big trees that are generated by Selva are unrealistic. The discrepancy between Selva and either the distanceindependent models or the forest at Paracou should also be replaced in the light of emergent properties (Maurer, 2005; Reuter et al., 2005). Model equations (growth equation, death, and recruitment probabilities) are fitted at the individual level, whereas model predictions are assessed at the stand level, through emergent characteristics such as the spatial pattern (quantified by Ripley’s K-function), the diameter distribution, and the relationship between spatial pattern and diameters (that is the competition pattern, as quantified by the g function). Understanding the emergence of these stand properties from individual equations can be done using aggregation techniques.

4.5.

Extensions

The first extension deals with the development of aggregation techniques to circumvent spatial interactions. Developments may come from the theory of marked point processes, and more precisely assessing the dependence between marks and locations of points (Stoyan and Penttinen, 2000; Schlather et al., 2004). In addition, as the present study was limited to the predicted stationary state, a second extension would be to generalize the present results to transient states. Preliminary results that simulate forest logging indicate that the discrepancies that are observed in the stationary state are maintained in a similar way during transients.

Acknowledgments We thank Pierre Auger for valuable comments at the beginning of this work and an anonymous reviewer for helpful remarks.

217

references

Acevedo, M.F., Urban, D.L., Ablan, M., 1995. Transition and gap models of forest dynamics. Ecol. Appl. 5, 1040–1055. Acevedo, M.F., Urban, D.L., Shugart, H.H., 1996. Models of forest dynamics based on roles of tree species. Ecol. Model. 87, 267–284. Berec, L., 2002. Techniques of spatially explicit individual-based models: construction, simulation, and mean-field analysis. Ecol. Model. 150, 55–81. Bolker, B.M., Pacala, S.W., 1997. Using moment equations to understand stochastically driven spatial pattern formation in ecological systems. Theor. Pop. Biol. 52, 179–197. Bolker, B.M., Pacala, S.W., 1999. Spatial moment equations for plant competition: Understanding spatial strategies and the advantages of short dispersal. Am. Nat. 153, 575–602. Cressie, N., 1991. Statistics for Spatial Data. John Wiley & Sons, New York, p. 900. Deutschman, D.H., 1996. Scaling from trees to forests: the problem of relevant detail. Ph.D. Thesis. Cornell University, Ithaca, NY, 258 pp. Deutschman, D.H., Levin, S.A., Devine, C., Buttel, L.A., 1997. Scaling from trees to forests: analysis of a complex simulation model. Sci. Online 277. Dieckmann, U., Herben, T., Law, R., 1997. Spatio-temporal processes in plant communities. Interim Report IR-97-026, International Institute for Applied Systems Analysis, Laxenburg, Austria, 25 pp. Dieckmann, U., Law, R., Metz, J.A.J. (Eds.), 2000. The Geometry of Ecological Interactions—Simplifying Spatial Complexity. Cambridge University Press, Cambridge, UK, p. 580. Fairfield Smith, H., 1938. An empirical law describing heterogeneity in the yields of agricultural crops. J. Agric. Sci. 28, 1–23. Ford, E.D., Diggle, P.J., 1981. Competition for light in a plant monoculture modelled as a spatial stochastic process. Ann. Bot. 48, 481–500. Fulton, M.R., 1991. A computationally efficient forest succession model: design and initial tests. Forest Ecol. Manage. 42, 23–34. Garcı´a, O., 2003. Dimensionality reduction in growth models: an example. Forest Biometry, Model. Inform. Sci. 1, 1–15. Garcı´a Vidal, O., 1974. On mathematical stand models. Technical report, Management section, Forestry Division, Santiago, Chile, 22 pp. Gourlet-Fleury, S., 1997. Mode´lisation individuelle spatialement explicite de la dynamique d’un peuplement de foreˆt dense tropicale humide (dispositif de Paracou - Guyane franc¸aise). The`se de doctorat. Universite´ Claude Bernard-Lyon I, Lyon, France, 274 pp. Gourlet-Fleury, S., Cornu, G., Je´sel, S., Dessard, H., Jourget, J.G., Blanc, L., Picard, N., 2005. Using models for predicting recovery and assessing tree species vulnerability in logged tropical forests: a case study from French Guiana. Forest Ecol. Manage. 209, 69–85. Gourlet-Fleury, S., Guehl, J.M., Laroussinie, O. (Eds.), 2004. Ecology and Management of a Neotropical Rainforest. Lessons Drawn from Paracou, a Long-Term Experimental Research Site in French Guiana. Elsevier, Paris, p. 311. Gourlet-Fleury, S., Houllier, F., 2000. Modelling diameter increment in a lowland evergreen rain forest in French Guiana. Forest Ecol. Manage. 131, 269–289. Grimm, V., 1999. Ten years of individual-based modelling in ecology: what have we learned and what could we learn in the future? Ecol. Model. 115, 129–148. Hara, T., Wyszomirski, T., 1994. Competitive asymmetry reduces spatial effects on size–structure dynamics in plant populations. Ann. Bot. 73, 285–297.

218

ecological complexity 3 (2006) 209–218

Huston, M.A., DeAngelis, D.L., Post, W.M., 1988. New computer models unify ecological theory. BioScience 38, 682–691. Iwasa, Y., Andreasen, V., Levin, S.A., 1987. Aggregation in model ecosystems. I: Perfect aggregation. Ecol. Model. 37, 287–302. Law, R., Dieckmann, U., 1999. Moment approximations of individual-based models. Interim Report IR-99-043, International Institute for Applied Systems Analysis, Laxenburg, Austria, 15 pp. Levin, S.A., Grenfell, B., Hastings, A., Perelson, A.S., 1997. Mathematical and computational challenges in population biology and ecosystems science. Science 275, 334–343. Levin, S.A., Pacala, S.W., 1997. Theories of simplification and scaling of spatially distributed processes. In: Tilman, D., Kareiva, P. (Eds.), Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions. Princeton University Press, Princeton, NJ, pp. 271–295. Lischke, H., Lo¨ffler, T.J., Fischlin, A., 1998. Aggregation of individual trees and patches in forest succession models: capturing variability with height structured, random, spatial distributions. Theor. Pop. Biol. 54, 213–226. Maurer, B.A., 2005. Statistical mechanics of complex ecological aggregates. Ecol. Complex. 2, 71–85. Murray, A.G., 2001. The use of simple models in the design and calibration of a dynamic 2D model of a semi-enclosed Australian bay. Ecol. Model. 136, 15–30. Neeff, T., Biging, G.S., Dutra, L.V., Freitas, C.C., dos Santos, J.R., 2005. Modeling spatial tree patterns in the Tapajo´s forest using interferometric height. Rev. Bras. Cartogr. 57, 1–6. Pacala, S.W., Deutschman, D.H., 1995. Details that matter: the spatial distribution of individual trees maintains forest ecosystem function. Oikos 74, 357–365. Parrott, L., 2005. Quantifying the complexity of simulated spatiotemporal population dynamics. Ecol. Complex. 2, 175–184. Picard, N., Bar-Hen, A., Franc, A., 2001. Spatial pattern induced by asymmetric competition: a modelling approach. Nat. Resour. Model 14, 147–175. Picard, N., Franc, A., 2001. Aggregation of an individual-based space-dependent model of forest dynamics into distribution-based and space-independent models. Ecol. Model. 145, 69–84.

Picard, N., Franc, A., 2004. Approximating spatial interactions in a model of forest dynamics. Forest Biometry, Model. Inform. Sci. 1, 91–103. Picard, N., Gourlet-Fleury, S., Favrichon, V., 2004. Modelling the forest dynamics at Paracou: the contributions of five models. In: Gourlet-Fleury, S., Guehl, J.M., Laroussinie, O. (Eds.), Ecology and Management of a Neotropical Rainforest. Lessons Drawn from Paracou, a Long-Term Experimental Research Site in French Guiana. Elsevier, Paris, pp. 281–296. Rathbun, S.L., Cressie, N., 1994. A space-time survival point process for a longleaf pine forest in southern Georgia. J. Am. Statist. Assoc. 89, 1164–1174. Reuter, H., Ho¨lker, F., Middelhoff, U., Jopp, F., Eschenbach, C., Breckling, B., 2005. The concepts of emergent and collective properties in individual-based models—Summary and outlook of the Bornho¨ved case studies. Ecol. Model. 186, 489–501. Ritchie, M.W., Hann, D.W., 1997. Implications of disaggregation in forest growth and yield modeling. Forest Sci. 43, 223–233. Schlather, M., Ribeiro, P.J., Diggle, P.J., 2004. Detecting dependence between marks and locations of marked point processes. J. R. Stat. Soc. Ser. B 66, 79–93. Stoyan, D., Penttinen, A., 2000. Recent applications of point process methods in forestry statistics. Stat. Sci. 15, 61–78. Taylor, L.R., 1961. Aggregation, variance, and the mean. Nature 189, 732–735. Uchman´ski, J., Grimm, V., 1996. Individual-based modelling in ecology: what makes the difference? Trends Ecol. Evol. 11, 437–441. Van Nes, E.H., Marten Scheffer, M., 2005. A strategy to improve the contribution of complex simulation models to ecological theory. Ecol. Model. 185, 153–164. Wilson, W.G., 1996. Lotka’s game in predator-prey theory: linking populations to individuals. Theor. Pop. Biol. 50, 368–393. Wilson, W.G., 1998. Resolving discrepancies between deterministic population models and individual-based simulations. Am. Nat. 151, 116–134. Zeide, B., 1993. Analysis of growth equations. Forest Sci. 39, 594–616.