Relating the strength of density dependence and the spatial distribution of individuals

Spatial patterns in ecology contain useful information about underlying mechanisms and processes. Although there are many summary statistics used to quantify these spatial patterns, there are far fewer models that directly link explicit ecological mechanisms to observed patterns easily derived from available data. We present a model of intraspecific spatial aggregation that quantitatively relates static spatial patterning to negative density dependence. Individuals are placed according to the colonization rule consistent with the Maximum Entropy Theory of Ecology (METE), and die with probability proportional to their abundance raised to a power α, a parameter indicating the degree of density dependence. Our model shows quantitatively and generally that increasing density dependence randomizes spatial patterning. α = 1 recovers the strongly aggregated METE distribution that is consistent with many ecosystems empirically, and as α → 2 our prediction approaches the binomial distribution consistent with random placement. In between our model predicts less aggregation than METE, but more than random placement. We additionally relate our mechanistic parameter α to the statistical aggregation parameter k in the negative binomial distribution, giving it an ecological interpretation in the context of density dependence. We use our model to analyze two contrasting datasets, a 50 ha tropical forest and a 64 m2 serpentine grassland plot. For each dataset, we infer α for individual species as well as a community α parameter. We find that α is generally larger in the tightly packed forest than the sparse grassland, and the degree of density dependence increases at smaller scales. These results are consistent with current understanding in both ecosystems, and we infer this underlying density dependence using only empirical spatial patterns. Our model can easily be applied to other datasets where spatially explicit data are available.


INTRODUCTION
negative density dependence and more spatial aggregation with weaker density  First, we briefly review the METE and random placement predictions. In Methods, 92 we derive our new probability distribution and explain how we will compare it to data, as 93 well as introduce the datasets that we will be using for comparison. In Results, we look at where Z Π is a normalization factor, and λ Π is the Lagrange multiplier determined by the 101 constraint condition 102 ∑ n nΠ(n|A, A 0 , n 0 ) = n 0 A A 0 . (2) In the case of a bisection, A = A 0 /2 and the Π function simplifies to 103 Π(n|n 0 ) = 1 which is independent of n. This means that given n 0 individuals, any arrangement of 104 them on the two sides of a bisected plot or quadrat is just as likely as any other.
n R individuals on the right is 116 p L = n L + 1 n L + n R + 2 .
To make our notation consistent with that above, let the number on the left be n and 117 the total number to be n 0 . The probabilities of a new individual arriving on the left or on 118 the right are then: 119 p L (n|n 0 ) = n + 1 n 0 + 2 , p R (n|n 0 ) = n 0 − n + 1 n 0 + 2 . (4) If we place n 0 individuals using this rule, the resulting probability distribution is given by 120 Eq. 3.

121
Random placement 122 Another model for spatial ecology, perhaps the simplest, is the random placement model  Connell, 1971)).

138
In the case of a bisected plot, each death must be on the left or right and the 139 probability of a single death on the left, p D,L , or right, p D,R , is Deriving the Π distribution 141 Now that we have the colonization and death rules (Eqs. 4 and 6), we can derive the 142 general Π α (n|n 0 ) for bisections. At each step in the model, we will have one death population size, we implicitly assume some degree of negative density dependence.
However the strength of the density dependence is determined by the magnitude of α, 148 where stronger density dependence corresponds to larger α. We can solve for the 149 resulting steady state distribution where we reach an equilibrium in the spatial pattern.

150
There are several approaches for deriving the steady state solution for such a system.

151
Here, we equate the rates leaving and entering any individual state Π α (n|n 0 ). We take 152 the probability that we start with n individuals on the left, one on the right dies and then where C(n 0 , α) is the overall normalization that does not have a closed analytic form. In 160 the case that n 0 (α − 1) is large, an approximate form for the normalization is See Supporting information S1 for the details of this derivation.

164
Beyond the first bisection 165 Bisecting the plot more than once (into quadrants, then into 8 cells, etc.) when 166 comparing the model to data allows us to consider how aggregation changes with scale, as well as obtain multiple data points for individual species. We use the following 168 method when we discuss bisecting a single plot more than once. 169 We begin by bisecting the plot in half in one direction, then bisecting each of the 170 resulting plots in the opposite direction. We alternate this bisection pattern until we have  is preferred (see Supporting information S2). 187 We instead find the maximum likelihood α given the data. Inferring α using this  Table S1).

191
For determining α for individual species, we will require multiple bisections and the 192 9/27 sample size p will be roughly the number of cell pairings, p ≈ 2 b−1 . There will be fewer data points in practice as some cell pairings will be empty. 194 We can also define a community α, assuming each species follows the same death 195 rule with identical α. In this case, we will have a larger sample size. For a single and 4, we do not include these error bars as they are smaller than the data points.

205
Data used 206 We will compare our results to two contrasting datasets. First, we will use data from a

222
Comparison to METE and random placement 223 Figure 1 compares the bisection predictions for Π(n) from METE, random placement, 224 and our density dependent model for various α, at n 0 = 10 and 50. In general, our model 225 predicts that increasing negative density dependence (larger α) leads to more random 226 spatial patterning, and less density dependence (smaller α) leads to stronger aggregation. 227 We can relate our distribution directly to both the METE and random placement enough (Supporting information S4 shows this result analytically). For 1 < α < 2, we 232 vary continuously between METE and random placement. We can make the distribution 233 even more spatially aggregated than METE with α < 1 and even less than random 234 placement (overdispersed) with α > 2. 235 We can also relate this distribution to the commonly used conditional negative Note that this approximation holds for 1 ≤ α ≤ 2, which is the ecologically relevant 242 range for most species. See Supporting information S5 for the derivation.

243
Individual species 244 Since the Π function is defined on the species level, we can consider each species 245 separately and find the maximum likelihood α for each. To do this we have to go beyond 246 the first bisection to get multiple data points for the same species at smaller scales.

247
For the serpentine data, we exclude Eriogonum nudum from the following figures as and is higher overall at the BCI dataset, even though the absolute scale is much larger.

253
The spread in α is quite large, but this variation is expected considering the small number 254 of data points, especially for rarer species. Most species have an α between 1 and 2, 255 which is somewhere between the aggregation predicted by METE and random placement.

257
We can instead treat α as a community parameter, using each species as a single data 258 point to recover a community α. Figure 3 shows the direct comparison between our for each of the models are shown in Table 1. In our analysis, we consider α both as a separate parameter for each species (as in Fig. 2,   291 and as a community parameter where each species has the same α (as in Fig. 3 and 292 Fig. 4). The community α is harder to interpret ecologically, but we include it to allow 293 for comparisons with models with community level aggregation parameters (e.g. Conlisk

372
In the BCI data, we do not find any species with high α and high abundance 373 (Supporting information Fig. S4), but otherwise we do not see a trend in abundance. We 374 also find no evidence of a trend with species' mean dbh (Supporting information Fig. S5),

408
They find that increasing randomization decreases the predicted slope of the species-area 409 relationship at the same scale, and therefore upscaling METE will overpredict species 410 richness. In addition, they analyze the effect of variation in aggregation among species, 411 which slightly decreases the slope at small scales and increases the slope at larger scales.

412
This results in a species-area relationship that more closely resembles a power law.  Table 2. Comparison of the Akaike Information Criterion (AIC) for α defined at the individual species level and at the community level in both the serpentine and BCI data and across scales. The AIC is lower at the species level in all cases.   Figure 3. 95% contour intervals for the different predicted probability distributions overlayed with (3a) the serpentine dataset, and (3b) the BCI dataset. For our density dependent model with a community α, α = 1.0003 maximizes the log-likelihood for the serpentine dataset, and α = 1.12 maximizes log-likelihood for BCI dataset.