Layout 1 INTRODUCTION Food web models provide a toolkit to study the relationship between species (or trophic groups) and communities. In order to better understand this relationship, structural (Dunne et al., 2002) and dynamical (Jordán et al., 2008) sensitivity analysis can be performed on food web models. The global (i.e., community-wise) effects or local (i.e. species-specific) disturbances can be simulated and quantified in several ways. These models can help to better understand several actual threats to aquatic ecosystems, including the appearance of introduced and invasive species, overfishing and local extinction of species. In all these cases, there is a major change concerning one focal species and effects spread out in the interaction network, influencing directly and indirectly many other species. In order to assess the effects on the whole community, we need general, predictive modelling tools. This way, we can better understand keystone species, indirect relationships, functional diversity and the vulnerability of communities. A central question is how do the effects of single- species perturbations (deletions) cascade across the network? The community effects of individual species may depend on their network position (Pimm, 1980): this can be quantified by several network metrics (ranging from the number of neighbours to trophic height), quantifying how richly they are connected to others and where they are in the network. Our key question is how to compare structural and dynamical predictions on the role individual species play in food webs. METHODS We consider here a food web simulation model, based on differential equations, where the community response is measured in various ways, following the disturbance of particular species. The model is as simple as possible, so we can focus on the pure differences between response functions, i.e. how to best measure the multi-species effects of single-species perturbations. Our modelling framework is very general, applicable to ecological systems where trophic interactions dominate community dynamics (for example, open-water systems with a relatively low intensity of competition and non-trophic effects). Generating networks and a dynamical simulation model We investigated 1000 randomly generated networks. By using some restrictions, they are similar and therefore comparable with a previously analysed network (Móréh et al., 2018), which is based on a real food web (Jordán et al., 2009; Lin et al., 2004). The constraints and the common features were the following: all webs contain 15 trophic groups (3 producers and 12 consumers); there is constant number (36) of links between them; the number of top predators and trophic levels were maximized (maxtop ≤ 3, maxTL ≤ 4). This a simple model, suitable for dynamical analysis, representing basically a food web with low trophic resolution at lower levels and higher trophic resolution at higher levels. It remains for future studies to see how the aggregation procedure (i.e., describing microorganisms at higher details, see D’Alelio et al., 2016, 2019; Jordán et al., 2018) can possibly change the results. Beyond these structural constraints our other criterion was the dynamical stability of the networks, and their robustness against perturbations. Since the investigation of extinctions was out of our purposes, we looked for networks where all groups coexist during the whole ARTICLE Comparing community response indices in aquatic food web models Ágnes Móréh,1,2 Ferenc Jordán1,2,3* 1Balaton Limnological Institute, Centre for Ecological Research, Tihany, Hungary; 2Evolutionary Systems Research Group, Centre for Ecological Research, Tihany, Hungary; 3Stazione Zoologica Anton Dohrn, Naples, Italy ABSTRACT Aquatic ecosystems face several major challenges from the introduction and invasion of species, to overfishing. In order to better manage these situations, we need predictive models, where diverse scenarios can be simulated and tested. One key challenge to address is how to quantify the relationships between single-species disturbances and their multispecies effects. Mapping the spread of direct and indirect effects in food webs helps to link species to communities. Since food webs are complex networks of interactions, it is typically not easy to make predictions, so modelling and simulation may help to reveal general patterns. In food web simulations, one can quantify the effects of local perturbations on other species, i.e., community response. This may provide information about the relative importance of individual species and it is also useful to assess the vulnerability of the whole com- munity to local changes. However, community response can be measured in several ways and various response functions give dif- ferent results. In order to better understand their similarities and differences, we present a comparative study on a reasonable set of community response functions in food web simulations. These results contribute to build more predictive, multi-species models for systems-based conservation and management. Non -co mmerc ial us e o nly Aquatic food web models 95 simulation process, none of them extinct even in the course of perturbations. We modelled the dynamic behaviour of the networks in the same way as we did in Móréh et al. (2018). The dynamics can be described as follows: (eq. 1) where Bi means the biomass of species i, the consumption of the consumers can be characterised by a Holling type- III functional response (h = 2), where ωi is the relative consumption rate of species i when consuming ρ, B0 means the half-saturation density. The relative consumption rates are set to be proportional to the number of preys (ωi = 1/n). The increase function of the basal species (i ≤ 3) is described by the logistic growth model, where ri and Ki are the intrinsic growth rate and the carrying capacity, respectively. The external mortality rates (di) of the basal species are set to 0. The increase and decrease of the consumer species (i > 3) depend only on the intensity of consumption of their predators and their preys (ri = Ki = 0), and the external mortality rate, which is set to 0.1 for all consumers. Since our purpose was to focus on the impact of network topology on dynamics, so we didn’t model the consumption and conversion rates explicitly, but assumed the strength of a predator-prey link (ε) proportional to the number of preys (ε = 1/n). For the sake of the same cause, the parameters (even the mortality rates) were fixed, only the different distributions of the links between the groups determinate the topology, thus the dynamics, too. The process of dynamical simulations had 3 main parts. First, the integration of the system (Hindmarsh et al., 2005) until it settled to a fixed point; if any of the species is extinct, the integration was terminated and the adjacency-matrix was casted off. Second, if all species reached a fixed point and coexisted, we made a preliminary stability analysis on the system by arbitrary changing the biomass-values of all species. After the system returned to the original equilibrium after this pulse perturbation event, we made a systematic perturbation process on all consumers (the producers were part of the dynamical system but their community effects were not evaluated). We changed the mortality rate (di) of the species in question increasing it by 10%. This analysis can be considered as a press perturbation experiment (sensu Bender et al. (1984)). Third, if the system is robust against all perturbation events (there is no extinction), not only the adjacency- matrix of the network, but also the biomass-values registered before and after all perturbation events are hold and stored. Community response indices We are interested in the effect of perturbing species i on the ecosystem including it, thus we have to make a relationship between the original and the perturbed state of the system. However, only one species is perturbed, all others’ state changes; the whole system is settled to a new equilibrium state. In our modelling framework, we consider the biomass of species but also abundance could have been used (as another, measurable, quantitative trait). All species “answer” to the perturbation event by its own biomass-change; thus, the community response (CR) to the perturbation of the ith species (CRi) can be determined as the sum of all these answers. There is not a single formula used to determine the response of community to any kind of perturbations, but more different approaches exist and are used in parallel. The investigated types and their variants are summarized in Tab. 1. “Functional importance index” (FIi) is defined by Hurlbert (1997), and calculates the summation of the differences of biomass-values measured before and after perturbation. This fact discriminates this index from all others discussing below, since the difference gives information about the net changes in a numerical way, but gives less information about the order of magnitude of the changes. All indices discussed further calculate the rate of biomasses, which, however, shows the measure of the changes, not the exact numerical values of them. “Interaction strength index” (ISIi) used by Okey (2004) calculates simply the sum of the rates of biomass- values after and before perturbation. Thus, this index takes not only the measure of the changes into account, but inherently the direction of the changes, too. “Paine’s interaction strength” (PIi, (Paine, 1992)) can be derived from Okey’s interaction strength index by subtracting 1 from the rates calculated in ISIi. In fact, these indices are complements of each other: while the ISIi gives the exact rate of the biomass-values, the PIi shows the measure of changes. In addition, PIi enables to ignore the sign of changes by calculating absolute values (Tab. 1). The next two indices incorporate a further aspect, namely the concept of keystone species considering its core meaning: “A keystone species is one whose effect is large, and disproportionally large relative to its abundance” (Power et al., 1996). The index named “keystoneness” (KNi) is also derived from ISIi divided it by the percent of the system’s overall living biomass represented by group i before it was perturbed (Okey, 2004). The “community importance index” (CIi) also provides an index of the impact of a species relative to its biomass Non -co mmerc ial us e o nly Á. Móréh, and F. Jordán96 (Power et al., 1996). It calculates the sum of the opposite (sign) of the real value of the change of all species divided by the percent of the overall biomass of the groups represented by species i before it was perturbed. Although all discussed indices are based on the comparison of the biomass-values before and after the perturbation event, and the summation of the changes of each groups, there are two important points of view, that discriminate the four variants of them (columns in Tab. 1): i) taking the sign of the biomass-changes of every single species into account or not; ii) taking the self-effect of the perturbed species into account or not (this means feedback: if a particular species is disturbed, effects spread in the network and indirect chain effects can influence the disturbed species later again, looping back to it). An additional distinction between the different types of CR- indices is, whether they incorporate the original biomass of the perturbed species (considering the keystone- concept), or not. Structural network indices In order to quantify the position of individual species in food webs, we calculated 18 structural importance indices that are summarized in Tab. 2. For detailed description of the indices see Móréh et al. (2018). Although all of the 18 indices express the positional importance in some wise, they differ from each other, too. There are local indices not taking indirect effects into account (e.g., D), while there are meso-scale or non-local indices (e.g., centralities). The topological importance indices consider binary interactions, others (e.g., WI) can quantify weighted webs. While some of them (e.g., s, K) characterise directed networks, others do undirected ones (e.g., D). Correlation analysis We used the non-parametric Spearman’s rank correlation test to investigate the connection between the different CR-functions. Beyond that we studied the correlations of CRs both with the measures of self-effects and the different structural network indices. Its purpose was to investigate the consistency of the relationships between the indices and the different types and variants of CRs. Since we were interested in the strength of correlations, not the sign of them, we used the absolute values of the ρ-values for further analysis. We considered a connection strong if |ρ| ≥ 0.5. RESULTS Correlation between the different response functions In Tab. 3, we summarize not only the strengths of correlations (|ρ|) between the different response functions, but also between the measure of self-effects or original biomass-value of the perturbed species and the CRI-types. We considered the correlation stronger if |ρ| ≥ 0.5. Investigating the relationships between the different variants within the same type we can see that the indices calculated with or without self-effect (in grey in Tab. 3) Tab. 1. Summary of all investigated indices: the main types re shown in rows, while the different variants are collected in columns. Bi bef and Bi aft denote the biomass values before and after the perturbation event; %Bi bef is the percentage of the affecting consumer in the whole system before it was perturbed (%Bbef = (Bi/Σn 1Bj)x100,n = 15). The original form of each types used in the literature are highlighted with grey background; variants calculated with self-effect (i ∈ N, where N = {1...15}) are marked with star. sign ignored with sign References i ∉ N i ∈ N i ∉ N i ∈ N Functional importance index (FIi) |FIi|* FIi* Hurlbert (1997) Livi et al. (2011) Interaction strength index (ISIi) - - ISIi* Okey (2004) Paine’s interaction strength index (PIi) |PIi|* PIi* Paine (1992) Móréh et al. (2018) Keystoneness (KNi) - - KNi* Okey (2004) Community importance index (CIi) |CIi|* CIi* Power et al. (1996) Okey (2004) Non -co mmerc ial us e o nly Aquatic food web models 97 show stronger correlations in most cases (|FIi| vs. |FIi|* shows the strongest, FIi vs. FIi* shows the weakest correlation). In turn, indices taking the sign of changes into account or ignoring them are much less correlated (framed values in Tab. 3). Investigating the connections between the different types Tab. 2. Summary of the structural network indices. Index name Notation Reference Degree D Wassermann and Faust (1994) Weighted degree wD Centrality indices Wassermann and Faust (1994) Betweenness centrality BC Closeness centrality CC Keystone index and its components Jordán et al. (1999) Keystone index K Bottom-up effects Kbu Top-down effects Ktd Direct effects Kdir Indirect effects Kindir Status index and its components Harary (1959) Status s Contra-status s’ Delta-status Δs Positional importance based on indirect effects Jordán et al. (2003) Topological importance, max step = 1 TI1 Topological importance, max step = 3 TI3 Topological importance, max step = 5 TI5 Weighted topological importance max step = 1 WI1 Weighted topological importance max step = 3 WI3 Weighted topological importance max step = 5 WI5 Tab. 3. The strengths of Spearman-correlations (|ρ|) between the different CR-types and variants. Bold values mean higher degrees of correlations (|ρ| ≥ 0.5). Relationships between the variants calculated with or without self-effect within the same type are marked with grey background. Values of correlations between indices taking the sign of changes into account or ignoring them (within the same type) are framed. |FIi| |FIi|* FIi FIi* ISIi ISIi* |PIi| |PIi|* PIi PIi* KNi KNi* CIi CIi* |CIi| |CIi|* |FIi| 1 0.98 0.32 0.47 0.31 0.16 0.93 0.84 0.31 0.16 0.65 0.66 0.23 0.36 0.39 0.22 |FIi|* - 1 0.43 0.46 0.38 0.14 0.93 0.89 0.38 0.14 0.65 0.66 0.3 0.34 0.4 0.26 FIi - - 1 0.35 0.58 0.06 0.36 0.58 0.58 0.06 0.1 0.09 0.55 0.01 0.6 0.66 FIi* - - - 1 0.07 0.14 0.45 0.27 0.07 0.14 0.77 0.76 0.17 0.36 0.33 0.44 ISIi - - - - 1 0.67 0.27 0.39 1.0 0.67 0.09 0.1 0.94 0.57 0.26 0.3 ISIi* - - - - - 1 0.12 0.04 0.67 1.0 0.32 0.32 0.53 0.92 0.23 0.31 |PIi| - - - - - - 1 0.9 0.27 0.12 0.67 0.67 0.18 0.34 0.46 0.26 |PIi|* - - - - - - - 1 0.39 0.04 0.44 0.45 0.36 0.12 0.64 0.54 PIi - - - - - - - - 1 0.67 0.09 0.1 0.94 0.57 0.26 0.30 PIi* - - - - - - - - - 1 0.32 0.32 0.53 0.92 0.23 0.31 KNi - - - - - - - - - - 1 1.0 0.03 0.59 0.31 0.46 KNi* - - - - - - - - - - - 1 0.02 0.59 0.3 0.46 CIi - - - - - - - - - - - - 1 0.44 0.32 0.4 CIi* - - - - - - - - - - - - - 1 0.28 0.41 |CIi| - - - - - - - - - - - - - - 1 0.94 |CIi|* - - - - - - - - - - - - - - - 1 Self-effects |diffi| 0.74 0.84 0.68 0.28 0.5 0.06 0.74 0.86 0.5 0.06 0.46 0.47 0.45 0.19 0.42 0.41 diffi 0.74 0.84 0.68 0.28 0.5 0.06 0.74 0.86 0.5 0.06 0.46 0.47 0.45 0.19 0.42 0.41 ratei 0.42 0.53 0.84 0.16 0.53 0.14 0.41 0.71 0.53 0.14 0.09 0.08 0.57 0.16 0.68 0.79 |ratei-1| 0.42 0.53 0.84 0.16 0.53 0.14 0.41 0.71 0.53 0.14 0.09 0.08 0.57 0.16 0.68 0.79 ratei-1 0.42 0.53 0.84 0.16 0.53 0.14 0.41 0.71 0.53 0.14 0.09 0.08 0.57 0.16 0.68 0.79 %Bi bef 0.66 0.66 0.09 0.76 0.11 0.33 0.67 0.45 0.11 0.33 1.0 1.0 0.01 0.6 0.3 0.46 Non -co mmerc ial us e o nly Á. Móréh, and F. Jordán98 and variants of CRIs (indices possibly considering sign and or self-effect), we can see that the correlation can be stronger if the compared indices take sign into account inherently (ISIi, KNi) or are calculated without absolute value (PIi, CIi). In case of indices can be calculated with or without sign (FIi, PIi, CIi), the correlation is stronger if sign is ignored (e.g., |FIi| vs. |PIi|) regardless of ignoring self-effects or not. The question is whether an index uses the difference or the rate of biomass-values measured before and after the perturbation event. FIi-indices use the difference of these, thus, it is strongly dependent on the original values of not only the perturbed (Bi), but also the all other species’ biomass. FIi (calculated with sign and with self- effect) is and exception: its correlation is very weak with the original percentage of the perturbed species’ biomass. Using the rate of the biomasses makes the changes independent of the real biomass-values, these indices emphasizes the measure of the changes. |FIi|* and |PIi|* show the strongest correlations with any types of self- effects, while the CIi-type indices correlate stronger only with the rate-types of self-effects (excepting CIi*, which is, however, in stronger correlation with the original percentage of the perturbed species’ biomass). While ISIi shows a moderate correlation with the strength of self- effects, dividing it by %Bi bef deteriorates this connection by including the dependency of the original amount of the perturbed species (KNi). The difference between the variants within the same type can be enlightened also by comparing the 12,000 original values of each types (1000 networks, 12 perturbation events). Fig. 1A shows the difference of the indices calculated with and without self-effect. The largest impact of self-effect can be seen in cases of KNi and ISIi, while the CIi-indices show smaller difference in this regard. The variants with self-effect are almost always greater than the variants ignoring it, except FIi and PIi, where this difference is always negative (marked with grey boxes in Fig. 1). Investigating the differences between the variants calculated with or without taking sign into account, we get that the variants without sign are always greater than the variants calculated with absolute value (Fig. 1). It means that using sign always underestimates the strength of perturbation effects because of the changes with opposite sign cancelling each other out. Correlation between the different CR- and structural network indices The results of the comparison of the different CRI- values and the structural indices of each investigated networks are shown in Fig. 2. First of all, there is a conspicuous difference between two groups of network indices. While the centrality indices (CC, BC), the unweighted versions of topological importance indices (TIn), the keystone index (K) and its components referring the direct or indirect effects for a node (Kdir, Kindir) show very weak connection with almost all types of CR-indices, the status-indices (s, s’, ∆s), the weighted versions of topological importance indices (WIn) and the degree (wD), and the components of keystone index quantifying the bottom-up and the top-down effects of a node (Kbu, Ktd) show higher average correlation strength; in case of some particular indices the connection is quite strong. Examining the different response-indices more carefully we can see that the types ignoring not only the sign of changes but also the original percentage of the perturbed species (|FIi|, |FIi|*, |PIi|, |PIi|*) show the strongest correlations with all network indices belonging the group that show higher average correlation with the response-indices. In this regard, the self-effect of the Fig. 1. The results of the comparison of the variants of the dif- ferent CR-types. A) Difference of types calculated with self-ef- fect or ignoring it; the boxplots show the absolute values of the difference; the originally negative values are marked by grey colour. B) Difference of variants calculated with or without sign. All values are positive. Note, that KNi and ISIi have been calcu- lated with sign. Non -co mmerc ial us e o nly Aquatic food web models 99 perturbed species does not matter. Taking sign into account the strength of correlations is rather moderate, or weak in some cases, and the differences also are much smaller. Moreover, ignoring not only the keystone concept (%Bbef) but also the self-effects causes a bit stronger correlation. When sign is ignored, calculating with %Bbef (|CIi| and |CIi|*) can increase the correlation values to a bit higher level in case of some network indices referring also to some aspects of keystoneness (e.g. K, Kdir), but these strengths remain still rather weak. In contrast with that, in case of status indices the strength of correlations is much higher, albeit not as strong as ignoring %Bbef. All in all, in case of the structural indices showing stronger correlations, the strength is influenced much more by sign than self-effect. However, the difference of variants with or without self-effect is not as unambiguous as in case of the sign, the averages of correlation-values in cases where self-effects are ignored are higher. DISCUSSION Considering the powers and the limits of these different formulas of community response functions we can choose the one that suits our purposes best. If our study is confined only the investigation of a single food web, in respect of the way of comparison it does not matter whether we use a formula using the difference or the rate of the biomass-values. Although neither the exact numerical difference nor the rate of change can handle the two important point of view of a change together: a large difference doesn’t mean large rate of change if the original biomass is high enough, in the same way, small difference can give both small and large rate depending on the original values. Despite that in case of studying more networks it might be more useful comparing only the magnitude of changes between the different networks increasing the comparability of the results. Likewise, the inclusion of self-effect in the sum is Fig. 2. The strength of correlation (|ρ|) between the studied network indices and the different CR-types grouped into columns based on the three investigated point of view, whether the sign of changes, the self-effect of the perturbed species and the keystone concept are considered or not. The network indices (rows) are sorted by descending the average correlation strength between all CR-indices and the structural index in question. Non -co mmerc ial us e o nly Á. Móréh, and F. Jordán100 determined by the current question. In case of rate-type CR-functions the correlations between the two variants are very good. In case of functional importance indices there can be larger differences between the results calculated by the two ways because of the higher dependency of the original biomass-values of the perturbed species (that is supported by the higher correlations with the percentage of it; Tab. 3). If the investigation of the connections between the species’ structural positions and any changes in a trait of them (say the biomass) is among the purposes of a study, including or ignoring self-effects has no particular impact on the strength of correlations; although the values without self-effects are higher in general, the differences are very small. All in all, from this point of view, the question of ignoring sign is much more important than ignoring self- effect, or even more calculating the difference or the rate of biomasses. Calculating with sign underestimates the net measure of changes because of changes with opposite sign can cancel out each other. Theoretically, it is possible that while a corresponding CR-index gives the result of 0, the whole structure of the ecosystem changes. Since there are indices, that inherently take sign into account (ISIi, KNi), it is better to use a form with an equivalent “structure” (e.g., PIi instead of ISIi). This problem can also lead up to different conservation biology approaches. Comparing small negative and large positive effects, we can understand if it is better to minimize our impact on nature (small negative better than large positive) or to implement positive effects (large positive better than small negative). Choosing a response function that does not provide information about increase or decrease of population size, only about the change in population size is supported by a conservation philosophy suggesting that minimizing the human impact (size of change) might be preferred over trying to help natural systems (direction of change). CONCLUSIONS The community-wide effects of single-species perturbations can be modelled, simulated and measured. Mapping the spread of effects in food webs is possible from several approaches, including a general bottom-up vs top-down comparison (Stäbler et al. 2019), in single- node vs pairwise perturbation simulation studies (Móréh et al. 2018) or in single-node vs multi-node structural analyses (Jordán et al. 2019). Linking individual species to communities may also consider reproductive systems (D’Alelio 2017), trait-based aggregation protocols (Jordán et al. 2018) or the relationship between carbon fluxes and interaction strengths (Scotti et al. 2012). Better databases, including long-term ecological data (D’Alelio et al. 2016) provide a robust empirical basis for ecosystem modelling (D’Alelio et al. 2019) and all these approaches may contribute to multi-species, system-based conservation efforts and fisheries management. Some results presented in this paper contributes to more predictive modelling. Based on Tab. 3, we see, for example, that simulating the effects of overfishing does not need to consider self-effects (feedback) if the sign (good or bad effect) of the interaction is not considered, only its magnitude (big or small effect). If signs are considered, however, calculating self-effects can be important for predictability. Fig. 2 suggests that the best structural predictions can be achieved when neither effect sign nor initial biomass are considered. In this case, weighted and possibly also indirect structural indices are the best ones to predict systems dynamics. These hints help to offer predictive simulation models for fisheries in aquatic ecosystems. ACKNOWLEDGEMENTS We thank Sándor Imre Piross and Anett Endrédi for useful advices. The research was supported by the National Research, Development and Innovation Office (NKFIH) grants OTKA K 116071 and GINOP-2.3.2-15- 2016-00057. Corresponding author: jordan.ferenc@gmail.com Key words: Food web; keystone species; perturbation; community response. Received: 16 October 2019. Accepted: 23 December 2019. This work is licensed under a Creative Commons Attribution NonCommercial 4.0 License (CC BY-NC 4.0). ©Copyright: the Author(s), 2019 Licensee PAGEPress, Italy Advances in Oceanography and Limnology, 2019; 10:8621 DOI: 10.4081/aiol.2019.8621 REFERENCES Bender EA, Case TJ, Gilpin ME, 1984. Perturbation experiments in community ecology: theory and practice. Ecology 65:1-13. D’Alelio D, Libralato S, Wyatt T, Ribera d’Alcalà M. 2016. Ecological-network models link diversity, structure and function in the plankton food-web. Sci. Rep. 6:21806. D’Alelio D, Hay Mele B, Libralato S, Ribera d’Alcalà M, Jordán F. 2019. Rewiring and indirect effects underpin modularity reshuffling in a marine food web under environmental shifts. Ecol. Evol. 9:1-16. Non -co mmerc ial us e o nly Aquatic food web models 101 Dunne JA, Williams RJ, Martinez ND, 2002. Network structure and biodiversity loss in food webs: robustness increases with connectance. Ecol. Lett. 5:558-567. Harary F, 1959. Status and contrastatus. Sociometry 22:23. Hindmarsh AC, Brown PN, Grant KE, Lee KE, Serban R, Shumaker DE, Woodward CS, 2005. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Softw. 31:363-396. Hurlbert SH, 1997. Functional importance vs. keystoneness: Reformulating some questions in theoretical biocenology. Austral Ecol. 22:369-382. Jordán F, Takács-Sánta A, Molnár I, 1999. A reliability theoretical quest for keystones. Oikos 86:453-462. Jordán F, Liu W, van Veen F, 2003. Quantifying the importance of species and their interactions in a host-parasitoid community. Comm. Ecol. 4:79-88. Jordán F, Okey TA, Bauer B, Libralato S, 2008. Identifíing important species: Linking structure and function in ecological networks. Ecol. Mod. 216:75-80. Jordán F, Liu W, Mike Á, 2009. Trophic field overlap: A new approach to quantify keystone species. Ecol. Mod. 220:2899-2907. Jordán F, Endrédi A, Liu WC, D’Alelio D. 2018. Aggregating a plankton food web: mathematical versus biological approaches. Mathematics 6:336. Jordán F, Pereira J, Ortiz M. 2019. Mesoscale network properties in ecological system models. Curr. Opin. Syst. Biol. 13:122- 128. Lin HJ, Shao KT, Hwang JS, Lo WT, Cheng IJ, Lee LH, 2004. A trophic model for Kuosheng Bay in Northern Taiwan. J. Mar. Sci. Technol. 12:424-432. Livi CM, Jordán F, Lecca P, Okey TA, 2011. Identifying key species in ecosystems with stochastic sensitivity analysis. Ecol. Mod. 222:2542-2551. Móréh Á, Endrédi A, Jordán F, 2018. Additivity of pairwise perturbations in food webs: topological effects. J. Theor. Biol. 448:112-121. Okey TA, 2004. Shifted community states in four marine ecosystems: some potential mechanisms. PhD Thesis, University of British Columbia, Vancouver. Paine RT, 1992. Food-web analysis through field measurement of per capita interaction strength. Nature 355:73-75. Pimm SL, 1980. Food web design and the effect of species deletion. Oikos 35:139-149. Power ME, Tilman D, Estes JA, Menge BA, Bond WJ, Mills LS, Daily G, Castilla JC, Lubchenco J, Paine RT, 1996. Challenges in the quest for keystones. BioScience 46:609-620. Scotti M, Gjata N, Livi CM, Jordán F. 2012. Dynamical effects of weak trophic interactions in a stochastic food web simulation. Community Ecol. 13:230-237. Stäbler M, Kempf A, Smout S, Temming A, 2019. Sensitivity of multispecies maximum sustainable yields to trends in the top (marine mammals) and bottom (primary production) compartments of the southern North Sea foodweb. PLoS One 14:e0210882. Wassermann S, Faust K, 1994. Social network analysis: Methods and applications. Cambridge University Press, Cambridge: 825 pp. Non -co mmerc ial us e o nly