Interaction modifications lead to greater robustness than pairwise non-trophic effects in food webs

Considerable emphasis has been placed recently on the importance of incorporating non-trophic effects in to our understanding of ecological networks. Interaction modifications are well established as generating strong non-trophic impacts by modulating the strength of inter-specific interactions. For simplicity and comparison with direct interactions within a network context, the consequences of interaction modifications have often been described as direct pairwise interactions. The consequences of this assumption have not been examined in non-equilibrium settings where unexpected consequences of interaction modifications are most likely. To test the distinct dynamic nature of these ‘higher-order’ effects we directly compare, using dynamic simulations, the robustness to extinctions under perturbation of systems where interaction modifications are either explicitly modelled or represented by corresponding equivalent pairwise non-trophic interactions. Full, multi-species representations of interaction modifications resulted in a greater robustness to extinctions compared to equivalent pairwise effects. Explanations for this increased stability despite apparent greater dynamic complexity can be found in additional routes for dynamic feedbacks. Furthermore, interaction modifications changed the relative vulnerability of species to extinction from those trophically connected close to the perturbed species towards those receiving a large number of modifications. Future empirical and theoretical research into non-trophic effects should distinguish interaction modifications from direct pairwise effects in order to maximise information about the system dynamics. Interaction modifications have the potential to shift expectations of species vulnerability based exclusively on trophic networks.

1. Considerable emphasis has been placed recently on the importance of incorporating non-trophic 11 effects in to our understanding of ecological networks. Interaction modifications are well 12 established as generating strong non-trophic impacts by modulating the strength of inter-13 specific interactions. 14 2. For simplicity and comparison with direct interactions within a network context, the 15 consequences of interaction modifications have often been described as direct pairwise 16 interactions. The consequences of this assumption have not been examined in non-equilibrium 17 settings where unexpected consequences of interaction modifications are most likely. 18 3. To test the distinct dynamic nature of these 'higher-order' effects we directly compare, using 19 dynamic simulations, the robustness to extinctions under perturbation of systems where 20 interaction modifications are either explicitly modelled or represented by corresponding 21 equivalent pairwise non-trophic interactions. 22 4. Full, multi-species representations of interaction modifications resulted in a greater robustness 23 to extinctions compared to equivalent pairwise effects. Explanations for this increased stability 24 despite apparent greater dynamic complexity can be found in additional routes for dynamic 25 feedbacks. Furthermore, interaction modifications changed the relative vulnerability of species 26 to extinction from those trophically connected close to the perturbed species towards those 27 receiving a large number of modifications. 28 5. Future empirical and theoretical research into non-trophic effects should distinguish interaction 29 modifications from direct pairwise effects in order to maximise information about the system 30 dynamics. Interaction modifications have the potential to shift expectations of species 31 vulnerability based exclusively on trophic networks. 32

Introduction 36
There is a building appreciation that to improve our understanding of the population dynamics of 37 ecological communities it is necessary to move beyond studies that focus on a single interaction 38 process at a time (Kéfi et al., 2012;Levine, Bascompte, Adler, & Allesina, 2017). Trophic interaction 39 modifications (Terry, Morris, & Bonsall, 2017;Wootton, 1993) are class of higher-order process 40 where a consumer-resource interaction is modulated by additional species. Examples include 41 associational defences (P. Barbosa et al., 2009), fear effects (Sih, Englund, & Wooster, 1998), certain 42 impacts of ecosystem engineers (Sanders et al., 2014) and foraging choices (Abrams, 2010). It has 43 been empirically demonstrated that many strong non-trophic effects are caused by such processes 44 (Ohgushi, Schmitz, & Holt, 2012;Preisser, Bolnick, & Benard, 2005;Werner & Peacor, 2003) with 45 large implications for community structure and dynamics (M. Barbosa,Fernandes,Lewis,& Morris,46 2017; Donohue et al., 2017;Matassa, Ewanchuk, & Trussell, 2018; van Veen, van Holland, & Godfray, 47 2005). 48 Approaches to understanding interaction modifications often try to distil the inherently multi-49 species process into a pairwise non-trophic effect (or 'trait-mediated indirect effect') from the 50 modifier species onto one or both of the recipient species (Okuyama & Bolker, 2007). This allows the 51 direct comparison of non-trophic and trophic interaction strengths (Preisser et al., 2005) and 52 network structure (Pilosof, Porter, Pascual, & Kéfi, 2017) but is a representation of a different class 53 of dynamic process (Terry et al., 2017). While this simplification is fully justified in studies of 54 communities at equilibrium (Grilli, Barabás, Michalska-Smith, & Allesina, 2017), the consequences of 55 this assumption in a transient, fluctuating or heavily perturbed system has yet to be fully explored. 56 Previous studies introducing interaction modifications to trophic networks (Arditi, Michalski, & 57 Hirzel, 2005;Goudard & Loreau, 2008;Lin & Sutherland, 2014) have demonstrated their potential 58 impact on the dynamics of ecosystems, but not whether this is attributable to the higher-order 59 nature of interaction modifications, or rather shifts in connectance and interaction strength. 60 One such important case is the dynamics of ecological systems in the face of species removal, where 61 there is the potential for secondary extinctions and eventually the collapse of the ecosystem (Dunne, 62 Williams, & Martinez, 2002). This aspect of stability, often described as 'robustness', is an important 63 question both from the perspective of managing anthropogenic change and in terms of 64 understanding the fundamental stability of ecological communities. Since empirically testing how 65 whole communities respond to extinctions is difficult or impossible (although see Sanders et al. 66 (2018)), a number of studies have attempted to determine the properties that make ecological 67 communities robust through simulation (Dunne & Williams, 2009;Säterberg, Sellman, & Ebenman, 68 2013). However, incorporating the acknowledged flexibility of ecological networks is a perennial 69 challenge for such studies (Montoya, Pimm, & Solé, 2006). 70 The impact on robustness in ecological networks of one specific subset of interaction modifications, 71 those caused by the flexible foraging in response to resource availability, has been examined by a 72 number of studies (Valdovinos, Ramos-Jiliberto, Garay-Narváez, Urbani, & Dunne, 2010). Approaches 73 have included topological rewiring (Gilljam, Curtsdotter, & Ebenman, 2015;Kaiser-Bunbury, Muff, 74 Memmott, Müller, & Caflisch, 2010;Staniczenko, Lewis, Jones, & Reed-Tsochas, 2010;Thierry et al., 75 2011), multi-species functional responses (Uchida & Drossel, 2007) and adaptive dynamics models 76 (Kondoh, 2003). These models showed that the additional dynamic process impacted robustness, 77 but only addressed a restricted subset of interaction modifications caused by predator switching. 78 However, consumption rates are influenced by more than the just the prey available to them. 79 Interaction modifications can also be caused by the threats to the consumer (Suraci,Clinchy,Dill,80 Roberts, & Zanette, 2016), facilitation (Bruno, Stachowicz, & Bertness, 2003), associational 81 susceptibility (Underwood, Inouye, & Hambäck, 2014) or mutualistic defence (Holland, Ness, Boyle, 82 & Bronstein, 2005) amongst others (Ohgushi et al., 2012). This introduces considerable additional 83 dynamic connectance to ecological communities, yet studies of generic interaction modifications in 84 large networks are limited (Arditi et al., 2005;Bairey, Kelsic, & Kishony, 2016;Garay-Narváez & 85 Ramos-Jiliberto, 2009;Goudard & Loreau, 2008;Lin & Sutherland, 2014), with most theoretical 86 analyses of interaction modifications focussing on small community units (Abrams, 2010;Bolker, 87 Holyoak, Křivan, Rowe, & Schmitz, 2003;Holt & Barfield, 2013). 88 Calls to incorporate the full panoply of non-trophic effects into our understanding of ecological 89 networks have built substantially in recent years (Fontaine et al., 2011;Ings et al., 2009;Kéfi et al., 90 2012; Levine et al., 2017;Ohgushi et al., 2012;Olff et al., 2009) and the first empirical inventories are 91 being established (Kéfi et al., 2015;Kéfi, Miele, Wieters, Navarrete, & Berlow, 2016). Theoretical 92 analyses can play a significant role in guiding the empirical creation of networks, identifying the 93 information necessary to best understand these systems. Here we demonstrate the distinctive 94 nature of interaction modifications compared to pairwise non-trophic effects through a direct 95 standardised comparison of their impacts on the robustness of large artificial networks. We then 96 examine the role of interaction modifications in determining the consequences of species loss and 97 species level vulnerability to secondary extinction. 98

Methods 99
We conducted our analyses using model food webs where system dynamics are derived from 100 metabolic scaling relationships. As detailed below, interaction modifications were introduced to 101 communities at an initial equilibrium before mortality was added to a single species at a time and 102 the resultant extinctions after the model was integrated to a new equilibrium examined. 103

Bio-energetic Model 104
The change in biomass density, , of each species in the community was modelled using a simple 105 Lotka-Volterra type model with a linear (Holling Type I) functional response and logistic intrinsic 106 growth rates, parameterised using body-mass relationships: 107 Each species was assigned a body mass ( ) drawn from a distribution based on their trophic level 109 (see SI 1) which was then used to calculate further parameters using quarter-power body-mass 110 scaling laws (Yodzis & Innes, 1992). Relative intrinsic growth or metabolic loss rates, , were set at 1 111 for all producers (trophic level 1) and = −0.1 −0.25 for each consumer (trophic level ≥2). 112 Consumer-specific consumption rates were set at = −0.25 , where the generality term was 113 set to 1/ , the number of resources of each consumer . Assimilation efficiencies, , for each 114 trophic interaction were drawn from a uniform distribution, ~ (0.05: 0.15). Carrying capacities, 115 , were drawn from ~(1: 10) for producers and ~10 (2:3) for consumers. Initially, external 116 mortality, , was set to 0 and modification terms, , set to 1. 117 Trophic topology and population densities of starting communities were determined as follows. 118 Initial trophic topologies were generated using the niche model (Williams & Martinez, 2000), with 35 119 species and a connectance of 0.14, removing cannibalistic interactions. Each population density was 120 initially set to 10 and these were then numerically integrated to a stable equilibrium (see SI 4 for 121 criteria). Only connected communities with at least 18 persisting species were retained and used as 122 starting communities for the robustness tests. Properties of the starting communities are described 123 in SI 2. 124

Specification of TIMs 125
Trophic interaction modifications were introduced to this community through modification terms, 126 , that specify the impact of modifier species on the consumption of species by species , with 127 the constraint that ≠ ≠ . These are positive numbers that multiply the attack rate as a function 128 of the divergence of the biomass density of a modifying species, , from its starting equilibrium 129 value * . Where there are multiple modifications of the same interaction these were assumed to 130 combine multiplicatively = ∏ (Golubski & Abrams, 2011;Goudard & Loreau, 2008). The 131 TIMs were distributed at a defined rate of occurrence such that each potential modification 132 (combination of trophic interaction and third species) had an equal chance of existing. 133 When < 1, the interaction is weakened and when > 1 the interaction is strengthened. In 134 our model can not be negative, preventing the reversal of the direction of the trophic 135 interaction. Where an increase in the modifier species leads to an increase in the strength of the 136 interaction, we describe it as a facilitating TIM. This can be said to be beneficial to the consumer and 137 detrimental to the resource, in the sense of the immediate impact from an increase in the modifier 138 population. We term the reverse situation an interfering TIM. See Figure 1 for diagrammatic 139 explanation of these terms. 140

Alternative representation of TIMs as Pairwise Non-Trophic Effects 146
As an alternative to modelling the impact of trophic interaction modifications in full (using 'higher 147 order' terms), equivalent pairwise non-trophic effects (NTEs) can be derived. These match the 148 impact of the full TIM model, the only distinction being the NTE from a modifier to a trophic 149 interactor is no longer dependent on the biomass of the other member of the trophic pair ( Figure 2). 150 To maintain parity with the TIM model, this was done by first partitioning the process into trophic 151 and non-trophic components, then fixing the value of the other trophic interactor to the density at 152 the original equilibrium: = * . Full steps of the derivation are in SI 6. A trophic interaction 153 onto a resource influenced by a modifier can be partitioned from: 154 (2) 155 to: 156 The corresponding terms for the interactions onto the consumer, , are: 158 As an example, a facilitating TIM, when increased, would lead to a greater than 1. This would 160 lead to a negative NTE on the resource (independent of ) and a positive NTE on the consumer 161 (independent of ). This partitioning process is only straightforward when each trophic interaction 162 is modified by at most one modifier species because of the synergistic relationship between multiple 163 modifiers assumed by our model. 164

Functional Form of TIM 173
To represent the relationship between the density of the modifier species and the modification of 174 the interaction we used a Gompertz sigmoidal curve parameterised to control features of ecological 175 relevance, detailed in full in SI 3. This function links the divergence of from its start point, 176 log 10 ( / * ), and three control parameters, τ (the proportional change on a logarithmic scale from 177 the equilibrium to a threshold where the interaction is most greatly impacted by the change in 178 modifier), α (the direction and maximum rate of change in the modified interaction as the modifier 179 density increases) and σ (the range of magnitudes over which spans), depicted in Figure 3. A 180 positive α denotes a facilitatory modification and a negative α an interference modification. A key 181 attribute of our parameterisation is that when = * , = 1. This maintains the original trophic 182 interaction strength when the modifier is at its starting equilibrium point -in effect these strengths 183 are assumed to already incorporate changes to the interaction at the equilibrium. For each robustness test, the external mortality rate ( ) of a single species (to which we will refer 196 for convenience as the 'targeted' species, although the mortality could be attributable to a range of 197 non-directed processes) was then set to 1, the system numerically integrated to a new steady state 198 and the status of each species in the resultant community assessed. Species were considered 199 'extinct' if their biomass fell permanently below 10 −4 (several orders of magnitude below the 200 starting values of most of the populations, SI 2), 'functionally extinct' if their biomass fell to below 201 1/10 th of their starting value and considered to have 'exploded' if the final density was over 10 times 202 the starting value. This was repeated targeting each species in turn, for each community, to give a 203 total of 3736 completed robustness tests (93.8% successful integrations, SI 4). All analyses were 204 carried out in R v.3.5.0 (R Core Team, 2018) using the deSolve numerical integration package 205 . 206 We then repeated the robustness tests after conversion of the TIMs into pairwise NTE form, and in a 207 third case without any TIMs. In each case the underlying trophic interactions and distribution of 208 TIMs was identical. 209

Test 2. TIMs and distribution of extinctions 210
To examine the relationship between the distribution of TIMs and consequent secondary 211 extinctions, the robustness tests of the full TIM case as described above was repeated with a higher 212 occurrence rate of TIMs (0.08) and a relaxation of the assumption that only one modification can 213 affect each trophic interaction. Results are reported for parameters drawn from ~(−3: 3), 214 ~(0.1: 3), ~ (−1: 1). Further tests, with different TIM occurrence rates and distributions of 215 shape parameters, reached qualitatively similar results. Properties of the community and the 216 relationship between the 'targeted' species and the extinct species were then calculated. 217 The trophic distance was calculated as the shortest route, in terms of the number of trophic links, 218 between the targeted species and each secondarily extinct species. This was compared to the 219 distribution of trophic distances between targeted and secondarily extinct species in the tests 220 without TIMs and to the baseline distribution of trophic distances between all species in the 221 networks. The number and class (whether they are initially beneficial or detrimental for the focal 222 species) of the NTEs affecting each extinct species, both just from the targeted species and from all 223 species, was compared to the expected number of NTEs upon each species given its trophic 224 connectance. 225

Results 226
The introduction of TIMs greatly increased the number of extinctions (Figure 4)

248
Compared to the expected number of NTEs from TIMs of a particular sign for each species with a 249 given connectance (mean = 3.79), secondarily extinct species tended to be recipients of both more 250 beneficial (mean = 4.53) and more detrimental (mean = 4.02) NTEs. When counting just TIMs from 251 the targeted species, secondarily extinct species tended to have more beneficial (0.289) and fewer 252 detrimental (0.178) NTEs than the expected number (0.207) (t-tests paired by species ID, all 253 p<0.0001, n =3576). 254 The targeting of species that caused more TIMs induced more extinctions (Figure 6a

265
Discussion 266 Here we have examined the ecological robustness of networks to demonstrate that interaction 267 modifications have distinct dynamic effects. In ecological systems, efforts to include non-trophic 268 effects into food webs (Fontaine et al., 2011;Pilosof et al., 2017) should take into account the 269 higher-order nature of interaction modifications. Without quantitative information of the 270 topological and strength distribution of TIMs it is not yet possible to precisely calculate their impacts 271 relative to trophic interactions. Nevertheless, it is clear that they have the potential to change our 272 expectation of both the relative vulnerability of species to extinction and the impact the loss of 273 certain species will have (Donohue et al., 2017). 274 We found that despite an increase in the number of 'moving parts' the full TIM model was more 275 robust than the exclusively pairwise model. Generally, unstructured dynamic models of community 276 stability suggest that a higher connectivity is destabilising (May, 1972;Uchida & Drossel, 2007) 277 making our result somewhat surprising, but there are a number of explanations. The overall level of 278 dynamic connectance (fraction of realised pairwise links) between our two representations of TIMs 279 is the same. In part, our result can be attributed to the tendency for an overall decline in population 280 densities under the perturbation. When the strength (and disruptive influence) of the resultant NTE 281 is dependent on the product of two species, in general the overall strength of NTEs will decline 282 faster than when the strength is dependent on just one species. Furthermore, if the density of the 283 other interactor falls to zero, the NTE under the full TIM model will also fall to zero, but with the 284 pairwise model the NTE is maintained (even though the mechanism is no longer extant). Hence, as 285 species become extinct the effective dynamic connectance of the full TIM model declines faster than 286 in the pairwise case. 287 There are also differences in the prevalence and rapidity of more specific feedback loops, although 288 despite our relatively simple model this factor is challenging to trace in these complex systems 289 (Dambacher & Ramos-Jiliberto, 2007). Feedbacks modulating the strength of the NTEs are heavily 290 dependent on the relative speed and strength of multiple feedback loops, which are derived from 291 individual species properties. Nonetheless, as the impact of the more complex full TIM model 292 derived NTEs are dependent on more species, in general, it appears that mitigating feedback loops 293 are more immediate in the full TIM case. The effects caused by an interaction modification (however 294 modelled) will cause an immediate increase in one affected species and decrease in the other. In the 295 full TIM model, where the strength of the interaction modification is dependent on × , the 296 change to the NTE will be mitigated since and move in opposite directions. In the pairwise case 297 the strength of the NTE on each species will rapidly diverge. That said, algebraic solutions of the 298 consumer-resource module show that in the longer term an increase in an attack rate can lead to a 299 reduction in both of the interacting species equilibrium values (Morin, 2011), highlighting the 300 multiple timescales over which these processes operate. 301 At least in this simple model, it appears that TIMs determine extinction risk in a somewhat 302 predictable way. The extinction risk is shifted from those species that are closely trophically 303 connected to perturbed species, to be more evenly distributed throughout the network. Overall, the 304 species that went extinct tended to be those that were affected by more TIMs. When considering 305 NTEs from all species, this was true of both 'beneficial' and 'detrimental' NTEs. Those species that 306 directly benefit from the presence of other species are clearly affected by extinction cascades. 307 However, those species receiving an above-expectation number of detrimental NTEs were also at 308 higher risk of extinction. It would be expected that these species would benefit from the reduction 309 or removal of their inhibitors. Note that the correlation between the increased species level 310 connectance and the number of modifications received was taken into account by comparing 311 expected number of TIMs on a per-species basis. 312 There are two complementary explanations for this impact of detrimental NTEs. Firstly, in a complex 313 model certain species do increase their density, at least for a period of time, which would have 314 negative effect on the focal species. Initial responses of the populations to a small reduction in the 315 population of a single species were as commonly positive as negative (51.1% of non-neutral growth-316 rate responses to a 1% density reduction in a random species were positive in the TIM models, SI 5). 317 The low number of species maintaining large increases at the end of the process (what we defined 318 as population explosions) shows that such increases were relatively transient. When considering 319 only TIMs from the targeted species, which always decrease in density, extinct species received 320 fewer than expected detrimental NTEs, supporting this direct mechanism. Secondly, in the simple 321 model of population dynamics we used, an increase in the attack rate of a consumer (which would 322 be classed as 'beneficial') eventually leads to a reduction in the equilibrium density of the consumer 323 since the population of resource is pushed down to a level that cannot sustain as high a consumer 324 population. 325 Much previous work has shown that complex networks are stabilised by consistent patterns in key 326 parameters, which can be derived from body-mass scaling rules (Brose, Williams, & Martinez, 2006;327 Otto, Rall, & Brose, 2007). In our model this source of stability declines, as these patterns are 328 disrupted by interaction modifications that effectively push each value out of the allometrically 329 specified range in both directions. However, attack rates in real systems do take place in the context 330 of other species and significant disruptions to pairwise interaction strengths derived from laboratory 331 experiments caused by interaction modifications have been empirically observed (Jonsson,332 Kaartinen, Jonsson, & Bommarco, 2018). Our results suggest that although species that rely on 333 others are indeed more sensitive, the overall level of dynamic linkage is as important and that both 334 potential 'beneficial' and 'detrimental' processes should be considered. 335 Our choice of function to represent interaction modifications, although mathematically complex in 336 form, offers certain advantages compared to previous linear (Arditi et al., 2005;Bairey et al., 2016), 337 exponential (Goudard & Loreau, 2008;Lin & Sutherland, 2014) or rational function (Sanders et al.,338 2014) models. In particular, it has ability to directly and independently control salient features of the 339 function with clear ecological relevance (distance to threshold, maximum rate of change, range 340 between maximum and minimum). The dependence on the relative divergence from an equilibrial 341 starting point rather than the absolute value of the density of the modifier is also likely to be more 342 representative of ecological processes and considerably more straightforward to parameterise in an 343 ecologically meaningful manner. 344 Nevertheless, this is still a highly simplistic model. The interaction modifications included in this 345 study were introduced at random, in the sense that each potential modification had an equal chance 346 of existing. Considerable stabilising structuring has been observed within non-trophic effects 347 networks (Kéfi et al., 2015;Kéfi, Miele, et al., 2016) but there is not yet sufficient data to be able to 348 determine whether there are consistent features across ecosystems. Analyses of the stability of 349 systems at equilibrium demonstrate the key role of the distribution of interaction modifications 350 (Terry, Bonsall, & Morris, 2018) and the potential for emergent structure from the combination of 351 trophic and non-trophic interaction networks. The distribution of interaction modifications will have 352 important consequences for structural network properties such as trophic coherence, modularity 353 and nestedness that have been demonstrated to affect stability (Johnson, Domínguez-García, 354 Donetti, & Muñoz, 2014;Stouffer & Bascompte, 2011;Thébault & Fontaine, 2010). 355 Interaction modifications do not solely affect consumption rates and there is considerable scope to 356 develop and test theoretical models of additional processes, such as non-trophic interaction 357 modifications and modification of other aspects of consumer functional responses (Kéfi et al., 2012). 358 The assumption of multiplicative combination of multiple TIMs in a synergistic manner in our second 359 analysis is also questionable. While empirical patterns have not been quantified, it is a reasonable 360 estimate that many modification effects act antagonistically to each other (Golubski & Abrams, 361 2011). For instance, the presence of a second fear-inducing predator may well have less effect than 362 the first, dampening the impact of a change in either modifier population. 363 Further opportunities are also open to introduce specific accounting for the timescale of changes 364 and the size of perturbations. It is possible that higher-order interactions are stabilising against small 365 perturbations, which may make the system as a whole more susceptible to large impacts, such as 366 the extinction of certain species (Levine et al., 2017). TIMs also have the potential to create the 367 necessary positive feedback structures to maintain alterative stable states (Holt & Barfield, 2013;368 Kéfi, Holmgren, & Scheffer, 2016). As yet, the prevalence of these features is largely unknown. The 369 speed of interaction modifications themselves can also vary. While many TIMs are behaviourally 370 mediated and occur essentially instantaneously, others are due to accumulated environmental 371 changes (Sanders et al., 2014) or evolutionary processes (Benkman, Siepielski, & Smith, 2012) and 372 operate at somewhat slower timescales. 373

Conclusion 374
In summary, interaction modifications are potent forces that introduce distinct dynamics to 375 ecological networks. This distinctive nature of interaction modifications is of relevance for dynamic 376 systems in many fields that make use of networks (Strogatz, 2001). The complexity of networks is 377 more than the product of connectance and the number of interacting units and our work highlights 378 the importance of transient dynamics in the collapse of communities (Hastings, 2001). 379 Despite long-standing calls for the inclusion of non-trophic effects into the mainstream of ecological 380 network science that has been long dominated by food webs (Ings et al., 2009), and the publication 381 of the first empirical community level non-trophic network (Kéfi et al., 2015), there remains a great 382 number of significant unknowns about the role of non-trophic effects at the network scale. Our work 383 shows that maintaining interaction modifications as distinct processes within empirical and 384 theoretical networks, rather than as pairwise non-trophic effects (Grilli et al., 2017), will enable us to 385 have a more complete understanding of the system dynamics and allow better predictions of 386 community responses to perturbations. This need not require significant additional effort on the 387 part of the original investigator, but would be challenging for others to retroactively discern from 388 published pairwise interaction networks. Analyses of network robustness are used extensively to 389 understand anthropogenic impacts on natural communities (Evans, Pocock, & Memmott, 2013;390 Kaiser-Bunbury et al., 2017;Säterberg et al., 2013); as the development and analysis of 391 comprehensive interaction networks expands (Kéfi, Miele, et al., 2016)

Supplementary Information
Code and simulation data are available on the Open Science Framework repository. https://osf.io/w83br/ DOI: 10.17605/OSF.IO/W83BR

Method for determining body masses from trophic level
Body masses are assigned based on trophic level in the original, pre-integration to equilibrium, communities. First, trophic levels calculated on using the pathway method of Levine (1980), rounded into integer trophic levels. These trophic levels were then used to assign relative body masses specified on a logarithmic scale (base 10). We used distributions with a lower variance than empirical distributions to aid the stability of our systems by distributing body-masses in a relatively hierarchical manner with trophic level. i.
Producers were assigned a body mass of 0.
ii. Primary consumer (Trophic level 2)) body-masses were drawn from a normal distribution: (0.65, = 0.5), truncated at 0 to remove the possibility that consumers are smaller than herbivores. iii.
Higher-level consumer (Tropic level 3+) body masses were specified based on the sum of repeated draws from a normal distribution for each trophic level above 2 (h = TL-2), standardising the variance

Information about starting populations
The baseline communities are available in R list format as part of the online supplementary material. Summary statistics are shown in the Figure S2. Eighteen species was chosen as a lower cut off-point. B: Distribution of 3983 species across trophic levels, as calculated using the unweighted-pathway method. 8 species that had a calculated trophic level above 4 are not shown. C: Relationship between trophic level of a species and its initial equilibrium biomass, showing a general negative relationship. 67 populations has an initial biomass below 0.01 and are not shown for clarity. D: Initial Equilibrium biomass as a fraction of each species carrying capacity. Consumer equilibrium biomass was considerable below carrying capacity in all cases, while producer biomass varied from near capacity to well below.

Derivation of Gompertz Control Parameters
The basic form of the Gompertz equation is: In order to maintain the sign of , we define it based on the logarthim of the Gompertz function which can be negative log 10 ( ) = ( , ) The core results that we wish to achieve are: i. The core Gompertz equation ( − − ) ranges from 0-1, therefore the criterion iv) can be achieved by setting a= σ.
Criterion ii) can be simply met as max ( ) = , which can be rearranged to give = .
Criterion iii) can met from solving ''( ) = 0 for x, which is = ( )/ . This can then be rearranged to give = , where c is given by the previous result.
Criterion i) can be simply met, since when = 0, − = 1, we can find = − , where b and a are set from previous results.
Finally to address v) and fix the intrinsic bias of the core equation towards either positive or negative modifications depending on the sign of the modification, we randomly rotate half of the functions around the point (0,0) by setting = { −1 +1 with equal probability.

Criteria for establishing Equilibriums
Finding of Initial Equilibria.
Each model generated as described in the main text was initiated with a biomass density of 10. All models are available as R functions in the online code appendix. Models were integrated using the deSolve library within R , using the default 'lsoda' solver which switches automatically between stiff and non-stiff methods as necessary for the problem. Integration was continued with regular checks whether finishing conditions were met. If the biomass density of a species fell below a density of 10^-10 it was removed (its density was set to zero). After each extinction those consumers left without any remaining resources were also set as extinct. If extinctions resulted in the splitting of the community into multiple sub-webs the largest was kept and all the species in the smaller web were set as being extinct.
Models were run for a minimum of 400 time steps, continued until the community was considered to have reached an equilibrium. This was determined by the mean absolute ratio between the mean population of each species over the previous 200 time unit period and the next was less than 1%. The density of surviving at the last time point was then used as starting values with a multiple rootfinding algorithm from the R package rootSolve (Soetaert 2009) using the Newton-Raphson method, with the requirement for positive roots. This approach allowed a faster and accurate identification of true equilibria for slow to converge populations. In a handful of cases this root finding process determined the population biomass equilibrium to be below 10^-10. In these cases the model integration was declared a failure.

Robustness Analysis
The community models including external mortality and potentially TIMs were generated as described in the main text, with starting values as determined above. Models were integrated as above, in blocks of 200 time units. If the mean absolute ratio between the mean population of each species over the previous 200 time unit period and the next was less than 1%, the community was considered to have sufficiently stabilised. Consumers without any resources were immediately removed but trophically disconnected subunits were retained. If no acceptable equilibrium point had been reached after 100'000 time units or the integrator tolerances were exceeded the run was declared a failure. For statistical tests only complete pairwise comparisons were used. Failures were dominated by integrator tolerance issues, and were more common with the TIM model:

Response of populations to drop in population being as likely positive as negative.
To test the initial direction of response to perturbation we examined the direction of movement of populations when the biomasses were set the starting value with the exception of the 'targeted species' which was set to 99% of its initial value. This was done across all 200 tested communities, for each species. Of 75960 responses, 61.8% were neutral, 18.6 per negative and 19.5% were positive.

Partition of Modified Trophic Interaction into Trophic and Non-Trophic Components
For case where i is the focus, j's are the predators of i, l's are the prey of i and k's are the modifier, the original impact term is: First expand as: Then partition into a direct trophic effect and a non-trophic effect: Make the assumption that the in the non-trophic term is fixed at * : This approach is only valid where there is only a single TIM per interaction. Multiple TIMs would require considerable further complex terms to account for synergistic effects of TIMs. A similar process, fixing can then be followed to find the impact on the other species.