SPATIAL PATTERN INDUCED BY ASYMMETRIC COMPETITION: A MODELING APPROACH

ABSTRACT. The paper addresses the question: how does asymmetric competition for light affect the spatial pattern of trees? It is based on an individual‐based spatially explicit model of forest dynamics, whose growth equations are derived from gap models. The model is calibrated on a stand of natural rainforest in French Guiana, where the tree pattern exhibits regularity at short distances (< 10 m) and clustering at medium distances (∼ 30 m). The model reproduces the regularity but not the clustering. As mortality and recruitment have been modeled so as to favor a random pattern, we conclude that regularity emerges from the asymmetric competition in the growth submodel. Also the scale at which regularity appears is linked to the range of interactions between trees.

diameter or height distribution or classes; tree models follow the trajectory of every single tree. Among the latter, a distinction is generally made between distance-dependent and distance-independent models (Bruce and Wensel [1988]). It relates to how interactions between trees are considered: in distance-dependent models, the location of every tree is known, so that a neighborhood of interacting trees around a subject tree can be explicitly stated. In distance-independent models, on the contrary, there is no spatial information at the tree level, so that interactions between trees result from an averaging procedure over space. Distance-independent models are generally simpler and require less calculation than distance-dependent models do. One question then is to know whether distance-dependent interactions between trees can be degraded into distance-independent processes without loss of realism (Deutschman [1996], Levin et al. [1997], Pacala and Deutschman [1995]). This matter in fact goes well beyond the scope of forest models since it concerns the whole field of spatio-temporal dynamics (Bascompte and Solé [1995], Czárán and Bartha [1992]).
Consider on one hand the spatial pattern formed by trees (Moeur [1993]), and on the other hand a variable of interest that results from growth processes, for instance wood volume. For a distanceindependent model, there is no connection between the two. In a distance-dependent model on the contrary, there is reciprocal action of one on the other. For a stand manager, the question of interest would be: how does the coupling between spatial pattern and growth processes modify wood volume prediction? In this paper, however, we address the counterpart question: how does this coupling affect spatial pattern? More specifically, we are concerned with the incidence of competition on spatial pattern (Brisson and Reynolds [1997], Hanus et al. [1998], Pielou [1962]).
In the era of modern computing, it is possible to track numerically complicated models with several interacting processes (Levin et al. [1997]). However to have access to analysis and understanding, we have favored simplicity. Hence we restrict to one aspect of the previous question, that is: how does asymmetric competition for light affect the spatial pattern of trees?
Our contribution to this question relies on a simple model, that will be described in the first section. Trees are considered as cylindrical boles topped by a flat disk of foliage. The competition pressure exerted on a subject tree is a function of disk overlaps with higher trees. The growth model is calibrated on 12 stands of natural rain forest in French Guiana. Independently of growth data, these stands exhibit regularity between trees on short distances (< 10 m) and clustering on medium distances (∼ 30 m). Our results and discussion then concentrate on the following validation of our approach: whether the spatial patterns simulated by the model correspond to the ones observed on the field or not.

Methods.
2.1 Site description. The data of this study come from Paracou site, located 40 km west from Kourou,French Guiana (5 • 15'N,52 • 55'W). The climate is equatorial with two seasons, an average annual rainfall of 3160 mm, and a dry season that lasts three months. The relief consists of small hills (less than 50 m high) separated by wet areas, with medium slopes (30% maximum).
On this site 12 plots of natural rain forests were set in 1984 by the CIRAD-Forêt. Each plot is a 6.25 ha square with a buffer zone 50 m wide. From 1986 to 1988, the plots underwent three silvicultural treatments: three plots were logged for hardwood, three plots were logged for hardwood plus selective clearing, three plots were logged for both hardwood and fuelwood plus selective clearing, and the last three plots were left untouched as controls. On each plot, the circumference of every tree with a DBH (diameter at breast height) ≥ 10 cm is measured with a precision of 0.5 cm. Its spatial coordinates (± 50 cm) and its species are noted too. Measurements have been carried annually from 1984 to 1995, and once every two years since. Height measurements on a sample of trees have been achieved in 1997 too. Other stand characteristics such as contour lines have also been collected. All data are gathered into a data base coupled to a GIS. A more precise description of Paracou trials is given by Favrichon [1998].

Model construction.
Model equations are derived from forska-type gap models (Leemans [1991], Prentice and Leemans [1990], Prentice et al. [1993], Desanker and Prentice [1994], Lindner et al. [1997]). Gap models (see Shugart [1984], Shugart et al. 1992] for a review) were initially designed to model forest succession over long time laps. They have also been applied to forest management issues (Shugart et al. [1980], Waldrop et al. [1986], Mohren et al. [1991]). A recent development of the gap model forska (Lindner et al. [1997]) is even concerned with modeling tree growth in pure even-aged beech stands. As gap models are applied to more precise issues, their validation is made more acute: while former gap models are tested for their ability to reproduce specific composition (Botkin et al. [1972], Shugart and West [1977]), the latest are tested for their ability to reproduce diameter distribution of every species within the stand (Waldrop et al. [1986], Prentice and Leemans [1990], Prentice et al. [1993]).
As gap models simulate the trajectory of every single tree, they are individual-based models. However tree interactions are not modeled with reference to a spatially explicit tree neighborhood; stand characteristics such as vertical light distribution, stand density, etc., are considered homogeneous over a plot whose size corresponds to the size of the crown of a dominant tree (Prentice and Leemans [1990]). Then every tree interacts with that homogeneous (at the plot level) field, so that it is not possible, for instance, to identify the trees that compete with a designed tree. Hence, gap models rely on a stand description at the tree level and a description of interactions at the plot level.
In our study we aim to see how tree competition can induce spatial pattern. So we need a description of interactions at the tree level, such as in the model sortie (Pacala et al. [1993(Pacala et al. [ ], [1996). A tree is described by its diameter, its height, and its spatial coordinates on the plane. Interactions between trees are summed up by an additional tree variable, the normalized leaf area index. As in gap models, forest dynamics is split into three components that are growth, mortality, and recruitment. We actually focus on growth. Mortality and recruitment are reduced to the simplest expression.
Although the study site comprises more than 200 species, species diversity was not considered (see discussion for a justification of this choice).
2.2.1 Growth submodel. For a tree labeled i, let D i be its diameter, H i its total height, and q i ∈ R 2 its spatial coordinates on the ground. When no confusion is possible, we shall omit the label i. The starting point is the growth equation of forska (Prentice and Leemans [1990], Leemans [1991]): where red is a reduction factor that describes competition for soil resources, s is the vertical leaf distribution within the tree from its lowest branch at height B to its total height H, P (z) is a factor that depends on the available light at the point of coordinates z ∈ R 3 in space, and g and r are parameters. As we are only interested with competition for light, we set red = 1.
The vertical leaf distribution is expressed as where La is the total leaf area and σ is a distribution on [B, H] that describes the shape of the crown. Usually the distributions found are the uniform distribution σ(z) = (H − B) −1 as in forska, or the Dirac distribution σ(z) = δ(z = H) as in jabowa-foret like gap-models (Botkin et al. [1972], Shugart and West [1977]). One question is to choose the precise form of σ. However the calculuses that follow (up to equation (11)) aim at showing that, with a few assumptions, equations are in fact independent of σ.
According to the pipe theory (Mäkelä [1977]), the total leaf area La is proportional to the cross-sectional area of sapwood S: Then assuming that S is proportional to D b , a coefficient C exists such that Equation (4) has been used in various models, including gap models, with b ranging from 1.93 (Kohyama [1989]) to 2.9 (Shugart and West [1977]). Values close to 2 are mostly found (Kohyama [1989]). Hence, as in the gap models jabowa (Botkin et al. [1972]) and forska (Prentice and Leemans [1990]), we have assumed b = 2. Equation (4) can be derived from other considerations: let R be crown radius; a fractal dimension a of crown can be defined so that La is proportional to R a (Mäkelä [1997], Zeide and Pfeifer [1991]). Then, assuming that crown radius is linked to diameter by a simple proportionality coefficient (Bossel and Krieger [1991]), Values of a found in the literature (Mäkelä [1997], Zeide and Pfeifer [1991]) range from 2.13 to 2.76.
The term −r z in equation (1) represents the loss due to respiration and, more generally, wood maintenance. Let H be the mean height: H = zσ(z) dz (for the Dirac distribution, H = H). Then using equations (2) and (3), total loss s(z)rz dz can be rewritten as rcSH.
Assuming that there is a form factor η such that H/H = η, total loss appears to be proportional to HS, which is a rough estimate of the volume of sapwood. So the term −rz in equation (1) expresses a proportional relationship between maintenance losses and the volume of sapwood. Using equation (4), total loss can finally be rewritten as θHD b , where θ is a proportionality coefficient (this expression will be used later in equation (14)).
The term Q = s(z) gP (q, z) dz in equation (1) represents volume gain due to photosynthesis: g is the potential gain when light intensity is not limiting, whereas P is a reduction factor that integrates competition for light. In forska, P (q, z) is calculated as ϕ(I(q, z)), where I(z) is light intensity at the point of coordinates z in space, and ϕ(I) is the photosynthesis response curve of a leaf. The function ϕ has a hyperbolic form. This dependence is observed for a leaf over a few minutes, but when scaling up to a crown during a year, the dependence between average photosynthesis rate and average light intensity tends to be more linear (McCrady and Jokela [1998], Ruimy et al. [1995]). Suppose then that P (z) is proportional to I(z). Light intensity is exponentially decreasing according to a Beer-Lambert law with absorption coefficient k, so that where I 0 is light intensity above the canopy and LAI(x, z) is the leaf area index at height z above the point of coordinates x ∈ R 2 on the ground. The leaf area index can be interpreted as the number of the leaf layer that stands above height z. Its calculation, for a tree labeled i, supposes that its competing neighborhood is known: where l ij ≥ 0 are weights that define the competing neighborhood: l ji > 0 if the tree labeled j casts a shadow on the tree labeled i. Competition for light is asymmetric, meaning l ij > 0 ⇒ l ji = 0. The first term in integral (7) represents self-shading and is noted LAI self i , whereas the second term represents competition and is noted LAI comp i . Suppose that tree crowns do not penetrate each other and that light beams are vertically directed, so that Using equations (2), (4) and (5) and the fact that [B,+∞[ σ(z) dz = 1, LAI self (B) appears to be a constant C/πρ 2 that we call LAI max . Then using equations (2), (4), (6), (8) and (9), the volume gain Q = s(z) gP (q, z) dz can be rewritten as For the Dirac distribution, Q self = exp[−kLAI max ]. For any density function, we can make the transformation from z to u = H z σ(h) dh, which yields Equations (10) and (11) show that the dependence of Q on D and LAI comp is independent of the choice of σ. Hence we shall not have to specify a particular form to σ.
In our model, the competitive pressure exerted by a tree labeled j on a tree labeled i is supposed to be proportional to where I(p) is the indicator function of proposition p, (I(p) = 1 if p is true, = 0, otherwise) and ν (d, u, v) stands for the area of the intersection of two disks of radius u and v a distance d apart. It can be checked that (l ji ) has the asymmetric property. Following equation (10), we also re-express the volume gain Q i of tree i as and P is a decreasing function such that P (0) = 1 and P (+∞) = 0. The precise form of P will be exposed later.
We call L the normalized leaf area index. It can be interpreted as the number of trees that cast a shadow over the subject tree. This variable can also be considered as a distance-dependent competition index (see Biging and Dobbertin [1992], [1995]) for a review on competition indices, that would be classified as an influence-zone overlap index. The major difference with indices such as Gerrard's competition quotient or Bella's competitive influence zone overlap results from the asymmetric term I(H j > H i ).
Equation (1) predicts volume increment. To separate diameter and height increments, we need a further relation. In forska, this is done by assuming a relationship between height and diameter. We have preferred to model volume increment and diameter increment separately. Given our assumptions, equation (1) can be rewritten as Differential equations are transposed into increment equations with a time step of one year. Then, knowing annual diameter increment ∆D, annual height increment ∆H can be derived as follows: Diameter increment is predicted empirically from D and L. Data show that, knowing D and L, the conditional distribution of diameter increment is well adjusted by a modified Weibull law with two parameters defined as follows: a random variable X follows the modified Weibull law with parameters (α, β) if X β follows an exponential law with parameter α. Diameter increment is then predicted by where α and β linearly depend on L and D, respectively, and Γ is for the gamma function.
Function P in equations (14) and (15) is also adjusted to data. First, we define the potential diameter increment ∆D pot at level τ as the (1 − τ )% quantile of the conditional distribution of diameter increment knowing D and L. For our modified Weibull distribution with parameters (α, β), it gives Potential diameter increment can be interpreted as the increment obtained when all photosynthetic gains are invested into diameter growth. Setting ∆H = 0 in equation (15) yields Eliminating ∆D pot from equations (17) and (18), and considering small trees for which D → 0, where β 0 is the value of β when D = 0. Remembering that P (0) = 1, the final expression for P is: where α 0 is the value of α when L = 0. As α is a linear function of L (see equation (21) below) the decrease of P with L is slower than an exponential decrease (as given by equation (6)). Although the exponential model, which finds its foundation in the analogy with the Beer-Lambert law, has been widely used in gap and process-based models (Bossel and Krieger [1991], Mäkelä [1997]), studies suggest that light extinction in canopies would be slower indeed (Chen et al. [1993]).

2.2.2
Mortality submodel. Let f (D, t) be the diameter distribution at time t. Supposing that diameter growth G and mortality m depend on diameter and time only, and neglecting diffusion due to the variability in diameter increments, it can be shown (Kohyama [1989]) that In the stationary state, ∂f /∂t = 0, mortality can thus be estimated as In the stationary state, diameters are distributed according to an exponential law. The mean diameter growth G is estimated from equation (16) by empirically relating L to D. Mortality is modeled as a stochastic process: at each time step and for every tree, a random number between 0 and 1 is drawn. When it is less than m(D), the tree dies.
2.2.3 Recruitment submodel. The number of trees is kept constant by making equal the number of dead with the number of recruited trees every year. Saplings are located at random on the plot. Their initial diameter is D 0 and their initial height H 0 .

Parameter estimation.
To estimate the dependence of α and β on L and D, a data set with the following variables was extracted from Paracou trials: DBH in 1988 (D 88 ), annual diameter increment between 1988 and 1991 (∆D 88 ), and spatial coordinates. As diameter increment is small on average compared to the precision of diameter measurements, it was smoothed over three consecutive years. More precisely, ∆D 88 is the slope of the linear regression of D i over i, (i = 1988 to 1991). We used the period 1988 91 that follows the silvicultural treatments to get the greatest variability in ∆D 88 . In evaluating ∆D 88 , palm trees, trees with an irregular circular shape, and trees that are not present over the whole period 1988 91 are disregarded. The normalized leaf area index is calculated from diameter and spatial coordinates (equations (12) and (13)), assuming that the relation of order on heights is the same as the relation of order on diameters (this hypothesis is not too strong: on a subsample of 591 trees whose height has been measured, the normalized LAI has been calculated from equation (12) with or without replacement of I(H j > H i ) by I(D j > D i ); Pearson's correlation coefficient between both values exceeds 0.95). To handle side effects, trees which are less than 10 m away from a border are also disregarded. The data set finally contains 23843 trees. After numerous experimentation, the selected model is: Parameters are estimated by minimizing the sum of the squared differences between ∆D 88 and diameter increment as predicted by equation (16).
An empirical relationship between mean L and diameter is derived from the same data set, but with restriction to control plots (7892 trees). An exponential function fits the data: Parameters are estimated by minimizing the sum of squared differences between observed normalized LAI and their value as predicted by equation (23).
To characterize diameter distribution in the stationary state, we used DBH measurements in 1984, that is to say before any silvicultural treatment. This data set comprises 46476 trees. Diameter distribution is an exponential law with parameter λ. Then mortality can be calculated from equation (20) (22) and (23), respectively. The calculus of m is done numerically (see Figure 1).
The crown radius / diameter ratio ρ is predicted from measurements made in 1990 of 182 individuals in one of the plot that underwent the first treatment. Crown areas S were estimated from aerial photographs, and ρ is estimated by linear regression between S/π and D.
Estimation of γ and θ relies on the measurements of DBH and height for 591 individuals in a control plot, in 1997. Tree heights were measured with a Bitterlich relascope. We proceeded as follows: suppose γ and θ are known; replacing L with L in equation (16) yields a differential equation on D. Also replacing L with L in equation (15) yields a differential equation on H. These two equations can be solved numerically. A height versus diameter function H = φ(D) is thus obtained, that can be compared with the measurements. Hence γ and θ are estimated so as to minimize the sum of squared difference between measured heights, and heights predicted by φ (see Figure 2).
To initialize saplings, we take D 0 = 1 cm and H 0 = φ(D 0 ) = 2.8 m. This is less than the minimum diameter for census at Paracou site (D min = 10 cm), which supposes that the equations adjusted for D ≥ D min can be extrapolated to D ≤ D min . This is done so that the random positioning of saplings does not affect too heavily the spatial pattern of trees of DBH ≥ 10 cm during the growth from 1 to 10 cm, competition has time to express. The total number of trees on the plot has to be increased accordingly. A density of 3840 trees of DBH ≥ 10 cm on a plot together with an exponential diameter distribution (equation (24)) yields a total number of 9445 trees of DBH ≥ 1 cm on a plot. This number can be compared with a sapling density, as estimated in 1986 on all plots with a sampling rate of 2.3%. The measured density of saplings (H > 1.50 m and D < 10 cm) is 4059 on a plot, which yields a total density of about 7900 trees more than 1.50 m high on a plot. Thus 9445 trees of DBH ≥ 1 cm on a plot may be an over-estimation.

Model runs.
Model equations are summed up in Table 1. They describe forest dynamics in the stationary state. Mortality causes the diameter distribution to become exponential when sufficient time elapses from the initial distribution. To minimize the effects of the initial spatial pattern, the model has to be run for a long time. A first initial state is obtained by locating trees at random, drawing their diameter from an exponential law with parameter λ, and setting H = φ(D). Then the model is run for 1000 years. We checked that this duration was long enough to reach the stationary state, without any relevance to the underlying biological system (the joint distribution of diameters and heights on the stand reaches a stable density in less than 300 years). The final state that is obtained is used as the initial state for further simulations.
Results are based on the repetition of 100 independent runs that last 500 years each.

Characterizing spatial patterns.
A spatial pattern is considered as one realization of a homogeneous and isotropic point process (see Cressie [1991] for a presentation of point process theory). Then it can be characterized through Ripley's K function (Penttinen et al. [1992]). Roughly speaking, NK(r) is the mean number of further points within a distance no more than r from a randomly chosen point. For a homogeneous Poisson point process, K(r) = πr 2 . When the point process generates regular (respectively, attractive) spatial patterns, the K function is less (respectively, more) than πr 2 for at least some r.
We use Ripley's estimator of the K function. To assess departure from randomness, we consider the position of the estimated K function with respect to an envelope based on 100 realizations of a homogeneous Poisson process with the same intensity.

Results.
3.1 Parameter estimates. Table 2 contains a list of definitions of all parameters and their estimates with 95% confidence limits. The correlation coefficient between observed crown radius and the values predicted by equation (5) is just correct (0.60). Yet data show that a proportionality relationship between R and D is the most reasonable. In particular an allometric function does not provide a better fit. Table  3 compares the value of ρ that we have found (0.09 m.cm −1 ) with some published values. Our value falls within the range of those obtained in other forest formations.   Leemans [1990] (see equation (6)  To compare our values of γ and θ with those found in the literature, γ and θ first have to be related to the parameters that are more commonly used in gap models. Following the calculus of Prentice and Leemans [1990], we consider the case of optimal growth (L = 0) in equation (14): where κ is the initial slope of height versus diameter), it becomes: since P (0) = 1. Considering that dH/dD ≈ H/D ≈ κ, we finally obtain: where ∆D max and ∆H max are maximum annual increment in diameter and height of small trees.
Comparing equation (14) with equation (1) (using equation (4)) yields θ ≈ Cr. Table 3 sums up a few values of γ and θ that can be derived from the literature on gap models. Compilation of sylvics (Prentice and Helmisaari [1991]) and process-based models (Bossel and Krieger [1991], Mäkelä [1997]) provide other values. Our estimates of γ and θ are within the scope of usual values.

Control run.
Model checking is based on height versus diameter plots. Figure 2 compares measurements with model. Although the general trend is well reproduced, it appears that height variability is less in simulations than in reality. Figure 3 shows height and diameter distributions. There is good agreement between simulated and measured diameter distribution. As for height, the mode of the distribution is slightly higher in simulations than in reality. and the prototype grain is a disk whose radius follows an exponential law. The agreement between simulations and reality is quite good. It may be noticed that the highest values of normalized LAI are closed to 3. This remark may be linked to conclusions of botanists who worked in French Guiana (Birnbaum [1997], Lescure [1978], Oldeman [1974]), and who observed that mature forest was vertically structured into three main canopy layers.
Finally, the density of trees larger than 10 cm in DBH is 4394±88 on a simulated plot. That is 13% more than the mean density on real plots. As diameter distribution is well simulated for D > 10 cm (Figure 3), it implies that simulations grow too few trees of DBH less than 10 cm, and/or the total number of trees larger than 1 cm is over-estimated. This had already been noticed in Section 2.3. Figure 5 shows Ripley's K function for stand number 1 of Paracou trials, together with the envelopes of 100 simulations of a homogeneous Poisson process with the same intensity, and the envelopes of 100 model runs. Consider first the right plot where the pattern is made from all trees of DBH ≥ 10 cm (that is to say, all the   the stairsteps line represents the estimate for the stand, plain lines are the envelopes of 100 simulations of a homogeneous Poisson process with the same intensity as the stand, and the shaded area is the envelope of 100 repetitions of the model. We have represented K(r)/π − r rather than K(r), so that positive, respectively negative, values indicate clustering, respectively regularity. The left plot is limited to trees larger than 22 cm, whereas the right plot takes into account all the trees on the stand. The underlying strip indicates departure from randomness (black indicates regularity, grey indicates clustering).

Spatial patterns.
trees that have been measured). On distances less than 10 m, the curve is under the envelope of Poisson processes, which indicates regularity. On distances between 20 and 60 m, the curves are above the envelope, which denotes clustering. To conserve space, the curves for the other 11 plots are not presented here. However, they all present significant regularity on short distances (less than 10 m). Clustering is generally present around 30 m, but its intensity varies from one plot to another. On one plot (number 7) clustering does not appear before 50 m, and on another (number 8) there is no clustering at all. For r > 50 m, there is considerable variation in spatial pattern from one plot to another.
On the left plot, trees of DBH ≥ 22 cm are extracted and we consider the spatial pattern they form. Repulsion is still present under 10 m but clustering has disappeared. Consider the spatial pattern defined by all trees larger than a diameter given on the vertical axis. Its Ripley's K(r) function is summarized by the corresponding horizontal strip, as in Figure 5: black indicates regularity, grey indicates clustering. In the right plot (model runs), regularity or clustering is indicated if at least 20 out of 100 simulations show it.
However, when considering bigger trees (DBH ≥ 22 cm, left plot), Ripley's K function of stand number 1 lies within the envelope of model runs on distances less than 20 m. Hence the spatio-temporal process defined by the model is compatible with regularity as observed on real tree patterns for trees of DBH ≥ 22 cm. However it cannot account for clustering on medium distances.
In the same way as the pattern of trees with DBH ≥ 10 or 22 cm has been studied, we can consider the pattern of trees with DBH ≥ D min for any D min . Figure 6 summarizes the results for D min ranging from 10 to 70 cm. Figure 5 represents qualitative rather than quantitative results, since it assesses whether there is regularity, randomness or clustering without comparing Ripley's K function between simulated and real patterns. It shows that the model is able to qualitatively reproduce the regularity on real plots, even when scaling up from small to bigger trees. But it does not account for clustering.

Modeling choices.
The empirical modeling of diameter increment is not successful: the Pearson correlation coefficient between observed and predicted increment is 0.24 (Table 2), which means that less than 6% of the variability of ∆D is explained by the model. There may be two reasons: (1) the model, as defined by equations (16), (21) and (22), is not relevant; (2) DBH and the normalized leaf area index are not the right explicative variables of diameter increment.
First, the fit of the conditional distribution of ∆D knowing D and L to the modified Weibull distribution (equation (16)) is quite good, and problems may not come from this choice. On the contrary, the dependence of the parameters of the modified Weibull distribution on D and L is not well summarized by equations (21) and (22). Better relationships could be proposed, in particular if α is made dependent on D, but they introduce new parameters for a small increase in the percentage of explained variability. In the same way, we could describe the conditional distribution of ∆D by a Burr distribution, which generalizes the Weibull distribution (Lindsay et al. [1996]), but that would increase the number of parameters too.
To assess the relevance of the model (defined by (16), (21), (22)), we have compared its predictions with those of an additive model with loess smoother. Estimates are obtained through a backfitting algorithm. The additive model yields a Pearson correlation coefficient between predicted and observed values hardly better than the former (0.26 as compared to 0.24).
Secondly, as the number of individuals and the variability in ∆D are considerable, the DBH and the normalized leaf area index do significantly contribute to explain ∆D. However these are not the best explicative variables. We can compare with selva, a model of forest dynamics built on the same study site of Paracou by Gourlet-Fleury [1997]. In selva, 23% of ∆D variability is explained from three explicative variables (and 7 parameters): DBH, the number N of trees larger than the individual within 30 m from it, and the variation ∆N of N . In fact diameter and normalized leaf area index may account for local interactions only. So if the explanatory power of our model was to be really estimated, we should be able to separate within the diameter increment which part is due to local interactions and which part is due to other factors. However our goal is not to build a precise empirical model of diameter growth, but to produce a realistic process of competition for light.
We have preferred a simple and explanatory model to a precise one. This is also the reason why species diversity has deliberately (and quite surprisingly for gap models) been disregarded. Taking species diversity into account could increase model precision: in selva, the part of explained variability reaches 42% when distinguishing 15 groups of species (65 parameters). There is a debate on whether the high species richness in tropical rainforest corresponds to a high diversity in ecological niches, or whether it can be reduced to a considerably smaller number of functional groups (Denslow [1987], Hubbell and Foster [1986]). Without taking part in it, we believe that the behavior of the system has to be understood globally, even if crudely, before dividing it into smaller part. Hence taking into account species diversity will be handled in later work.
Although individual characteristics are badly predicted, stand characteristics such as DBH, height or L distributions are correctly modeled (see Figures 3 and 4).
In summary, our approach was to start with functional equations, such as these of gap models, then to calibrate them on field data. Eventually, the set of equations that are gathered in Table 1 define a spatio-temporal process, which can be considered independently of the underlying biological system. Then the question of interest is: can this spatio-temporal process account for the spatial patterns of trees as observed in the field? 4.2 Model behavior. The modeling of diameter increment was based on data collected just after the silvicultural treatments in order to get the maximum of contrast in L and ∆D. Mortality submodel however supposes stationarity. Hence our model cannot account for transient phases, such as evolution after clear-cutting.
As we aim at assessing the role of asymmetric competition on spatial pattern, mortality and recruitment submodels have been designed so as to be as neutral as possible to spatial pattern. Mortality depends on diameter only. Saplings are located at random. So recruitment tends to drive the spatial pattern towards randomness and cannot in any case account for regularity. Imagine now that the growth submodel is distance-independent: then the spatial pattern would remain random (although this result can be obtained with a little thinking, it has been checked on simulations where distance-dependent equations (12) and (13) have been replaced with distance-independent equation (23)). As regularity is observed, it must come from the growth submodel which is distance-dependent.
As diameter is exponentially distributed, small trees represent the main part of the population. Then, when considering all trees on the plot, their spatial pattern may mask the spatial pattern of bigger trees who have experienced competition for a longer time. This may explain in Figure 5 why repulsion on simulated stands is less marked when trees of DBH ≥ 10 cm are taken into account (right plot), than when trees of DBH ≥ 22 cm only are considered (left plot). As repulsion results from the growth process, some time is necessary for competition to express itself before repulsion becomes apparent.
It may seem paradoxical that the growth process, which acts on localized trees, affects their spatial pattern. Of course, when growing, trees do not move. "Movements" are generated by the deaths of trees, immediately replaced with saplings. We can justify the observed regularity by the following argument: a tree close to a big tree grows more slowly than an isolated tree of the same size. Then it stays a longer time in the same diameter class. As mortality rate is a function of diameter only, to reach the same size an isolated and a competing tree will undergo different cumulative survival probabilities. It implies that to observe one tree of a given size near a big tree, we have to start from a bigger contingent than to observe one tree of the same size in an open spot. As the number of trees of a given size is constrained by diameter distribution, the over-mortality of competing trees generates regularity.
The previous argument stands if the stand is made of open and dense spots, that is to say, if the environment is heterogeneous. Suppose now that light conditions are the same everywhere: then, despite its being distance-dependent, the growth submodel would behave like a distance-independent process. So the pattern would be random. For a random spatial pattern, the local density varies more intensively than for a regular pattern. In other words, open and dense spots are more contrasted in a random pattern than in a regular pattern. Then randomness tends to favor regularity and vice versa, until a kind of balance is reached.
What was not a priori obvious was whether, with the parameters calibrated on field data, the balance would tip to detectable regularity or not. Among the set of parameters (see Table 2) a few directly influence competition: α 1 (through equation (21)) and β 0 (through equation (19)) fix the strength of competition, whereas the crown radius / diameter ratio ρ quantifies the range of interactions between trees. More precisely, ρ is the only parameter on which the value of the normalized LAI depends (see Table 1). Thus we also tried to assess the link between the greatest distance of repulsion, as it appears on Figures 5 or 6 and the range of interactions between trees. Figure 7 shows that the distance of repulsion is increasing with ρ. Hence the spatial pattern on short distances reflects the crown overlap.
Eventually, our model can be considered as a spatio-temporal process that accounts for regularity but does not explain clustering on medium distances. Point processes that reproduce the whole pattern of trees have been proposed in the literature (Hanus et al. [1998], Moeur [1997]).
Yet these point processes are built so as to obtain the required spatial pattern, without any consideration to the biological processes that drive the growth of trees. Our model on the contrary is based primarily on the competition for light. As it does not account for the observed clustering on medium distances, let's now discuss other possible causes of departure from randomness.

Possible causes of a nonrandom spatial pattern.
We have shown that introducing an asymmetric reducer can create realistic regularity. Other processes are likely to cause departure from randomness, in particular: (1) competition for soil resources, (2) habitat heterogeneity, (3) seed dispersion.
First, competition for soil resources could easily be added to our model by defining the factor red in equation (1). Analogously to gap models, we could set where V max is the maximum volume bearable by the habitat and U is the distance of influence of a tree when searching for water and nutriments. This index is symmetric since c ij > 0 ⇒ c ji > 0. In the model selva (Gourlet-Fleury [1997]), competition is globally described through the number N i of trees larger than individual i within 30 m from it. We can write N i = j =i c ji , where c ji = I ( q i − q j < 30) I (D j > D i ).
This index is asymmetric. The interesting point is that selva generates clustering at 30 m. Hence the asymmetry of the competition index does not necessarily imply regularity! Secondly, spatial heterogeneity of the habitat may cause departure from randomness, provided that the character that is varying deeply affects the development of the plants. Then the spatial pattern of plants would reflect the spatial pattern of the environmental character. It could be incorporated into the model through a spatially explicit reduction factor red = f (x).
The two former points concern the growth submodel. Recruitment submodel too is likely to act upon the spatial pattern. For instance, a probability of seed dispersion that decreases with the distance from the parent tree would yield clustering.
The final question is how these different processes interact. If each process causes departure from randomness at a different scale, then it should be possible to distinguish each effect. At Paracou, in particular, it is likely that regularity on short distances comes from competition for light, whereas clustering on medium distances comes from soil heterogeneity or seed dispersion.