Analytical solution of a model for complex food webs

We investigate numerically and analytically a recently proposed model for food webs [Nature {\bf 404}, 180 (2000)] in the limit of large web sizes and sparse interaction matrices. We obtain analytical expressions for several quantities with ecological interest, in particular the probability distributions for the number of prey and the number of predators. We find that these distributions have fast-decaying exponential and Gaussian tails, respectively. We also find that our analytical expressions are robust to changes in the details of the model.

In ecosystems, species are connected through intricate trophic relationships [1,2] defining complex networks [3,4,5,6], the so-called food webs. Understanding the structure and mechanisms underlying the formation of these complex webs is of great importance in ecology [7,8,9,10]. In particular, food web structure provides insights into the behavior of ecosystems under perturbations such as the introduction of new species or the extinction of existing species. The nonlinear response of the elements composing the network leads to possibly catastrophic effects for even small perturbations [11].
Recently, Williams and Martinez have proposed an elegant model of food webs -the "niche" model-that just with a few ingredients successfully predicts key structural properties of the most comprehensive food webs in the literature [1]. Numerical simulations of the niche model predict values for many quantities typically used to characterize empirical food webs that are in agreement with measured values for seven webs in a variety of environments, including freshwater habitats, marine-freshwater interfaces and terrestrial environments.
We investigate the niche model from a theoretical perspective. We study analytically and numerically the behavior of key quantities for sparse food webs, i.e., webs with L ≪ S 2 , where L is the number of trophic interactions between species and S is the number of species in the web. This is the limit of interest in ecology because (i) for most food webs reported in the literature the directed connectance, defined as C = L/S 2 takes small values, and (ii) it corresponds to the limit of large web sizes S [9,10]. We calculate the probability distributions of number of prey and of number of predators and find that for C ≪ 1 they depend only on one parameter of the model-the average number z of trophic links in the network. These distributions give valuable information about the structure of the network [12] and enable us to calculate other interesting quantities such as the fraction of "top," "intermediate," and "basal" species, and the standard deviation of the "vulnerability" and "generality" of the species in the food web [1]. Our results provide compact patterns that describe the structure of the food webs generated by the niche model. These patterns could not have been predicted from the numerical simulations reported in Ref. [1] and may be of practical and fundamental importance for the study of empirical food webs. Moreover, we test our analytical predictions with empirical food webs and find agreement.
Next, we define the model. Consider an ecosystem with S species and L trophic interactions between these species. These species and interactions define a network with S nodes and L directed links. In the model, one first randomly assigns S species to "trophic niches" n i mapped into the interval [0,1]. For convenience, we will assume that the species are ordered according to their niche number, i.e., n 1 < n 2 < ... < n S .
A species i is characterized by its niche parameter n i and by its list of prey. Prey are chosen for all species according to the following rule: a species i preys on all species j with niche parameters n j inside a segment of length r i centered in a position chosen randomly inside the interval [r i /2, n i ], with r i = xn i and 0 ≤ x ≤ 1 a random variable with probability density function The values of parameters b and S determine the average connectivity z ≡ 2L/S of the food web and the directed connectance C = L/S 2 [1,13]. One can also express the average number of prey per species as Sr, where the bar indicates an average over an ensemble of food webs. It then follows that the connectivity is z = 2Sr, the number of directed links is L = S 2 r, and the connectance is C = r. One can also obtain these expressions in terms of b using the equality r = In the niche model, isolated species-that is, species with no prey or predators-are eliminated and species with the same list of prey and predators -that is, trophically-identical species-are "merged" [14].
First, we address the statistics of the number of prey. For large S, the number of prey of a species i is k i = Sr i , so that the probability distribution p prey is given directly by the distribution of r. Specifically, p prey (k) = p(r)/S. The cumulative probability P (r ′ > r) = 1 r dr ′ p(r ′ ) is the area of the region R of the n − x diagram bounded by lines x = 1, n = 1, and the hyperbole r = nx where p n (n) = 1 is the probability density function of n.
The integration of (2) gives rise to a function involving hypergeometric functions [15]. To obtain a more "physical" solution, one can derivate (2) twice to obtain dp(r) dr In the limit C ≪ 1, one has b ≫ 1 (see [14]), so that p x (x) ≃ be −bx , and the term in the right-hand side vanishes exponentially, indicating the p(r) and P (r ′ > r) have exponentially decaying tails [16].
To obtain a simpler analytical solution for p(r) than given by the hypergeometric functions, we approximate p x in the entire x-range by an exponential. We expect the results to be the same for x ≪ 1 [14] because p x takes non-vanishing values only for small x. Under this approximation, the integration of (3) yields where is the exponentialintegral function [15].
The probability distribution p prey (k) is obtained from (4) making the substitutions r = k/S and b = S/z, yielding We compare in Figs. 1(a) and (b) the predictions of (4) with numerical simulations. We find close agreement between our analytical expression and the numerical results. In particular, they show an exponential decay for large k. The deviations observed for small values of k are due to the fact that k j = Sr j is an average value implying that it is a good approximation only when the fluctuations of k j are small, which is no longer true for small k. Next, we address the statistics of the number of predators. Note that for r ≪ 1 [14], the predators of species i have, to first approximation, niche values n j > n i and that the segment r j is placed with equal probability in the interval [0, n j ]. Therefore, the probability for a species j to prey on i is r j /n j = x j n j /n j = x j , implying that the average probability for the species with n j > n i to prey on species i is x.
If we assume that S ≫ 1, the number of predators of i is the result of S − i independent "coin-throws" with probability x of being a predator and probability 1 − x of not being a predator, implying that the probability of species i having m predators is given by the binomial distribution. It then follows that the distribution of number of predators for a general species is the average over the different binomials . As predicted, we find a regime where the distribution is uniform followed by a Gaussian decay. We test our analytical predictions with empirical data [1] for (e) pprey(k) and (f) p pred (m) for Bridge Brook (solid line) and St. Martin (broken line).
In the limit of interest, S ≫ 1, x ≪ 1, and Sx = z, one can approximate the binomial distribution by a Poisson and the sum by an integral where γ is the "incomplete gamma function" [15,17]. For m < z/2, the function γ is approximately constant, while it decays with a Gaussian tail for m ≈ z. In Fig. 1(c) and (d), we compare the predictions of (7) with numerical solutions and find good agreement. In Figs. 1e,f, we compare our analytical predictions, Eqs. (5)-(7), with data from recent food webs: Bridge Brook (S = 25, z = 8.6) and St. Martin Island (S = 42, z = 9.8). We find that the distributions of number of prey is well approximated by the data and that the distributions of number of predators are "noisy" but still show the expected cutoff for m ≈ z and is approximately constant for m < z as predicted by Eq. (7). This agreement is remarkable since the webs analyzed are quite small so one might not expect the theoretical expressions to hold.
Next, we evaluate the fraction of top T , intermediate I and basal B species. As the names indicate, top species have no predators and basal species have no prey, while intermediate species are those with both prey and predators. The fraction of intermediate species is just I = 1 − (T + B). The fraction T of top species is, by definition, Note that a similar result is obtained if one calculates the sum (6) for m = 0. Since typically 5 < z < 20, Eq. (8) can be approximated simply as T = 1/z. To calculate the fraction B of basal species, we note that a species has no prey only if its range r falls in a region with no species [18]. In the limit of large S, the probability density for finding an empty interval of length δ is Se −Sδ , as predicted by the canonical distribution [19]. Thus, the probability of finding a species-free segment of length larger than r is e −Sr , which gives the probability for a species of range r not to prey on other species. Using (4), it follows that the average probability is: In the model [1], isolated species are eliminated, so they are not counted towards top or basal species. To correct the estimates (8)- (9) for this effect, we remove the isolated species. We estimate the number of isolated species to first order by assuming that having no prey is statistically independent of having no predators, implying that the fraction of isolated species is just the product of the fractions of top and basal species. This assumption does not take into account the possibility that a species with no prey is likely to have a low niche value n and hence it has a high probability to have predators. Nonetheless, this simple approximation provides an upper bound for the number of isolated species which leads to a lower bound on T and B, In Fig. 2, we compare our analytical predictions for the fraction of top and basal species with numerical simulations of the model. As expected, Eqs. (8)-(10) provide bounds for the numerical results. Finally, we calculate the standard deviations of the vulnerability and generality of the species in food webs generated according to the model. The vulnerability of a prey is defined as its number m of predators, and the generality of a predator as its number k of prey. Following Ref. [1], we define the normalized standard deviations of the vulnerability as σ 2 V = m 2 /m 2 −1 and of the generality as σ 2 G = k 2 /k 2 − 1. By definition, one has m = k = z/2 for both cases. To evaluate σ V , we first calculate m 2 . Equation (7) yields m 2 = z 2 /3 + z/2, so that We next calculate σ G , for any value of C, by direct evaluation of k 2 . If S ≫ 1, the number of prey of a species having a range r is k = Sr, and we find that k 2 /k 2 = r 2 /r 2 = 8(b + 1)/[3(b + 2)], implying that For C ≪ 1, σ G becomes a constant with value 5/3, a result that can also be obtained from (5). We show in Fig. 3 the results for our analytical expressions (11) and (12) and compare them with results from numerical simulations of the niche model. We have also studied the robustness of our predictions to changes in the particular formulation of the details of the model. The nature of approximations used in the derivations of the expressions for the distributions of the number of prey and predators, Eqs. (5)-(7), allow us to conclude that: (i) the distribution of the number of predators does not depend on the specific form of p(x). The only requirement is that the connectance C = x/2 tends to zero under some limit, so that z = SC remains finite when S tends to infinity; and (ii) the distribution of the number of prey depends on the functional form of p(x), but Eq. (7) will still be is obtained for all p(x) decaying exponentially as C tends to zero. Thus, it appears that our findings are robust under quite general conditions, a result that is not possible to obtain without an analytic treatment of the problem.  (12). Note that as for Fig. 2, removing isolated species leads to slightly less good agreement with the simulation results for σV . However, the removal of isolated species does not appear to be a factor in the deviations found for σG. The reason why σG underestimates the simulation results at small z values relates to the fact that kj = Srj is a good approximation only when the fluctuations of kj are small, which is no longer true for small k.
Our results are also of interest for a number of other reasons. First, we demonstrate for the first time that the distributions of number of prey and number of predators have different functional forms. Second, we show that both distributions have characteristic scales, i.e., both distributions have well defined means and standard deviations as S increases to infinity. Third, we find that the functional forms of the distributions of number of prey and number of predators depend only on the average connectivity z, and agree with empirical data. This result is rather surprising in face of the complexity of the empirical and model food webs. Finally, we show that other quantities of biological interest also depend exclusively z.