text:
# Diversity patterns from ecological models at dynamical equilibrium

## I Introduction

One of the most beautiful instances of statistical laws in biology is the species-area law, which relates the area $`A`$ of a habitat and the number of different species $`S`$ coexisting there. In its qualitative form, this law has already been stated by Alexander von Humboldt in the 19th-century: Larger areas harbor more species than smaller ones (see Rosenzweig, 1999). The most commonly used quantitative relation between number of species and area has the form of a power-law,

$$SA^z.$$
(1)
The exponent $`z`$ depends on the geographical characteristics of the region under consideration and on the taxon considered (see, however, He & Legendre, 1996 for a different view).

Despite the large number of studies on biodiversity patterns, their relationship to underlying population dynamic models has hardly been explored until now. This is the scope of the present paper. We obtain scaling laws for biodiversity as a function of external control parameters like the immigration rate and the total amount of abiotic resources. Relating these external parameters to the area of our model island, we obtain a species area relationship that compares favorably simplifito field data.

In this approach, biodiversity is established and mantained by an immigration flux, much like in the phenomenological theory of island biogeography (MacArthur & Wilson, 1963). At long times, the system always evolves to a statistically stationary state where immigrations and local extinctions balance. It is then characterized by a constant turnover of species, but global quantities do not change on the average, and ecological variables like the number of species, the number of trophic levels, the number of links per species or the species abundances attain a time-independent distribution. This state can be called a dynamical equilibrium (MacArthur & Wilson, 1963). We emphasize, however, that a biosystem in this state is strongly driven: its biodiversity depends on a nonzero flux of immigrations. It cannot be described by a given fixed point of population dynamics. The static description of ecosystems, very common in classical population dynamics, seems thus inadequate to model biodiversity.

A large number of field studies support the functional relation given by Eq. (1), with exponents which range from $`0.13`$ to $`0.18`$ for nested areas in the mainland, in the interval $`0.250.33`$ for groups of nearby islands, go up to $`0.50.8`$ for archipelagoes, and range in $`0.170.72`$ for habitat islands (Rosenzweig, 1995, Begon et al., 1998).

Searching for an explanation of insular biodiversity patterns, MacArthur and Wilson (MacArthur & Wilson, 1963) proposed in the sixties an equilibrium theory of island biogeography. According to this theory, the number of species on islands is the result of a dynamical balance between the arrival of new species (immigration) and local extinction. Many field studies and statistical analysis of data have been carried out to define the applicability limits of the theory. They range from island defaunation experiments (Simberloff & Wilson, 1969) and subsequent analysis (Simberloff, 1969, Heatwole & Levins, 1972) to the study of insular biodiversity patterns (Gilpin & Diamond, 1976 ) and of fluctuations of the number of species (Gilpin & Diamond, 1980, 1981; Manne et al., 1998). On the other hand, the lack of explanatory power of MacArthur and Wilson’s theory has been criticized (Williamson, 1989; Whittaker, 1992), and corrections hav been proposed (Simberloff, 1969). The major shortcoming in our view is that the approach is not explicitly founded on ecological dynamics.

Many authors have investigated the species-area relationship without relying on the equilibrium theory of island biogeography. A classical model was proposed by Preston (Preston, 1962), who derived a power-law relationship from the assumption that the abundance of species is characterized by a log-normal distribution. This assumption however has been questioned, since field study support broader distributions. Recently, Harte and collaborators obtained a power-law relationship between species diversity and area from the hypothesis that the spatial distribution of individuals is self similar with respect to the operation of cutting a small area from a larger one (Harte et al., 1999). Thus their model applies to continental nested subareas. An analytical relation between the scale-invariant spatial distribution and the species distribution was subsequently derived (Banavar et al., 1999). Still, the hypothesis of a self-similar distribution of abundances is a strong assumption, even if it seems supported by some field data.

Several recent models combine immigration or speciation and ecological dynamics, trying to overcome the major shortcoming of the previous approaches. Wissel modelled an ecosystem of similar species (Wissel, 1992). He combined the effects of environmental and demographic stochasticity together with interspecies competition thus obtaining a power law species-area relationship. Durrett and Levin proposed a model where speciation is coupled to contact dynamics to mimic ecological processes, obtaining nearly power law species area relationships (Durrett & Levin, 1996). Caldarelli et al. (1998) and Drossel et al. (2000) coupled population dynamics to speciation and immigration processes to simulate changes in biodiversity over evolutionary times. Loreau and Moquet modelled immigrations of plant species from a large pool to an island (Loreau & Mouquet, 1999). They included explicit competition for space in the framework of the equilibrium theory. A recent model for species turnover has reproduced power law distributions for species abundances and species lifetimes as observed in field studies, providing moreover a mechanism for switching from power-law to log-normal distributions as parameters are changed (Solé et al., 2000). The combination of a diffusion mechanism coupled to spatial noise (and leading to local extinctions) generates a power law relationship between number of species and area with exponent $`z0.25`$, even if interaction among species is not explicitly represented (Pelletier, 1999).

Another class of models which combine immigrations and population dynamics is that of species assembly (Post & Pimm, 1983; Drake, 1990; Case, 1990; Morton & Law, 1997; Happel & Stadler, 1998; Schreiber & Gutierrez, 1998). In these models, a community is constructed through local immigrations from a regional species pool. After every immigration, the new community is tested for persistence (i.e., the property that no species gets extinct even in the limit of infinite time). Imposing the condition of persistence bounds these models to the limit of very rare immigrations, while in our model the immigration rate plays the role of the essential control parameter.

The classical view of ecosystems represents them as static ensembles of species at a stable fixed point of population dynamics. In this framework, a key result is constituted by May’s theorem (May, 1972). He showed that a large ecological system formed by $`S`$ species with a total number of $`C`$ random connections with average value zero and dispersion $`\sigma `$ will have, with probability one, no stable fixed point if the variance of the interactions verifies $`\sigma ^2>C/S`$. This result, which has been slightly corrected in more recent years (Cohen & Newman, 1985), sets an upper limit to the amount of complexity allowed by the stability condition, and breaks the traditional view of complexity as a force increasing (static) stability.

The representation of ecosystems as static entities was already challenged by MacArthur and Wilson’s theory. According to them ecosystems are ever changing ensembles of species, evolving in a stationary state where the average number of species remains constant in time. In the same spirit, we take a dynamical view of ecological processes. We model immigrations to an insular ecosystem as a flux of species coming from a continent. Starting from an empty island, new species arrive at random and build the model food webs. This is a possible way of assembling the ecological community (see also Drake, 1990). Population dynamics leads ultimately to species extinctions, but the ecosystem stays in general far from any fixed point of the dynamics.

After a transient period, a statistically stationary state where extinctions balance new immigrations is reached. As will be seen, this stationary state is rather far from static equilibrium: In fact, if the immigration rate is switched off, we observe that several species go rapidly extinct and the ecosystem reaches, more and more slowly as more species go extinct, a static fixed point with a very small number of species. Although we used different models to describe immigrations and population dynamics, we observed that the main statistical features of the equilibrium state are the same for different descriptions. This makes us confident that our results are robust and fairly independent of the details of the modellization.

The relationship between the time scale of the external driving (the immigration rate) and that of the internal ecological dynamics assumes in this context a key importance. If the typical time scale for immigrations is very large compared to the time scale for equilibration of the ecological dynamics, then a new fixed point is reached after every new arrival. In this case, immigrations act only as proposal of new species, but do not shape the ecosystem, whose properties will be independent of the immigration rate. If, on the other hand, the immigration of new species is very fast compared to the ecological dynamics, the ecosystem will be determined only by immigrations and will resemble a random assemblage of species. At intermediate time scales non trivial equilibria emerge and both immigrations and ecological dynamics play an important role. In this regime, the average number of species in the stationary state, $`S`$, increases as a power law of the immigration rate $`I`$,

$$SI^a,$$
(2)
with values for the exponent typically in a narrow range, $`0.4<a1`$, depending on the system parameters.

We shall relate Eq. (2) to the species-area relationship in Sec. IV B, where we assume that the immigration rate is proportional to the linear size of the island (MacArthur & Wilson, 1963). The corresponding exponent $`z`$ depends on parameter values, but is typically comprised in a narrow range, from $`z=0.52`$ to $`0.56`$. These values are very close to those observed for biodiversity in whole archipelagos, when one single source of immigrants is considered. Since our model does not consider spatial structure, it is implicitely assumed that our “island” is either isolated or corresponds to a whole archipelago. In this case, our results show a remarkable agreement with experimental data. It is thus tempting to speculate that the effect of area on biodiversity is largely due to the immigration rate, and that the exponent $`z`$ observed for islands of the same archipelago is smaller because also interchanges of species have to be considered. We shall argue that an approximately power law species-area relationship could be a generic feature arising from the ecological dynamics and the existence of a statistical stationary state, therefore independent of the details of the system.

The immigration rate influences also other properties of the food webs: The number of trophic levels increases with $`I`$, consistently with the observation that food chain length is positively correlated to the size of the ecosystem (Schoener, 1989) and with a recent simulation (Spencer, 1997). Also the number of links per species and the total biomass change with the immigration rate. The distribution of species abundances has a power-law shape, with an exponent close to $`1`$ and slowly decreasing with the immigration rate. The first result is in agreement with field observations (Pielou, 1969), considerations based on the theory of multiplicative processes (Kerner, 1957; Sornette, 1998; Biham et al., 1998) and results from the simulations of a similar model (Solé et al., 2000), and seems to be rather general. The distribution of species lifetimes has also, in an intermediate range, almost power law shape, with an exponent close to $`2`$ slowly decreasing with the immigration rate, as it should be expected. Also in this case this is in agreement with field observations (Keitt & Marquet, 1996; Keitt & Stanley, 1998) and with the results of the model by (Solé et al., 2000).

The paper is organized as follows. In Sec. II we introduce an Individual Based Model of ecological dynamics based on stochastic dynamics. In Sec. III we present a formulation of ecological dynamics based on continuous deterministic models. Since we are interested in the comparison between these two description levels, the results will be presented together in Sec. IV. We conclude with an overall discussion.

## II An individual based model

Recently, population ecology has started to use individual based (or individual oriented) models (IBM) as a complementary tool in the study of ecological dynamics (Łomnicki, 1999; Grimm et al., 1999). One of the main interests of such approach is that it allows the explicit modelization of individual characteristics, like the age of the individuals in a population (influencing the time of breeding or the moment at which they die), or the energy that they store and require to move and survive (Bascompte et al., 1997). Most IBM studies refer to concrete problems where a few species of known characteristics interact to produce a well defined behaviour or pattern, which the IBM should recover or predict (Fahse et al., 1998; Spencer, 1997). Another interest of IBM is in what has been termed virtual ecology: The comparison between real data and simulated data obtained from a system where realistic restrictions have been considered might allow the design of better protocols for recruitment and observation (Berger et al., 1999; Hall & Halle, 1999).

Simulations of very large systems with many individuals and/or many species have not been undertaken until recently because of computational limitations. Thus the IBM approach was restricted to few species in relatively small lattices representing real space, with one to few individuals per lattice site. More ambitious problems, like the relation between theoretical results for deterministic continuous models and their IBM counterpart, were addressed only recently (Keitt, 1997). Some authors derived time-continuous models from the more basic description of the flow of energy between constituents (Svirezhev, 1997) or among individuals (Wilson, 1998; Solé et al., 1999). This is indeed a very relevant point. One would expect that the coarse-grained higher-level description represented by deterministic models captures the essential features of lower-level individual-based models. This is in fact the phylosophy behind our approach: In the IBM, as well as in the higher-level models to be introduced in the forthcoming sections, we study the predictions of the model from a statistical point of view, ignoring details that will necessarily be different in different models.

### A Ecological dynamics

Consider a large area $`A_{\mathrm{cont}}`$ on which a maximum number $`M_h`$ of basal species coexist and up to $`M_a`$ animal species compete for resources. The ecological interactions in this community will be defined through a matrix with entries $`C(m_i,m_j)`$. Depending on the values of $`C(m_i,m_j)`$ and $`C(m_j,m_i)`$ we will determine the trophic relationship between individual $`i`$ of species $`m_i`$ and individual $`j`$ of species $`m_j`$, as we shall explain later. We have considered two possible algorithms to determine the non-zero elements of the interaction matrix. Our first election corresponds to the cascade model (Cohen et al., 1990), which returns a network with topological properties comparable to those of real ecosystems. In this case, the distinction between basal and animal species automatically arises from the ecological relationships given by the interaction matrix. We define $`M=M_h+M_a`$ to be the total number of species in the system. If the number $`m_i=1,\mathrm{}M`$ specifies a peching order for feeding, the algorithm works as follows: Any species $`m_i`$ can feed only on species $`m_j`$ which is lower in the order, that is, $`m_j<m_i`$. This avoids the formation of loops. A link to any of the potential prey species is established with probability $`\mathrm{}/M`$. If the value of $`\mathrm{}`$ is fixed (according to real observations) to be around four, this model returns the correct proportions of basal, intermediate, and top species, a maximum number of levels typically around ten, and a distribution of the number of predators per prey which agrees with field observations (Cohen et al., 1990).

A second possibility for the interaction matrix consists in randomly assigning $`\mathrm{}`$ preys to each of the $`M_a`$ animal species. This would correspond to a disordered situation where no processes have acted in order to select the topology of the ecological network. In this case, and only for implementation purposes, basal species occupy positions $`m_i=1,\mathrm{}M_h`$, and animal species occupy $`m_i=M_h+1,\mathrm{},(M_h+M_a)`$. The interaction matrix has the form

$$\left(\begin{array}{cc}0& 0\\ & \\ C(m_a,m_h)& C(m_a,m_a)\end{array}\right)$$

where $`m_h`$ and $`m_a`$ indicate basal species and animals, respectively. The statistical properties of the system do not depend on the form chosen for the interaction matrix. As we will see, the relevant quantities take the same form in the cascade model case (CM) and in the random matrix case (RM). For both algorithms, the values of the matrix elements are randomly chosen from a uniform distribution in $`[0,E]`$, where $`E[20200]`$ (see Table I). The value of the matrix coefficients is proportional to the energy gained by individual $`i`$ when feeding on individual $`j`$ and represents a sort of assimilation efficiency (see below).

In determining the matrix $`C`$, we have essentially defined a structured ecosystem in a very large area with many species. This is what we consider to be the continent, which will be the source of immigration of propagules to an island of area $`N_h`$. This last quantity can be understood as the maximum number of patches covered by grass, for instance, and acts as a limiting value (together with the basal growth, to be defined) for the number of animals that will inhabit the island.

As time proceeds, and once we properly define the immigration mechanism, we will have a number $`n_h(t)`$ of individuals in basal species present on the island and a number $`n_a(t)`$ of individuals belonging to animal species. The total number of individuals in a wide sense (say patches of grass plus animals) is $`n(t)=n_h(t)+n_a(t)`$. Each individual is characterized by an energy $`e(i)`$. Individuals reproduce provided their value of $`e(i)`$ is large enough. Basal species increase their energy at a constant rate. Animals dissipate energy as time elapses, and increase the value of $`e(i)`$ through predation, which happens stochastically. At each time-step the following rules are implemented:

1. Pair formation. At each time step, we randomly form $`n(t)/2`$ pairs of individuals, independently of their specific affiliation. If $`n(t)`$ is odd, one individual remains without partner. This rule can be understood as a mean field picture of a space-explicit approach. Different possibilities are i) $`hh`$ pair: the two grass patches are not consumed by animals and keep their energy, ii) $`ah`$ pair: if the matrix element $`C(m_a,m_h)`$ is positive (meaning that the individual $`a`$ feeds on $`h`$), predation is possible, iii) $`aa`$ pair, allowing predation between different animal species depending on the matrix coefficients.
2. Predation and Feeding. Either of the individuals in each pair $`(i,j)`$ can feed on its partner, according to the ecological relations defined in the matrix $`C(m_i,m_j)`$. Predation happens when $`C(m_i,m_j)0`$ and $`C(m_j,m_i)=0`$. In this case,

$$\begin{array}{ccc}e(i)\hfill & \hfill & e(i)+C(m_i,m_j)\times \frac{e(j)}{E_{\mathrm{rep}}}\hfill \\ e(j)\hfill & \hfill & 0,\hfill \end{array}$$

where $`E_{\mathrm{rep}}`$ is an energy scale related to reproduction, that will be defined below. The energy received is proportional to the matrix element, but also to the energy stored in the predated individual. In this sense, to eat a new born is not equivalent to eating an adult close to its reproductive energy (which fixes the maximum energy). Furthermore, an individual with a total energy $`e(i)>E+E_{rep}`$ cannot further increase the value of $`e(i)`$. In addition, if the value of the fraction $`e(j)/E_{rep}`$ is larger than unity the rule is modified as $`e(i)e(i)+C(m_i,m_j)`$.

If both $`C(m_i,m_j)`$ and $`C(m_j,m_i)`$ are non-zero (or both zero), no interaction takes place.
3. Basal growth. Every individual belonging to a basal species increases its energy at each time step by a net amount $`b`$,

$$e(i)e(i)+b\mathrm{if}m_i[1,M_h].$$
4. Dissipation. At each time step, and for each of the alive animals, $`e(i)e(i)d`$, where $`d`$ defines the dissipation rate. It takes the same value for all species in our model.
5. Reproduction. If $`e(i)E_{rep}`$, the individual $`i`$ is allowed to reproduce. In the case of basal species, the new individual is introduced into the system provided there is place, i.e. if $`n_h(t)<N_h`$. The individual which reproduces loses an amount of energy $`\delta `$,

$$\begin{array}{ccc}e(i)\hfill & \hfill & e(i)\delta \hfill \\ e(k)\hfill & \hfill & \delta ,\hfill \end{array}$$

where $`e(k)`$ is the energy of the new born.
6. Death. An individual can die for three different reasons: If its energy reaches zero, if it is eaten by a predator, and with a fixed probability $`p_d`$ per time step.

Table I resumes the parameters of the model and the approximate range of values used in our simulations. We will present results for some representative cases. No qualitative differences were observed for comparable sets of parameters.

### B Modelling immigrations

The initial quenched matrix $`C(m_i,m_j)`$ can be thought of as the pool of species in the continent, where a very large area (with its resources) allows the coexistence of all possible species, in our case $`M=M_h+M_a`$. An island has a finite area $`N_h`$ and harbors only a subset $`M_{\mathrm{isl}}`$ of $`M`$.

The immigration flux $`I`$ can take values from the set $`\{\mathrm{}\mathrm{\hspace{0.33em}1}/4,\mathrm{\hspace{0.33em}1}/3,\mathrm{\hspace{0.33em}1}/2,\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}3},\mathrm{}\}`$ only. If $`I=q1`$, then $`q`$ new individuals randomly chosen from any of the possible $`M`$ species in the pool arrive to the island at each time step. If $`I=1/q<1`$, one individual is introduced every $`q`$ time steps. Other situations, which imply a less smooth flux are excluded in the following.

Our simulations show that the average number of species coexisting on the island depends very strongly on the vertical transmission of resources, as is well known to happen in real ecosystems (Rosenzweig, 1995). High dissipation relatively to basal growth ($`d`$ close to $`b`$) turns into few species on the island. For $`db`$ (a factor of 2 or 3 may suffice) the average number of species coexisting when $`I`$ is large enough approaches the maximum number $`M`$.

With the addition of a constant flux of species from the continent to the island, the system is poised to a state of dynamical equilibrium, where the number of species that disappear due to the ecological interactions or to demographic stochasticity is balanced by the new incoming species. The immigration rate might produce a rescue effect for species with few individuals, close to extinction, and at the same time includes in a natural way one form of environmental stochasticity.

Thus, the incoming flux of individuals from the continent, the immigration rate, becomes our main variable. By changing its intensity, we can calculate the average number of species $`S`$ present in the statistically stable regime on an area $`AN_h`$. Moreover, assuming a relation between the immigration rate $`I`$ and the area $`A`$, we shall derive the species-area law resulting from the ecological dynamics of the IBM and compare it with field measurements.

## III Deterministic continuous model

In this section we present the deterministic continuous models of population dynamics adopted in our simulations. All individuals in each species are grouped together and we represent them through a single dynamic variable, the density of biomass (or abundance) of species $`i`$ at time $`t`$, $`N_i(t)`$.

### A Ecological dynamics

The density of biomass of the species $`i=1,\mathrm{},S`$ evolve through a system of differential equations,

$$\frac{dN_i}{dt}=\alpha _iN_i(t)\beta _i\left(N_i(t)\right)^2+\underset{j}{}\gamma _{ij}N_j(t)N_i(t),$$
(3)
determining the growth rate of the biomasses as a function of the abundances of all species in the ecosystem.

Species with biomass less than a predefined threshold value $`N_c`$ go extinct and are eliminated from the system. This mimics the effect of demographic stochasticity and the fact that species are made of discrete entities.

The term $`\alpha _i`$ stands for the dissipation of energy following from the biological activity of the members of species $`i`$ (movement, extraction of nutrients, basal metabolism), as well as the death rate of individuals, and corresponds to the quantities $`d`$ and $`p_d`$ of the IBM. The term $`\beta _iN_i(t)`$ is known as self-damping. It expresses a negative feedback of $`N_i`$ on its own growth rate, which has been shown in some circumstances to be necessary in order to stabilize the model. The terms $`\gamma _{ij}N_j`$ represent the biomass transfered per unit time from species $`j`$ to species $`i`$ if the sign is positive, and from species $`i`$ to species $`j`$ if it is negative, thus modelling prey-predator interactions. They are the counterpart of the matrix $`C(m_i,m_j)`$ in the IBM.

Energy flows into the system through the coupling of basal species to external resources, which are formally represented as an additional “species” $`N_0`$ whose equation will be specified below. Terms of the form $`\gamma _{i0}N_0`$ are thus equivalent to the parameter $`b`$ in the IBM. However modelled, external resources introduce in the system an energy scale $`R`$ which limits its total biomass.

The quantity $`g_{ij}\gamma _{ij}N_j`$, equal to the energy transfer from prey species $`j`$ to predator species $`i`$ per unit of predator, is called the predator’s functional response to prey $`j`$. We studied two different variants of the continuous model, with different functional responses and different equations for the resources.

* Model A. Generalized Lotka-Volterra equations with constant $`\gamma _{ij}`$’s randomly drawn from a uniform distribution, $`\gamma _{ij}[0,\gamma _{\mathrm{max}}]`$.

In order to represent competition among basal species, we introduce a fictitious dynamics for the resources $`N_0(t)`$, modelling them through an equation of the same kind as (3),

$$\frac{1}{N_0}\frac{dN_0}{dt}=\gamma _0R+\underset{j}{}\gamma _{0j}N_j,$$
(4)
where the constants $`\gamma _{0j}`$ have all negative signs and we assume that at least one basal species is present.

The predator’s functional response $`g_{ij}=\gamma _{ij}N_j`$ is proportional to prey biomass and belongs to the more general category of prey dependent functional responses.

There is to observe that in Lotka-Volterra equations, the quantity $`\alpha /\gamma _{\mathrm{max}}`$ introduces an energy scale in the ecosystem, aside to the other energy scales $`R`$, $`N_c`$ and $`\alpha /\beta `$. One can then expect that different regimes are present for different relative values of these energy scales, and this is indeed what we observe. Not all these regimes are biologically meaningful. For instance, there is a regime, corresponding to small values of the dimensionless parameter $`u=\gamma _{\mathrm{max}}N_c/\alpha `$, where dissipation dominates and only basal species can survive in the long run. We shall describe shortly in the Appendix this garden regime. More details on the regimes of Lotka-Volterra equations will be given in a forthcoming work.
* Model B. Ratio dependent functional response.

Another possibility is that the energy scale of $`1/\gamma _{ij}`$ is determined by the biomasses $`N_i`$ and $`N_j`$. This choice is known as ratio dependent predator response (Arditi & Ginzburg, 1989), since the functional response $`\gamma _{ij}N_j`$ depends on the ratio between the prey biomass and the predator biomass. In this case the $`\gamma _{ij}`$’s are not anymore constant, but they are inversely proportional to some linear combination of $`N_i`$ and $`N_j`$. The simplest possibility is that the prey $`j`$ has a unique predator $`i`$, in which case the functional response is given by

$$g_{ij}(N_i,N_j)\gamma _{ij}(N_i,N_j)N_j=\left(a+bN_i/N_j\right)^1.$$
(5)
In the case of several predators for the prey $`j`$, different generalizations of Eq. (5) have been proposed (Arditi & Michalski, 1995; Schreiber & Gutierrez, 1997; Drossel et al., 2000). We adopt our own generalization, which reads

$$\gamma _{ij}=\frac{\alpha _{ij}^{}c_{ij}}{b_jN_j+_{kP(j)}c_{kj}N_k},$$
(6)
where $`i`$ is the predator and $`j`$ the prey. As in model (A), we assume $`\gamma _{ij}=\gamma _{ji}`$. Here $`c_{ij}`$ and $`b_j`$ are dimensionless coefficients, $`\alpha _{ij}^{}`$ expresses the rate at which a single individual of species $`i`$, in absence of competition, consumes a corresponding quantity of biomass from species $`j`$, and $`P(j)`$ indicates the set of predators of species $`j`$.

The above equations, unlike model (A), explicitly represent the competition among predators of the same prey. It is then possible to model external resouces as a constant flux of energy available to basal species,

$$N_0(t)R.$$
(7)

Ratio dependent and prey dependent functional responses have been supported and criticized in several papers (see Abrams & Ginzburg, 2000 for a recent review). We do not want to enter such a debate here. Any functional response is just a crude representation of a much more complicated situation, in which spatial distributions of individuals, foraging strategies and mating behaviour are involved. The point raised in this paper is that, although model (A) and model (B) may have different scaling properties, the scaling behaviour of biodiversity is robust with respect to changes in the functional response, in an appropriate range of parameters. Indeed, we observe that model (B) gives results qualitatively similar to those of model (A) in a range of $`u=\gamma _{\mathrm{max}}N_c/\alpha `$ where $`S`$ scales as $`\mathrm{log}(R/N_c)`$.

### B Immigration and ecological parameters

At time $`t=0`$ no species is present on the island. New species arrive one after another, at fixed intervals of time, $`T_{\mathrm{mig}}`$. Between successive arrivals, population dynamics equations are integrated and species may go extinct.

For every new species, the ecological parameters are chosen at random and kept fixed until the species goes extinct. This means that new species are not related to species already on the island, that is, the continental pool is considered infinite respect to the number of species on the island.

New species have no predators on the island, and a number of preys $`\mathrm{}`$ randomly extracted between one and $`\mathrm{}_{\mathrm{max}}`$ (in most simulations we used either $`\mathrm{}_{\mathrm{max}}=4`$ or $`\mathrm{}_{\mathrm{max}}=8`$). The $`\mathrm{}`$ preys are extracted with uniform probability among the $`S(t)`$ existing species, regarding the external resources $`N_0`$ as a normal prey. This operation defines the ecological network.

For every link, the interaction strengths $`\gamma _{ij}`$ are extracted from an uniform distribution in $`[0,\gamma _{\mathrm{max}}]`$ in the case of model (A). In the case of model (B), the parameters $`c_{ij}`$ are extracted uniformly in $`[0,c_{\mathrm{max}}]`$ and the remaining parameters are fixed at $`\alpha _{ij}^{}1`$, $`b_j1`$. In both cases, $`i`$ is the predator and $`j`$ is the prey, and we then make the assumption that the interaction strengths are antisymmetric: $`\gamma _{ji}=\gamma _{ij}`$. We also studied the case of reduced efficiency, $`\gamma _{ij}=\eta \gamma _{ji}`$, with $`i`$ predator, $`j`$ prey and $`0.5<\eta <1`$, without observing any qualitative difference. In case of model (A), the parameter $`\gamma _0`$, proportional to the growth rate of the resources $`N_0`$, is set to $`\gamma _0=\gamma _{\mathrm{max}}`$. As a simplification, the dissipation parameters are the same for all species $`\alpha _i=\alpha `$, $`\beta _i=\beta `$.

Colonizing species arrive with very small populations. This assures that they are rapidly eliminated if they are not fit for the island ecosystem. We used initial values $`N_{\mathrm{ini}}=N_c`$ and $`4N_c`$. Reducing the initial size reduces spurious effects due to species with very short permanence in the system, and has only a very small influence on the statistical patterns described later.

### C Discussion of the modelling choices

The choice that new species have preys but not predators on the island has to be justified. This rule is aimed at forbidding the formation of ecological loops. Moreover, we want that newcomers have considerable chances of surviving, otherwise the rate of arrival of species with non vanishing permanence time would have large fluctuations, increasing the fluctuations of all ecological variables. With this rule, the resulting ecological networks resemble very much those obtained with the application of the cascade algorithm, used in the IBM model.

Our simulations considered several representations of the ecological dynamics at the individual and at the population level. Although our consistent results let us believe that the models capture generic properties of ecological networks, there are a number of alternative (and equally plausible) dynamical rules or generalizations of our rules that would be interesting to consider. We shortly discuss some of them. This is of course a strong simplification. On the other hand, extracting at random for each species parameters $`\alpha _i`$ and $`\beta _i`$ would lead the average value of these parameters to decrease towards zero, due to the advantage conferred by small dissipation rates and self-damping. In order to avoid this effect, we believe that it is necessary to model a trade-off between dissipation, $`\alpha _i`$, and predation efficiency $`\gamma _{ij}`$, in such a way that the latter is an increasing function of the former. Another critical point is that all pairs of coefficients $`\gamma _{ij}`$ and $`\gamma _{ji}`$ have, in our model, opposite sign, so that symbiotic relationships are not represented.

Concerning the network structure, it is certainly a big simplification to build links independently one of each other. For a more realistic model one needs some measure of the distance among species in some multidimensional space. It would then be possible, for instance, to extract at random the first prey of a new species, and then to extract the remaining preys with a probability depending on the distance from the first prey. Another simplification, in the framework of the continuous model, is that a new species has only preys and not predators. This was imposed in order to avoid the formation of ecological loops, and is very similar to the procedure adopted in the cascade model for network formation. It is thus conforting that ecological networks constructed both using the cascade model and random matrices containing loops (in the IBM framework) returned the same qualitative behaviour.

For very long times, coevolution of species would become relevant. This might be taking into account by making appropriate modifications in the ecological network or in the interaction parameters. Discarding coevolution is justified if the time scale of the simulation is much shorter than the time after which at least one species in the ecosystem mutates. Nevertheless, the latter time scale is expected to decrease as the number of species in the ecosystem, $`S(t)`$, increases, until a point where it is not anymore possible to neglect coevolution. Such a situation is worth considering and will be the subject of future work.

## IV Statistical features of the stationary state

For every set of observations, we shall present both results obtained with the continuous model and results obtained with the IBM, when available. In fact the two descriptions produce the same qualitative behavior, even if simulations are much faster for the continuous model, so that it is possible to simulate larger systems and to obtain better statistics.

When we extract at random an ecological network with a high number of species (up to 1000) and a low number of links per species (for instance four), the ecological dynamics leads to the extinction of most of the species, until only very few ones are represented in the system. This result does not seem to depend on the way in which the links have been extracted (either using random matrices or through the cascade model), on the parameter values, or on the kind of ecological dynamics represented (individual based or continuous, with constant $`\gamma _{ij}`$ or with ratio dependent response functions).

In presence of a constant flux of immigrant species, the ecosystem, initially empty, increases very fast in diversity until it reaches a number of species which remains on the average stationary in time, although characterized by large fluctuations. This process is illustrated in Fig. 1. The system differs significantly from a static network: In fact, if we stop immigrations we notice an abrupt decrease in the number of species until a fixed point of much lower diversity is reached (see Fig. 1).

Choosing as time unity the quantity $`1/\alpha `$ and as biomass unity the external resources $`R`$, the dynamical equations can be written as a function of four dimensionless parameters. For model (A) they are:

$$\gamma _{\mathrm{max}}^{}=\gamma _{\mathrm{max}}R/\alpha ,N_c^{}=R/N_c,\beta ^{}=\beta R/\alpha ,T_{\mathrm{mig}}^{}=\alpha T_{\mathrm{mig}}.$$
(8)
For model (B) the first parameter is substituted by

$$c_{\mathrm{max}}^{}=c_{\mathrm{max}}/\alpha .$$
(9)
Together with $`\mathrm{}_{\mathrm{max}}`$, the maximal number of preys for a colonizing species, these parameters determine the model ecosystem.

We observe turnover of species in the system, and even a complete change in the species composition of the island in the course of time. In the IBM model, where we use a fixed continental pool of species, the presence of different basal species determines the intermediate and top species allowed by the subnetwork on the island. Due to stochastic effects, we observe an alternation between different basal species (often after a long time interval), and consequently a complete renewal of the island ecosystem (see Fig. 2). This picture agrees qualitatively with experiments on island repopulation (Simberloff & Wilson, 1969; Heatwole & Levins, 1972; Simberloff, 1969), where, after defaunation, a different specific composition was obtained.

### A Species diversity and relation among time-scales

When a stationary state is achieved, we observe that, for some range of parameters, the average number of species $`S`$ increases as a power law of the immigration rate $`I=1/T_{\mathrm{mig}}`$,

$$SS_0+cI^a,$$
(10)
with a constant $`S_0`$ usually very small.

We define a new exponent $`b`$, that we call competition exponent, as

$$b=\frac{1}{a}1.$$
(11)
In the stationary state, the average extinction rate per species is $`e=1/(T_{\mathrm{mig}}S)`$ and increases with the number of species at equilibrium as $`S^b`$ (for $`SS_0`$ when only the immigration rate increases), whence the name of competition exponent:

$$eS^b.$$
(12)
Since the exponent $`b`$ is larger than zero, the larger the number of species, the smaller the time scale for the extinction of a single species. The fact that the extinction rate increases with the number of species has been postulated in the theory of island biogeography. However, we find this result not as a phenomenological law, but as a generic feature of the dynamics of random ecological networks. We shall first present results from the IBM and then compare them with those from the continuous models, Eq. (3).

The results of the IBM model are summarized in Fig. 3. Four curves of the average number of species $`S`$ as a function of the immigration rate $`I`$ are plotted. Curves qualitatively similar to ours were obtained in other simple models for island colonization (Rosenzweig, 1995; Loreau & Mouquet, 1999), but the functional relationship between $`I`$ and $`S`$ was not investigated.

All curves in Fig. 3 show two plateaus corresponding to a i) low diversity regime (small immigration rate), and ii) disordered species composition (large immigration rate), which are linked by a transition region with power-law shape described by Eq. (10). The effective exponent obtained from a power-law fit is $`a0.75`$. This intermediate regime is the most interesting one, since here both the immigration rate and the ecological organization play a relevant role in setting the average number of species at equilibrium. In the lower plateau, the dynamics is dominated by the internal ecological processes, while in the higher plateau the fast arrival of new individuals controls the species composition on the island. This latter situation is analogous to the one observed in the defaunation experiment reported in (Simberloff & Wilson, 1969; Simberloff, 1969).

The results for the continuous model are completely similar. The regime between the two plateaus is represented in Fig. 4. Notice that in this case we do not have a fixed continental pool, or, in other words, the parameter $`M`$ has to be interpreted as infinite. Some of the curves show a curvature in log-log plots, which can be eliminated introducing $`S_0`$ as additional fit parameter and plotting ($`SS_0`$) versus $`I`$. The effect of $`S_0`$ is thus to reduce the effective exponent $`a`$ as the immigration flux is reduced.

The observed exponents range from $`a=0.42`$ to $`a=1`$. Interestingly, the case with $`a=1`$, corresponding to a competition exponent $`b=0`$ (not represented in Fig. 4), refers to a case where basal species were not in competition, since we used constant parameters $`\gamma _{ij}`$ and constant resources $`N_0(t)R`$. In all other cases the exponent $`b`$ was positive.

We now discuss shortly the behavior of biodiversity with the parameters of the ecological equations. Keeping fixed the other dimensionless parameters given in Eqs. (8) and (9), biodiversity increases logarithmically with the resources $`R/N_c`$:

$$SA_1+A_2\mathrm{log}\left(R/N_c\right).$$
(13)
This result holds for the IBM and for the continuous models, but for model (A) it is only valid in some range of parameter values. In fact, model (A) at fixed $`\gamma _{max}^{}`$ is found, for large $`R/N_c`$, in the garden regime, where only basal species survive, and the scaling behavior is different there (see Appendix). The exponent $`a`$, defined in Eq. (10), changes only very weakly with $`R/N_c`$.

Biodiversity also increases with the maximal transfer rate, either $`c_{\mathrm{max}}^{}`$ for the case of Eq. (6) or $`\gamma _{\mathrm{max}}^{}`$ for model (A), when all other parameters in Eq. (8) are kept constant. In the first case, the number of species tends to zero as $`c_{\mathrm{max}}^{}`$ approaches unity (for $`N_{\mathrm{ini}}=N_c`$). In the second case, at small $`\gamma _{\mathrm{max}}^{}/N_c^{}`$ we reach the “garden regime” (see Appendix) where the number of species is almost independent of $`T_{\mathrm{mig}}`$. Both limits can be interpreted as corresponding to $`a=0`$, thus $`b=\mathrm{}`$. The exponent $`a`$ then increases slowly with $`\gamma _{\mathrm{max}}`$.

The effect of the parameter $`\beta ^{}`$, when other parameters in Eqs. (8) and (9) are fixed, depends on the immigration rate. While biodiversity increases for increasing $`\beta ^{}`$ at small immigration, the opposite happens if immigration is large. Thus, $`IS`$ curves relative to different values of $`\beta ^{}`$ should cross at some point. This behavior reflects on the fact that the exponent $`a`$ decreases with increasing $`\beta ^{}`$, while $`S_0`$ increases. As in the case of $`\gamma _{\mathrm{max}}^{}`$, the decrease of $`a`$ can be explained by the fact that, at larger $`\beta ^{}`$, the probability that a colonizer has a positive growth rate becomes smaller. The positive effect on $`S_0`$, on the other hand, is due to the fact that the larger $`\beta ^{}`$, the more likely that two predators feeding on the same prey species can coexist.

We also measured the whole stationary distribution of the number of species $`s`$, $`P(s)`$. We plot in Fig. 5a the ratio between the variance $`V_s=s^2s^2`$ and the mean $`S=s`$ as a function of the mean number of species $`S`$. The measure is much noisier than that of the average value, and it is strongly affected from possible sampling errors. The variance is typically lower than the mean for small $`s`$, but at some point increases faster than the mean, becoming larger than it at large $`s`$. Since for a Poissonian distribution the variance and the mean are equal, the distribution $`P(s)`$ is narrower than a Poissonian one for small $`s`$ and broader than it for large $`s`$ (see Fig. 5b, where for each curve a Poissonian distribution with the same mean value has been represented as a dashed line). Notice however that the comparison with a Poissonian distribution is very good for some range of parameters.

### B Species area relationship

In order to investigate the dependence between biodiversity and area, we have to fix the relationship between area and immigration rate $`I`$. Usually, a positive correlation is expected even though the actual dependence may vary with the species considered. We restrict our study to the assumption that the immigration rate is proportional to the linear size of the island (MacArthur & Wilson, 1963):

$$I=1/T_{\mathrm{mig}}=I_0+kA^{1/2}.$$
(14)
We also include a constant $`I_0`$ to take into account that, for islands in an archipelago, the immigration rate depends much more on the geometry of the archipelago and on its distance from the mainland than on the value of the area. For isolated islands and whole archipelagos, on the other hand, there is usually a unique source of immigrants from the mainland, and the effect of area on immigration is expected to be important. In view of this situation, in the first part of the discussion the constant $`I_0`$ will be neglected. The parameter $`k`$ is related to the distance from the mainland.

In the framework of the IBM, the other parameter influenced by area is the number of patches $`N_h`$, which is taken to be proportional to area: $`N_h=A`$, with an appropriate choice of units. Thus we simulated systems with different values of $`N_h`$, varying the immigration rate as above. Our main result is that we always obtain a power-law dependence of the number of species with the area, with typical values of $`z`$ in the interval $`0.60.8`$, as it is observed for the case of archipelagos, to which our immigration model should apply. Taking into account additional sources of immigrations, like closeby islands, is expected to reduce the dependence of the immigration rate on area, and thus to cause a decrease of the effective $`z`$ values, making them more similar to the values observed in groups of neighboring islands.

For the IBM, we represent in Fig. 6 the number of species as a function of area, with immigration rate $`Ik\sqrt{A}`$ and $`N_hA`$, for different values of the continental pool $`M`$. It can be seen that the species area relationship bends for large areas, apparently tending to an asymptote. Increasing the species pool $`M`$ increases the asymptote, but leaves the value of the effective exponent unchanged. Parameter values are $`E=150`$, $`E_{rep}=200`$, $`d=1`$, $`\delta =20`$, $`p_d=0.002`$ and $`k=0.353`$.

We now come to the species-area relationship in the continuous formulation of the ecological dynamics. We should first discuss how the parameters of the ecological equations, Eqs. (8) and (9), depend on area. The variables $`N_i`$ of ecological equations have the meaning of spatial densities of individuals. Thus the equations Eq. (3) are invariant with respect to changes of area. There is however another important determinant of the ecological dynamics: The threshold $`N_c`$ below which extinctions happen. Two different cases have to be considered:

1. $`N_c`$ independent of area, in other words there is a critical density below which the species go extinct, as it happens for the Allee’s effect (Allee et al., 1949). Such a situation is expected, for instance, if the individual of the species are uniformly dispersed in the area $`A`$ so that, below the critical density $`N_c`$, they can not find mating partners.

In this case, the ecological dynamics is invariant under changes of the area and, in particular, the extinction rate does not depend on $`A`$. Assuming that the immigration rate increases with area as in Eq. (14), and that the area is much larger than $`(I_0/k)^2`$, we find

$$SA^z,z=a/2,$$
(15)
where $`a`$ is the exponent in Eq. (10).
2. Extinctions depend on the absolute number of individuals, $`N_iA`$ and the extinction threshold is thus

$$N_c1/A.$$
(16)
1. First we consider this case together with $`I=I_0`$ independent of area (corresponding to $`k=0`$ in Eq. (14)). From equation (13), we obtain then that the number of species increases logarithmically with area:

    $$SA_1+A_2\mathrm{log}(A).$$
    (17)
    This relation is indeed observed for birds in the central islands of the Solomon archipelago, which are all very close to at least one other large central islands (Diamond & Mayr, 1976). For such islands, the immigration rate can be expected to be rather independent of area.
    2. Extinctions depend on the absolute number of individuals and the immigration rate increases with island size $`(A(I_0/k)^2)`$,

    $$N_c1/A,I=kA^{1/2}.$$
    (18)
    In this case, we can not rely on previous results, and we have to perform new simulations, scaling the parameters as in Eq. (18).

We note that, in both cases 2(a) and 2(b), model (A) can become problematic at large area. In fact, as area increases the coupling constant $`\gamma _{\mathrm{max}}N_c/\alpha 1/A`$ becomes smaller, and the system becomes dominated by dissipative effects. The result is that, unless $`\alpha =0`$, Eq. (29) in the Appendix would not be satisfied as area increases, and the system would meet a “garden regime”, in which only basal species survive (see Appendix). It is not surprising that Lotka-Volterra equations go in trouble for very low densities: in fact, they are analogous to equations of chemical kynetics, and when the density of the species involved becomes too small, the assumptions on which they are founded break down. To avoid such a problem, we used the ratio dependent model (B) together with Eq. (18). Nevertheless, our numerical study shows that also model (A) provides comparable results in a suitable range of parameters.

Results are plotted in Fig. 7, and yield an approximately power law species area relationship, for large enough values of the parameter $`\alpha k`$. In the log-log plot, the curves $`S(A)`$ show a negative curvature, that can be eliminated introducing a new parameter $`A_0`$,

$$Sc(k)\left(AA_0(k)\right)^{z(k)}.$$
(19)
Both the effective exponent $`z(k)`$ and the limit area $`A_0(k)`$ increase slowly with the immigration parameter $`k`$. For the curves in Fig. 7, the exponent $`z(k)`$ ranges from $`z=0.49`$ at $`k=2`$ to $`z=0.56`$ at $`k=50`$.

We notice however that the scaling form Eq. (19) is only approximate, that a scaling of the form $`A^z\mathrm{log}(A)`$ would probably be more adequate, and that the immigration rate is a better scaling variable than area. We extract an exponent from our approximately power law species area relationship only for the sake of comparison with field data.

Our model of immigrations applies to whole archipelagos or to isolated islands, because we consider a unique source of immigrations. It is remarkable that the exponent $`z`$ observed for a set of Pacific archipelagos has the value $`z=0.54`$ (Adler, 1992; Rosenzweig, 1995), in very good agreement with the results of our simulations. Real data are shown for comparison in Fig. 7b.

It is striking that models at different description levels, as the IBM and the continuous model (B), yield very similar species area relationships. We compared other statistical features of the individual based and the continuous models, finding that they are qualitatively very similar (this holds for the Lotka-Volterra model as long as the garden regime and the opposite low dissipation regime are avoided). In the following we will present a complete analysis of the stationary state for data obtained from the IBM and from model (B).

### C Distribution of species abundances

We measured distributions of species abundances, defined as the probability density $`p(N)`$ of species with $`N`$ individuals (or total biomass equal to $`N`$ in the case of the continuous model), both for the IBM and for the continuous model. We observe a good qualitative agreement of the results in both approaches.

In the framework of the IBM, we measured the distribution of species abundances for three values of the immigration rate corresponding to the three regimes in Fig. 3 (slow, intermediate and fast driving). Fig. 8 represents the frequency $`p(N)`$ with which species formed by $`N`$ individuals were recorded. Different curves refer to different immigration rates (other parameters as in Fig. 3). All curves show an initial fast decaying part corresponding to species that go extinct almost immediately after arriving to the island. Since these species do not find preys to feed on, their initial energy decays exponentially and they die out of starvation. The relevant part of these distributions results from species which play a role in the ecological network. This part shows a power-law decay of the form

$$p(N)N^\xi ,$$
(20)
with $`\xi 1`$. Finally, the external resources set the value of $`N`$ at which an exponential cut-off appears. Our results are in good agreement with field measures of diversity, many of which also return a power-law distribution of species abundances with an exponent in the range $`11.25`$ (Pielou, 1969; Solé et al., 2000).

The same results are obtained in the framework of the continuous model. In this case, however, we observe that the exponent $`\xi `$ increases slowly with immigration rate, tending to $`\xi 1`$ in the limit $`I0`$. The maximum value that we found in our simulations is $`\xi 1.25`$, still compatible with observational data. A sample of results is reported in Fig. 9a. In Fig. 9b the decrease of the exponent $`\xi `$ with the immigration rate is shown.

### D Lifetime distribution

The distribution of lifetimes of species is shown in Fig. 10 in a log-log plot, for several values of the immigration rate and of other parameters. After an initial part where the distribution is almost uniform, corresponding to species with very short permanence time, we observe an approximately power law decay of the probability density for a range of at least one and half decade

$$p(\tau )\tau ^\eta ,$$
(21)
when using different parameter values, we found values of the effective exponent $`\eta `$ between 2.1 and 2.8.

The average lifetime $`\tau `$ in the equilibrium state is related to biodiversity through the relation

$$\tau =T_{\mathrm{mig}}ST_{\mathrm{mig}}^{b/(1+b)},$$
(22)
which follows from its definition and from the properties of the stationary state. Here $`T_{\mathrm{mig}}`$ represents the average time between succesful immigrations which contribute to the island biodiversity. Thus the average lifetime decreases with the immigration rate $`1/T_{\mathrm{mig}}`$ and, consistently, the value of the exponent increases, as it is shown in Fig. 10b.

Our results compare qualitatively well with observed patterns (Keitt & Marquet, 1996; Keitt & Stanley, 1998). It was in fact observed that the time of permanence of birds in local patches follow a distribution approximately of the form (21) with effective exponent $`\eta =1.6`$, indeed smaller than the typical values found in our simulations. A result which compares better to this last value has been found, using a model without explicit ecological dynamics, in (Solé et al., 2000).

### E Network organization

The structure of the ecological network does change, even if very slowly, with changing immigration rate. We have examined in particular the number of trophic levels, the number of links per species and the total biomass.

We define the trophic level of a species as the minimal number of links connecting it to resources. In all our simulations the number of trophic levels varies between four and ten. It shows a tendency to increase with immigration rate, as it is illustrated in Fig. 11.

The average number of links per species, counted as average number of preys, is shown in Fig. 12 as a function of the immigration rate. It changes very slowly (logarithmically) and, in some cases, in a non monotonic way (for most curves we only observe either the increasing or the decreasing part). A similar pattern is observed as a function of the resources $`R`$ (see Fig. 12b). Thus, as a function of the number of species, the number of links per species behave non monotonically. It also depends weakly on the maximum number of links allowed when the new species is added to the ecosystem, $`\mathrm{}_{\mathrm{max}}`$.

The total biomass also increases approximately as a power law of the immigration rate, as shown in Fig. 13. The exponent ranges from 0.15 to 0.58.

## V Relationship to MacArthur & Wilson’s theory

We have already seen in the previous section that quantitative biodiversity patterns can be derived from a balance between external driving through immigrations and the intrinsic population dynamics of the ecosystem. Here we want to relate these results to existing phenomenological approaches, in particular MacArthur and Wilson’s (1963, 1967) theory of island biogeography.

In an ensemble average (or, equivalently, in an average over long times), the response of the system to a constant immigration flux $`I`$ reduces to a stationary extinction $`E`$. This can be expressed as a function of the avergage number of species $`S`$ in the system, $`E=E(S)`$. Of course, this function depends also on the model paramenters and on qualitative features of the immigration flux. Since the system is in a stationary state, immigration and extinction have to compensate each other on average:

$$IE(S)=0.$$
(23)
This balance is indeed as postulated by MacArthur and Wilson. The immigration flux measures the average number of new species arriving per unit of time. The functional form of the ‘extinction curve’ $`E(S)`$ can now be obtained from the underlying population dynamics. As explained above, we find

$$E(S)=E_0S^{1+b},$$
(24)
in the scaling regime where the number of species increases as a power-law of the immigration rate. The exponent $`b`$ was called competition exponent and has been introduced in Eqs. (11) and (12). In fact, from the stationary solution of Eq. (23), we recover Eq. (10)

$$S=\left(\frac{I}{E_0}\right)^{\frac{1}{1+b}}.$$
(25)
This equation allows to derive the exponent $`z`$ of the species-area relationship. Indeed, if we assume that

$$IA^s,E_0A^ϵ,$$
(26)
we find

$$SA^z,\text{with}z=\frac{s+ϵ}{1+b}.$$
(27)
For the population models, we assumed $`s=1/2`$ (the immigration rate is considered proportional to the linear size of the island), and we obtained $`ϵ=0`$ (in fact, the number of species increases as the logarithm of $`R/N_c`$ at fixed $`T_{\mathrm{mig}}`$, both in the IBM and in the continuous model (B)).

We remark that here, as in the explicit population models, we are assuming that the only source of immigrants is a continent far apart. The exponents that we find should then be compared to the exponents observed for isolated islands or archipelagos, while the exponent computed among islands of the same archipelago is expected to be lower, due to a reduced dependence of immigration rate on area. Another point to remark is that the immigration rate $`I`$ measures the flux of new species arriving on the island. If this flux is assumed to originate from a ‘continent’ of $`M`$ species, $`I`$ can be related to the total immigration flux $`I_0`$ by correcting for the immigrations already present on the island. The simplest ansatz is $`I(S)=I_0(1S/M)`$ (MacArthur and Wilson, 1963, 1967). Expressed in terms of $`I_0`$, the average number of species $`S`$ reaches a saturation value of order $`M`$. (The pool of immigrant species is indeed finite in our IBM, but infinite in the continuum models. Hence, this correction becomes important in comparing results of the different models).

## VI Summary and Conclusions

We have presented a study of biodiversity in an insular ecosystem at the individual and the population level. Our interest has been focused on the statistical properties of the dynamical stationary state and on the scaling relations between the system variables. Instead of describing detailed situations in which some particular species and their exact interactions with their known preys and predators are included, we let ecological networks self-organize through random assemblage of species, ecological dynamics and possible extinctions.

Our main result is that, in a broad range of parameters, biodiversity scales approximately as a power law of the immigration rate. The value of the exponent varies slightly when the parameters of the models are changed, but the qualitative features of the stationary state are quite robust.

The behavior of biodiversity with immigration rate allows to derive a species area relationship with a power law shape, if we assume that the immigration rate is proportional to the linear size of the island. Such a model of immigration considers as unique source of diversity a flux of species from a continent far apart, thus it should be compared to observations relative to isolated islands or to whole archipelagos. The agreement is in this case rather good: the observed value of the effective exponent $`z`$ on archipelagos is $`z=0.54`$ (Adler, 1992; Rosenzweig, 1995), while we typically get, with the continuous model, values between $`0.52`$ and $`0.56`$ and, with the individual based model, values between $`0.6`$ and $`0.8`$. Thus, the comparison of our two description levels points out to species-area law of the type (1) as a generic feature of a broad set of ecological models with random interactions.

Our models reproduce qualitatively other features observed in real ecosystems. We observe a power law distribution of population abundances, i.e. the number $`p(N)`$ of species with $`N`$ individuals approximately decreases as $`p(N)N^\xi `$. This is expected to be a general consequence of the multiplicative nature of population dynamics equations. The exponent $`\xi `$ found in our simulations is close to unity, in favourable comparison with field data, and increases with the immigration rate.

We also observe a broad distribution of the time of permanence of species in the system $`\tau `$, as it has been observed in the field (Keitt & Marquet, 1996; Keitt & Stanley, 1998) and in a related model (Solé et al., 2000). The average permanence time is proportional to the number of species and inversely proportional to the immigration rate, $`\tau =ST_{\mathrm{mig}}`$, so that it decreases with the immigration rate. The fact that it is observed, both in field studies and in models, that its distribution is broad, could help to reconcile the apparent dichotomy between fugitive species and permanent species: These groups of species could correspond to the two extreme cases of a unique distribution of permanence times (Schoener & Spiller, 1987). The approximate power law shape of the distribution of times of permanence in the island is reminiscent of the analogous distribution of the lifetime of genera in the fossil record, which is approximately given by $`P(\tau )\tau ^\eta `$, with $`\eta 2`$, close to what is observed in our model for very small immigration rates and also close to ecological observations. It is tempting to speculate that this similarity points out at similar mechanisms acting on the time scales of ecosystem dynamics as on the timescales of macroevolution.

The number of trophic levels in the food web is also strongly influenced by the immigration rate. We typically find from four to ten trophic levels, depending on parameters, and with a tendency for the number of levels to increase with immigration rate. Hints to the correlation between immigration rate and number of levels can be found in the fact that the length of food chains appears to be positively correlated to habitat area, although the data are quite poor (Schoener, 1989; Spencer, 1997). Our results suggest that one of the factors limiting the length of food chains is the immigration or speciation rate. Notice that in our model no other limitations to the length of food chains exist: energetic considerations would limit the number of levels to a value $`\mathrm{log}(R/N_c)`$, much larger than the one observed.

An important result of our study is that the observed statistical patterns are rather robust with respect to changes in the dynamical rules of the model. One example is the representation of space in the IBM model. Although one could think that in this case explicit space is needed, we modelled space only in an effective way, increasing the immigration rate and the resources with the area. We believe that this effective approach captures the main features of the behavior of biodiversity with area, even if important issues, like for instance the presence of many different habitats, are not represented in the model.

As Pimm poses it, (…) it is pointless to try to justify models’ equations biologically – their assumptions are almost bound to be wrong. (…) The concern should not be whether the assumptions are wrong (they are!), but whether it matters that they are wrong. (Pimm, 1991). It seems that the statistical laws and the scaling relationships that we observed are generic properties of complex ecosystems, that is an unavoidable results of a minimal set of rules governing population dynamics and immigrations. Thus the strategy is to look for the simplest set of rules which appear sensible and which allow to derive the observed statistical patterns of biodiversity.

## Acknowledgments

Discussions with David Alonso, Lloyd Demetrius, Barbara Drossel, Lorenz Fahse and Martin Rost are gratefully acknowledged.

## Appendix. The garden regime

In the case of model (A) of continuous dynamics, we observed that only basal species could survive on the long run for parameter values in a certain range. In this case, observing the system on a very long time scale, we did not see any stationary state (although a stationary state was reached for a much smaller system), but we saw a number of species steadily and slowly increasing in time (see Fig. 14).

We call such a regime the garden regime, since predators are absent. Its statistical properties are peculiar: the distribution of biomasses is narrow and peaked at very low values, and the distribution of lifetimes is bimodal, with a high peak for very short lifetimes corresponding to non-basal species, and a shorter one for large lifetimes, corresponding to basal species. (Even if at any moment there are many more basal species than predators, the number of predators passing through the system and almost immediately going extinct is much larger than the corresponding number of basal species.) Finally, the distribution of the transfer rate $`\gamma _{0i}`$ from the external resources to species $`i`$ is strongly peaked close to the maximum allowed value $`\gamma _{\mathrm{max}}`$. All these features can be easily rationalized as follows.

Let us consider a situation where $`S`$ basal species coexist feeding on the abiotic resources $`N_0`$. Our calculations (unpublished) show that at the static fixed point of the corresponding ecological equations, Eqs. (3) and (4), all biomasses are positive if and only if all differences in the coefficients $`\gamma _{0i}`$ are smaller than $`B/S`$, where $`B=\gamma _0\beta R/\alpha `$. Thus, as $`S`$ increases, the coupling constants between the basal species and the environment deviate from the initial distribution and become more and more similar. Notice that, if $`\beta =0`$, the coefficients $`\gamma _{0i}`$’s should be exactly equal to guarantee coexistence (this expresses in this context the “principle of competitive exclusion”). This conclusion does not vary qualitatively if one considers Eq. (6) instead of the model with constant coefficients $`\gamma _{ij}`$. Thus, in a system without predators and high immigration of basal species, we would expect to find many basal species with very similar biomasses, all of order $`R/S`$ ($`S`$ should be smaller than $`R/N_c`$) and coefficients $`\gamma _{0i}`$ very close one to each other in value.

Let us now consider the arrival of a predator to such a system, considering first Lotka-Volterra equations (model (A)). The rate of growth of the predator cannot be larger than $`r_{\mathrm{max}}`$ given by

$$r_{\mathrm{max}}=\mathrm{}_{\mathrm{max}}\gamma _{\mathrm{max}}\frac{R}{S}\alpha \beta N_c.$$
(28)
Since $`R/S`$ is larger than $`N_c`$, we find that, for large $`S`$, a non-basal species can colonize only if

$$\gamma _{\mathrm{max}}\frac{1}{\mathrm{}_{\mathrm{max}}}\left(\frac{\alpha }{N_c}+\beta \right).$$
(29)
Every time we observed in our simulations ecosystems mainly composed of basal species and whose biodiversity has a time behavior similar to that in Fig. 14, the above condition was not fulfilled.

In model (B), Eq. (6), $`r_{\mathrm{max}}`$ equals $`\mathrm{}_{\mathrm{max}}\alpha ^{}c_{\mathrm{max}}\alpha \beta N_c`$. This quantity must be larger than zero, otherwise no species would survive. Thus the “garden regime” that we described in this appendix can not be found in simulations of ratio dependent functional responses, unless we decide to use different parameters for basal and non-basal species.


format: text

text:
# Diffraction-limited bispectrum speckle interferometry and speckle polarimetry of the young bipolar outflow source S140 IRS1

## 1 Introduction

The formation process of massive stars is still not well understood (cf. Garay & Lizano 1999 for a recent review). The evolutionary timescales for massive stars are very short, as they start to burn hydrogen while still accreting material from the surrounding protostellar cloud. Young massive stars also affect their environment strongly by driving very powerful winds, molecular outflows, and jets, and also by emitting intense UV radiation. An especially interesting aspect of massive star formation is the observational result that many massive protostars drive very energetic outflows that are often more massive than the central protostar (cf. Churchwell 1997). To reach a better understanding of the physics of massive star formation, observations with high spatial resolution are crucial in order to disentangle the numerous different physical processes taking place simultaneously. This is the motivation for our high-resolution study of the massive protostellar object S140 IRS1, well known for driving a massive bipolar molecular outflow.

S140 is an HII region at the south-east edge of the L1204 dark cloud and part of a cloud complex located at the edge of a prominent infrared emission ring, known as the Cepheus ring. This ring is probably the result of a supernova explosion and stellar winds from massive stars in the open cluster NGC 7160, close to the center of the ring (cf. Kun et al. 1987.)

About $`1^{}`$ north-east of the S140 HII region, Rouan et al. (1977) detected strong far-infrared emission at a position at which no objects could be seen in visible light. Harvey et al. (1978) found that the spectral energy distribution of this infrared source, called S140 IRS, strongly increases between 2 $`\mu `$m and 100 $`\mu `$m and estimated an infrared luminosity of $`2\times 10^4L_{}`$ (for a presumed distance of 1 kpc). Beichman et al. (1979) carried out 20 $`\mu `$m observations of this region and were able to resolve the infrared emission into three individual sources. At 20 $`\mu `$m, the dominant source IRS1 is 7 – 10 times brighter than the two other sources IRS2 and IRS3. The luminosity of IRS1 was estimated to be $`L5\times 10^3L_{}`$ (Lester et al. 1986), suggesting it to be a deeply embedded ($`A_V30`$ mag; Harker et al. 1997) early B-type star with a mass of about $`10M_{}`$. However, these estimates of the stellar parameters appear to be quite uncertain. A determination of the extinction towards IRS1 based on 3 $`\mu `$m ice-band spectroscopy by Brooke et al. (1996) yielded a column density of $`N_\mathrm{H}=1.4\times 10^{23}\mathrm{cm}^2`$ corresponding to an extinction of $`A_V70`$ mag. This would suggest a somewhat higher luminosity and thus also a higher mass for IRS1.

A strong molecular CO outflow in the S140 IRS region was first detected by Blair et al. (1978). The total outflow mass was estimated to be $`64M_{}`$ (Bally & Lada 1983), i.e. about 6 times more than the mass of the central protostar. Minchin et al. (1993) studied the CO molecular line emission in the S140 region with the JCMT and found a bipolar outflow morphology. S140 IRS1 lies just in the middle between the blue- and red-shifted outflow lobes; since the other two infrared sources IRS2 and IRS3 are clearly not located on the outflow axis, IRS1 can reliably be assumed to be the source of this outflow. The position angle of the outflow axis is about $`160\mathrm{°}`$. Since the blue- and red-shifted outflow lobes overlap strongly, the outflow axis is believed to be rather close to the line of sight.

In previous studies, near-infrared images of S140 were obtained for example by Harker et al. (1997) and Yao et al. (1998). Near-infrared polarization maps were reported by Lenzen (1987), Whitney et al. (1997), and Yao et al. (1998). These observations, however, were seeing-limited and therefore did not have the resolution required to study the inner environment ($`1000`$ AU) of the central source S140 IRS1. In order to get a better insight into the nature of this interesting object, we have carried out bispectrum speckle interferometry and speckle polarimetry of S140 IRS1.

## 2 Observations

The speckle interferograms were obtained with the 6 m telescope at the Special Astrophysical Observatory (SAO) in Russia on 19 September 1999. The data were recorded through a $`K^{}`$-filter with a central wavelength of 2.165 $`\mu `$m and a bandwidth of 0.328 $`\mu `$m. 5552 speckle interferograms of S140 IRS1 and 5736 speckle interferograms of the reference stars HIP 110410 and HIP 110498 were taken with our Hawaii array speckle camera. The exposure time per frame was 150 ms, the pixel size was 27.0 mas and seeing was $`1^{\prime \prime }`$. The images were reconstructed using the bispectrum speckle interferometry method (Weigelt 1977; Lohmann et al. 1983; Weigelt et al. 1991). The object power spectrum was determined with the speckle interferometry method (Labeyrie 1970). The bispectrum of each frame consisted of 100 million elements. The resulting image (Fig. 1, upper right) has a diffraction-limited resolution of 76 mas. In order to increase the SNR we also reconstructed an image with a reduced resolution of 130 mas (60% of the diffraction limit), which is shown in Fig. 1, upper left.

The speckle polarimetry observations were carried out on 13 June 1998 using the SAO 6 m telescope, our two-beam polarimeter, and a NICMOS 3 camera. The data recording parameters were as follows: 2800 polarized object speckle frames and 2800 reference star (HIP 109981) frames recorded through a $`2.11\mu `$m-filter with central wavelength / FWHM bandwidth of 2.110 $`\mu `$m / 0.192 $`\mu `$m, 2400 polarized object speckle frames and 2400 reference star (HIP 110498) frames recorded through a $`K^{}`$-filter with central wavelength/FWHM bandwidth of 2.165 $`\mu `$m / 0.328 $`\mu `$m, exposure time per frame 150 ms, pixel size 30.1 mas, field of view $`7.7^{\prime \prime }\times 7.7^{\prime \prime }`$ (256 $`\times `$ 256 pixels), and seeing $`1^{\prime \prime }`$. Our 2-beam speckle polarimeter has the following standard design: achromatic K-band collimator objective to obtain a parallel beam behind the f/4 primary focus of the SAO 6 m telescope, rotatable $`\lambda /2`$ plate for rotating the polarization direction, Wollaston prism to obtain two polarized images separated by $`2.9^{\prime \prime }`$ on the detector, and a second achromatic K-band objective behind the $`\lambda /2`$ plate and Wollaston prism (more technical details of the speckle polarimeter will be described elsewhere). Polarized speckle interferograms were recorded at the four $`\lambda /2`$ plate positions $`0\mathrm{°}`$, $`22.5\mathrm{°}`$, $`45\mathrm{°}`$, and $`67.5\mathrm{°}`$. For each of these four data sets and the corresponding reference star data sets, polarized images were reconstructed using the bispectrum speckle interferometry method. In each of the four reconstructed images two images of the object were obtained which were polarized parallel and perpendicular to the separation vector of the two images. From these polarized images polarization maps were reconstructed using the standard technique (see e.g. Whitney et al. 1997 and references therein). In order to improve the SNR of the images, the resolution was reduced to 280 mas. In Fig. 1 (lower left) polarization vectors are plotted at the positions where the SNR of the degree of polarization (derived by splitting the data sets) is better than 5 and the error of the polarization vectors is smaller than $`10\mathrm{°}`$.

## 3 Results

### 3.1 Morphology in the $`K`$-band image

In Fig. 1 we present our reconstructed $`K`$-band images of S140 IRS1. Our images show a wealth of previously unseen details. The central source IRS1 has a point-like appearance. A series of diffraction rings are clearly visible around IRS1. While the first ring is perfectly closed, the following diffraction rings are broken. A detailed investigation of the data indicates that the central source is marginally resolved with a FWHM Gaussian diameter of approximately 20 mas (corresponding to a physical size of $`20`$ AU at the assumed distance of 900 pc) and is probably elongated in a north-south direction.

Our images show a bright extended feature of diffuse emission pointing from IRS1 towards the south-east. This emission has a very clumpy structure with several prominent knots as indicated in our contour image. Clump A is located on the axis defined by the general shape of the emission structure. Three other clumps, denoted as B, C, and D, can be seen to the right or left of the axis. Additionally, a rather extended feature of emission, denoted as G, is seen south-east of IRS1. Nevertheless, the general shape of the diffuse emission follows a well defined direction. The line connecting the central maximum with clump A has a position angle of $`154\mathrm{°}\pm 3\mathrm{°}`$ (measured counter-clock wise from north). The two features denoted as E and F seem to trace a structure that is oriented roughly perpendicular to the outflow direction.

### 3.2 Polarization maps

Our polarization maps show that the light from the elongated feature is strongly polarized, with degrees of polarization ranging up to 50%. Our two polarization maps, which were taken in slightly different filters (see above), agree very well. If we draw lines perpendicular to each polarization vector in our maps, nearly all these lines cross within less than $`0.3^{\prime \prime }`$ of the central intensity peak. Thus, our maps are very consistent with a centro-symmetric polarization pattern around the central source. Comparison of our polarization maps with theoretical simulated polarization maps (e.g. Fischer et al. 1996) suggests that the light from the elongated feature is scattered light, originating from the central source.

## 4 Interpretation of the observed structures

The near-infrared emission is dominated by the central intensity peak, that contains $`35\%`$ of the total flux in our $`5^{\prime \prime }\times 5^{\prime \prime }`$ image. This peak probably marks the location of the central protostar S140 IRS1. An alternative explanation might be that the protostar is too deeply embedded to be directly visible in our image, and the intensity peak might be light from the protostar that is scattered at the walls of the outflow cavity. However, the observed polarization pattern which is centro-symmetric around the intensity peak strongly suggests that the location of the illuminating source actually is at the intensity peak of our image. The observed ratio of diffuse versus direct light together with typical optical properties of circumstellar dust (e.g. Preibisch et al. 1993) suggests an extinction of $`A_K4`$ mag ($`A_V35`$ mag) along the line-of-sight to the protostar, i.e., the dominant intensity peak seems to be a highly extincted image of the central protostar. The possibly elongated structure of the central source might be caused by bright emission from the innermost part of the larger extended structure, or might indicate that the central object is a binary star.

The most striking feature in our image is the bright extended and very clumpy structure pointing from IRS1 to the south-east. Interestingly, the position angle of this structure of $`154\mathrm{°}\pm 3\mathrm{°}`$ matches the direction of the blue-shifted CO outflow very well: the 12CO map of Minchin et al. (1993) gives a position angle of $`155\mathrm{°}\pm 15\mathrm{°}`$ for the blue-shifted outflow lobe. This strongly suggests that the elongated feature in our image is related to the blue-shifted CO outflow. The axis of the CO outflow is known to be closer to the line-of-sight than to the plane of the sky (Minchin et al. 1993), i.e. the south-eastern lobe of the outflow is oriented roughly towards us, while the north-western lobe is pointing away from us. This suggests the following interpretation of the observed features:

The central object IRS1 is deeply embedded in a dense circumstellar envelope or perhaps a thick circumstellar disk. The outflow has cleared a cavity in the circumstellar material and what we see as the bright extended structure south-east of IRS1, is light from the central protostar that is scattered by dust grains at the inner wall of this outflow cavity into our direction. The very clumpy appearance of the emission suggests that the surface of the cavity wall is not smooth but highly structured. This might be the signature of the violent interaction between the outflowing material and the circumstellar envelope.

Our image traces the diffuse emission out to at least about $`2^{\prime \prime }`$ away from the central source. This suggests that the size of the circumstellar envelope or disk around IRS1 is at least $`2000`$ AU. The asymmetric general shape of our image, i.e. the fact that we do not see a counter lobe north-west of IRS1, is easily explained as a geometrical effect. The red-shifted outflow component, pointing in the north-west direction, has probably cleared a similar cavity. That cavity, however, is not visible in our near-infrared images because it is pointing away from us and the light predominantly escapes in a direction away from us.

## 5 Conclusions

Our high-resolution data provide new insight into the circumstellar environment of S140 IRS1 and enable us to study the structure of the envelope around the central protostar at scales of $`100`$ AU for the first time. Our images reveal bright emission that most probably is reflected light from the inner walls of a partially evacuated cavity in the circumstellar material around S140 IRS1. Due to the high spatial resolution of our data we even can see details in the structure of the cavity walls, which appear to have a highly inhomogeneous and clumpy surface. The fact that the orientation of the evacuated cavity perfectly agrees with the direction of the molecular outflow strongly suggests that this cavity has been carved out by the strong molecular outflow from S140 IRS1. This demonstrates the close relationship between the structure of the dense inner circumstellar material around the central protostar IRS1 at scales of a few hundred AU, and the large scale molecular outflow at scales of several 10 000 AU.

###### Acknowledgements.

We thank the referee for her helpful comments.