Stochastic model for the species abundance problem

Stochastic model for the species abundance problem in an ecological community Simone Pigolotti,1,2 Alessandro Flammini,3,2 and Amos Maritan3,4 1International School for Advanced Studies (SISSA-ISAS), Trieste, Italy 2Istituto Nazionale per la Fisica della Materia (INFM), UR-Trieste, Trieste, Italy 3Dipartimento di Fisica, Università di Padova, v. Marzolo 8, Padova, Italy 4Istituto Nazionale per la Fisica della Materia (INFM), UR-Padova, Padova, Italy (Received 3 March 2004; revised manuscript received 22 April 2004; published 30 July 2004) We propose a model based on coupled multiplicative stochastic processes to understand the dynamics of competing species in an ecosystem. This process can be conveniently described by a Fokker-Planck equation. We provide an analytical expression for the marginalized stationary distribution. Our solution is found in excellent agreement with numerical simulations and compares rather well with observational data from tropical forests. DOI: 10.1103/PhysRevE.70.011916 PACS number(s): 87.23.Cc I. INTRODUCTION One of the most widespread quantities employed in ecology to describe the biodiversity in a given ecosystem is the distribution of species abundance. In operational terms it can be defined as the histogram of the number of species (in a well-defined temporal and geographical context)consisting of a generic number of individuals, or, from a more theoretical perspective, as the probability that a generic species is composed by a certain number of individuals. Data collected in different locations suggests that the relative species abundance distributions show a certain degree of similarity [1]. To elucidate the causes that determine the shapes of these distributions and therefore their similarity is a problem of the utmost importance and not only of theoretical nature: to understand the motives that influences the relative rarity or commonness of different species can be of great help in determining policies for the conservation of the endangered ones. The first studies on this subject can be dated back to the 1940s and are due to Fisher [2]and Preston [3]. Their works were focused on finding distributions that could fit well a particular data set in an empirical way. In particular, Preston [3]argued that the probability of finding species with a certain number of of individuals xshould be lognormal distributed, while Fisher [2]proposed a function of the form e−ax/x, with aⰆ1, the so-called Fisher log series. Later, MacArthur [1]pointed out that similar distributions are found in very different ecosystems, suggesting that the shape of such distributions is to a large extent determined by very basic, general, and ecosystem-independent mechanisms. This in turn hinted at the possibility of predicting the shape of such distributions with simple and general models, without taking into account too many specific details of the ecosystem under consideration. Several models have been proposed that spoused this view [4–6]. Many of them restrict to modeling a single ecological community, a collection of similar species that feed on the same pool of resources in a local area. This definition implies that species belonging to the same community interact mainly in a competitive way: in particular, there are no prey-predator relationships among them. The particular case of a single ecological community can be framed in the wider context of a neutrality hypothesis. The concept of neutrality was introduced in the framework of a biomolecular evolution theory by Kimura [7], and then extended to other fields of biology. In the words of Hubbell [4], an ecological theory can be considered neutral when “…treats organisms in the community as essentially identical in their per capita probabilities of giving birth, dying, migrating and speciating. This neutrality is defined at the individual level, not the species level….” The question whether there exist ecological communities satisfying this assumption is still rather controversial [8], therefore it is crucial to understand what the consequences are of this zero-order hypothesis [9]. In the context of a neutral hypothesis, it is reasonable to describe the number of offspring to which any given individual gives place to as a stochastic variable. As a consequence, the number of individuals in a species at a given time can be regarded as a multiplicative random process. Here we present a model aimed at reproducing the features of species abundance distributions under a minimal set of assumptions: neutrality and the possibility to describe the birth process as a multiplicative random processes. The model translates in a Fokker-Planck equation for the species abundance distribution and is amenable to an analytical treatment. The solutions found are compared with the experimental data that we have. The shape of these solutions depend on one parameter and give, in the two limiting cases, both a lognormal-like curve and the Fisher log series. The paper is organized as follows. In the second section we will present the model and comment on the assumptions made. In the third, we will take the continuum time limit of our model and will provide an analytical solution for the marginalized stationary probability distribution function (PDF). In the fourth section we will analyze what kind of connection there is between our model and the neutral theory [4,6]. In the last two sections we compare our results with the experimental data and comment on them. II. THE MODEL Let us consider an ecological community consisting of a fixed number sof species. According to the MacArthur and PHYSICAL REVIEW E 70, 011916 (2004) 1539-3755/2004/70(1)/011916(5)/$22.50 ©2004 The American Physical Society70 011916-1 Wilson theory of island biogeography [10], the number of species in a community approaches a dynamical equilibrium between immigration, speciation, and extinction. We assume that we can neglect the fluctuations around this equilibrium value; in our model, when a species become extinct, it is immediately replaced by another one. We also assume that the net effect of the competitive interaction between species in the community is just to keep the total number of individuals in the community fixed—the resources available are enough to support just Nindividuals across all the species. This last assumption implies that the populations of the species undergo a zero-sum dynamics. This hypothesis is well confirmed by experimental data [3,10]; at the end of Sec. III we will show that relaxing these constraints does lead to similar conclusions in the large Nlimit. We introduce the s variables xi t, representing the population of the ith specie at (discrete)time t, with the condition: 兺 i=1 s xi t=N∀t. Let P共␭兲be the probability that an individual in the community has ␭offspring during one time step. Here neutrality plays a key role: our assumption implies that P共␭兲is the same for all individuals. The population of the ith species evolves according to the following equation: xi t+1 =N兺k=1 关xi t兴␭k,i t+b 兺j=1 s 共 兺k=1 关xi t兴␭k,i t+b 兲 共1兲 where [] means the integer part. We are assuming that the existence of species with a noninteger number of individuals is not too drastic. This might lead to roundoff problems only for rare species. At each time step (generation)we just sum the number of offspring of every individual belonging to that species and then add a small quantity b. This quantity becomes relevant only for small xi, and this describes the behavior of species near their extinction threshold. We are assuming that the net effect of extinctions, immigration, and speciation can be modeled in a simple way with this term, whose effect is to force the xi’s to be greater than zero. Indeed, for b=0, our system admits an absorbing state with only one xiequal to Nand the others equal to 0, the so-called monodominance [4]. Notice that species are only coupled through the denominator, which simply preserves the normalization condition. The number of individuals of each species will be typically large, so we apply the central limit theorem to the sum of random variables in this equation, obtaining the following model: xi t+1 =N␭ ¯ xi t+ ␴ 冑xi t ␰ i t+b 兺j=1 s共␭ ¯ xj t+ ␴ 冑xj t ␰ j t+b兲 共2兲 where ␭ ¯ and ␴ are the mean value and the root mean square deviation of the distribution P共␭兲, and the ␰ ’s are uncorrelated Gaussian variables with zero-mean and unit variance. It is worth noting the relation between our model and the multiplicative process introduced by Kesten in [11]. Kesten studied random multiplicative processes of the form Xt+1 =␭tXt+bt, where Xtis the variable and both ␭and bare random variables. He found that, depending on the mean value of ␭and on the boundary conditions, one retrieves a lognormal or a power-law regime. Models for ecology and economics based on this kind of processes were proposed by Sornette [12]and Solomon [13]. In our model the number of individuals of different species can be thought as following coupled Kesten-like processes. The coupling is a consequence of the constrain that keeps fixed to Nthe number of individuals in the community and that is enforced in Eq. (1) by the factor Nand by the denominator. III. THE CONTINUUM LIMIT In order to obtain some analytical result, we take the continuous time limit of this model, by introducing the time interval dt in the following way: ␭→1+␭dt b→bdt ␴ → ␴ dt.共3兲 By means of this substitution, our model becomes xi t+dt =xi t+dt共␭ ¯ xi t+ ␴ 冑xi t ␰ i t+b兲 1+dt N兺j=1 s共␭ ¯ xj t+ ␴ 冑xj t ␰ j t+b兲 .共4兲 Expanding the denominator and using the fact that 兺jxj =N, we get the Langevin equation x ˙i=fi共x兲+兺 j=1 s Bij共x兲 ␰ j共5兲 where fi共xi兲=b 冉 1− s Nxi 冊 Bij共xគ兲= 冉 ␦ ij −xi N 冊 冑xj.共6兲 The Fokker-Planck equation [14]associated with this Langevin equation is P ˙共xគ,t兲=−兺 i=1 s ⳵ i 再 −fiP共xគ,t兲+D兺 j ⳵ j关gji共xគ兲P共xគ,t兲兴 冎 共7兲 with D= ␴ 2/2 and gij共xគ兲=gji共xគ兲=兺 kBikBjk = 冉 ␦ ij −xj N 冊 xi.共8兲 We search for a solution of this equation satisfying detailed balance, i.e., Pstfi=D兺j ⳵ j共gijPst兲. Defining the marginalized PDF PIGOLOTTI, FLAMMINI, AND MARITAN PHYSICAL REVIEW E 70, 011916 (2004) 011916-2 p共x兲= 冕 0 ⬁ 兿 j⫽idxjPst共xគ兲共9兲 we can easily obtain an equation for p共x兲 b 冉 1−sx N 冊 p共x兲=Dd dx 冋 冉 x−x2 N 冊 p共x兲 册 .共10兲 This equation can be easily solved, giving: p共x兲⬀x ␤ −1 冉 1− x N 冊 ␤ 共s−1/N兲−1 ␤ =b D.共11兲 Notice that this distribution correctly shows the monodominance behavior ␦ 共0兲or ␦ 共N兲in the limit ␤ →0. Finally, if we fix ␮ = ␤ s/N, in the limit for N→⬁we obtain p共x兲= ␮ ␤ ⌫共 ␤ 兲x1− ␤ e− ␮ x.共12兲 In Fig. 1 we plot simulation of the stationary PDF for various value of the parameter ␤ , and check the validity of (12). Instead of having a system of stochastic differential equations, it is possible to take into account the interaction of a species with the ecosystem in an averaged way. Let us consider the Langevin equation: x ˙共t兲=b+␭ ¯ x− ␥ x+D冑x ␰ 共13兲 where the parameter ␥ takes into account the effect of competition. In order to have normalizable solutions, we have to require that ␥ ⬎␭ ¯ . When this condition holds, it is straightforward to show that the stationary pdf satisfying detailed balance is the same as (12), with ␮ =−共␭ ¯ − ␥ 兲/D. Notice that in this case, the detailed balance solution is exact; it is also remarkable that the stationary distribution (12)can be achieved without fixing neither the number of species, nor the number of individuals. IV. CONNECTION WITH THE NEUTRAL THEORY An interesting question is whether there is some relationship between our model and the neutral theory [4], as formulated by Volkov et al. [6](see also [5]). More precisely, one could ask if our model arises from the continuum limit of a master equation similar to the one proposed for the neutral theory. Let us write the master equation for the generic birth and death process P ˙共x兲=d共x+1兲P共x+1兲+b共x−1兲P共x−1兲 −关d共x兲+b共x兲兴P共x兲共14兲 for x⬎1. We set d共1兲=0 to avoid that a species disappears without being replaced by another one [6]. Equation (14)can FIG. 2. Fit of Barro Colorado Island (BCI)and Pasoh species abundance data, Preston plot [4]. Comparison between our solution and lognormal. Fitted value of the parameters of our distribution are ␤ =0.23 and ␮ =0.010 for the BCI; ␤ =.37 and ␮ =0.015 for Pasoh. In absence of an objective estimate of the error bars on the observational data, both our result and the lognormal give a reasonable fit. FIG. 1. Simulation of marginalized stationary pdf for various values of the diffusion coefficient D, compared with theoretical curves. For all curves b=1,␭ ¯ =1s=100, N=109. Curves are binned linearly with binning size ␦ x=104. Notice that as Dincreases the curve approaches the Fisher log series. STOCHASTIC MODEL FOR THE SPECIES ABUNDANCE…PHYSICAL REVIEW E 70, 011916 (2004) 011916-3 be converted into a Fokker-Planck equation assuming that x is a continuous variable and that b共x兲and d共x兲and P共x兲are smooth enough that we can expand them into a Taylor series. Thus, for example, d共x+1兲P共x+1兲−d共x兲P共x兲= ⳵ ⳵ x关d共x兲P共x兲兴 +1 2 ⳵ 2 ⳵ x2关d共x兲P共x兲兴. Taking the first-order expansion, we obtain from (14)the following Fokker-Planck equation: P ˙共x兲=− ⳵ ⳵ xJ共x兲共15兲 with −J共x兲=关d共x兲−b共x兲兴P共x兲+1 2 ⳵ ⳵ x兵关d共x兲+b共x兲兴P共x兲其 +... . J共x兲has the meaning of a probability current and we can write the general form of the stationary solution satisfying detailed balance as a function of d共x兲and b共x兲setting J共x兲=0 p共x兲=p共1兲 d共x兲+p共x兲exp 冋 − 冕 1 xd共x⬘兲−b共x⬘兲 d共x⬘兲+b共x⬘兲dx⬘ 册 .共16兲 One can easily check that our stationary PDF (12)is recovered, provided the following particular choice of the birth and death coefficients: d共x兲=1+ ␮ 2x−b b共x兲=1− ␮ 2x+b.共17兲 This choice [6]implies that there is balance between immigration and emigration in each species. The more general case, in which this equilibrium does not hold, is treated in [15]. V. COMPARISON WITH EXPERIMENTAL DATA Among the most reliable data on single-trophic species distribution of species abundance is the tropical forest census [16]. In order to make a coarse-graining, a Preston plot is used; data are collected via a logarithmic binning in base 2, and species at the edge between two consecutive binning are equally divided between them. Since we have a continuous probability density, we compared the histogram with the integral over the bins of the distribution with the experimental data and made a least-square fit of the parameters ␤ and ␮ , plus the normalization. We found a good agreement of our predicted curve with the histogram; in Fig. 2 the comparison between our solution and the lognormal is shown. Notice that the two distributions have the same number of fitted parameter. It would be interesting to compare our distribution with data collected form other kind of ecosystems and to try to clarify the dependence of our free parameter ␤ from ecological quantities, such as the immigration pressure, the speciation rate, and the extinction threshold. VI. DISCUSSION AND PERSPECTIVES The model we introduce admits a family of stationary pdf depending on the parameter ␤ . This parameter fully determines the shape of the distribution: for ␤ Ⰶ1, one recovers the Fisher log series; while for ␤ large, one obtains a lognormal-like distribution. As we already pointed out, both of these distributions are well known in the population biology literature as possible candidates to be the “right” distributions found in nature. There is some analogy between our model and the Kesten process. Indeed, also the Kesten process admits two different regimes: one lognormal and one with a power-law tail. The main differences is that in our case the multiplicative random process is applied to the square root of the variables rather than to the variable itself. As a consequence, in the Kesten case, the exponent of the power-law tail of the stationary distribution is always greater than one, while the small ␤ regime of our system is characterized by a power-law tail over many decades, with an exponent that is always less than 1. The cutoff due to the conserved number of individuals ensures the normalization of these long-tailed distributions. It is remarkable that our distribution is the same as found in studies made by Kerner in the 1950s [17]on the invariant measure in a system of Lotka-Volterra equations with purely asymmetric couplings. In that work, the interactions are only of predator-prey type, and the system is deterministic, while we are considering a stochastic system with purely competitive coupling. The discovery of the same distribution in such different models suggests that there might exist some deeper and more general mechanism determining the statistical behavior of the ecosystems, regardless of the type of interactions among species. ACKNOWLEDGMENTS This work is a byproduct of many discussions with J. Banavar and I. Volkov. [1]R. MacArthur, Proc. Natl. Acad. Sci. U.S.A. 43, 293 (1957). [2]R. A. Fisher, A. S. Corbet, and C. B. Williams, J. Anim. Ecol. 12,42(1943). [3]F. W. Preston, Ecology 29, 254 (1948) [4]S. P. Hubbell, The Unified Neutral Theory of Biodiversity and Biogeography (Princeton University Press, Princeton, NJ, 2001). [5]A. McKane, D. Alonso, and R. V. Solé, Phys. Rev. E 62, 8466 (2000); e-print physics/0305022. [6]I. Volkov, J. R. Banavar, S. P. Hubbel, and A. Maritan, Nature (London)424, 1035 (2003). [7]M. Kimura, The Neutral Theory of Molecular Evolution (CamPIGOLOTTI, FLAMMINI, AND MARITAN PHYSICAL REVIEW E 70, 011916 (2004) 011916-4 bridge, University Press, Cambridge, UK, 1983). [8]D. Tilman, Plant Strategies and the Dynamics and Structure of Plant Communities (Princeton University Press, Princeton, NJ, 1988). [9]J. Harte, Nature (London)424, 1006 (2003). [10]R. MacArthur and E. O. Wilson, The Theory of Island Biogeography (Princeton University Press, Princeton, NJ, 1967). [11]H. Kesten, Acta Math. 131, 207 (1973). [12]D. Sornette, Physica A 250, 295 (1998). [13]S. Solomon, e-print cond-mat/9901250. [14]C. W. Gardiner, Handbook of Stochastic Methods (SpringerVerlag, Berlin, 1985). [15]I. Volkov, J. R. Banavar, S. P. Hubbel, and A. Maritan, preprint. [16]R. Condit, S. P. Hubbell, and R. B. Foster, J. Trop. Ecol. 12, 231 (1996). [17]E. H. Kerner, Bull. Math. Biophys. 19, 121 (1957); Bull. Math. Biophys. 21, 217 (1959). STOCHASTIC MODEL FOR THE SPECIES ABUNDANCE…PHYSICAL REVIEW E 70, 011916 (2004) 011916-5