## Abstract

Patterns of species interactions affect the dynamics of food webs. An important component of species interactions that is rarely considered with respect to food webs is the strengths of interactions, which may affect both structure and dynamics. In natural systems, these strengths are variable, and can be quantified as probability distributions. We examined how variation in strengths of interactions can be described hierarchically, and how this variation impacts the structure of species interactions in predator–prey networks, both of which are important components of ecological food webs. The stable isotope ratios of predator and prey species may be particularly useful for quantifying this variability, and we show how these data can be used to build probabilistic predator–prey networks. Moreover, the distribution of variation in strengths among interactions can be estimated from a limited number of observations. This distribution informs network structure, especially the key role of dietary specialization, which may be useful for predicting structural properties in systems that are difficult to observe. Finally, using three mammalian predator–prey networks (two African and one Canadian) quantified from stable isotope data, we show that exclusion of link-strength variability results in biased estimates of nestedness and modularity within food webs, whereas the inclusion of body size constraints only marginally increases the predictive accuracy of the isotope-based network. We find that modularity is the consequence of strong link-strengths in both African systems, while nestedness is not significantly present in any of the three predator–prey networks.

## 1. Introduction

Ecosystem dynamics are strongly sensitive to the structural organization of trophic interactions among species. Food webs are often characterized by qualitative networks, where only the presence/absence and direction of interactions (links) are established among species (nodes) [1]. Quantitative networks include link-strengths, which establish the relative importance of each link (figure 1*a*). Although the potential importance of link-strengths has long been recognized [2], this information is rarely included in empirical food-web analyses [3]. Both empirical and theoretical works suggest that variance in strength among links in a network has a large impact on structure and, by extension, dynamics [4,5]. Because species interactions vary over space and time, and are estimated imprecisely, the strengths of individual links are probabilistic (where a given link-strength has an associated probability; figure 1*b*). However, it is not clear how to integrate this information into ecological networks, or how link-strength variation impacts structure. This is important because the structure of species interactions ultimately affects ecosystem dynamics [6,7].

Specific patterns of interaction are often observed in ecological networks [8–10]. In general, such patterns arise when physical or behavioural limitations direct interactions between species or groups of species. For example, spatial partitioning of vegetation on the landscape may tend to attract particular groups of herbivores and, by extension, their predators, and this may lead to resource partitioning dictated primarily by habitat heterogeneity [11]. Alternatively, some of these structural patterns of interaction may contain information regarding the role of physiological factors, such as body size, in shaping the organization of ecological communities [12]. Body size constraints are often realized by the inability of predators to consume prey owing to differences in their respective body sizes [12–14], thereby placing strong limitations on the structure of food webs.

Body size constraints may affect two structural features of antagonistic networks that have potentially important dynamical implications [9,15]: nestedness, which arises when predators have incompletely overlapping diets, and modularity (compartmentalization), where some predators exhibit strong overlap in their prey assemblages but show little or no overlap with other predator groups [9,16]. Nestedness is often observed in mutualistic networks [17] and in some antagonistic networks [10], and may occur when organisms follow hierarchical rules of interaction [14], including those generated by body size constraints [13]. Theoretical food-web models that assign links hierarchically and that allow some prey overlap among predators reproduce many statistical features of empirical food webs [18]. Recent theoretical work has shown nested interactions to be unstable in antagonistic networks [9,15]. In contrast, modularity may stabilize food webs [9,19] by minimizing extinction cascades [20]. However, analysis of the incidence and impact of nestedness and modularity has been confined to qualitative networks. The presence of nestedness and modularity, as a function of link-strength, is of key importance because weak and strong links may differently affect the organization of ecological networks [21].

Although link-strength variability adds substantial complexity to the description of ecological networks, it is necessary to understand whether, and to what extent, variation in link-strengths impacts network structure. To assess the impact of link-strengths on the organization of ecological networks, we introduce a method by which stable isotope ratios of predator and prey species can be used to estimate the relative magnitudes of species interactions in a probabilistic framework. An isotope-based approach permits analysis of food webs that are difficult to observe and/or too large to manipulate. This may be useful for understanding how characteristics of ecological networks change over long periods of time or across large spatial ranges. Here, we show that link-strength variability can be assessed as a hierarchy of distributions, such that ecological variation and measurement uncertainty are quantified on multiple scales within a network, and that these distributions can be estimated from incomplete data. We then show how stable isotope ratios can be used to build predator–prey networks, and use three mammalian communities quantified from stable isotope data to assess whether link-strength distributions predict nestedness and modularity in these communities, and the effect of body size constraints on these structures. Our results show that even limited knowledge of the variance in strength among links provides information regarding the topological properties of three empirical predator–prey networks, and that structural organization changes as a function of link-strength across networks. In addition, although body size limitations constrain species interactions, our results indicate that the inclusion of such constraints only marginally improves only predictions of network structure in these communities.

## 2. Methods

Links that denote interactions between species within predator–prey networks can be described qualitatively (i.e. by the presence or absence of links) and quantitatively (i.e. by links that are weighted with single values typically normalized to range between 0 and 1; figure 1*a*). Here, we expand the quantitative approach to describe links probabilistically, such that link-strengths for a given pairwise interaction have associated probabilities. We define link-strength in terms of the proportional biomass flow from prey to a consumer (cf. Bascompte *et al.* [8]). Thus, links connecting a single consumer to its prey denote the proportional contribution of biomass of prey to a predator and sum to one (figure 1*b*). Furthermore, link-strength probability distributions are described hierarchically, such that variance among links and variance within links are considered separately.

### 2.1. Link-strength variability in predator–prey networks

Quantitative ecological networks have been shown to have long-tailed link-strength distributions [8,22,23]. This distribution describes the variation of strength *among* links in a system (hereafter network-level interaction distribution, NLID; figure 1*c*). The long-tailed nature of this distribution reveals that most interactions are relatively weak and some are strong [23]. The specific shape of these distributions can have a significant impact on the stability of theoretical ecosystem models [24,25], though its influence on structure is not well understood.

Link-strengths for any observed ecological network are the estimates of the interactions among species within a community, and therefore all such networks are inherently probabilistic. The variation of strength *within* links in an antagonistic network reflects the temporal and spatial variation in prey consumption, dietary variation among individuals and measurement uncertainty. Within-link variability is typically considered only in theoretical investigations or in controlled laboratory settings. We distinguish such distributions by referring to them as pairwise interaction distributions (PIDs; figure 1*d*), which are unique for each link in the system. We note that the NLID and the set of PIDs describing links in a network are hierarchically related, such that the mean and variance of the NLID is ultimately determined by the moments of all PIDs in a system. In general, a PID between an interacting consumer and resource is difficult to measure; we approach this problem by using stable isotope ratios to estimate link-strengths.

Ratios of stable isotopes (typically of carbon and nitrogen) can be used to accurately quantify trophic relationships between consumers and prey [26–28]. Estimating trophic interactions with stable isotopes is a particularly advantageous approach for systems that are difficult to observe, such as when animals are elusive or primarily nocturnal. Moreover, this isotope-based approach enables an accurate quantification of link-strength probabilities for a given trophic interaction, thereby accounting for both ecological variability and measurement error. The carbon and nitrogen isotope ratios of a consumer, with respect to potential prey, characterize the *isotopic niche space* of the consumer [27,28]. It has been shown that the isotopic niche is comparable to the traditional concept of a dietary niche [27]. Posterior probability distributions of the proportional contributions of each prey to a consumer's diet were numerically estimated with a Bayesian isotope mixing model [28], and represent the PIDs of the network. These proportional contribution estimates incorporate the uncertainty inherent to isotopic approaches, including analytical error, isotopic variability in prey and consumers, poorly constrained fractionations between trophic positions [29], isotopic similarity of multiple prey [30] and underdetermined or non-unique solutions [28]. PIDs calculated from a mixing model are each marginal distributions of a single large multivariate (Dirichlet) distribution. As such, each marginal PID is constrained such that the link-strengths connecting a given consumer to all of its prey must sum to one.

Posterior distributions of Bayesian isotope mixing models are challenging to determine analytically, but can be estimated numerically using importance sampling approaches [28]. Numerical approximations of the PIDs and NLID of a food web can be assembled directly from mixing model output. A set of potential contribution-to-diet values are calculated, where each element of the set represents a ‘potential’ combination of link-strengths between the consumer *j* and its *n* potential prey and sums to one. For our study, we generated 1000 such combinations, and organized them as rows of a matrix **W**_{j}, where each column represents different potential contributions of prey *k* (out of *n* available prey). A PID for a given consumer *j* linked to a prey resource *k* is then described by a vector of 1000 potential ‘prey contribution values’, **w**_{jk}, with each prey (*k*) represented by a column in the matrix **W**_{j} (see figure 1*e* for an illustration of this matrix). Overall, a set of randomly drawn rows for each matrix **W**_{j}, across *m* consumers in the system, represents a set of link-strengths connecting every consumer and resource in the predator–prey network, and is equivalent to a single quantitative network, where each link has a single link-strength. By extension, the full set of matrices **W**_{j} for *j* = 1, … ,*m* for each of *m* consumers in the system, represented as , describes the probabilistic network, or alternatively, an ensemble of potential quantitative networks. The NLID for the probabilistic network can be assembled by uniformly sampling across all PIDs (**w**_{jk}) and across every matrix (**W**_{j}), such that the row vectors describing a potential set of pairwise interaction strengths linking predator *j* to all prey are drawn randomly across all consumers in 𝕎. Accordingly, the NLID is the frequency distribution of all link-strengths in a probabilistic food web. To evaluate the structural properties of the probabilistic predator–prey networks, we considered each potential quantitative network in the ensemble independently, such that a range of values were calculated for a given structural property.

This framework permits a static representation of link-strength distributions, whereas they are likely time-dependent in nature, varying with season or the age structures of populations. Because we are using stable isotope ratios from hair and bone tissues, the estimated link-strength distributions represent interactions averaged over months and years, respectively. This time-averaging quality of isotope-derived link-strengths averages out potential short-term trophic relationships. If prey-switching was a common strategy for a particular consumer, our mixing-model approach would generate link-strength distributions with either high variability or multi-modality, depending on the number of resources being interchanged. A dynamic approach, whereby link-strength distributions are time-dependent, could be assessed using isotopic techniques by serially sampling consumer tissues (e.g. sea otter whiskers [31]).

### 2.2. Structure as a function of link-strengths

To quantify structural properties as a function of link-strength, we developed a cut-off algorithm to sequentially measure structure, as links with low contribution-to-diet values are successively ‘trimmed’ from the system [21]. This procedure was carried out for every potential (1000) quantitative network in the probabilistic ensemble. This cut-off algorithm iteratively eliminates the weakest links in a system, allowing us to explore the structure of both the entire and just the core resources used by consumers. Given a potential quantitative network within an ensemble, a link exists if and only if it has a link-strength >*i*, where *i* is the cut-off value. The structural properties of each qualitative network are then analysed for a set of cut-off values *i* (i.e. ). For each measured network property, this algorithm quantifies five values (one for each cut-off), describing how structural properties change as we move from describing the entire ecological network to just the strongly connected prey contributing to predator's diets.

## 3. Analysis and discussion

To determine the extent that link-strengths impact the structural organization of species interactions, we first provide a baseline by which (i) the effects of link-strength distributions on network structure can be quantified and (ii) ecosystems of varying species richness and predator/prey ratios can be compared. The organization of species interactions dictates the extent to which consumer species share their resources. We tested whether networks generated by randomly sampling link-strengths from the NLID reproduce the patterns of interaction observed in empirical predator–prey networks. Model systems drawn randomly from the NLID may retain the distribution of link-strengths of the original ecological network, but any correlation among link-strengths will disappear (figure 2*b,c*). Model networks drawn from the NLID do not have correlations between the link-strength distributions that may characterize predator–prey interactions. If there are strong ecological trade-offs driving species interactions, we would expect empirical networks to have strong correlations between link-strength distributions. For example, a weak interaction between a predator and a prey may exist owing to a strong interaction between the same prey and a different predator, leading to niche partitioning. However, networks drawn randomly from the NLID will not conserve this pattern. These NLID-derived model networks provide a benchmark to study the role of ecological factors other than those generating the distribution of link-strengths in shaping predator–prey networks. Accordingly, we use the NLID to determine the extent that ecological factors influence patterns of interaction by examining (i) whether the shape of the NLID is predictive of network structure and (ii) whether and to what extent the NLID can be estimated when the majority of interactions are either unknown or unobserved. We then use models drawn from the NLID to determine the influence of link-strengths on the patterns of interaction measured for three mammalian predator–prey networks.

### 3.1. Structural implications of the network-level interaction distribution

We found that networks with link-strengths drawn randomly from the NLID share general structural similarity with the corresponding empirical predator–prey networks (see the electronic supplementary material, appendix S1). Moreover, the specific shape of the NLID, as determined by the first two moments of the distribution, was a strong predictor of whether consumers had generalist or specialist diets. Because link-strength distributions of networks generated from the NLID cannot be correlated, knowledge of both the NLID and species-specific PIDs can be used to determine the extent that trade-offs in the strength of trophic interactions contribute to network structure. The shape of the NLID is a function of pairwise interactions in ecological networks, and knowledge of these interactions is often fragmentary. Therefore, the usefulness of this distribution is contingent on adequate sampling of interactions. In §3.2, we determine the effect that unobserved interactions have on an investigator's ability to accurately parametrize the NLID. Accurately estimating the NLID from incomplete data may enable predictions of network structure for systems that are not well known or difficult to observe.

### 3.2. Estimating the network-level interaction distribution from limited observations

Most ecosystems are difficult to observe and have an unknown number of interacting species. While the true NLID of a food web may be unknowable, we performed a numerical simulation to determine the accuracy of NLID estimation by subsampling observed interactions. The distribution of strengths among links (the true NLID) was assembled from a series of PIDs that were generated from random independent beta-distributions, which incorporate the constraint that link-strengths vary between 0 and 1; the frequency distribution for a random variable *X* = *x* that is beta-distributed is given by
where, for example, . We randomly sampled the two shape parameters from uniform ranges and to construct a continuum of food webs varying in size from 2 to 500 links without a specified structure, thereby relaxing the constraint that individual PIDs are dependent on other PIDs in the network. NLID assembly was then carried out, as described already, for link-strength vectors (**w**_{jk}) randomly drawn from each beta-distribution. In this case, link-strengths are not constrained to sum to one for a given consumer because we do not assume that all of a consumer's interactions are measured. If link-strengths are constrained to sum to one, the structure becomes implicit, and it is possible that different patterns of interaction may impact estimates of NLID shape; however, we isolate our simulation to the more general scenario where network structure is assumed to be unknown.

We estimated the accuracy by which subsamples of link-strengths from a food web of a given size can accurately estimate the true NLID with the following approach. We (i) randomly subsampled from the available PIDs, (ii) constructed an estimated NLID from the subsampled PIDs, and (iii) measured the absolute difference of the mean (mean error) and standard deviation (s.d. error) of the true and estimated NLIDs to assess the similarity between the two distributions. For food webs with 100, 200, … , 500 links, modelled results were fit to an exponential function to determine the rate at which the observed mean and standard deviation of the estimated NLID converged to the known moments of the true NLID as sampling effort increased (see figure 3 and electronic supplementary material, table S2).

In a quantitative or probabilistic food web, the shape of the NLID may constrain network structure and impact potential dynamics. Our numerical simulation revealed that the error of NLID estimates decreased exponentially as the proportion of sampled PIDs increased (see figure 3 and electronic supplementary material, table S2). Thus, sampling a small proportion of PIDs will yield accurate estimates of the NLID, particularly for large interaction networks. The exponential nature of this relationship indicates that, even for smaller systems, a relatively modest sampling effort will yield fairly accurate estimates of the true NLID. We suggest that an analytical determination of this relationship would be a worthwhile pursuit, but is beyond the scope of this paper. Because the shape of the NLID can be estimated when the network is incompletely sampled, it is useful to explore whether the NLID contains information regarding the topological organization of link-strengths.

Our results suggest that NLID is predictive of general structural features of a food web, and can be estimated from limited observations. Moreover, because the shape of the NLID is related to consumer specialization, knowledge of the NLID can be used to predict the likely distribution of consumer specialists in communities where many interactions cannot be measured. Knowledge of food webs increases incrementally as data are gathered [1], and these features of the NLID may be useful for generating Bayesian prior probabilities for structural properties that can, in turn, be updated as more information becomes available [32]. In §3.3, we use model networks drawn randomly from the NLID as a baseline by which nestedness and modularity of three empirical predator–prey networks derived from stable isotope data can be compared. We then assess the influence of link-strengths on the structural organization of these communities.

### 3.3. Case studies: ecological networks from Saskatchewan, Amboseli and Lake Naivasha

We examine three mammalian predator–prey networks using stable isotopes to quantify link-strength distributions (see electronic supplementary material, appendix S2 for treatment of isotopic data). One of these isotopic datasets is from a mammalian community in Saskatchewan, Canada (nine predators and eight prey; electronic supplementary material, figure S1), an inland boreal forest ecosystem interspersed with isolated rivers and wetlands [33]. Two others are from East Africa: Lake Naivasha (three predators and eight prey; electronic supplementary material, figure S2) [34], and Amboseli (four predators and 10 prey; electronic supplementary material, figure S3) [35], both savannah–woodland mosaics spanning high-altitude montane forests to low-altitude grassland-dominated environments across southern Kenya. By restricting our analysis to mammalian predator–prey networks, which are components of ecological networks, we can more accurately quantify mass flow between species. The importance of the large-bodied mammalian component of food webs has been recently supported by a growing body of evidence linking extinction of these species with large-scale ecosystem reorganization [36]. This approach contrasts with the investigation of species-rich food webs where link-strengths are often necessarily limited to qualitative descriptions. Furthermore, we limit our analysis of each community to locally abundant species more than 4 kg, thereby focusing on interactions among larger consumers and their potential prey. To incorporate body size constrains, forbidden interactions were built into mixing models for each empirical system based on body size differences between consumers and prey, following observations by [13] (see the electronic supplementary material, figures S2–S4). In short, if a link is forbidden, there is no link between the two species and the corresponding link-strength is zero. To assess the role of these constraints in structuring the empirical predator–prey networks, we compared two probabilistic network models: one that does not include body size constraints (model 1), and one that retains the same number and structure of forbidden links as the respective empirical network (model 2). Pairwise interaction distributions for these systems were quantified with the mixing model MixSIR [28]. The strengths of non-forbidden links for both model networks were randomly drawn from the corresponding NLIDs of empirical predator–prey networks. In model 1, all predators interact with all available prey, whereas model 2 shares the same qualitative link structure as the empirical network.

We quantified nestedness (𝒩) and modularity (ℳ) for each of the three empirical predator–prey networks, because these properties are thought to have important ramifications for the dynamics of ecological systems [9,10,20]. To quantify nestedness, we use the nestedness based on overlap and decreased fill metric [37], where a value of 0 denotes an un-nested network, and a value of 1 denotes a fully nested network. Modularity has traditionally been quantified as the local density of links in a network [38], although other metrics exist [39]. Because our predator–prey networks have two trophic levels, we estimate the average local link density with a bipartite clustering coefficient (see Zhang *et al.* [40]; electronic supplementary material, appendix S1). If the average local link density becomes greater than the overall link density, the network is modular and ; if the average local link density becomes less than the overall link density, the network is non-modular and (cf. Araújo *et al.* [41]).

Departures in nestedness or modularity estimates of the empirical networks from predictions generated by models 1 and 2 would indicate that link-strengths are constrained by different physical or biological factors. These factors include specialization on alternative prey, competition among consumers leading to low resource overlap, or spatial heterogeneity that affects prey consumption. Because some structural properties, such as nestedness [42], correlate with food-web size [43], we report the relative difference of structural property values between empirical and model network ensembles. These relative structural metrics are expressed as , where refers to the value of either nestedness () or modularity () at cut-off *i*. A value of ‘0’ indicates no difference between the empirical predator–prey network and corresponding model (i.e. high similarity, see electronic supplementary material, appendix S2); if the empirical network has a higher value than that of the model network; if the empirical network has a lower value than that of the model network.

Body size limitations are included as forbidden interactions because there is empirical evidence supporting the role of body size in constraining interactions in mammalian food webs [13]. Therefore, a model network with this forbidden link structure imposed *a priori* (model 2) has the same qualitative structure as that of the corresponding empirical predator–prey network; however, this is not guaranteed for the quantitative structure. A comparison of and for each system reveals that models without or with body size constraints (models 1 and 2, respectively) lead to similar structural predictions across cut-off values. Model 2 incorporates body size constraints and generates networks that are generally more similar to empirical networks (see the electronic supplementary material, table S3), showing that the inclusion of body size constraints does increase the predictive accuracy of network models. However, the addition of body size information does not lead to or values significantly different from those without body size constraints (model 1). This result is important because it shows that the organization of forbidden links imposed by body size constraints is not the primary driver of structure in these predator–prey networks, despite being the only determinant of structure when link-strengths are not considered.

Our structural analysis of the Saskatchewan system shows that it is not significantly nested (figure 4*a,d*). Modularity was marginally higher than expected for low cut-off values, and there was strong similarity between empirical and model network modularity for high cut-off values. These relationships were maintained when smaller predators were eliminated from the system (such that the species’ body size ranges were more comparable to those of the African predator–prey networks; electronic supplementary material, figure S6). In fact, our results reveal that and are strongly inversely related across cut-off values for all systems (particularly for model 2; see electronic supplementary material, figure S7), a relationship that has been observed in previous studies [42]. Our analyses show that this is also observed across cut-off values within probabilistic predator–prey networks. The structural properties of both African systems were qualitatively similar. Across all cut-off values, Amboseli and Lake Naivasha showed less nestedness than predicted by both model networks for cut-off values of 0.2–0.4 (figure 4*b,c*). As such, Amboseli and Naivasha are characterized by groups of consumers that show a low degree of overlap in prey use, particularly for prey forming the core of their diets (figure 4*e*,*f*).

A comparison of and across empirical networks has three important implications. First, model networks with interactions drawn from the empirical NLID (models 1 and 2) do not necessarily result in accurate predictions of nestedness and modularity across cut-off values, and accuracy is lowest for African systems (figure 4). For both African systems, structural differences between model and empirical networks are generally highest for cut-off values of 0.2–0.4, showing that stronger links are more modular than predicted. However, when the weakest links are considered, predictive accuracy is generally increased such that and values are closer to zero. If these systems were investigated without considering link-strengths, the minimal structural differences (low and values) would imply that model networks drawn from the NLID are wholly predictive of structure. Our analysis of higher cut-off values reveals such a conclusion to be false, highlighting the importance of considering the effects of link-strengths when quantifying the structure of species interactions. It is possible that weak and strong links differ in the processes that generate their patterns of interaction [21]. Moreover, the structural differences between weak and strong links are likely to influence ecosystem dynamics, and this deserves additional attention from both a theoretical and an empirical perspective.

Second, the strongest link-strengths of the Saskatchewan system are less modular, whereas both African systems exhibit higher modularity than predicted by our model networks. None of the three empirical predator–prey networks are highly nested. The similarity of the two African predator–prey networks relative to Saskatchewan may reflect landscape-driven constraints on species interactions. Recent work has shown the Serengeti food web to be compartmentalized with respect to spatial guilds, where local habitat mosaics result in compartments of tightly interacting species [11]. Our results lend general support to this finding, and we show that modularity is emphasized for stronger links in the network. Furthermore, our analysis suggests that systems that are more spatially homogeneous, such as the coniferous woodlands of Saskatchewan, may be less modular and more similar to ecological networks predicted from the NLID. Future work could use our approach to investigate the relevance of spatial organization in structuring trophic interactions across landscapes.

Third, the limited improvement in structural predictions when body size constraints are included (model 2) suggests that forbidden links do influence link-strength distributions in these systems, but not significantly (figure 4; electronic supplementary material, table S3). Foraging constraints due to body size limitations are often implicated in driving structure in terrestrial mammalian food webs [13]. When only the presence/absence of links are considered, nested structures are typically observed. However, when link-strengths are considered, models with body size constraints do not result in structural predictions statistically different from models without these constraints (electronic supplementary material, table S3). We do not suggest that body size is unimportant; however, our results show that it may not be the most significant contributor of quantitative (weighted) structure.

A detailed analysis of food webs where both the strength and variability of links can be quantified will enable a more complete understanding of the primary drivers of food-web structure. Importantly, the construction of predator–prey networks from the flow of biomass given by the isotopic distributions of interacting species is a framework that permits analysis of systems that are difficult to observe, and perhaps more importantly is extendable to historical or palaeontological communities, thereby increasing the temporal range over which community structure can be measured. In the future, an exploration of trophic interaction structure across space and time will be vital for understanding its origin and function, as well as its sensitivity to environmental or climatic perturbations.

## Acknowledgments

We thank J. Bascompte, D. Costa, N. J. Dominy, J. A. Estes, T. Gross, J. B. Hopkins, S. Kim, M. Mangel, J. W. Moore, S. D. Newsome, M. M. Pires, L. Rudolf, A. O. Shelton and P. V. Wheatley for many helpful comments and discussions. This work was funded by an Institute of Geophysics and Planetary Physics research grant to J.D.Y., a UC-Santa Cruz Committee-On-Research grant to P.L.K., a National Science Foundation Graduate Student Fellowship (NSF-GRF) to J.D.Y. and a Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) grant to P.R.G.

- Received June 15, 2012.
- Accepted June 29, 2012.

- This journal is © 2012 The Royal Society