Indirect Energy Flows in Niche Model Food Webs: Effects of Size and Connectance

Indirect interactions between species have long been of interest to ecologists. One such interaction type takes place when energy or materials flow via one or more intermediate species between two species with a direct predator-prey relationship. Previous work has shown that, although each such flow is small, their great number makes them important in ecosystems. A new network analysis method, dynamic environ approximation, was used to quantify the fraction of energy flowing from prey to predator over paths of length greater than 1 (flow indirectness or FI) in a commonly studied food web model. Web structure was created using the niche model and dynamics followed the Yodzis-Innes model. The effect of food web size (10 to 40 species) and connectance (0.1 to 0.48) on FI was examined. For each of 250 model realizations run for each pair of size and connectance values, the FI of every predator-prey interaction in the model was computed and then averaged over the whole network. A classification and regression tree (CART) analysis was then used to find the best predictors of FI. The mean FI of the model food webs is 0.092, with a standard deviation of 0.0279. It tends to increase with system size but peaks at intermediate connectance levels. Of 27 potential predictor variables, only five (mean path length, dominant eigenvalue of the adjacency matrix, connectance, mean trophic level and fraction of species belonging to intermediate trophic levels) were selected by the CART algorithm as best accounting for variation in the data; mean path length and the dominant eigenvalue of the adjacency matrix were dominant.


Introduction
Food webs are icons of complexity, depicting intricate networks of feeding interactions. Since food webs can be studied both from the point of view of population dynamics and that of matter and energy flows, they bridge community and ecosystem ecology. Moreover, their study has led to insights that apply to other complex systems [1][2][3].
Examining food webs reveals a wide variety of indirect interactions, such as indirect matter and energy flows, trophic cascades, apparent competition, indirect mutualism and commensalism, and exploitative competition [4]. Indirect flows take place when energy or nutrients move own) whose niche value lies within a specified range. For species i with niche value n i , the feeding range has width r i and can be centered anywhere in the interval [r i /2, n i ] [16].
The n-species Yodzis-Innes model uses consumer functional responses and the scaling of metabolic rate with body size [23][24][25] to add realism to a simple model of trophic dynamics. (Since, as described below, the analysis used to quantify flow indirectness requires a conservative currency, the model's state variables were taken to be the total energy content of each species.) Including a functional response that saturates at high prey density improves model realism by acknowledging the fact that there is a limit to how much food an individual can consume. The use of scaling relationships helps incorporate biologically reasonable sets of parameter values into a theoretical model. The model, which employs variables and parameters whose dimensions and values are listed in Table 1, is described below. In keeping with environ analysis convention, we consider energy to flow from column j to row i, not the other way around, as is the convention in dynamic food web modeling.
In the absence of consumers, producer j grows logistically at rate r j B j 1 À B j K j , where B j is the total energy content (or population biomass) of species j, r j is its maximum growth rate and K j is the environment's carrying capacity for species j. To obtain the rate at which species j is eaten by species i, we reason as follows. The rate of consumption of j by i is proportional to the population size of i, B i . The quantity y ij is the maximum rate at which species i can consume species j, divided by i's metabolic rate, x i . Multiplying this quantity by x i gives x i y ij , the maximum per-capita consumption rate for i preying on j. The functional response, F ij (B), gives the consumption rate as a fraction of this maximum rate, yielding x i y ij F ij (B)B i for the actual rate. However, the predator does not ingest and assimilate all the prey it captures, so its consumption rate must increase to compensate for this. Dividing the previously obtained rate by the predator's efficiency, e ij , accomplishes this, giving the expression We now turn to the functional response. Following [26], a sigmoidal (Holling Type III) functional response with predator interference [27] was chosen, in part because it stabilizes the dynamics of food web models [28]. In this model, the consumer's search rate is proportional to prey abundance raised to the non-negative power q [29], and consumers of a given species interfere with each other with strength c [27]. As a result, Following [26], the values q = 1 and c = 1 were used. (Table 1) This results in relatively high predator interference and a pronounced sigmoid functional response. The overall differential equation for producer species j is: Consumers of species i lose energy to metabolism at rate x i B i , gain it from prey item j at rate x i B i y ij F ij (B), and lose it to consumers of species k at rate x k y ki B k F ki (B)/e ki [19]. Overall, we have: Table 1 summarizes the model's parameters and their values. To parametrize the model, empirically documented relationships between trophic level and body mass [30] were used to assign body masses to species in the model. Following [31], the (generally non-integer) trophic level of each species was computed as the mean of two quantities: (1) the integer distance between the target species and the closest basal species (those that do not prey on any other species); and (2) Levine's [32] generally non-integer flow-based trophic position, computed under the assumption that predators receive equal fractions of their diet from all prey species. The equal flows assumption allows flow-based trophic positions to be assigned to species in a purely topological web. The expression for flowbased trophic position is: where TL i is the trophic level of species i, S is the total number of species in the food web, and p ji is the fraction of species i's diet provided by species j. The mean of this quantity and distance from a basal species was used because it can be computed from topological information and provides a close approximation to the true flow-based trophic position in food webs for which flow data is available [31]. Species were then assigned metabolic rates using the 3/ 4-power scaling relationship between metabolic rate and body size [23][24][25].
A new flow-based dynamic network analysis method called dynamic environ approximation (DEA [33]) was used to compute FI. The basic logic of DEA is as follows. If a food web has adjacency matrix A, then A k gives the number of paths of length k between each pair of species and P m k¼1 A k gives the total number of paths of length m or less between each such pair [34,35]. If the structure of the network changed over time, then the number of paths would be given by the product series DEA uses a related product series of matrices describing energy flows in the food web to trace the flows through the system. The flow matrix is then normalized by the total outflow from the donor species to create a matrix, G, of nondimensional flow intensities for each integer time step. Then, for a window of m time steps, we have integral flow is the matrix transforming inputs received at time t into eventual flows throughout the system, summing over all times up to t + m. Flow indirectness is then computed as FI ij (t) = (N ij (t) − G ij (t))/N ij (t) [33]. This method was used to find FI for each interaction in the web and the entries of the FI(t) matrix were then averaged.

Results
Over the full range of parameter values, the mean flow indirectness of the model food webs was 0.092, with a standard deviation of 0.0279. It increased with system size but peaked at intermediate connectance levels, resulting in the pattern seen in

Determinants of flow indirectness
A classification and regression tree (CART) analysis performed in R [36] with the package rpart [37] was used to explore which aspects of web structure were most strongly correlated with mean FI. The algorithm performed a split only when doing so increased R 2 by at least 0.01. Also, to avoid overfitting, 10 cross-validations were performed at each step.
Out of 27 potential ecological and graph-theoretic predictor variables, only five (mean path length, dominant eigenvalue of the adjacency matrix, connectance, mean trophic level and fraction of species belonging to intermediate trophic levels) were selected by the CART algorithm as best accounting for variation in the data. Two of these, mean path length and dominant eigenvalue of the adjacency matrix, were dominant. (Fig 3) The CART model accounted for 82.2% of the variation in FI. A full list and explanation of the potential predictor variables used is given in the Supporting Information. Fig 4 provides a more detailed look at the relationship between flow indirectness and these predictor variables.
Despite predicting the observed FI values very well, the CART model does not explain the curvilinear relationship between connectance and FI. A quadratic model was therefore used to test the relationship between connectance and path length. The model had an R 2 of 0.9639 and the ΔAIC between this model and a linear model was -240817. (In both models, the intercept was forced through zero, as no other value is possible.) Marks seagrass web, which has a much larger mean path length for its size and connectance than any of the other modeled or empirical webs, and the St. Martin Island web, whose dominant eigenvalue is zero because the web contains no cycles. Predictions of the flow indirectness values of the empirical webs were made by following the CART tree developed using the simulated webs. (Table 2) Freshwater food webs had consistently lower predicted FI values than terrestrial and marine webs. This is due to the freshwater webs' lower mean path lengths, which appear to simply be a result of their smaller sizes. Sizenormalized mean path length did not vary systematically with ecosystem type.

Discussion
The flow indirectness of a food web provides an indication of the extent to which pairwise interactions are mediated by the network in which they are embedded. When FI is high, a significant portion of the effect of an experimental manipulation of the abundance of one species on that of another will be determined by the rest of the food web and may not generalize to other systems.
We found that the mean flow indirectness over the full range of size and connectance values was 0.092, with various combinations of topological variables giving FI values as low as 0.01 or as high as 0.33 (Fig 3). The increase of FI with web size (Fig 2) is particularly important, as it indicates that, in real food webs, a substantial proportion of energy can be expected to travel over indirect paths (Table 2). Due to dissipation, each individual path carries little energy, but the number of paths makes up for this [6]. Cycling plays an important role in pathway proliferation, so it makes sense that the dominant eigenvalue of the network, a measure of the amount of cycling in the network [44,45], is strongly positively related to FI.
The quadratic relationship between connectance and path length explains the curvilinear relationship seen in Fig 2. It may be that increasing connectance first results in longer paths by connecting more species, but then the increasing number of links short-circuits long paths, resulting in a lower mean path length and correspondingly lower flow indirectness.
It will be useful to examine the effects of other system attributes on the flow indirectness and find out whether this quantity is linked to the vulnerability of food webs to species loss. [46] found that, for sixteen empirical topological food webs, vulnerability to cascading extinctions in the face of species loss was negatively correlated with connectance and uncorrelated with the prevalence of omnivory (in spite of the correlation between omnivory and connectance). However, the webs examined in that study had connectances ranging from 0.026 to 0.315-values falling within the range in which flow indirectness increases with connectance. (Fig 2) It would be instructive to examine the effects of species loss on model webs with higher connectance values, to see whether the positive relationship between connectance and robustness continues to hold and whether omnivory becomes a more important determinant of robustness as flow indirectness declines.
Two major frameworks exist for studying networks of trophic interactions and the movement of energy within ecosystems: those of community and ecosystem ecology. When food webs are studied from a community ecology perspective, the emphasis is on individual species and their population dynamics. Such webs are as detailed as possible but often omit parts of the biota at the study location, especially decomposers and detritivores in terrestrial systems. (The Coachella Valley, CA food web of [47] is a prominent exception.) By contrast, the ecosystem framework uses comprehensive, usually highly aggregated models that focus on the movement of energy and nutrients. Researchers working within these two frameworks have ignored each other's work to a remarkable extent.
Environ analysis [8,10,11] is a set of conceptual and mathematical tools for analyzing networks of stocks and flows. It has traditionally been applied to phenomenological models of real ecosystems. This is both a strength and a weakness, in that the analysis stays close to reality but is tied to a relatively small number of models that are usually highly aggregated. In particular, the six-compartment intertidal oyster reef model of [48] may be the Drosophila of environ analysis because of the number of techniques and hypotheses that have been demonstrated and tested using it [13,33,[49][50][51]. This study is not the first to apply environ analysis to a large synthetic model. [9] created models of ecosystems by assigning species to one of six functional groups: primary producers, herbivores, carnivores, omnivores, detrital feeders and detritus. Each functional group, including detritus, contained the same number of species, ranging from five to one hundred. Biologically plausible intergroup interactions were then randomly assigned. The model used linear donor-controlled dynamics with randomly selected coefficients. Thus, this model possessed some realism with regard to functional groups, very little with regard to network structure, and almost none (except in the case of flows to detritus) with regard to dynamics.
The advantage of theoretical models such as the Yodzis-Innes model is that they describe causal relationships between species. Compartment models, on the other hand, are typically phenomenological, "bookkeeping" models. When dynamical assumptions such as donor control are added to these models, they are typically very simple and lack biological justification. The relatively detailed causal assumptions and parameter constraints of the Yodzis-Innes model may be criticized as being overly complex and unrealistic, but they are probably less wrong than linear models with donor control, which assumes that the amount of prey eaten by a predator species depends only on the prey's population size. However, research on such dynamically simple models has produced insights into ecosystem function and network properties.
Working within the framework of linear steady-state models of conservative energy and matter flows and storages, [49] identified six network characteristics that directly increase flow indirectness in compartment models: number of compartments, connectance, storage, cycling, feedback and magnitude of direct flows. Most of the model webs in this study had much lower mean flow indirectness values than the ecosystem models examined in previous work [9,14]. This is likely due to the fact that niche model webs, unlike the models studied before, do not include detritus or detritivores. Therefore, they contain substantially less cycling than ecosystem-oriented models. Since cycling greatly increases the fraction of model currency traversing indirect paths [49], its absence must reduce flow indirectness.
This omission is significant because detritivory is a nearly universal feature of real food webs [52][53][54]. Thus, the current results almost certainly underestimate the true importance of indirect flows in natural food webs and future work should attempt to include detritus and detritivores. A logical next research step would be to use a version of the niche model modified to include detritus and detritivory, such as that of [55].
Other important topics for future research include the sensitivity of the results described here to parameter values and model assumptions and the examination of energy cycling in model webs with and without detritus [6,13]. In particular, the standard niche model uses niche values that are uniformly distributed between 0 and 1. However, the niche value has been found to be correlated with body size [16,21,22]. Therefore, the distribution of niche values should be derived from body size distributions observed in nature. The allometric diet breadth model [56] approaches this idea but relies on previously specified body size data, although this could be randomly generated. A simpler approach would retain all the assumptions of the niche model but use a more realistic distribution of niche values.
The results reported here helps bridge contemporary food web ecology and systems ecology, while providing a new way of looking at ecosystem complexity. It is also possible to apply dynamic environ approximation to non-trophic stock and flow networks such as dispersal networks [57] and human systems such as roads and economies, and doing so may provide insights into their function.

Web Construction and Simulation
The goal of this study was to explore the importance of indirectness in a commonly studied theoretical food web model, the niche model [16]. This model was selected because it is frequently studied and reproduces many features of real food webs with a fair degree of accuracy [16]. In this model, each species has a niche value, n i , a feeding range width, r i , that can be interpreted as the fraction of possible niche values that can be consumed by species i, and a feeding range center. The niche value for each species is drawn from a uniform distribution ranging from 0 to 1. Range centers are drawn from a uniform distribution ranging from r i /2 to n i . A uniform distribution is used both for its simplicity and to reflect the hypothesis that niche values in real ecosystems are roughly uniformly distributed.
Species' range widths are generated by drawing values from a beta distribution whose mean is twice the connectance of the web and multiplying them by the species' niche value. (Specifically, the niche model uses a beta distribution with α = 1 and b ¼ 1À2C 2C , where C = 2L/ S 2 , the connectance of the web, S is the number of species, and L is the number of links.) Since the expected value of this distribution is 2C and that of n i is 0.5, this procedure results in range width having an expected value of C. Because niche values are uniformly distributed on the [0, 1] interval and a consumer's feeding range width is the fraction of this interval that contains potential prey, the fraction of species a given consumer preys on is approximately its range width. This gives the food web the desired connectance. Each species is assumed to prey on all species within its range, including itself, and a food web directed adjacency matrix is assembled [16].
Candidate webs generated by this method were tested to ensure that they had at least one producer and consisted of only one set of connected species, termed a component in graph theory [35]. For the latter test, the Laplacian matrix, L, was used. This matrix is defined as the difference between the degree matrix D, which has the degree of the graph's nodes on the diagonal and zeros elsewhere, and the undirected adjacency matrix A from which self-loops are excluded, making a ii = 0. The equation for the Laplacian matrix is then L = D − A. The number of times 0 appears as an eigenvalue of the Laplacian is the number of components in the graph [58,59].
If a web passed these tests, trophic levels were assigned to each species as the mean of distance from the closest basal species and the flow-based trophic position method (Eq 4), computed under the assumption that predators receive equal fractions of their diet from all prey species [31,32]. (The equal flows assumption allows trophic levels to be computed for a purely topological web.) Taking the mean of these two methods used provides a good approximation to real trophic levels in quantitative food webs and ecosystem models [31]. Trophic levels were then used to assign body sizes as 10 T i −1 where T i is the trophic level of species i [30], and massspecific metabolic rates were assigned using 3/4-power scaling [23][24][25]. Initial abundances were drawn from a uniform distribution ranging from 0.5 to 1, ensuring that the simulation results were not artifacts of a particular set of initial conditions and that all species were initially present at ecologically significant levels. The simulation was then run for 1000 time steps using fourth-order Runge-Kutta integration with a step size of 0.01, after which time a steady state had been reached or closely approximated. At that point, any extinct species were removed and the simulation run for 1000 more time steps. In order to avoid transient dynamics, only this second run was analyzed with DEA.
The effect of food web size (10 to 50 species) and connectance (0.1 to 0.48, in increments of 0.02) on flow indirectness was examined. Because the niche model is stochastic, 250 model realizations were generated and simulated for each pair of size and connectance values.

Dynamic Environ Approximation
In the standard Yodzis-Innes model, the amount of energy gained by a predator in a predation event is less than the amount lost by the prey. The boundary inputs and outputs required to balance the system's energy budget are not explicitly tracked. Therefore, in order to create the conservative flow matrix required for environ analysis, producer growth was conceptualized as an input to the system, while uneaten or unassimilated food and metabolic losses were conceptualized as outputs. For each integer time step, a flow matrix consisting of the second terms of Eq 3 (with negative outflows from each compartment on the diagonal) was set and used to compute the throughflow-normalized flow matrix G (g ii = 0, g ij ¼ f ij t out j , where f ij is the energy flow from j to i and t out j is the total outflow from j). This was then used to perform DEA with a window size of 20, which previous work indicated is typically enough to capture all relevant dynamics [14,33]. As a large majority of simulation runs had reached or nearly reached a constant steady state, the N matrix was only computed for one starting time. Flow indirectness (FI) was then calculated as N − G and the mean for each web was computed. In 62 runs out of 205,000, FI values larger than 1 or less than 0 were obtained; it was concluded that the integration step size was too large for the dynamics of these runs and they were excluded from further analysis. FI values for diagonal entries, which represent cycles linking a species to itself, were taken to be 0.
All simulation and analysis code is given in Supporting Information S1 Source Code.
Supporting Information S1 Text. Listing of variables used in CART analysis. (PDF) S1 Source Code. C++ code for food web simulation and environ analysis. (ZIP)