# UNIFYING ECOLOGY ACROSS SCALES: PROGRESS, CHALLENGES AND OPPORTUNITIES

EDITED BY : Mary I. O'Connor, Diego Barneche, Julie Messier and Angelica L. Gonzalez PUBLISHED IN : Frontiers in Ecology and Evolution

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88966-292-0 DOI 10.3389/978-2-88966-292-0

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# UNIFYING ECOLOGY ACROSS SCALES: PROGRESS, CHALLENGES AND OPPORTUNITIES

Topic Editors: Mary I. O'Connor, University of British Columbia, Canada Diego Barneche, University of Exeter, United Kingdom Julie Messier, University of Waterloo, Canada Angelica L. Gonzalez, The State University of New Jersey, Camden, United States

Citation: O'Connor, M. I., Barneche, D., Messier, J., Gonzalez, A. L., eds. (2020). Unifying Ecology Across Scales: Progress, Challenges and Opportunities. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88966-292-0

# Table of Contents

*05 Editorial: Unifying Ecology Across Scales: Progress, Challenges and Opportunities* Mary I. O'Connor, Diego R. Barneche, Angélica L. González and Julie Messier *08 General Form for Interaction Measures and Framework for Deriving Higher-Order Emergent Effects* Elif Tekin, Pamela J. Yeh and Van M. Savage *20 The Impact of Spatial and Temporal Dimensions of Disturbances on Ecosystem Stability* Yuval R. Zelnik, Jean-François Arnoldi and Michel Loreau *30 The Sensitivity of Subsurface Microbes to Ocean Warming Accentuates Future Declines in Particulate Carbon Export* Emma Louise Cavan, Stephanie A. Henson and Philip W. Boyd *40 On the Interpretation of the Normalization Constant in the Scaling Equation* Karl J. Niklas and Sean T. Hammond *51 "Unifying" the Concept of Resource Use Efficiency in Ecology* Dorothee Hodapp, Helmut Hillebrand and Maren Striebel *65 Eco-Evolutionary Origins of Diverse Abundance, Biomass, and Trophic Structures in Food Webs* Jean P. Gibert and Justin D. Yeakel *76 Temperature Modifies Consumer-Resource Interaction Strength Through Its Effects on Biological Rates and Body Mass* Azenor Bideault, Michel Loreau and Dominique Gravel *87 Scaling Contagious Disturbance: A Spatially-Implicit Dynamic Model* Tempest D. McCabe and Michael C. Dietze *104 Integrating Multiple Data Types to Connect Ecological Theory and Data Among Levels* Jian D. L. Yen, Zeb Tonkin, Jarod Lyon, Wayne Koster, Adrian Kitchingman, Kasey Stamation and Peter A. Vesk *111 When Landscape Ecology Meets Physiology: Effects of Habitat Fragmentation on Resource Allocation Trade-Offs* Yaron Ziv and Goggy Davidowitz *119 The Metabolic Basis of Fat Tail Distributions in Populations and Community Fluctuations* Angel M. Segura and Gonzalo Perera *127 Effects of Climate Warming on Consumer-Resource Interactions: A Latitudinal Perspective* Priyanga Amarasekare *142 Interaction Dimensionality Scales Up to Generate Bimodal Consumer-Resource Size-Ratio Distributions in Ecological Communities*

Samraat Pawar, Anthony I. Dell, Tianyun Lin, Daniel J. Wieczynski and Van M. Savage

Neo D. Martinez

# Editorial: Unifying Ecology Across Scales: Progress, Challenges and Opportunities

Mary I. O'Connor 1 \*, Diego R. Barneche2,3, Angélica L. González 4 and Julie Messier 5

1 Department of Zoology and Biodiversity Research Centre, University of British Columbia, Vancouver, BC, Canada, 2 Australian Institute of Marine Science, Crawley, WA, Australia, 3 Oceans Institute, University of Western Australia, Crawley, WA, Australia, 4 Department of Biology, Center for Computational and Integrative Biology, Rutgers University, Camden, NJ, United States, 5 Biology Department, University of Waterloo, Waterloo, ON, Canada

Keywords: ecology, global change, models, allometry, information theory, metabolic theory of ecology (MTE), population dynamics

### **Editorial on the Research Topic**

**Unifying Ecology Across Scales: Progress, Challenges and Opportunities**

### THE NEED TO UNIFY ECOLOGY: A GRAND CHALLENGE FOR A CHANGING WORLD

Ecology is the science of how living systems grow, change, and persist. Although this is not the definition presented in most textbooks, this is the central theme of this scientific discipline as it is practiced in the current era of rapid global change. Change comes in many forms, from natural succession of communities to biological invasions affecting patterns of biodiversity to the collapse of fisheries. Understanding and forecasting natural change as well as mitigating undesirable anthropogenic change is one of the grand challenges we face in the twenty-first century.

 It has become clear that although we have vast and multifaceted ecological knowledge, it has yet to solidify into a coherent body of science. For example, over decades in the latter half of the twentieth century researchers, journals and even entire academic departments focused on specific branches of ecological thinking as narrowly defined as behavioral ecology, population ecology, or community ecology. Currently though, it has become clear that to understand our changing world and our place in it, examining the ecological changes afoot from a fragmented and narrow disciplinary perspective is insufficient. Population dynamics, ecosystem functions, individual behavior, and other aspects of living systems are deeply connected, and we cannot project changes in one without understanding how they are related to other processes across scales of space and time and levels of biological organization.

### CONTRIBUTIONS IN THIS SPECIAL TOPIC

In recent decades, efforts to unify ecological understanding across scales and levels of organization have flourished and disciplinary barriers have fallen (see, for example, Brown et al., 2004; Scheiner and Willig, 2008; Harte, 2011; Marquet et al., 2014). In this special topic, we highlight recent progress in the unification of ecological sciences. Papers in this issue can be organized in three groups. Some papers take a very broad view on the current state of ecology, proposing frameworks and approaches to organize understanding and guide new research (Fulton et al.; Hodapp et al.;

#### Edited and reviewed by:

Peter Convey, British Antarctic Survey (BAS), United Kingdom

> \*Correspondence: Mary I. O'Connor<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 25 September 2020 Accepted: 09 October 2020 Published: 30 October 2020

#### Citation:

O'Connor MI, Barneche DR, González AL and Messier J (2020) Editorial: Unifying Ecology Across Scales: Progress, Challenges and Opportunities. Front. Ecol. Evol. 8:610459. doi: 10.3389/fevo.2020.610459 Kempes et al.; Martinez; O'Connor et al.). A second group provides concrete tools and approaches to unifying ecological data and models, through specific and quantitative frameworks (McCabe and Dietze; Newman et al.; Niklas and Hammond; Tekin et al.; Yen et al.). A third group provides specific analyses that reveal new insights into how dynamics in natural systems are linked across scales and levels of organization (Amarasekare; Bideault et al.; Cavan et al.; Gibert and Yeakel; Pawar et al.; Segura and Perera; Zelnik et al.; Ziv and Davidowitz). Together, this issue is full of big ideas complemented by specific analyses, providing an exciting and welcome step in our ongoing effort for a unified science of ecology.

Fulton et al. offer a sweeping view of ecology, moving through what they identify as a transitional period now to a future in which we will have developed a more coherent ecological science through better models. This paper is a mustread, providing an impressive and inspiring perspective on the successes of synthesis in ecology through models. Among the many successes they highlight is the recognition that we understand that human activities and evolutionary change both occur within ecological time scales, and "this new reality drives home that ecology lacks a universally solid foundation from which to extend existing theories and modeling approaches to easily absorb the many interacting components and scales." Fulton et al. identify immediate challenges as the need to embrace new modeling approaches that span more dimensions and allow more dynamic and evolving parameterizations and structures to allow models to incorporate new sources of data. They emphasize that these advances require not only technical advances in modeling, but also coordinated advances in ecological theory and concepts, as well as multiscale data streams. To continue the ecological tradition of using models to explain and predict, now in a rapidly changing world, they argue for continued investment in an approach that unifies disparate ideas while also allowing multiple approaches to be employed to solve the problems we face.

Other contributions propose frameworks to unify concepts across scales, disciplines, and levels of organization. Consistent with Fulton et al.'s call for greater attention to a universally solid foundation for ecological science, O'Connor et al. propose strengthening the foundation by integrating concepts of information and information theory into ecology's fundamental principles. They review major concepts in information theory, and highlight how these fundamental principles of science have not been fully explored in ecological understanding. Their paper provides a set of ecological principles grounded firmly in a science that includes information processing as fundamental to living systems, along with energy flows and material cycles. Similarly, with firm grounding in these first principles, Kempes et al. provide scaffolding that explains how evolution must occur within the boundaries of the physical world, and what these physical constraints imply for all that we still do not know about evolution. Martinez extends a unifying framework further with allometric networks, and like several other proposed frameworks in this issue, Martinez includes humans and socioecological systems in his unifying scheme. Hodapp et al. provide additional structure through emphasizing resource use efficiency as a unifying concept across diverse living systems, which may differ in the identity and number of limiting resources, trophic levels, and ecological, temporal and spatial scales.

Inspired by the broad roadmap laid out by contributions like those described above, we might then ask, "how do we proceed and implement a unifying framework?". This is a hard question. Thoughtful and constructive answers are provided in every paper in this edition. Yen et al. identify the statistical challenge of integrating data—often collected specifically at one level of biological organization—across multiple levels of organization. They provide an approach integrating multiple data types and parameterize more complex, process-explicit models, providing a path forward to testing multi-scale theory using data spanning many organizational levels. Niklas and Hammond's contribution dives into understanding body form and function in generalized terms, and renders the allometric scaling model more biologically relevant and useful by exposing the value and meaning of the often-overlooked scaling normalization constant. Newman et al. identify and explain specific challenges that we face when trying to understand emergent patterns of complex ecological systems at landscape scales: coarse-graining, the middle number problem of the domain of data complexity, and non-stationarity issues when predicting future ecological states from models based on adjustable parameters. McCabe and Dietze offer a solution to another landscape-scale problem—scaling contagious disturbance. They provide an approach that allows modeling of such events in the context of climate change, thereby linking community-scale disturbance events with broader scale climate change and feedbacks. Cavan et al. also work with climate change models and finer scale ecological dynamics by using the metabolic theory of ecology (MTE) to link rates of carbon flux in microbial systems to global climate feedbacks.

Tekin et al.'s contribution provides an approach that could potentially also serve as an important part of a unifying framework by demonstrating that many if not all measures of pairwise interactions—present in a diverse array of complex systems across very different scientific disciplines—can be derived from and analyzed in the context of a single mathematical framework. They build upon this finding to introduce consistent measures of higher order interactions. In this way they provide a practical tool and an insightful guide to thinking about complexity that is likely to benefit unification efforts across ecological (and other scientific) disciplines.

The third set of papers in this issue dives deeper into ecological and evolutionary dynamical models. These contributions implement unifying concepts by integrating ideas from diverse areas of thought into a single framework for population dynamics, and then analyzing the models to understand possible patterns and processes in nature. Gibert and Yeakel demonstrate how evolutionary change within the context of ecological dynamics can provide unexpected links between individual traits and the structure of an emergent food web. Ziv and Davidowitz link landscape-scale features (patch size and fragmentation) to evolutionary outcomes (selection for life-history traits increasing survival) and population dynamics (likelihood of population extinction). The need for strong quantitative general predictions of fluctuations of species abundance are addressed by Segura and Perera through deriving a theoretical model that links metabolism to abundance fluctuations in ecological populations and communities. Amarasekare carefully considers how physiological effects of temperature on the rates that drive consumer resource dynamics can lead to very different consequences of warming and climate change in different thermal regimes (e.g., tropics vs. temperature). Tackling the same question from a different angle, Bideault et al. find that the way temperature dependences are distributed among the processes that determine consumer-resource interactions and outcomes can substantially affect the ultimate outcome of warming. Pawar et al. build on an immense database of consumer-resource interactions and related traits to observe that body size distributions differ for consumers foraging in two-dimensional and three-dimensional habitats. They provide a theoretical explanation for this pattern, and use their joint mathematical/empirical approach to offer new insights about how consumer resource interactions are constrained and how they flourish in different environments. Zelnik et al. provide a theoretical analysis of the effects of different spatial and temporal properties of disturbance on ecosystem stability. Their modeling approach provides insights into the tight link between temporal and spatial dimensions on ecosystem responses to disturbance. By integrating Earth System Models and satellitederived temperature data, and using principles from MTE, Cavan et al. estimate changes in global exports of marine carbon. Their findings suggest a main role of temperature sensitivity shaping carbon outputs and call for the need of incorporating organismal temperature dependence into biogeochemical models. Each of these contributions is exciting, novel and potentially transformative to how we think about ecological dynamics and the structures they produce. These are the bricks and mortar of a unified science, providing specific and clear

### REFERENCES

new ideas and observations to our body of knowledge. Many work at the intersection between population dynamics and the unifying approach of MTE, which is based on general principles of how ecological rates of energy flux and material cycling scale with body size and temperature (Brown et al., 2004).

The contributions in this issue advance current themes in unifying efforts, specifically those employing MTE (Bideault et al.; Cavan et al.) or its core allometric scaling model (Kempes et al.; Martinez; Niklas and Hammond), those integrating data and models (Fulton et al.; Newman et al.; Yen et al.), bridging between local scale dynamics and climate change models (Cavan et al.; Fulton et al.; McCabe and Dietze), and those focusing on resource use and trade-offs (Hodapp et al.; Ziv and Davidowicz). Newer themes in unifying ecological science are also emerging, in the form of emphasis on disturbance and its multi-scale consequences (McCabe and Dietze; Segura and Perera; Zelnik et al.), on trait based perspectives on consumer-resource interactions and interaction strengths (Amarasekare; Bideault et al.; Pawar et al.; Tekin et al.), eco-evolutionary feedbacks (Gilbert and Yeakel; Kempes et al.) and on information as a unifying principle for mainstream ecological thinking (O'Connor et al.). Together, the full issue opens doors and guides our steps through them in our collective efforts to better understand our changing biosphere through ecological science.

## AUTHOR CONTRIBUTIONS

All authors created the vision for the special issue, served as editors, and discussed the main concept for the editorial manuscript. MO'C wrote the manuscript. All authors reviewed and edited it.

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 O'Connor, Barneche, González and Messier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# General Form for Interaction Measures and Framework for Deriving Higher-Order Emergent Effects

#### Elif Tekin1,2, Pamela J. Yeh1,3 and Van M. Savage1,2,3 \*

1 Department of Ecology and Evolutionary Biology, University of California, Los Angeles, Los Angeles, CA, United States, 2 Department of Biomathematics, David Geffen School of Medicine, University of California, Los Angeles, Los Angeles, CA, United States, 3 Santa Fe Institute, Santa Fe, NM, United States

Interactions are ubiquitous and have been extensively studied in many ecological, evolutionary, and physiological systems. A variety of measures—ANOVA, covariance, epistatic additivity, mutual information, joint cumulants, Bliss independence—exist that compute interactions across fields. However, these are not discussed and derived within a single, general framework. This missing framework likely contributes to the confusion about proper formulations and interpretations of higher-order interactions. Intriguingly, despite higher-order interactions having received little attention, they have been recently discovered to be highly prevalent and to likely impact the dynamics of complex biological systems. Here, we introduce a single, explicit mathematical framework that simultaneously encompasses all of these measures of pairwise interactions. The generality and simplicity of this framework allows us to establish a rigorous method for deriving higher-order interaction measures based on any of the pairwise interactions listed above. These generalized higher-order interaction measures enable the exploration of emergent phenomena across systems such as multiple predator effects, gene epistasis, and environmental stressors. These results provide a mechanistic basis to better account for how interactions affect biological systems. Our theoretical advance provides a foundation for understanding multi-component interactions in complex systems such as evolving populations within ecosystems or communities.

## Keywords: complex biological systems, emergent patterns, higher-order interactions, ecological interactions,

### INTRODUCTION

biodiversity

Because of their key role in understanding the dynamics of complex biological, physical, and social systems, there is a long and rich history of studying interactions and their consequences (Wootton, 1993; Billick and Case, 1994; Darling and Côté, 2008; Mihaila et al., 2010; Hamilton, 2011; Toprak et al., 2013; Barrios-O'neill et al., 2014; Foucquier and Guedj, 2015; Palmer et al., 2015; Podgornaia and Laub, 2015; Nishikawa and Motter, 2016; Shi, 2016). These approaches have often been complemented and enhanced by network theory that has led to important advances in prediction of patterns (Segrè et al., 2005; Yeh et al., 2006; Braun and Shah, 2015). Studies of interactions and networks have heavily and almost solely focused on two-component interactions.

#### Edited by:

Mary I. O'Connor, University of British Columbia, Canada

#### Reviewed by:

Marcelo Michel Rivadeneira, Centro de Estudios Avanzados en Zonas Áridas (CEAZA), Chile Juan Ernesto Guevara Andino, Field Museum of Natural History, United States

> \*Correspondence: Van M. Savage<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 11 May 2018 Accepted: 02 October 2018 Published: 23 October 2018

#### Citation:

Tekin E, Yeh PJ and Savage VM (2018) General Form for Interaction Measures and Framework for Deriving Higher-Order Emergent Effects. Front. Ecol. Evol. 6:166. doi: 10.3389/fevo.2018.00166 As a result, we have a thorough understanding of pairwise interactions and the various interaction categorization measures defined across distinct subject areas such as statistical physics, gene networks, and prey-predator systems. Yet despite this focus on interactions at the pairwise level, a general and comprehensive framework for two-way interaction categorizations has not been established. Moreover, this lack of a framework makes it particularly challenging to address and incorporate higherorder interactions that must be scaled up from any specific pairwise interaction model. This also likely contributes to the paucity of studies or higher-order interactions—defined here as interactions among three or more components—in the literature. Indeed, many studies have hypothesized or either implicitly or explicitly assumed that higher-order interactions are extremely rare and/or insignificant (May, 1972; Van Belle, 2011; Wang et al., 2015; Wootton and Stouffer, 2016). In contrast, recent studies have provided evidence that there is a large amount of higher-order interactions, suggesting a critical need to include higher-order effects to better understand complex systems and alterations of ecosystem processes (Weinreich et al., 2013; Taylor and Ehrenreich, 2015a,b; Beppler et al., 2016; Tekin et al., 2016; Levine et al., 2017; Mayfield and Stouffer, 2017).

Constructing an integrated theoretical framework that encapsulates all the pairwise interaction measures and rigorously constructs higher-order interaction measures by building up from a generalized formula for pairwise interactions would be extremely valuable for gaining insights into understanding complex systems. Here, we provide a perspicuous path to this general theory of interaction measures at pairwise and higherorder levels. We aim consequently to enhance the understanding of interactions via the comparison of different interaction measures and via the search for higher-order interactions, potentially providing profound insights into complex systems research.

In **Table 1** we give several of the most common choices for interaction measures that have been introduced and utilized in diverse fields. Diverse examples of interaction measures include (i) covariance for calculating the joint variability of multiple random variables (Rice, 2003), (ii) mutual information (MI) defined to quantify the information gained about one random variable through knowledge of the other random variable (Cover and Thomas, 2012), (iii) joint cumulants in statistical physics, also known as an n-point correlation or Ursell functions in quantum field theory (Kendall and Stuart, 1969), (iv) statistical measures of analysis of variance (ANOVA) for detecting the differences of means and variances among different groups (Cohen, 2008), (v) additive and multiplicative models introduced for quantifying the interactions among multiple predators in their ability to affect the survival of a prey population (**Figure 1**) (Sih et al., 1998), and (vi) a commonly used measure of Bliss or epistatic Independence for identification of drug-drug and gene-gene interactions (Bliss, 1939). Although these measures are all grounded on the unifying objective of quantifying the dependence or interaction between different components (referred to as objects, variables, or factors), a general theory for deriving these measures does not exist. Indeed, some of the interaction measures mentioned above are strongly associated with each other (see Similarities across different interaction measures), hence introduction of a general framework of interactions becomes much more needed.

Identifying higher-order interactions is essential because the behavior of complex systems can be unpredictable due to interactions among the pairwise subsets of system components (Case and Bender, 1981; Billick and Case, 1994). However, generalizing 2-way (pairwise) interaction measures to three or more component systems is often challenging, often done incorrectly or incompletely, and requires a heuristic process of learning the complex calculations or directly implementing software packages with implicit choices and computations. Consequently, measures for higher-order interactions are often not known or not well-defined (Tekin et al., 2017). All in all, a simplified yet rigorous approach for studying pairwise and higher-order interactions is vitally important for creating a basis to explore the patterns and consequences of emergent phenomena in a wide range of systems.

One important point to consider is that proceeding from lower-order to higher-order interaction measures requires the resolution of ambiguities concerning what exactly higher-order interactions represent (Foucquier and Guedj, 2015). In this respect, two recent studies on higher-order drug interactions with a major focus on three-way interactions—in terms of data representation and analysis—shed light by distinguishing between net and emergent higher-order interactions (Beppler et al., 2016; Tekin et al., 2016). A net interaction, which is the more commonly considered type of higher-order interaction, refers to the total interaction that results from effects at all levels. In contrast, emergent interactions arise only when all of the components (or objects) or some subset of components are combined. Similarly, a more comprehensive and general approach is needed to evaluate higher-order interaction effects in a wide range of fields. From an ecological perspective, when conservation management is the goal for a prey population, understanding the actual interactions among all predators can be critical. We cannot assume that more predators will increase prey risk, and similarly, without proper emergent interactions analysis, we do not know which predator removal or decrease will ultimately increase a prey's number. Therefore, a generalized approach for quantifying higher-order interactions must directly distinguish between net and emergent interactions and incorporate effects resulting from different levels of interaction [where the level corresponds to the number of components in any given subset of component, ranging from one component (level 1) up to the whole system (level N)]. This higher-order interaction framework would require data for responses for each subsets of components—singles, pairs, triples, etc.

In this paper, we first establish a general framework for analyzing two-way interactions and show how our framework can be used to derive the many interaction measures described above. Indeed, we show how all the measures mentioned above (such as mutual information, ANOVA, multiplicative models, etc., see **Table 1**) fit within this framework. Next, we present a conceptual and theoretical advance for building higher-order interaction measures that starts from a general framework of

TABLE 1 | Summary table of two-way interaction measures with our generalized formula TXY = F (X ◦Y) − (F (X ◦ I) △F (Y ◦ I)).

pairwise interactions and the definition of no-interaction at the pairwise level. In doing so, we show how different interaction measures are generated by different choices of mathematical operations and functionals, potentially allowing exploration of all possible interaction measures via all possible choices of operations and associated identity elements. Importantly, our general theory of interactions facilitates the understanding of connections between different measures that are frequently used in different research areas and will thus be a useful guide for studies of interactions in many distinct fields.

### Similarities Across Different Interaction Measures

Despite the fact that there are a wide range of choices for pairwise interaction metrics (see **Table 1**), interactions are defined generically based on a prediction or expectation that signifies no interaction and corresponds to each respective measure having a value of 0. Given the no-interaction expectation—equivalently referred to as additive, multiplicative, or independence expectations—the classification of interactions is evaluated based on deviations from this no-interaction case. When the combined effect is sufficiently greater or weaker than the expected effect of no-interaction, the interaction is classified as positive (synergy or cooperation) or negative (antagonism or interference), respectively (see **Figure 1** for a toy example and an ecological case study). Interpreting the magnitude of interaction often requires rescaling (normalizing), using methods such as those developed by Segrè et al. (2005) and Tekin et al. (2016), and then testing the significance by various statistical analysis methods such ast-tests, bootstrapping experiments, and Bayesian statistics (Jakulin, 2005).

Importantly, some interaction measures (as briefly overviewed in the Introduction and presented in **Table 1**) are strongly associated with each other in terms of their construction and baseline expectation of no-interaction. For instance, the underlying additive model of ANOVA is strongly linked with the additive model of Multiple Predator Effects (MPE) studies (Sih et al., 1998), where the expected effect of combination is sum of their individual effects. On the other hand, there is a correspondence between the multiplicative model of MPEs—the expected effect of combination is product of their individual effects—and the Bliss Independence model of drug interactions and epistasis (Beppler et al., 2016). These correspondences suggest that further exploration and potential applications of interaction measures to different systems could be fruitful. Our general framework uncovers these similarities and hence proves to be very useful in understanding the interaction measures and their similarities as well as their applications to different systems.

When proceeding from simpler forms of interactions (i.e., pairwise interactions) to higher-order interactions (i.e., threeway, four-way, N-way interactions), two different types of interaction classifications arise, namely net and emergent interactions (Beppler et al., 2016). Net interactions measure whether any effective interaction exists at all, whereas emergent interactions measure whether there are interactions beyond what is expected from the "sum" of the lower-order parts. Inherently, characterizing higher-order interactions necessitates data for all the system outcomes (i.e., responses) in the presence of all the subsets of components. For example, for the characterization of three-way interaction, one would need responses in the presence of each single component alone, responses under all pairwise component combinations, and responses under all three components. Overall, better understanding of how higher-order interaction measures translate into these two different types of higher-order interactions is vital to better predict dynamics of complex system cooperations.

### METHODS

### Generalized Form of Interaction Measures General Theory of Two-Way Interactions

Here, we introduce a mathematical framework that generalizes and unifies interaction measure formulations for two

FIGURE 1 | Example representation and ecological case study of 2-way and 3-way interactions for consumer-resource systems. A toy example of an interaction study design for multiple predator effects (MPE) is given to represent 2-way and 3-way interactions, where mortality rates of prey population (small fish) are dependent on the consuming ability of certain predator or combination of predators. Schematic of no predator, single predator, pairwise combination, and triple combination of fish predator populations are illustrated in 8 different panels, where the pairwise panels of XY, XZ, YZ are chosen to represent the three distinct types of interactions: synergy, antagonism, and no-interaction by the multiplicative model with the rescaling introduced in Segrè et al. (2005) and Tekin et al. (2016). In particular, synergistic interaction corresponds to the case when the mortality rate with two predators is sufficiently larger than the expected mortality rate of predators based on single-predator effects. Conversely, predators can act antagonistically, meaning that the mortality rate in the presence of combined predators is sufficiently less than the expected mortality rate based on single-predator effects. In the latter scenario, predators do not interact, hence expected prey mortality is the same as actual prey mortality with both predators in the environment. Classification of interactions is typically done by evaluating the effects of the combined components relative to the effects of independent single components or lower-order combination effects. Given the raw data (see table), we employ the analysis of variance (ANOVA) as well as multiplicative model of interactions and present p-values and interaction measure calculations (T, see section Methods for the definition), respectively. For 3-way combinations, two distinct types of interactions are characterized by the multiplicative model, where net interaction represents the overall effect and emergent interaction represents the three-way combination effect that is beyond any pairwise interactions. The same interaction study design has been performed for an ecological case study by Coors and De Meester (2008), where stressors are a predator (X), parasite infection (Y), and pesticide exposure by carbanyl (Z) on a water flea Daphnia magna survival. A MATLAB code package that has the toy data with the implementation of the interaction measure calculations are included as part of our paper (see Data Sheet 1).

components (equivalently referred to as factors or objects). In a broad context, interactions are determined based on the evaluation of the combined effect relative to a null expectation of no-interaction based on the single-component effects. Hence, for defining a generalized formula for interactions, it is essential to introduce notation, properties of operations for combining components, and formulations for assessing the effects of combined or single components. Here, we define an algebraic operation ◦ (associative and commutative) that combines two components X and Y and that has identity element I, and a real-valued functional, F, (a function of functions) and a second algebraic operation △ (associative and commutative) that combines this functional with different arguments (**Figure 2**). Then, we express the definition of nointeraction as F (X ◦Y) = (F (X ◦ I) △F (Y ◦ I)), and define the generalized measure for quantifying interactions between two components as

$$\mathcal{Z}\_{XY} = \mathcal{F}\left(X^{\circ}Y\right) - \left(\mathcal{F}\left(X^{\circ}I\right) \bigtriangleup \mathcal{F}\left(Y^{\circ}I\right)\right) \tag{1}$$

The identity element in the functional form of the interaction measure is provided in the generalized formula to explicitly make the point that the absence of another component should yield no interaction. In other words, combining any component (X) with the identity element (Y = I) should be equivalent to no interaction, corresponding to I = 0 for addition and to I = 1 for multiplication. Moreover, the arithmetic operations are restricted to either addition or multiplication to preserve the associativity and commutativity properties of combining more than two components. This preservation property is because there is no time component in our framework that corresponds to introducing different components to the system. In other words, the ordering of different components in the interaction definition does not matter, hence TXY = TYX for two components and TXYZ = TYXZ = TZXY for three components (see next section). By definition, TXY = 0 generally means there is no interaction or dependence between two components, whereas deviation of TXY from zero suggests an interaction. Different choices of operations and functional F correspond to different choices of interaction measures as derived below. Furthermore, the sign and magnitude of TXY indicates the type and strength of interaction after proper rescaling (Segrè et al., 2005; Tekin et al., 2016). Rescaling has a firm conceptual foundation that has been proven useful in many areas and is defined relative to important baselines, such as a maximally synergistic scenario or special cases of antagonism, to appropriately reflect the strength of any interaction (Segrè et al., 2005; Sanjuán and Elena, 2006; Tekin et al., 2016). Similarly, two-way and higher-order interaction measures need a rescaling method (a normalization procedure) to distinguish different interaction categorizations.

### General Theory of Higher-Order Interactions

As we proceed from two-way to higher-order interactions, there are two distinct types of interactions: (1) net interactions and (2) emergent interactions (**Figure 2**). The first form measures whether an effective or net interaction exists at all, hence defined as an effect of a combination that is different than predicted from solely the effects of non-interacting single components. On the other hand, the second form measures whether there are emergent interactions beyond what is expected from the lower-order parts of the whole combination. Below, by following a similar notation as for two-way interactions, we introduce the functional forms of generalized higher-order interaction measures. For that, we first note that the ordering of the components (or objects) does not matter as we assume the absence of a time component in our framework. Therefore, the interaction measures are symmetric with respect to each component—hence X (Y ◦Z) is equivalent to X ◦Y ◦Z or X ◦Z ◦Y and vice versa—because algebraic operations of ◦ and △ hold commutativity and associativity properties.

### **Net interactions**

Net interactions are measured by a direct extension of the generalized formula for two components (**Figure 2**) (Sanjuán and Elena, 2006; Staude et al., 2010; Beppler et al., 2016). For example, when three components are completely independent or do not interact at any level, then the functional with the combined components as the argument is equal to the combination of functionals of single components under the operation △. Therefore, adding a third component (Z) into the generalized version of the two-way interaction measure (Equation 1), the generalized version of the three-way net interaction measure becomes

$$\mathcal{Z}\_{\text{XYZ, NET}} = \mathcal{F}\left(\mathbf{X}^{\complement}\mathbf{Y}^{\complement}\mathbf{Z}\right) - \left(\mathcal{F}\left(\mathbf{X}^{\complement}\mathbf{I}\right)\triangle\mathcal{F}\left(\mathbf{Y}^{\complement}\mathbf{I}\right)\triangle\mathcal{F}\left(\mathbf{Z}^{\complement}\mathbf{I}\right)\right) \tag{2}$$

Next, we derive the net interaction formulation for a combination of an arbitrary number, N, of components. Representing the set of components by X1, X2, · · · , XN, and assuming there is no interaction between component X1 and the rest of the components, we have

$$\mathcal{F}\left(X\_1^\complement X\_2^\complement \cdots \curvearrowright X\_N\right) = \mathcal{F}\left(X\_1^\complement I\right) \triangle \mathcal{F}(X\_2^\complement \cdots \curvearrowright X\_N^\complement I)$$

By induction or equivalently applying the same realization repeatedly to the right side of the above equation yields the N-way version of the generalized formula as given by

$$\mathcal{F}\left(\mathbf{X}\_1^{\odot}\mathbf{X}\_2^{\odot}\cdots\mathbf{X}\_N^{\odot}\right) = \mathcal{F}\left(\mathbf{X}\_1^{\odot}I\right)\Delta\mathcal{F}\left(\mathbf{X}\_2^{\odot}I\right)\Delta\dots\Delta\mathcal{F}\left(\mathbf{X}\_N^{\odot}I\right) \tag{3}$$

Using a big △ notation to denote the application of the operation △ to a sequence of numbers—similar to sigma (or pi) notation for summation (or multiplication)—our generalized formula for quantifying net N-way interaction is given by

$$\mathfrak{T}\_{X\_1 X\_2 \dots X\_N, \text{ NET }} = \mathcal{F} \left( X\_1^{\odot} X\_2^{\odot} \dots \stackrel{\circ}{\ } X\_N \right) - \underset{\text{i}}{\stackrel{\Delta}{\ }} \mathcal{F} \left( X\_i^{\odot} I \right)$$

### **Emergent interactions**

Emergent interactions are defined as interactions that exist even after the exclusion of all interaction effects that are due to lower-order parts and may contribute to a combined overall (net) interaction (**Figure 2**). In the drug interaction and epistasis terminology, emergent interactions are concretely defined in our recent work and clearly contrasted with the definition of net interaction (Beppler et al., 2016; Tekin et al., 2016). Specifically, an emergent interaction is an interaction that exists beyond the effects that are due to all lower-order parts, whereas a net interaction is any effect beyond only the individual component effects. In an article by Staude et al. (2010) on covariance,

FIGURE 2 | Generalized form of interaction diagrams and measures. Two-way and three-way interaction diagrams and measures (T) over distinct interaction types, general functionals (F), and operations (, △) are given. Empirically, classification of pairwise interaction requires measurements from each single component alone as well as measurements from each pairwise combination of components. On the other hand, classification of three-way interactions requires measurements from all possible subsets: measurement under each single-component alone, measurement from all pairwise subsets, and measurement from triple-component combination. Theoretically, for two-way interactions, the definitions of net and emergent interactions are identical. On the other hand, net and emergent interactions are not the same at the three-way level because the net three-way combination subsumes the pairwise combination effects of each pair (XY, XZ, YZ). Therefore, the emergent 3-way interaction (the diagram on the bottom far left) is expressed by subtracting all lower-order contributions from the net 3-way interaction. The interaction classification is obtained from the sign of the interaction measure, T, where T = 0 represents no-interaction, T < 0 represents synergy (positive interaction), and T > 0 represents antagonism (negative interaction). See Figure 1 for a toy example of consumer-resource systems that depicts different interaction classes: synergy, no-interaction, and antagonism.

this concept is referred to as the dependence in the higherorder combination that is not embedded through the lowerorder correlations. In information theory, Jakulin et al. (Jakulin, 2005) referred to this idea similarly and stated that interaction information among different attributes can be gained by the information that is not present in any subset. By definition, at the two-way level, an emergent interaction is identical to the net interaction because there is nothing from which to emerge except the single-component effects. On the other hand, at the next level, the emergent three-way interaction measure captures the interaction that does not originate from the pairwise parts, and hence emerges only in the presence of all three components. For more than three components, it matters what level of lowerorder interactions needs to be excluded. There can be different choices based on the specific system or the main purpose of the interaction identification, as revealed in more detail below and in the **Appendix A** in Supplementary Material.

As a start, to measure the emergent three-way interaction, we calculate how much of the three-way net interaction originates from each of the isolated pairs and subtract that from the net (or total) three-way interaction. In this case, the contribution that comes solely from the pair X, Y represents the situation that only X, Y interacts within the three-way combination of X, Y, and Z. In other words, the remaining component Z within the combination does not interact with the pairwise part of X and Y or either X or Y individually. Employing the two-way interaction formula of generalized interaction measure (Equation 1) and the associative properties of operations, this is equivalent to saying that F (X ◦Y ◦Z) = F ((X ◦Y ◦ ) ◦Z) = F ((X ◦Y) ◦ I) △F (Z ◦ I). Hence, the two-way combination of components X, Y contributes to the entire (net) three-way interaction (Equation 2) by an amount of F ((X ◦Y) ◦ I) △F (Z ◦ I) −F (X ◦ I) △F (Y ◦ I) △F (Z ◦ I), which is equal to F (Z ◦ I) △ TXY when △ is multiplication and is equal to TXY when △ is addition. Similarly deriving the contribution coming solely from the interaction of components Y and Z and also from X and Z, the generalized formula of the emergent three-way interaction is given by disentangling the sum of the pairwise parts from the three-way interaction. Therefore, given the multiplication operation of △, the emergent 3-way interaction among components X, Y, and Z is

$$\mathfrak{T}\_{XYZ,\text{ EMERGINT}} = \mathfrak{T}\_{XYZ,\text{ NET}} - \mathcal{F}\left(X^{\circ}I\right) \bigtriangleup \mathfrak{T}\_{YZ} - \mathcal{F}\left(Y^{\circ}I\right)$$

$$\Delta \mathfrak{T}\_{XZ} - \mathcal{F}\left(Z^{\circ}I\right) \bigtriangleup \mathfrak{T}\_{XY} \tag{4}$$

On the other hand, the weighting parameters (e.g., F (X ◦ I)) disappear when △ is the addition operation, so an emergent interaction is given by

$$\mathfrak{T}\_{XYZ, \text{ \textquotedblleft EMGRENT}} = \mathfrak{T}\_{XYZ, \text{NET}} - \mathfrak{T}\_{YZ} - \mathfrak{T}\_{XZ} - \mathfrak{T}\_{XY} \tag{5}$$

These expressions mean that when some interaction exists but all the pairs are independent (or not interacting), then the emergent interaction is equal to the net interaction because any interaction can only be coming from the combined effects of all three drugs for this special case. When pairwise interactions exist and do not cancel, the net and emergent interaction measures will virtually always differ. Notice also that all interaction measures are symmetric with respect to each component. Moreover, when a single component is independent (or not interacting) with the others, the emergent three-way measure vanishes, as desired.

For N components (X1, X2 ... XN), the emergent N-way interaction formula can be derived similarly by starting from the definition of the generalized formula for net N-way interactions and subtracting the appropriate lower-order effects. The derivation for generalized formula of emergent interactions with four components is described in detail in **Appendix A** in Supplementary Material.

### RESULTS

### Construction of Specific Two-Way Interaction Measures From General Theory

We will now derive how our generalized framework of two-way interactions leads to numerous interaction measures introduced in studying generic multi-component systems (**Table 1**). For this, we start by defining X and Y as random variables and will show that the generalized formula reduces to the covariance formulation and mutual information by the choice of different algebraic operations of and △. Given that F is the expectation operator, and both and △ are the operation of multiplication, TXY (Equation 2) becomes covariance and hence measures the joint variability between X and Y, as

$$\mathcal{Z}\_{XY} = E\left(X \ast Y\right) - E\left(X \ast 1\right) \ast E\left(Y \ast 1\right) = E\left(XY\right) - E\left(X\right)E\left(Y\right)$$

Here, we note that when X and Y are independent, their covariance is zero, i.e., TXY = 0. Notably, the two-way covariance measure is analogous to concepts in theoretical physics such as joint cumulants or Ursell functions (Kendall and Stuart, 1969) and the Isserlis formula in statistics (Isserlis, 1918). However, the higher-order interaction formulations of these concepts differ from the higher-order covariance when there are more than three components, as derived in **Appendix A** in Supplementary Material.

Following the covariance derivation using our general theory of two-way interactions, we now show that choosing appropriate functional and algebraic operations leads to the mutual information measure of interactions. Hence, we define F to be negative entropy, △ to be addition, and ◦ to be addition (i.e., I = 0). Here, the ◦ operation represents the combination of two components and hence defines joint entropy under the functional F. For these choices, the generalized formula yields the mutual information

$$\mathfrak{T}\_{XY} = -H\left(X, Y\right) + H\left(X\right) + H\left(Y\right)$$

Here, associating the random variables X and Y with attributes, the interaction formulation TXY measures the correlation between the attributes based on concepts of information theory. By definition, when X and Y have no interaction, their mutual information vanishes because no information is gained about one attribute through the other attribute (Jakulin, 2005).

Proceeding to cases for which the arguments X and Y represent experimental data rather than standard random variables, we discuss how interactions are measured by the formula TXY by evaluating the effects of combined components inferred from system responses (F) such as growth, survival, or mortality rates (**Figure 1**). These experimental data types and response functions are used in many important and frequently-used measures of interactions, such as ANOVA, epistatic additivity, Bliss Independence, and multiple predator effects measures.

We begin by reviewing how ANOVA (Analysis of Variance) translates into a measure for the prediction of interactions and then deriving its correspondence with our general interaction formula. In the case of two-component systems, ANOVA interactions are quantified by a 2 by 2 factorial design with four different observations—corresponding to no components (0), component X alone, component Y alone, and components X and Y together. In this case, two-way ANOVA tests whether the presence of Y affects the impact of X according to a measure based on linear deviations. When the effects of observations are denoted by S, no interaction is described as S (X) −S (0) = S (X + Y) −S (Y). Note that this equation can be re-expressed in the form S (X + Y) −S (0) = S (X) −S (0) + S (Y) −S (0). Therefore, redefining the response measurement to be F (X):=S (X) −S (0), the interaction is quantified by the significance of deviation from the additive effects of two components, i.e., F (X + Y) −F (X) −F (Y). This corresponds to a general formula of interactions, TXY, with ◦ and △ being addition operations

$$\mathfrak{T}\_{XY} = \mathcal{F}\left(X + Y\right) - \mathcal{F}\left(X\right) - \mathcal{F}\left(Y\right)$$

Regarding the use of ANOVA for interaction classifications, work by Billick and Case (1994) and Wootton (1994) demonstrate that the results of ANOVA can be misleading if the transformation of data is not carefully done or when the empirical system does not align with the assumptions of ANOVA.

Following the similar component types and response measurements, our general interaction measure produces another widely-used model of interactions, namely the multiplicative-risk model. This measure has been introduced by Soluk and Collins (1988) for quantifying multiple predator effects (MPE), and it overcomes the ambiguities of the additive prediction by ANOVA when the single components have large impacts on the response function (F). In such cases, the additive expectation can never be achieved because the predicted combined effect is stronger than is feasible or possible for the two together (for example you cannot have a stronger effect than complete killing, see **Appendix B** in Supplementary Material for details). The multiplicative model formulates no-interaction as F (X + Y) = F (X) F (Y), based on the idea of proportional decreases in population size. It is thus much more akin to Bliss Independence discussed below. Importantly, the multiplicative model of interaction can be written as a general interaction formula, where ◦ and △ are addition and multiplication operations, respectively. Hence,

$$\begin{aligned} \mathfrak{T}\_{XY} &= \mathcal{F}(X+Y) - \left(\mathcal{F}(X+0) \* \mathcal{F}(Y+0)\right) \\ &= \mathcal{F}\left(X+Y\right) - \mathcal{F}\left(X\right)\mathcal{F}\left(Y\right) \end{aligned}$$

associates the null model of no-interaction based on the product of response measurements (F), which is prey survival rates for the MPE system (Sih et al., 1998; Barrios-O'neill et al., 2014). The subscript notation (such as FXY) is still typically used to express response measurements under different combinations of predators, or similar component types in different settings, even though X and Y are actually being summed together in this case. Hence, the equation above takes a form TXY = FXY − FXFY. As explained in **Appendix B** in Supplementary Material, this multiplicative model is preferred over additive models that do not always yield plausible expectations of combination effects. As a further note, this multiplicative model of two-way interaction can be tested by merely applying twoway ANOVA on log-transformed data because testing for log F (X + Y) = log F (X) + log F (Y) is equivalent to testing for F (X + Y) = F (X) F (Y) by properties of logarithms. However, the multiplicative model at higher-order levels is not simply equivalent to the log-transformed ANOVA as detailed below. Note that following the same logic of additive and multiplicative models, where the interaction effect signifies the deviation from the null model of additivity or multiplicativity, alternative measures to two-way ANOVA [namely Hedge's d and log response ratio (LnRR)] have been proposed to determine interactions based on the deviation from the corresponding model expectation standardized to the effect size (Gurevitch et al., 2000; Crain et al., 2008; Côté et al., 2016). Similar correction terms for the effect sizes can be added to net and emergent higher-order interaction measures for meta-analysis studies.

Intriguingly, there is a strong correspondence between the multiplicative model introduced for MPEs and the frequently used drug-interaction model of Bliss Independence. This indicates that the Bliss measure can also be expressed by utilizing the general interaction formulation we propose here, i.e., TXY = FXY − FXFY. In drug interaction studies, response measurements, F, are given as growth rates of pathogens in the presence of single or multiple drug environments relative to the no-drug environment. Relative growth is commonly referred to as relative fitness and is indeed analogous to the survival rates of prey species in predator environments (Beppler et al., 2016). Bliss Independence categorizes interactions based on deviations from a null expectation (no-interaction case) that the percent decrease of the pathogen growth rate in the presence of drug X is not affected by the presence of the other drug Y. Analytically, the deviation from no-interaction instance is measured by F (X + Y) −F (X) F (Y) or FXY−FXFY, following a subscript notation as in the multiplicative model.

### Building Higher-Order Specific Interaction Measures From General Theory

In this section, we discuss the widely-used pairwise interaction measures described above and summarized in **Table 1**, and we derive the three-way interaction measures based on the corresponding no-interaction criterions and our general theory. For clarity, we follow the same order of measures as in the section on the construction of two-way interactions. Therefore, we start by reviewing the standard interaction measures when the components are random variables. In probability theory, the extensions of probabilistic measures of moments to multivariable systems are generally referred as higher-order cross moments. For example, adapting the generalized formula of covariance (see **Table 1**) to three variables yields a 3rd-order cross moment known as co-skewness that is used as a risk management measure in finance or in social problems (Cox L. A. Jr., 2009; Miller, 2013). To be more explicit, the net interaction measure (Equation 3) is given by TXYZ, NET = E (XYZ)−E (X) E (Y) E(Z) (Stratonovich, 1967; Staude et al., 2010), which vanishes when X, Y, and Z are mutually independent with each other and each pairwise product. Next, substituting the generalized formula of two-way interactions (i.e., covariance) and the net three-way interaction, the emergent interaction (Equation 4) becomes

$$\begin{aligned} \mathfrak{T}\_{XYZ, \text{ }\text{EMRGENT}} &= E\left(\text{XYZ}\right) - E\left(\text{X}\right)E\left(\text{YZ}\right) - E\left(\text{Y}\right)E\left(\text{XZ}\right) \\ &- E\left(\text{Z}\right)E\left(\text{XY}\right) + 2E\left(\text{X}\right)E\left(\text{Y}\right)E\left(\text{Z}\right) \end{aligned}$$

which is equivalent to the three-way covariance [i.e., E ( (X − E (X)) (Y − E (Y)) (Z − E (Z)))] or the non-normalized formula of co-skewness as discussed above. Notably, this threeway measure is equivalent to the Isserlis formula (Isserlis, 1918), 3-point connected correlation, or equivalently the Ursell function with three random variables (Kendall and Stuart, 1969; Staude et al., 2010). However, this does not hold true for more than three components (**Appendix A** in Supplementary Material).

Continuing with the component types as random variables, we now show that the general formula of higher-order interactions reduces to two interaction measures of paramount importance in information theory. Defining the functional F as the negative of entropy and following the algebraic operations that yield mutual information (see **Table 1**), the net interaction formulation coincides with the total correlation measure as given by

$$\mathfrak{T}\_{\text{XYZ,NET}} = -H\left(X, Y, Z\right) + H\left(X\right) + H\left(Y\right) + H(Z)$$

This measures the total dependence among three attributes X, Y, and Z (Watanabe, 1960; Jakulin, 2005). By incorporating the pairwise mutual information into the emergent interaction formulation when △ is the addition operation (Equation 5), we attain the multiple mutual information measure defined by McGill (McGill, 1954; Sun Han, 1980; Jakulin, 2005).

$$\begin{aligned} \mathfrak{T}\_{XYZ, \text{ EMERGENT}} &= -H(\mathcal{X}, Y, \mathcal{Z}) + H(\mathcal{X}, Y) + H(\mathcal{X}, \mathcal{Z}) \\ &+ H(Y, \mathcal{Z}) - H(\mathcal{X}) - H(Y) - H(\mathcal{Z}) \end{aligned}$$

Thus far, we have introduced the commonly used interaction and dependence measures from probability theory. Next, we move onto the standard measure of ANOVA, which analyzes the effects of multiple components via the factorial experiment design. As in the two-way interactions, we simply use the underlying additive model of three-way ANOVA and refer to it as just an additive model (**Table 1**). When the response function F is defined relative to the control (no component is present), the threeway net interaction formula of the additive model is given by F (X + Y + Z) −F (X) −F (Y) −F (Z). Now, substituting each pairwise contribution in the generalized formula of emergent interactions (**Figure 2**), the emergent interaction measure is identical to the underlying model of three-way ANOVA.

$$\begin{aligned} \mathfrak{T}\_{XYZ, \text{ EMERGENT}} &= \mathcal{F}(X+Y+Z) - \mathcal{F}(X+Y) - \mathcal{F}(X+Z) \\ &- \mathcal{F}\left(Y+Z\right) + \mathcal{F}\left(X\right) + \mathcal{F}\left(Y\right) + \mathcal{F}\left(Z\right) \end{aligned}$$

This correspondence is proven explicitly in the **Appendix C** in Supplementary Material by recognizing that three-way ANOVA determines how the two-way interaction between the components X and Y is affected by an addition of a third factor Z. This in-depth derivation of the underlying model of three-way ANOVA helps to compare the statistical measure of ANOVA with other interaction classification methods.

Due to the additive model's implausible predictions under certain cases, many MPE studies have pursued the implementation of the multiplicative model (**Appendix A** in Supplementary Material). Subsequently, we now derive the definition of three-way interactions under the multiplicative model of MPE studies or equivalently the Bliss Independence formula. This derivation is very similar to the foundation of the covariance for three random variables as △ is the multiplication operation in both concepts (**Table 1**). Defining the functional F as survival rate relative to the control, the net interaction measure becomes

$$\mathfrak{T}\_{XYZ,\ NET} = \mathcal{F}\left(X + Y + Z\right) - \mathcal{F}\left(X\right)\mathcal{F}\left(Y\right)\mathcal{F}\left(Z\right)$$

or TXYZ, NET = FXYZ − FXFYFZ using subscripts for denoting the combination of components. In this case, substituting the terms corresponding to each pairwise contribution, such as FXYFZ − FXFYFZ = FZFXY for XY combination, the emergent interaction becomes

$$\begin{aligned} \mathfrak{T}\_{XYZ,\text{ EMERGENT}} &= \mathcal{F}\_{XYZ} - \mathcal{F}\_X \mathcal{F}\_{YZ} - \mathcal{F}\_Y \mathcal{F}\_{XZ} - \mathcal{F}\_Z \mathcal{F}\_{XY} \\ &+ 2 \mathcal{F}\_X \mathcal{F}\_Y \mathcal{F}\_Z \end{aligned}$$

Using a fitness definition for F, TXYZ, NET and TXYZ, EMERGENT correspond to the net and emergent interaction formulas with three components in the fields of drug interactions and gene epistasis (Sanjuán and Elena, 2006; Beppler et al., 2016).

### DISCUSSION

Multi-level interactions play an important role in explaining the characteristics of systems that arise as collective behaviors or responses among different components. As an example in conservation biology, mammalian herbivory, invertebrate herbivory, and nutrient levels all interact to affect the survival, growth rates, and fecundity of rare plant species in a forest ecosystem (McGill et al., 2006; Savage et al., 2007; Webb et al., 2010; Dávalos et al., 2014; Enquist et al., 2015). Given this importance across diverse fields, many different interaction classification methodologies have been introduced. However, a unified theory behind these choices has been lacking, and as a result, higher-order interaction measures are often incorrectly derived and interpreted and thus higher-order interactions are not well understood. Moreover, the lack of a general theory behind the pairwise interaction measures makes it difficult to establish a common understanding of higher-order interactions and to design well-defined measures that go beyond pairwise interactions. In our study, we introduce a pioneering approach to resolve these uncertainties and for studying interactions in complex systems. We further propose a general theory of interaction measures to uncover interactions at pairwise and higher-order levels. Our framework provides an advance in such a way that specific cases and measures can be derived from the definition of what an interaction means—the deviation from an expectation that there is no interaction—and what response measurements are needed to evaluate for the interactions—the functional with different arguments. From that we establish a rigorous framework for extending and generalizing these concepts and measures to higher-order levels.

In doing so, we show the utility of our general framework by deriving pairwise interaction measures across different fields covariance, mutual information, cumulants, ANOVA, Bliss independence, and epistatic additivity (**Table 1**). Furthermore, we demonstrate that all pairwise interaction measures can be easily extended to characterize the net interaction—an overall effect resulting from all levels—and the emergent interaction—a measurement of an overall interaction effect relative to interactions at lower-order levels. We also discuss their specific representations, correspondences, and interpretations in a diverse set of study areas. Overall, our advance is of importance for studies of diverse complex systems that rely on a large number of interactions and emerge across multiple levels.

Despite a rich focus on pairwise interactions in the literature, recent studies have shown that higher-order interactions are pervasive in drug and gene systems (Taylor and Ehrenreich, 2015a,b; Tekin et al., 2016), and that the typical approach of focusing solely on pairwise interactions is most likely inadequate for successfully attaining important insights into research on complex systems (Weinreich et al., 2013; Beppler et al., 2016; Levine et al., 2017; Mayfield and Stouffer, 2017). Therefore, higher-order interactions (i.e., among any number of components) must be incorporated to fully grasp system dynamics. In this respect, our study provides a promising framework to uncover emergent phenomena in multi-component systems such as protein and gene interaction networks, multiple predator effects, food webs, interacting agents in economics, voting behaviors in political science and the cohesion dynamics of social groups.

As shown by our previous studies (Yeh et al., 2006; Tekin et al., 2016), empirical data of drug interactions separate into a trimodal distribution—clearly separating synergy, nointeraction, and antagonism—and clearly different than would arise from a random model of interactions because stochastic effects would have smoothed out this distribution and led to a more uniform distribution. Furthermore, in our previous study (Tekin et al., 2016), we assessed the statistical significance of interactions by testing against the null hypothesis of random interactions that would on average lead to a uniform distribution that ranges from −1 to 1 across the entire interaction scale. To quantify these comparisons, we used a Silverman test (Silverman, 1981) and found the higher-order emergent interactions are significantly different from random effects. These results reveal higher-order interactions are not simply due to random effects. Also worth noting is that uncertainty and error estimation due to measurements can be dealt with using bootstrapping strategies (Cruz-Loya et al., 2018).

Throughout the text, we also point out a well-founded concept of "rescaling" in physics that has been proven to enhance the categorization of pairwise and higher-order emergent effects (Segrè et al., 2005; Tekin et al., 2016). Without employing an appropriate rescaling approach, the magnitude of interactions is often not apparent and hence can lead to misleading results and erroneous insights. Therefore, appropriately understanding and rescaling the interaction formulations is an essential part of studying and understanding multi-component interactions in any system.

We further note that many systems have additional complexities besides just many interacting components or higher-order interactions. For example, in predator-prey systems, both predators and prey can evolve, and there are direct feedbacks in the interactions. Both of these effects can lead to drastic changes in multiple predator effects (MPEs) (Barrios-O'neill et al., 2014). In addition, pathogens can evolve in response to drug combinations in ways that are almost impossible to predict due to rugose fitness landscapes. To predict and understand the dynamics of pathogen populations thus requires further theory (Palmer et al., 2015). Our general

### REFERENCES

framework of interactions, with additional modifications, will help address some of the important questions about complexities in biological, physical, and social systems. Overall, understanding emergent features can help us to predict the dynamical consequences of complex interactions, including in questions of crucial significance to human and global health, such as combatting the evolution of resistance to antibiotics and mitigating detrimental impacts of climate change on the diversity and stability of food webs.

### AUTHOR CONTRIBUTIONS

ET, PY, and VS conceived of the project, designed the study, and wrote the manuscript. ET and VS developed new analytical tools. All authors gave final approval for publication.

### FUNDING

ET, PY, and VS received funding from James S. McDonnell Complex Systems Scholar Award, US National Science Foundation (ID: 1254159, http://nsf.gov/), UCLA Faculty Career Development Award, Hellman Foundation Award, and NIH/National Center for Advancing Translational Science (NCATS) UCLA CTSI Grant Number UL1TR001881.

### ACKNOWLEDGMENTS

We thank Nina Singh and Cynthia White for comments on the manuscript.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00166/full#supplementary-material

zooplankton and a larval amphibian. Ecotoxicol. Environ. Saf. 71, 209–218. doi: 10.1016/j.ecoenv.2007.08.007

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Tekin, Yeh and Savage. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Impact of Spatial and Temporal Dimensions of Disturbances on Ecosystem Stability

### Yuval R. Zelnik\*, Jean-François Arnoldi and Michel Loreau

Centre for Biodiversity Theory and Modelling, Theoretical and Experimental Ecology Station, CNRS and Paul Sabatier University, Moulis, France

Ecosystems constantly face disturbances which vary in their spatial and temporal features, yet little is known on how these features affect ecosystem recovery and persistence, i.e., ecosystem stability. We address this issue by considering three ecosystem models with different local dynamics, and ask how their stability properties depend on the spatial and temporal properties of disturbances. We measure the spatial dimension of disturbances by their spatial extent while controlling for their overall strength, and their temporal dimension by the average frequency of random disturbance events. Our models show that the return to equilibrium following a disturbance depends strongly on the disturbance's extent, due to rescue effects mediated by dispersal. We then reveal a direct relation between the temporal variability caused by repeated disturbances and the recovery from an isolated disturbance event. Although this could suggest a trivial dependency of ecosystem response on disturbance frequency, we find that this is true only up to a frequency threshold, which depends on both the disturbance spatial features and the ecosystem dynamics. Beyond this threshold the response changes qualitatively, displaying spatial clusters of disturbed regions, causing an increase in variability, and even a system-wide collapse for ecosystems with alternative stable states. Thus, spanning the spatial dimension of disturbances is a way to probe the underlying dynamics of an ecosystem. Furthermore, considering spatial and temporal dimensions of disturbances in conjunction is necessary to predict ecosystem responses with dramatic ecological consequences, such as regime shifts or population extinction.

### Edited by:

Diego Barneche, University of Sydney, Australia

### Reviewed by:

Bernardo R. Broitman, Centro de Estudios Avanzados en Zonas Áridas (CEAZA), Chile Mariana Bender, Universidade Federal de Santa Maria, Brazil

> \*Correspondence: Yuval R. Zelnik<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 29 August 2018 Accepted: 05 December 2018 Published: 21 December 2018

#### Citation:

Zelnik YR, Arnoldi J-F and Loreau M (2018) The Impact of Spatial and Temporal Dimensions of Disturbances on Ecosystem Stability. Front. Ecol. Evol. 6:224. doi: 10.3389/fevo.2018.00224 Keywords: localized disturbance, rescue effect, bistability, return time, variability, persistence

## 1. INTRODUCTION

Understanding the stability of ecosystems, i.e., their ability to recover and persist in the face of natural and anthropogenic disturbances, is of fundamental importance to ecology and conservation (May, 1973; Neubert and Caswell, 1997; Loreau and de Mazancourt, 2013). Ecosystems are spatially extended, comprised of multiple interacting communities in different locations, and therefore an important factor in understanding their stability is their spatial structure (Levin, 1992; Peterson et al., 1998; Wang and Loreau, 2016). However, while the influence of space on properties such as biodiversity and food web structure has been intensely investigated (Loreau et al., 2001; Chase and Leibold, 2002; Montoya and Sol, 2002; McCann et al., 2005), basic questions regarding spatial stability remain open. In particular, despite the fact that most disturbances (e.g., fires, pest outbreak, pollution runoff) are strongly heterogeneous in space, the impact of their spatial structure on stability is largely unknown. Similarly, their temporal dimension, e.g., their timespan or the frequency of their occurrence, is critical. Taken together, these dimensions span a vast space of possible disturbances that ecosystems can face (e.g., fires and storms). This, in part, explains why reaching a clear understanding of ecosystem stability has proven to be an extremely challenging endeavor.

Research on ecosystem stability has a long history in ecology, and numerous studies have investigated how various properties of disturbances affect ecosystem responses. The importance of spatial properties of disturbances, in particular, has been assessed by a few studies of regeneration dynamics under recurrent, spatially structured disturbances (Turner et al., 1993; Moloney and Levin, 1996; Fraterrigo and Rusak, 2008). These studies introduced the concept of landscape equilibrium and demonstrated how the spatial and temporal scales of disturbances can generate different stability patterns. A point not explicitly addressed in these studies, however, is the importance of rescue dynamics occurring at a regional scale when local recovery processes are too slow or fail altogether. This can occur in sufficiently connected ecosystems, following high-intensity disturbances (Foster et al., 1998; Fraterrigo and Rusak, 2008). In fact, recovery from a disturbance is a consequence of both local and regional processes. Local processes lead to recovery due to dynamics that are internal to local communities (e.g., birth and death of individuals), while regional processes lead to recovery by bringing in individuals from neighboring communities via dispersal (Turner, 1989; Leibold et al., 2004). These two processes mediate the large-scale system response to a disturbance, and their respective parts in this response is bound to strongly depend on the spatial connectivity of the system and, importantly, on the spatial structure of disturbances.

Recent work has made this relationship more explicit, by defining three distinct regimes of recovery from a single spatially heterogenous disturbance: Isolated, Rescue and Mixing (Zelnik et al., 2018). If a system is highly connected due to strong dispersal of organisms, then it is in the Mixing Regime, and the system's behavior at large scales is essentially an extended version of a local system (Durrett and Levin, 1994). At the other extreme, if dispersal is low and hence each site acts separately with its own local dynamics, then the system is in the Isolated Regime, and its large-scale behavior is an aggregation of many independent small systems (Tilman et al., 1998; Yachi and Loreau, 1999). In between these two extremes is the Rescue Regime, where systems with intermediate connectivity show large-scale rescue dynamics due to the interaction between limited dispersal and the system's behavior at the local scale (Peterson, 2000; Dai et al., 2013; Wang et al., 2017). For instance, in the study by Dai et al. (2013), a metapopulation of yeast exhibits a front structure which emerges due to interaction of dispersal with nonlinear local behavior of the yeast. A different example is found in the work of Wang et al. (2017), where the correlations between local bird populations, mediated by dispersal, leads to a spatial scaling law of the variability of populations across North America.

While the spatial structure of both system and disturbance plays no role in the Mixing regime, for weaker dispersal it does: in both the Isolated Regime and the Rescue Regime the spatial structure of the disturbance has significant effects as it can initiate qualitatively different responses that involve both local and regional processes (Zelnik et al., 2018). This is the case in an experimental study of a predator-prey protist system, in which local extinctions are met by rescue processes, which prevent synchronization of the regional metapopulation (Fox et al., 2017). We will therefore consider systems with intermediate dispersal, and focus on the effect of the spatial structure of disturbances as well as their temporal properties.

 If the disturbance is an isolated event, a natural measure to consider is the return time to the unperturbed state (May, 1973; Neubert and Caswell, 1997). On the other hand, in a regime of repeated disturbances (e.g., climatic events), measures of temporal variability are commonly used (Tilman et al., 2006). In the presence of alternative stable states, those repeated disturbances can cause a regime shift from one state to another. One well-known example is that of lake eutrophication (Carpenter, 2005) due to fertilizer runoff disturbances. Here the stability measure of interest is typically persistence, i.e., the probability that a system will remain in a desired state (Holling, 1973; Pimm, 1984). Importantly, these stability measures reflect not only the spatial and temporal properties of the disturbance, but also the dynamical features of the perturbed ecosystem. Exploring this interplay is the focus of our study, which we will address by considering three spatial ecosystem models with increasing nonlinear local dynamics, ranging from logistic growth to bistability. Under various perturbation scenarios we will measure their stability using return time, variability and persistence.

We begin by looking at the ecosystem's recovery following a single disturbance, and show that changing the spatial structure of the disturbance reveals two basic recovery trajectories: isolated and rescue. Isolated recovery trajectories reflect the local resilience of the system, while rescue trajectories involve spatial processes, and their dominance signals the failure of local processes. We thus argue that the relationship between spatial structure and recovery contains substantial information about the local dynamics of the system, both close to and far from equilibrium. We continue by exploring the temporal axis of disturbances, and demonstrate a direct link between return time (following an isolated disturbance event) and temporal variability (under a regime of repeated disturbances). We find that for low disturbance frequency patterns of variability do not contain additional information in comparison to the patterns of return time. However, past a frequency threshold (which depends on the system's internal dynamics) the variability patterns change. As we will argue, this signals the onset of a new dynamical regime driven by disturbances, which can lead to a regime shift—in our case a transition from a populated to a bare state (extinction).

Our work demonstrates that the spatial dimension of disturbances can be used to reveal information on the ecosystem's internal behavior. Furthermore, our results illustrate that the conjunction of the spatial and temporal properties of disturbances may lead to unforeseen dynamical responses, with drastic ecological consequences.

### 2. METHODS

### 2.1. Models

We assume for simplicity that the local community dynamics can be described by a single state variable N that represents the ecosystem's local biomass density. We study the dynamics in multiple locations in space using partial differential equations. We define three different models that differ in their local dynamics but have identical dispersal across space with linear diffusion. In all models the local biomass may reach a carrying capacity K, so that N = K (the populated state) is a stable steady state in all three models. An additional solution exists for N = 0 (the bare state), with its stability properties differing among models.

The first and simplest model (LG) describes local logistic growth coupled with dispersal:

$$N\_t = rN(1 - N/K) + d\nabla^2 N,\tag{1}$$

where Nt is the change in time of the local biomass and ∇ 2N is the second derivative in space of N (a diffusion term). Here r is the characteristic, local dynamical rate of growth, while the rate of spread by dispersal is governed by d. In this model the bare state N = 0 is an unstable solution. This is the classic model of population growth (Hall, 1988), shown to appropriately depict the dynamics of various biological systems, from the growth of unicellular organisms (Gause, 1934), to human populations (Marchetti et al., 1996).

The second model (AE) describes species dynamics with a strong Allee effect (Kramer et al., 2009), so that low biomass densities are not viable. Such dynamics have been found in a variety of species, ranging from the gypsy moth to woodland caribou (Kramer et al., 2009). The model reads

$$N\_l = rN(1 - N/K)(N/\alpha - 1) + d\nabla^2 N,\tag{2}$$

where α is the viability threshold, i.e. the minimal amount of biomass N that is necessary to allow positive growth. This model has two alternative stable states (N = 0, N = K) and a single unstable state (N = α), and we assume that 0 < α < K. This is the simplest model for dynamics with alternative stable states, a property that has been found in many ecosystems (Scheffer, 2009), such as lakes (Carpenter, 2005) and coral reefs (Nyström et al., 2000).

Finally, our third model (SR) describes dynamics with slow recovery following intense disturbances, and stands as an intermediate between the two previous models. It will help us to clarify the distinction between strong nonlinearity and bistability. Its main feature is that while there is only one stable equilibrium at N = K, far from this equilibrium the return rate is very slow compared with the return rate close to equilibrium. This could model succession dynamics, for which the recovery following strong disturbances (e.g., clearcutting) is very slow, as it involves the successive colonization by different species, and not simply the regrowth of the disturbed species (Uhl, 1987), or a weak Allee effect, a prevalent feature in population dynamics (Kramer et al., 2009). The model is:

$$N\_t = rN(1 - N/K)(N/K)^\prime + d\nabla^2 N,\tag{3}$$

where γ controls the nonlinearity of the dynamics, such that at high values of γ local recovery is very slow following highintensity disturbances.

For each model we can define a local potential (see left panels of **Figure 2**), such that its derivative with respect to N corresponds to the derivative of N with respect to time—i.e., the local dynamics. This means that the local dynamics follow the slope of this potential, so that the biomass density can be thought of as a ball moving from peaks to valleys in the landscape that the potential defines. In both the LG and SR models only one stable equilibrium exists, but the speed of return to the equilibrium may be much slower for low biomass density in the SR model. Two stable states exist in the AE model (the populated state and the bare state).

By rescaling time, space and biomass, we can effectively reduce the parameter space, and set r = 1, d = 1 and K = 1. Our results thus hold for any values of these three parameters. We set α = 0.4 to make sure that the AE model recovers from a single disturbance (see next subsection), and γ = 4 to make sure the return time far from equilibrium of the SR model is sufficiently slow. We focus on one-dimensional systems as they are simpler to analyze, but the qualitative results hold for other types of spatial structure such as twodimensional systems (see **Appendix D**). We use a system size of L = 500, which is large enough to allow for the spatial dynamics to manifest itself (so that the system is not in the Mixing Regime Zelnik et al., 2018), with periodic boundary conditions. For a clearer illustration, in **Figure 3** and **Figure S2** we show snapshots of a two-dimensional system of size 200 × 200.

### 2.2. The Spatial Dimension of Disturbances

We define a disturbance as a change in the state variable that is forced on the ecosystem. We consider a pulse disturbance occurring at a given time, with its full effect being applied at that time. This assumption is appropriate for the many types of disturbances that are faster than the dynamics of the ecosystem, and lends itself to mathematical analysis. We choose a disturbance that removes biomass (reduces N), so that a disturbance of strength s will reduce the overall biomass of the ecosystem by sK (but any negative values of N will be set to 0 for consistency). Once a disturbance takes place, the ecosystem may recover to its original state, or a regime shift can occur if the ecosystem is bistable. We are interested here in stability and recovery dynamics, and therefore focus on parameter values for which a single disturbance cannot lead to a regime shift.

Since a disturbance need not occur uniformly across space, we vary the spatial extent of the disturbance σ while keeping its overall strength s constant. A disturbance is performed by choosing its locus, and removing some biomass in a domain of size σ centered around the locus. We can vary the spatial extent from σ = 1 for a uniform disturbance across space, to σ = s for a localized disturbance.

To measure recovery we use the return time T defined as the time needed for the ecosystem to recover 90% of the biomass lost to the disturbance. While the choice of a threshold is arbitrary, its specific value has no significant effect on the results as long as it is not too close to either 0% or 100% (which roughly correspond to reactivity and asymptotic resilience, respectively Arnoldi et al., 2016). By avoiding these extreme values, we simply emphasize the role played by the overall recovery dynamics, rather than by the system's initial response or final convergence.

### 2.3. The Temporal Dimension of Disturbances

We consider a disturbance regime by repeatedly applying disturbances with a given average frequency f , over a time period τ . For simplicity we assume no correlation in space or in time, so that the time between disturbances is drawn from an exponential distribution with some average frequency (a Poisson process, see **Appendix B** for details), while the location of the disturbance's center is drawn from a uniform distribution.

We use two measures of stability for a system that is disturbed repeatedly, i.e., variability, which measures how far the system ventures from its average value, and persistence, which measures how likely it is to move to the basin of attraction of a different equilibrium. We define variability V as the variance in time of the total biomass of the system, given a regime of repeated disturbances. In order to neglect the effect of transients, we calculate V over the last 80% of the simulations, which last for 10, 000 time steps. We define the collapse probability C as the probability that the system will be in the bare state at the end of a simulation, such that C = 0 means no chance of a system collapse, while C = 1 means that a collapse is certain. We use a longer simulation time (100, 000 time steps) to calculate C since we are interested in predicting a collapse before it occurs. For each of these calculations we run 100 simulations with different randomizations of the location and time of disturbances.

### 3. RESULTS

### 3.1. Spanning the Spatial Dimension of Disturbances Reveals Local Ecosystem Dynamics

We begin by looking at the response of an ecosystem to a single disturbance with varying spatial extent σ. We focus on disturbances with a fixed overall strength s = s0 for simplicity and clarity, and relax this assumption in the discussion. Thus a global disturbance σ = 1 (**Figure 1**, right panels) occurs when N is decreased by s0K in the entire system, while a localized disturbance σ = s0 (**Figure 1**, left panels) occurs when N is set to zero in a domain of relative size s0.

The response to a disturbance can take two possible forms: isolated recovery due to local processes, and rescue recovery due to incoming biomass from outside the disturbed region. Isolated recovery dominates the system response when each site recovers without the aid of neighboring sites (**Figure 1**, right panels). In contrast, rescue recovery occurs when the disturbed region cannot recover without the rest of the system, or when the bulk of the recovery occurs due to such spatial dynamics (**Figure 1**, left panels).

The coupling of local dynamics and dispersal results in distinct recovery processes in the three models, as shown by

the trajectories in phase-space diagrams in the middle column of **Figure 2**. In these panels we unfold the recovery along two dimensions: the horizontal axis denotes the size of the disturbed region at a given time, while the vertical axis shows the biomass density in the disturbed region. Immediately after the disturbance, the system is along the dashed black curve, and it then changes over time until it enters the shaded region where it is considered to have recovered.

If a large part of the trajectory during recovery is horizontal, this means that the disturbed region is shrinking due to rescue recovery, which indicates a lack of local resilience, which would otherwise allow isolated recovery to take place. This behavior reflects the strong nonlinearity of local dynamics, which can be seen in the changing curvature of the local potential (**Figure 2**, left column). We can see that for the AE model (**Figure 2**, bottom row) the recovery is along a horizontal line for recovery scenarios with a sufficiently small spatial extent, so that regional processes bring about the recovery. In contrast, the recovery is entirely due to local processes in the LG model since the local dynamics are much faster here, while for the intermediate SR model a mixture of the two processes can be seen to take place.

These differences translate into markedly different values of the return time T (**Figure 2**, right panels). The vertical recovery trajectories that follow all disturbances in the LG model and large-sized disturbances in other models indicate isolated recovery, and hence small values of T. For the intermediate SR model localized disturbances lead to a larger contribution of rescue recovery, leading to a sigmoid shape of T as a function of disturbance extent σ. The AE model shows a similar behavior of larger T following localized disturbances, but the trend here shows a maximum for mid-sized disturbances. This occurs because in bistable systems, the most efficient way to perturb the system is to locally remove biomass just bellow the viability threshold, and then let the system collapse locally. Such a disturbance has an equivalent effect to that of a stronger disturbance that would remove all biomass over a larger region. The spatial recovery process will take longer to recover, thus giving larger return time values (see **Appendix A** for details).

FIGURE 2 | Contribution of isolated and rescue recovery as a function of disturbance spatial extent for the three models presented in the main text. The left column shows the local potentials defining local processes. Top row: Logistic (LG) model; Middle row: the highly non-linear SR model; Bottom row: bistable AE model. Middle column: isolated recovery on the y-axis, and rescue recovery on the x-axis. Black dashed line shows the equal disturbance strength used s = 0.2 for different disturbance extent σ. Blue lines are recovery trajectories, where recovery is considered complete when trajectories reach the gray shaded region. For the SR and AE model, as disturbances become more localized, a shift is observed from a dominant isolated recovery (upward trajectories) to a dominant rescue recovery (leftward trajectories), impacting return times (right column). The "x" marks in blue correspond to the different trajectories shown in middle columns. The green and magenta circles show initial states following two different disturbances (left and middle columns) and their associated return times (right column). The dotted line (bottom row) shows the local tipping point of the bistable AE model, beyond which local dynamics collapse to the bare state.

This explains the humped shape of return time as a function of disturbance extent. Since bistability is a sufficient condition for a hump-shaped relationship to occur, the latter could be used as an indicator of bistability. This illustrates the more general idea that considering the spatial dimension of disturbances can allow us to probe the local dynamics of a spatially extended ecosystem.

### 3.2. From Variability to Collapse Under Increasing Frequency of Disturbances

Natural ecosystems are constantly perturbed, leading us to consider a temporal dimension of disturbances, namely their average frequency. We therefore translate the results of the previous section on the response to a single disturbance (**Figure 3**, top) into an understanding of temporal variability under repeated disturbances (**Figure 3**, bottom). In fact, there is a direct link between the response to a single disturbance and temporal variability in response to repeated disturbances. Indeed, biomass fluctuations are the result of past disturbances, as they integrate short- to long-term responses of the ecosystem to individual disturbances (Arnoldi et al., 2016). Variability is a statistic of those fluctuations, and is therefore a function of both the integrated response to a single disturbance and the average frequency f of disturbances. More precisely, if g(t) traces the change in overall biomass through time following a pulse disturbance a time t = 0, then variability V can be expressed as V = f R ∞ 0 g 2 (t)dt (see Equation S15 in **Appendix B**). However, this identity assumes no interaction in space between the different disturbances, and therefore should not hold at high disturbance frequency.

As expected, at low frequency of disturbances the analytical approximation agrees with the numerical simulations quite well for all three models (**Figure 4**, second column). For higher frequencies (**Figure 4**, third column) where multiple disturbances often take place in the same time frame, we see a slight underestimation of the analytical approximation, although the general trend is well captured. Importantly, variability and return time show the same behavior. We see effects of regional processes on variability for more localized disturbances in both the SR and AE models, where the former shows a sigmoid shape while the latter has a hump shape, which is a consequence of the bistability in the AE model. We note that these trends hold in more general scenarios, such as disturbances with a random extent or following seasonal patterns (**Appendix D**).

At this point it would appear that the temporal dimension of disturbances f is not as informative on ecosystem behavior as the spatial dimension of disturbances σ. However, as f is increased further, a discrepancy between variability and its prediction based on recovery from a single disturbance starts to grow. This signals that the disturbances start to interact with each other, a phenomenon that is not captured by our approximation. Disturbances start to aggregate in space, which can substantially increase variability (**Appendix C**) due to large excursions toward low total biomass levels. For bistable systems such as the AE model, such excursions can lead to a collapse of the whole system. This is evident in the two last columns of **Figure 4**, in which we see, for the AE model, that the values of σ for which the

FIGURE 3 | Single and multiple disturbance regimes and the relationship between return time and variability. The left panels show time series of the overall biomass, while right panels are spatial snapshots of the corresponding time-series (red dots in the left panels). The response to a single disturbance is shown in the top left panel. We focus on two of its characteristics: return time T, and an integral measure of the transient g(t) (see main text). The response to multiple disturbances occurring randomly at an average frequency f is shown in the bottom left panel. It is summarized by its variability V (variance of overall biomass). In the limit of low f there is an inherent relationship between return time and variability in the sense that V can be approximated by f R ∞ 0 g 2 (t)dt. Simulations were made using the SR model with parameters values: s = 0.1, σ = 0.11, f = 0.025, and γ = 2. Random uniform noise was added in left panels to demonstrate how realistic time series might look like.

discrepancy of variability is highest precisely corresponds to the values of σ for which the collapse probability is most significant. Thus, at high frequencies, disturbances of similar strength but different spatial extent lead to dramatically different responses. This example highlights the fact that the combination of spatial and temporal dimensions of disturbances can have a drastic effect on ecosystem stability.

### 4. DISCUSSION

Investigating the role of the spatial and temporal dimensions of disturbances in ecosystem stability, we obtained four main results: (1) In comparison with a global disturbance, a localized one of the same strength can initiate a fundamentally different, and much slower, ecosystem response, especially when local dynamics are nonlinear. (2) The return time from a single disturbance and the temporal variability caused by repeated disturbances show the same trends, even for locally intense (and therefore nonlinear) disturbances. (3) The relationship between a system's response and the spatial extent of the disturbances it experiences reveals its underlying dynamics. For instance, a hump-shaped relationship between return time and the spatial extent of the disturbances may indicate bistability. (4) The correspondence between return time and variability breaks down for high disturbance frequencies. This discrepancy signals the

occurrence of spatial interactions between disturbed regions, which, in turn, may lead to a regime shift.

Although we considered simple spatially homogenous models, our results should apply to a wide range of ecosystems. Forests, savannah and shrublands might be good examples of ecosystems to which our models apply since disturbances such as fires and grazing occur frequently and are often localized, and the recovery of plant communities often follows complex succession dynamics driven by spatial processes (Adler et al., 2001; Turner, 2010; Staver and Levin, 2012). Our results, however, need not be restricted to such spatially homogeneous systems. Although we built our theory using spatially uniform models, this simplifying feature is not essential to our arguments, which only require a notion of locality. Therefore, our theory may also be relevant to less homogeneous ecosystems, such as mountain lake networks, coral reefs and riverine systems. Indeed, such ecosystems undergo different disturbances that are often strongly localized, and their dynamics may be sufficiently nonlinear (Knowlton, 1992; Campbell Grant et al., 2007; Forrest and Arnott, 2007).

Uniquely to our work, we considered systems locally pushed far from their equilibrium, and even to a different basin of attraction. In a marine ecosystem context, this could represent coral reefs (Nyström et al., 2000; Adjeroud et al., 2009) or rocky intertidal systems (Sousa, 1979; Paine and Levin, 1981), which frequently undergo intense disturbances (e.g., storm damage). These locally intense disturbances can allow rescue recovery, mediated by dispersal, to dominate the ecosystem response. In the case of the bistable (AE) model this glimpse outside the basin of attraction of the populated state is the direct cause of the hump-shaped trends of variability and return time as a function of disturbance extent. In fact, the front propagation that drives rescue recovery contains information about the ecosystem's basins of attractions, reflecting the existence of alternative stable states and the transient dynamics between them. Thus, by observing the ecosystem's response to localized disturbances, rescue recovery allows us to probe ecosystem dynamics far from equilibrium. For instance, comparing between different disturbed marine ecosystems may give further evidence that some have alternative states (e.g., coral reefs) while for others the dynamics show a succession process (e.g., rocky intertidal systems). This reasoning could be taken further by focusing on regions where rescue recovery takes place, e.g., analyzing the plant community structure at transition zones between grassland and forest in a savanna ecosystem (Augustine, 2003).

Right column shows the corresponding reconstruction of the return time curve, using 100 randomly chosen points (red asterisks) in the parameter space. The return time values are normalized by the disturbance strength s, while the normalized disturbance extent is defined as σ˜ = 1 − s/σ. Note that the hump (sigmoid) shape of the curve for the AE (SR) model are easily recognizable from these reconstructions.

Spanning the spatial dimension of disturbances could thus allow us to detect nonlinearities in ecosystem behavior, revealed by the increasing local intensity of disturbances (see **Figure 2**). One might expect that along the temporal dimension of disturbances, increasing their average frequency could also reveal nonlinear effects, since the ecosystem becomes more strongly disturbed. In fact, increasing frequency has only a trivial linear effect, as reflected by the relation we found between return time and variability (see **Figure 3**). Beyond some threshold, however, a response of a different kind emerges, due to spatial interactions between disturbed regions which aggregate in potentially largescale clusters. This causes a higher variability than expected and can, consequently, cause a global loss of persistence or a regime shift. Taking, once again, the example of corals reefs, we could ask how the impact of both natural and anthropogenic disturbances leads to a phase-shift from hard coral to fleshy algae dominance. A regime shift due to an aggregation of unrecovered regions would occur not as a typical tipping point due to loss of resilience (e.g., due to changing temperatures), but rather due to the crossing of a threshold for disturbance frequency. Importantly, in such a scenario the two dimensions, spatial and temporal, must be considered in conjunction. The threshold beyond which aggregation occurs depends strongly on the spatial extent of disturbances and hence the associated response is not a mere superposition of responses to single disturbances. In other words, this finding highlights and explains how the interplay between the spatial and temporal dimensions of disturbances can have drastic ecological consequences, such the loss of persistence. Since our findings are purely theoretical, it would be enlightening to elucidate the prevalence of this interplay in empirical systems that have undergone regime shifts (e.g., phase-shifts in coral reefs Nyström et al., 2000 or the desertification of the once green Sahara Ortiz et al., 2000).

As previously mentioned, in bistable systems the relationship between return time (as well as variability) and the spatial extent of disturbances is hump-shaped. This relation could be used as an indicator of bistability, assessed empirically by comparing time series of the same ecosystem in different regions with estimates of the intensity of single disturbances. Its implications for ecosystem management depend on the type of disturbances considered. Anthropogenic disturbances that are largely controlled, such as logging in forests (Chazdon, 2003) or large-scale fishing (Kaiser et al., 2006), can be better planned to avoid both an unpredictable yield due to high variability and an overall collapse. For many natural disturbances control is neither possible nor desired (e.g., fires in semi-arid ecosystems necessary for plant germination Wellington and Noble, 1985), but predicting their effects and the possibility of regime shifts is paramount (Kéfi et al., 2007).

In order to focus on the role of the spatial properties of disturbances and allow a clearer presentation, we conducted our analysis assuming disturbances of constant overall strength. It is straightforward to extend the analysis to more general settings, such as a random extent of disturbances and seasonal patterns (see **Appendix D** for details). It is particularly interesting to consider the case of different values of disturbance strength s. As shown in **Figure 5**, if we randomly choose a set of points with different values of strength s and extent σ, we can use these to reconstruct a normalized version of the dependency of the different stability measures on disturbance extent. Thus we can use the different phenomena described previously, such as a hump-shape relationship as an indicator of bistability, under more general conditions, thereby making our theory more empirically accessible.

Our work is a step toward a quantitative account of spatial and temporal dimensions of disturbances, and their interplay with local and regional ecosystem dynamics. This is an important goal in the context of global change. Disturbances are of increasing frequencies and occur at different scales (which is evident, e.g., in coral reefsJackson, 1991 and forests Turner et al., 1993), while the spatial structure of ecosystems themselves is altered by land use change, often causing fragmentation of the landscape (Harrison and Bruna, 1999). It is thus important to build a framework in which we can understand and predict the ecological impacts of this complex interplay.

### REFERENCES

### DATA AVAILABILITY STATEMENT

The results shown in this study are based on Matlab code that can be found in the Github repository of https://github.com/yzelnik/ stddes-scripts.

### AUTHOR CONTRIBUTIONS

YZ, J-FA, and ML designed the study. YZ and J-FA performed the research. YZ, J-FA,and ML wrote the manuscript.

### FUNDING

This work was supported by the TULIP Laboratory of Excellence (ANR-10-LABX-41) and by the BIOSTASES Advanced Grant, funded by the European Research Council under the European Union's Horizon 2020 research and innovation programme (grant agreement No 666971).

### ACKNOWLEDGMENTS

We wish to thank Bart Haegeman and Matthieu Barbier for helpful discussions and comments on a previous version of the manuscript. We would also like to thank two referees for their useful comments on the manuscript.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo. 2018.00224/full#supplementary-material

Jackson, J. B. (1991). Adaptation and diversity of reef corals. BioScience 475–482.

Holling, C. S. (1973). Resilience and stability of ecological systems. Annu. Rev. Ecol. Syst. 4, 1–23. doi: 10.1146/annurev.es.04.110173.000245

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Zelnik, Arnoldi and Loreau. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Sensitivity of Subsurface Microbes to Ocean Warming Accentuates Future Declines in Particulate Carbon Export

Emma Louise Cavan1 \*, Stephanie A. Henson2 and Philip W. Boyd1,3

1 Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, TAS, Australia, 2 National Oceanography Centre, Southampton, United Kingdom, 3 Antarctic Climate and Ecosystems CRC, University of Tasmania, Hobart, TAS, Australia

Under future warming Earth System Models (ESMs) project a decrease in the magnitude of downward particulate organic carbon (POC) export, suggesting the potential for carbon storage in the deep ocean will be reduced. Projections of POC export can also be quantified using an alternative physiologically-based approach, the Metabolic Theory of Ecology (MTE). MTE employs an activation energy (Ea) describing organismal metabolic sensitivity to temperature change, but does not consider changes in ocean chemistry or physics. Many ESMs incorporate temperature dependent functions, where rates (e.g., respiration) scale with temperature. Temperature sensitivity describes how temperature dependence varies across metabolic rates or species. ESMs acknowledge temperature sensitivity between rates (e.g., between heterotrophic and autotropic processes), but due to a lack of empirical data cannot parameterize for variation within rates, such as differences within species or biogeochemical provinces. Here we investigate how varying temperature sensitivity affects heterotrophic microbial respiration and hence future POC export. Using satellite-derived data and ESM temperature projections we applied microbial MTE, with varying temperature sensitivity, to estimates of global POC export. In line with observations from polar regions and the deep ocean we imposed an elevated temperature sensitivity (Ea = 1.0 eV) to cooler regions; firstly to the Southern Ocean (south of 40◦S) and secondly where temperature at 100 m depth <13◦C. Elsewhere in both these scenarios Ea was set to 0.7 eV (moderate sensitivity/classic MTE). Imposing high temperature sensitivity in cool regions resulted in projected declines in export of 17 ± 1% (< 40◦S) and 23 ± 1% (< 13◦C) by 2100 relative to the present day. Hence varying microbial temperature sensitivity resulted in at least 2-fold greater declines in POC export than suggested by classic MTE derived in this study (12 ± 1%, Ea = 0.7 eV globally) or ESMs (1–12%). The sparse observational data currently available suggests metabolic temperature sensitivity of organisms likely differs depending on the oceanic province they reside in. We advocate temperature sensitivity to be incorporated in biogeochemical models to improve projections of future carbon export, which could be currently underestimating the change in future POC export.

Keywords: particulate carbon export, microbes, warming, metabolic theory, activation energy

#### Edited by:

Diego Barneche, University of Sydney, Australia

#### Reviewed by:

C.-Elisa Schaum, Universität Hamburg, Germany Francisca C. García, University of Exeter, United Kingdom

> \*Correspondence: Emma Louise Cavan<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeography and Macroecology, a section of the journal Frontiers in Ecology and Evolution

Received: 18 October 2018 Accepted: 07 December 2018 Published: 07 January 2019

#### Citation:

Cavan EL, Henson SA and Boyd PW (2019) The Sensitivity of Subsurface Microbes to Ocean Warming Accentuates Future Declines in Particulate Carbon Export. Front. Ecol. Evol. 6:230. doi: 10.3389/fevo.2018.00230

### INTRODUCTION

The biological pump exports large amounts of carbon from the surface ocean to the deep, where it can be stored on climatically-relevant timescales helping to regulate atmospheric carbon dioxide levels (Volk and Hoffert, 1985; Falkowski et al., 1998). However future ocean warming threatens to decrease the amount of particulate organic carbon (POC) reaching the deep ocean (Laufkötter et al., 2016). The magnitude of deep ocean carbon storage by the biological pump is largely dependent on 3 factors; (1) the magnitude and size-partitioning of primary production taking place in the surface mixed layer, (2) the sinking rate of the particles formed initially from particle production, and (3) the organic remineralisation rate by organisms such as zooplankton and microbes that degrade POC (Buesseler and Boyd, 2009; Turner, 2015). Each of these parameters is a function of temperature, with warming increasing all three rates if other factors, such as nutrient supply, remain unchanged (López-Urrutia et al., 2006; Taucher and Oschlies, 2011; Iversen and Ploug, 2013).

Warming is projected to increase metabolic rates (Brown et al., 2004), which will increase both primary production (more carbon available to sink to deep ocean) and the metabolic rates of heterotrophs that consume phytoplankton and sinking POC (decreasing the carbon sink) (López-Urrutia et al., 2006; Taucher and Oschlies, 2011; Cavan and Boyd, 2018). In parallel warming-induced stratification and shoaling of the mixed layer will reduce nutrient inputs to the sunlit upper ocean (Bopp et al., 2001), expanding oligotrophic waters thus likely favoring smaller phytoplankton species and ultimately reducing primary production (Bopp et al., 2005). The decrease in nutrient supply to surface waters is typically considered to have the greatest effect on the magnitude of future primary production (Marañón et al., 2014) and thus overall primary production is expected to decline, although there is some disagreement amongst models (Laufkötter et al., 2015). One feedback from warming that will increase the carbon sink is the reduction in water viscosity, allowing particles to sink through the water column faster, escaping the upper ocean where remineralisation is most intense (Bach et al., 2012).

Export production is projected to decline by 1–12% depending on the Earth System Model (ESM) model used, due to declines in the magnitude and changes in the size-partitioning of primary production, and increased remineralisation of POC due to warming (Laufkötter et al., 2016). Remineralisation occurs throughout the water column and affects both the amount of POC exported and POC attenuation through the mesopelagic zone. Export or formation of detritus can either be parameterised by simple empirical algorithms, which may include primary production and temperature terms (Dunne et al., 2005; Henson et al., 2011; Britten et al., 2017) or in ESMs by phytoplankton aggregation and remineralisation parameterisations (Aumont et al., 2015).

Currently, there are a range of approaches to parameterise temperature-dependent terms in models. For example, in the biogeochemical model REcoM2, one theoretical exponential (Arrhenius type) relationship is used to describe all temperature-dependent terms; remineralisation of particulate and dissolved organic matter, silicon uptake, zooplankton respiration, and zooplankton grazing (Schourup-Kristensen et al., 2014). However, in the PISCES-v2 model, a constant Q10 of 1.9 from Eppley (1972), which correlated phytoplankton growth rates with temperature, is used to describe both phytoplankton growth rates and POC degradation (Aumont et al., 2015). A Q10 of 1.9 implies an increase in 10◦C will result in metabolic rates that are 1.9 times higher. None of the marine biogeochemical components of the ESMs account for adaptation, which may be an important response to ocean warming with knock-on effects for POC export. For example, a laboratory study has shown that after 100 generations (freshwater) phytoplankton can adapt to warming by down-regulating respiration relative to photosynthesis to maintain the carbon allocation efficiency needed for growth (Padfield et al., 2016). Although phytoplankton in nutrient-poor regions may exhibit a different adaptation response as nutrient limitation can supress the temperature dependence of phytoplankton (Marañón et al., 2014).

In some ESMs temperature dependence varies between different metabolic rates (e.g., heterotrophic processes have a higher rate of change to temperature than autotrophic processes) (López-Urrutia et al., 2006; Dunne, 2013) thus incorporating differences in temperature sensitivity are only captured at a coarse scale. At present these biogeochemical models do not acknowledge temperature sensitivity within rates i.e., differences in the response to temperature between species or between the same species living in a different biogeochemical province. However, there is unlikely to be one uniform metabolic response to warming globally.

Sensitivity to temperature can be quantified as the activation energy (Ea) of a metabolic reaction—the amount of energy needed for a chemical reaction to occur (Schoolfield et al., 1981; Clarke and Johnston, 1999; Yvon-Durocher et al., 2012). Ea can be estimated via the Metabolic Theory of Ecology (MTE), which describes how metabolic processes vary as a function of organismal mass and temperature (Brown et al., 2004):

$$
\ln\left(I\,M^{-0.75}\right) = \,-E\_a\left(\frac{1}{cT}\right) + \ln(i\_0) \tag{1}
$$

Where I is the metabolic rate normalized by mass (M), Ea is the activation energy, c is Boltzmann's constant (8.62 × 10−5 eV K−1 ), T is the temperature in Kelvin and finally i0 is a normalization constant. The latter can be used to describe how organisms or communities are adapted to living at different temperatures (Clarke, 2006). When computed in the context of MTE using Equation 1, the Ea is the slope between mass-normalized metabolism and temperature. MTE states all organisms have an activation energy of 0.6–0.7 eV (Gillooly et al., 2001; Brown et al., 2004), thus there is little natural variation in temperature sensitivity. Activation energies higher than this range suggest organisms are more sensitive to changes in temperature, as shown in Arctic, and Antarctic zooplankton [Ea = 1.29 eV, (Alcaraz, 2016)], mesopelagic heterotrophs [Ea = 0.9 eV, (Brewer and Peltzer, 2016)] and a mixed heterotrophic microbial community from the sub-Antarctic [Ea = 0.9 eV, (Cavan and Boyd, 2018)]. As primary production and respiration (remineralisation) are both metabolic rates, MTE can be used to describe how they might change with future warming in the global oceans, and thus their influence on POC export. Model experiments have confirmed that temperature sensitivities of metabolic rates have an important role in ecosystem function, by regulating the magnitude of primary production and respiration (Taucher and Oschlies, 2011).

If organisms do respond differently to warming based on the environmental conditions they currently encounter, then biogeochemical models are likely to be inaccurately projecting carbon sequestration under future warming scenarios. Therefore, the objective of this study is to use empirical algorithms, including MTE, satellite data, and ESM outputs to determine the effect of temperature sensitivity on future POC export projections by 2100. We varied the activation energy in line with observations to investigate how temperature sensitivity influences POC export. MTE only accounts for changes in metabolic rates due to temperature, and thus using this framework implies that nutrients, phytoplankton community composition, and ocean physical changes, such as stratification, are unaltered in the year 2100.

### METHODS

### Data

Annual satellite sea surface temperature (SST) and output from the Vertically Generalized Production Model (VGPM) for primary production for 2003–2016 (Oregon State University, 2017) were used as the baseline data for the beginning of this century. Coupled-Model Intercomparison Project-Phase 5 (CMIP5) ESM outputs over 1986–2005 (historical run) and under the IPCC Representative Concentration Pathway (RCP) 8.5 (business-as-usual run) for 2081–2100 were used to estimate the change in temperature at the sea surface and at 100 m depth and the change in export and primary production by the end of this century. See **Table 1** for a list of frequently used acronyms. The 8 models used were GFDL-ESM2G and GFDL-ESM2M (Dunne et al., 2013), HadGEM2-CC and HadGEM2-ES (Bellouin et al., 2011; Collins et al., 2011), IPSL-LR and IPSL-MR (Séférian et al., 2013), MPI-LR and MPI-MR (Jungclaus et al., 2013) (see **Supplementary Table 1** and Bopp et al. (2013) for a summary of model set-ups). These models were selected as they were the ones available in the CMIP5 (at the time of our analysis) that simulated export, primary production, and temperature.

### Satellite-Derived Climatologies