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

Abstract Considerable emphasis has been placed recently on the importance of incorporating non‐trophic effects into our understanding of ecological networks. Interaction modifications are well‐established as generating strong non‐trophic impacts by modulating the strength of interspecific 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 maximize information about the system dynamics. Interaction modifications have the potential to shift expectations of species vulnerability based exclusively on trophic networks.


| INTRODUC TI ON
There is a building appreciation that to improve our understanding of population dynamics within ecological communities, it is necessary to move beyond studies that focus on a single interaction process at a time (Kéfi et al., 2012;Levine, Bascompte, Adler, & Allesina, 2017).
Trophic interaction modifications (TIMs) (Terry, Morris, & Bonsall, 2017;Wootton, 1993) occur when a consumer-resource interaction is modulated by additional species. These are a class of higherorder processes since their effects are not fundamentally pairwise.
Approaches to understanding interaction modifications often try to distil the inherently multi-species process into a pairwise NTE (or "trait-mediated indirect effect") from the modifier species onto one or both recipient species (Okuyama & Bolker, 2007, Figure 1). This allows the direct comparison of non-trophic and trophic interaction strengths (Preisser et al., 2005) and network structure (Pilosof, Porter, Pascual, & Kéfi, 2017) but is a representation of a different class of dynamic process (Terry et al., 2017). This simplification can give valuable insights into communities at equilibrium (e.g. Grilli, Barabás, Michalska-Smith, & Allesina, 2017). However, the consequences of this assumption in a transient, fluctuating or heavily perturbed system have yet to be fully explored. Previous studies introducing interaction modifications to trophic networks (Arditi, Michalski, & Hirzel, 2005;Goudard & Loreau, 2008;Lin & Sutherland, 2014) have demonstrated their potential impact on the dynamics of ecosystems, but not whether this is attributable to the higher-order nature of interaction modifications, as opposed to shifts in connectance and interaction strength.
An important case is the dynamics of ecological systems in the face of species removal, where there is the potential for secondary extinctions and eventually the collapse of the ecosystem (Dunne, Williams, & Martinez, 2002). This aspect of stability, often described as "robustness," is important both from the perspective of managing anthropogenic change and in terms of understanding the fundamental stability of ecological communities. Since empirically testing how whole communities respond to extinctions can be difficult or impossible (although see Sanders, Thébault, Kehoe, and Frank van Veen (2018)), a number of studies have attempted to determine the properties that make ecological communities robust through simulation (Dunne & Williams, 2009;Säterberg, Sellman, & Ebenman, 2013). However, incorporating the acknowledged flexibility of ecological networks is a perennial challenge for such studies (Montoya, Pimm, & Solé, 2006).
Calls to incorporate the full panoply of NTEs into our understanding of ecological networks have built substantially in recent years (Fontaine et al., 2011;Ings et al., 2009;Kéfi et al., 2012;Levine et al., 2017;Ohgushi et al., 2012;Olff et al., 2009), and the first empirical inventories are being established (Kéfi et al., 2015;Kéfi, Miele, Wieters, Navarrete, & Berlow, 2016). Theoretical analyses can play a significant F I G U R E 1 Depiction of distinction between trophic interaction modifications (TIM) and consequent non-trophic effects (NTEs) when describing the impact of a modifier species on a consumerresource pair. In (a) the modifier species causes a facilitating TIM (blue dashed line) where an increase in the modifier leads to a beneficial NTE on the consumer and a detrimental NTE on the resource. Correspondingly, in (b) an interfering TIM causes a beneficial NTE on the resource and a detrimental NTE on the consumer role in motivating the empirical construction of networks, identifying the information necessary to best understand these systems. Here, we demonstrate the distinctive nature of interaction modifications compared to pairwise NTEs through a direct standardized comparison of their impacts on the robustness of large artificial networks.
We then examine the role of interaction modifications in determining species-level vulnerability to secondary extinction.

| MATERIAL S AND ME THODS
We conducted our analyses using model food webs where system dynamics are derived from metabolic scaling relationships. As detailed below, interaction modifications were introduced to a set of communities each at an initial equilibrium. Robustness to extinction was examined by introducing external mortality to a single species at a time and integrating the model to a new equilibrium. All analyses were carried out in r v.3.5.0 (R Core Team, 2018) using the deSolve numerical integration package (Soetaert, Petzoldt, & Setzer, 2010), and all code and data are available online.

| Bio-energetic model
The change in biomass density, B i , of each species in the community was modelled using a simple Lotka-Volterra type model with a linear (Holling type I) functional response and logistic intrinsic growth rates, parameterized using body-mass relationships: Each species was assigned a body mass (M i ) drawn from a distribution based on their trophic level (see Appendix S1: Section S1) which was then used to calculate further parameters using quarter-power body-mass scaling laws (Yodzis & Innes, 1992). Relative intrinsic growth or metabolic loss rates, r i , were set at 1 for all producers (trophic level 1) and r i = −0.1M −0.25 i for each consumer (trophic level ≥ 2).
Consumer-specific consumption rates were set at a ij = j M −0.25 j , where the generality term j was fixed at 1∕n, the number of resources of each consumer j. Assimilation efficiencies, e ij , for each trophic interaction were drawn from a uniform distribution centred around 0.1 (Moore & de Ruiter, 2012), e ij ∼  0.05,0.15 . Carrying capacities, K i , were drawn from K i ∼  1,10 for producers (to introduce a moderate degree of self-regulation) and K i ∼ 10  (2,3) for consumers (considerably higher than the starting populations and so introducing only a small amount of self-regulation). Initially, external mortality, m i , was set to 0 and modification terms, ij , set to 1.
The trophic topology and population densities of a set of 200 starting communities for the robustness tests were generated as follows. Initial trophic topologies were generated using the niche model (Williams & Martinez, 2000), with 35 species and a connectance of 0.14, removing cannibalistic interactions. Each population density was initially set to 10 and the system numerically integrated to a stable equilibrium (see Appendix S1: Section S4 for criteria). Only fully connected communities with at least 18 persisting species were retained. Properties of the starting communities are described in Appendix S1: Section S2.

| Specification of TIMs
Trophic interaction modifications were introduced through modification terms, ijk , that specify the impact of modifier species k on the consumption of species i by species j, where i , j and k are all different. These are positive numbers that multiply the attack rate as a function (detailed in the next section) of the divergence of the biomass density of a modifying species, B k , from its starting equilibrium value B * k . Consequently, when ijk is smaller than 1, the interaction is weakened and when ijk is greater than 1, the interaction is strengthened. Where multiple species modify the same interaction, these effects were assumed to be synergistic and combine multiplicatively, ij = ∏ k ijk (Golubski & Abrams, 2011;Goudard & Loreau, 2008). In our model, ijk cannot be negative, preventing the reversal of the direction of the trophic interaction.
Interaction modifications cause both positive and negative NTEs ( Figure 1). Where an increase in the modifier species leads to an increase in the strength of the interaction, we describe the modification as a facilitating TIM. This can be said to be beneficial to the consumer and detrimental to the resource, in the sense of the immediate impact from an increase in the modifier population. We term the reverse situation an interfering TIM.

| Functional form of TIM
To represent the relationship between the density of the modifier species and the modification of the interaction, we used a Gompertz sigmoidal curve parameterized to control features of ecological relevance, detailed in full in Appendix S1: Section S3. This function links the magnitude of divergence of B k from its start point B * k , as log 10 B k ∕B * k , and three control parameters ( Figure 2). These are (the maximum rate of change in the modified interaction as the modifier density increases), (the proportional change from B * k to reach the threshold point of maximum response) and (the range of magnitudes over which ijk spans). A positive α denotes a facilitatory modification and a negative α an interfering modification. A key attribute of our parameterization is that when B k = B * k , ijk = 1. This maintains the original trophic interaction strength when the modifier is at its starting equilibrium point-in effect, these strengths are assumed to already incorporate the effect of the modifier to the interaction at the equilibrium.

| Alternative representation of TIMs as pairwise non-trophic effects
As an alternative to modelling the impact of TIMs in full (using "higher-order" terms), equivalent pairwise effects can be derived.
(1) F I G U R E 3 Illustrative response surfaces showing the distinction between non-trophic effects caused by pairwise interactions (a) and those from a higher-order interaction modification process involving three species (b). Per capita effects are shown in the vertical axis on a logarithmic scale. In both cases, the response to the modifier is identical at the starting density of the trophic interaction partner of the focal species, shown by the black line. However, the density of the trophic interaction partner of the focal species (left-hand axis) does not affect the non-trophic effect in the pairwise case (a), but does with the full interaction modification model (b). Note also that as the interactor's density becomes small, the strength of the non-trophic effect rapidly declines Section S6. A trophic interaction affecting a resource i influenced by a modifier k can be partitioned from: to: The corresponding terms for the interactions affecting the consumer, j, are as follows: As an example, a facilitating TIM, when B k increased, would lead to a ijk greater than 1. This would lead to a negative NTE on the resource i (independent of B j ) and a positive NTE on the consumer j (independent of B i ). This partitioning process (and hence direct comparison between full and pairwise models) is only straightforward when each trophic interaction is modified by at most one modifier species because of the synergistic relationship between multiple modifiers assumed by our model.

| Test 1. Comparison of robustness between pairwise and higher-order models of TIMs
To compare the consequences of introducing TIMs by these two approaches, we conducted three sets of robustness tests using the same set of trophic networks. The first used the full TIM model, the second used TIMs that had been converted to pairwise form and a third case without any TIMs.
We randomly added TIMs to the set of initial communities such For each robustness test, the external mortality rate (m i ) of a single species (which we will refer to as the "targeted" species for convenience, although the mortality could be attributable to a range of non-directed processes) was then set to 1, the system numerically integrated to a new steady state and the status of each species in the resultant community assessed (Appendix S1: Section S4). Species were considered "extinct" if their biomass fell below 10 −4 (several orders of magnitude below the starting values of most of the populations, Appendix S1: Section S2), "functionally extinct" if their biomass fell to below 1/10th of their starting value and considered to have "exploded" if the final density was over 10 times the starting value. Robustness tests were repeated, targeting each species in turn for each community, to give a total of 3,736 completed tests (93.8% successful integrations, Appendix S1: Section S4).

| RE SULTS
The introduction of TIMs greatly increased the number of extinctions ( Figure 4). However, almost double the number of extinctions were observed when the interaction modifications were represented with the pairwise model compared to when they

| D ISCUSS I ON
Here, we have examined the ecological robustness of networks to demonstrate that interaction modifications have distinct dynamic effects. In ecological systems, efforts to include NTEs into food webs (Fontaine et al., 2011;Pilosof et al., 2017) should take into account the higher-order nature of interaction modifications. Without quantitative information of the topological and strength distribution of TIMs, it is not yet possible to precisely calculate their impacts relative to trophic interactions. Nevertheless, it is clear that they have the potential to change our expectation of both the relative vulnerability of species to extinction and the impact the loss of certain species will have (Donohue et al., 2017).
We found that when TIMs were fully represented our model systems were notably more robust than under an exclusively pairwise model, yet the overall number of "direct" relationships between species through all types of interaction is the same between our two representations of TIMs. While higher trophic connectance can increase robustness to secondary extinctions (Dunne & Williams, 2009), increased connectance due to TIMs has the opposite effect in our model. However, when isolated from connectance, the higherorder nature of trophic interactions increases robustness. This aligns with the results of Bairey et al. (2016) who found that higher-order interactions could increase the persistence and local stability of F I G U R E 4 Boxplot showing that robustness was significantly lower when (TIMs) were modelled as pairwise interactions rather than directly as full interaction modifications. Both cases produced significantly more extinctions than the 'No TIM' case (paired t tests, p < .0001, n = 3,942) and connectance in higher-order systems is a challenge (Golubski, Westlund, Vandermeer, & Pascual, 2016). Our results highlight the need to develop how the complexity of higher-order systems can be best discussed and analysed, since differences between these models would not be picked up under classic ecological complexity measures (e.g. May, 1972).
In part, our result can be attributed to the tendency for an overall decline in population densities after perturbation (Appendix S1: Section S7). In the full TIM model, the strength (and disruptive influence) of a modification is dependent on the product of two species and so (in the context of generally declining species densities) will decline faster than when the strength is dependent on just one species. Weaker NTEs will reduce the potential for cascading effects from species removal, analogous to the impact of trophic link strength (Zhao et al., 2016). Furthermore, if the density of either consumer or resource falls to zero, any TIM impact on the species it is trophically linked to would also fall to zero under the full TIM model. However, with the pairwise model, the NTE is maintained even though the mechanism is no longer extant. Hence, as species become extinct, the effective connectance of the full TIM model declines faster than in the pairwise case. These mechanisms would only be evident in non-equilibrial analyses.
There are also differences between the full TIM and the pair- At least in this simple model, it appears that TIMs disperse relative extinction risk from those species that are closely trophically connected to the perturbed species, to be more evenly distributed throughout the network. Trophic interaction modifications greatly increase network connectivity-for instance, our inclusion of randomly distributed TIMs in our second analysis brought the mean overall path length between species down from 3.0 to 1.6 and halved the mean network breadth from 4.2 to 2.1. While it is likely that in real ecological systems, many TIMs (for example, those caused by consumer-avoidance responses) link trophically close species, others such as those caused by ecosystem engineer species may well create the long links that cause this short-circuit effect.
Overall, the species that went extinct tended to be those that were affected by more TIMs. When considering NTEs from all species, this was true of both "beneficial" and "detrimental" NTEs. Those species that directly benefit from the presence of other species are clearly affected by extinction cascades. However, those species receiving an above-expectation number of detrimental NTEs were also at higher risk of extinction despite the prospect that these species would benefit from the on-average reduction or removal of their inhibitors. The correlation between the increased species-level trophic connectance and the number of modifications received was accounted for by comparing expected number of NTEs on a per-species basis. In part, the vulnerability of species subject to detrimental NTEs can be explained by temporary increases in modifier species.
Initial responses to a small reduction in the population of a single species were as commonly positive as negative (51.1% of non-zero growth-rate responses to a 1% density reduction in a random species were positive in the TIM models, Appendix S1: Section S5). The low number of species maintaining large increases at the end of the process (what we defined as population explosions) shows that such increases were relatively transient. When considering only TIMs from the targeted species, which always decrease in density, extinct species received fewer than expected detrimental NTEs, supporting this direct explanation. Our results suggest that, although species facilitated by other species are indeed more sensitive to extinction, it is the overall number of relationships with other species that is critical. Both apparently "beneficial" and "detrimental" processes should be considered on an equal footing.
Much previous work has shown that complex networks are stabilized by consistent patterns in key parameters, which can be derived from body-mass scaling rules (Brose, Williams, & Martinez, 2006;Otto, Rall, & Brose, 2007). In our model, this source of stability declines, as these patterns are disrupted by interaction modifications that effectively push each attack rate value out of the allometrically specified range in both directions. Tracking modification dynamics during simulations is complex, but it is clear that attack rates shifted considerably in our model. To take one illustration, the extinction of a modifier species in our second test would cause a median attack rate change of a factor of 5.62, but with a long tail of larger modifications. Attack rates in real systems take place in the context of other species, and significant disruptions to pairwise interaction strengths derived from laboratory experiments caused by interaction modifications have been empirically observed (Jonsson, Kaartinen, Jonsson, & Bommarco, 2018). The extent of attack rate and interaction strengths variability is an important qualifier to observed allometric patterns.
Our choice of a Gompertz function to represent interaction modifications, although mathematically complex in form, offers certain advantages compared to previous linear (Arditi et al., 2005;Bairey et al., 2016), exponential (Goudard & Loreau, 2008;Lin & Sutherland, 2014) or rational function (Sanders et al., 2014)  Nevertheless, this is still a highly simplistic model. The interaction modifications included in this study were introduced at random, in the sense that each potential modification had an equal chance of existing.
Considerable stabilizing structuring has been observed within a rocky shore NTE network (Kéfi et al., 2015;Kéfi, Miele, et al., 2016) (Kalinkat et al., 2013) and may act to limit the potential impact of increases in trophic interaction rates mediated by TIMs by placing an upper limit on the overall intake rates. There is considerable scope to compare the impact of modifications to different aspects of nonlinear functional responses beyond overall consumption rate (Kéfi et al., 2012). Secondly, further work is needed to explore the consequences of nonlinear combination of multiple TIMs, beyond the synergistic assumptions used here. It is a reasonable estimate that many modification effects act antagonistically to each other (Golubski & Abrams, 2011). For instance, the presence of a second fear-inducing predator may well have less effect than the first, dampening the impact of a change in either modifier population. However, the empirical base to parameterize any such analyses is currently very small.
Further opportunities for future work include introducing specific accounting for the time-scale of changes and the size of perturbations. It is possible that higher-order interactions are stabilizing against small perturbations, which may make the system as a whole more susceptible to large impacts, such as the extinction of certain species (Levine et al., 2017). Trophic interaction modifications also have the potential to create the necessary positive feedback structures to maintain alternative stable states (Holt & Barfield, 2013;Kéfi, Holmgren, & Scheffer, 2016). As yet, the prevalence of these features is largely unknown. The speed of interaction modifications themselves can also vary. While many TIMs are behaviourally mediated and occur essentially instantaneously, others are due to accumulated environmental changes (Sanders et al., 2014) or evolutionary processes (Benkman, Siepielski, & Smith, 2012) and operate at somewhat slower time-scales.

| CON CLUS IONS
In summary, interaction modifications are potent forces that introduce distinct dynamics to ecological networks. This distinctive nature of interaction modifications is of relevance for dynamic systems in many fields that make use of networks (Strogatz, 2001) since our work shows that the complexity of networks is more than the product of connectance and the number of interacting units. Despite long-standing calls for the inclusion of NTEs into the mainstream of ecological network science that has been long dominated by food webs (Ings et al., 2009), and the publication of the first empirical community level non-trophic network (Kéfi et al., 2015), there remains a great number of significant unknowns about the role of NTEs at the network scale. Our work shows that maintaining interaction modifications as distinct processes within empirical and theoretical networks, rather than as pairwise NTEs (Grilli et al., 2017), will enable a more complete understanding of the system dynamics and allow better predictions of community responses to perturbations.
When documenting non-trophic interactions, identifying processes as interaction modifications need not necessarily require significant additional effort on the part of the original investigator, but would be challenging for others to retroactively discern from published pairwise interaction networks. Analyses of network robustness are used extensively to understand anthropogenic impacts on natural communities (Evans, Pocock, & Memmott, 2013;Kaiser-Bunbury et al., 2017;Säterberg et al., 2013); as the development and analysis of comprehensive interaction networks expand (Kéfi, Miele, et al., 2016), we must incorporate interaction modifications appropriately.

DATA AVA I L A B I L I T Y S TAT E M E N T
All R code used in this study and the generated simulation results are available on the Open Science Framework: https ://doi. org/10.17605/ OSF.IO/W83BR (Terry, 2019) or github.com/jcdterry/ TIMs_and_Robustness.