The Minimum Environmental Perturbation Principle: A New Perspective on Niche Theory

Contemporary niche theory is a powerful conceptual framework for understanding how organisms interact with each other and with their shared environment. Here we show that a large segment of niche theory is equivalent to a Minimum Environmental Perturbation Principle (MEPP): ecosystems self-organize into a state that minimizes the impact of organisms on their environment. Different choices of environmental dynamics naturally give rise to distinct dissimilarity measures for quantifying environmental impact. The MEPP allows for the analysis of ecosystems with large numbers of species and environmental factors and provides a new avenue for analyzing ecological invasions. The MEPP also rigorously connects ecological bistability with the existence of multiple minima in a statistical-physics inspired landscapes. We show that the presence of environmental feedbacks where organisms can produce new resources in addition to depleting them violates the global MEPP. However, even in the presence of such feedbacks, a weaker, local version of the MEPP still applies in a limited region of resource space.

The concept of a "niche" has long been central to ecological theory. Over the past fifty years, this concept has been significantly clarified and refined, beginning with the pioneering work of Levins, MacArthur and Tilman, and more recently consolidated by Chase and Leibold [1][2][3]. As illustrated in Figure 1(a), this framework highlights two distinct aspects of an organism's niche: its requirements for survival and its impact on the environment. The state of the environment is mathematically represented by a vector R in an abstract space, whose components R α (α = 1, 2 . . . M ) correspond to concentrations of resources, populations of predators, stress intensity, or any other factors that affect an organism's growth rate [4]. In this space, the "requirement niche" of each species i (i = 1, 2 . . . S) can be encoded by a zeronet-growth isocline (ZNGI). This line or surface separates the environmental states where the growth rate g i (R) is positive from the states where it is negative. The "impact niche" is represented by an impact vector q i (R), which specifies the magnitude and direction of the environmental change induced by a single individual of the species. The intrinsic dynamics of the environment are contained in the supply vector h(R), which could represent an external supply of abiotic nutrients, regrowth of biotic resources like grasses or algae, natural death of predators, or any other environmental process that is independent of the organisms whose niches are being analyzed.
For simple environments with just two relevant environmental factors, one can graphically investigate all the possible outcomes of the ecological dynamics by drawing the ZNGI's, impact vectors and supply vectors. As the dimension of R grows, the space becomes harder to visualize, and the number of possibilities grows exponentially.
The high-dimensional situation is especially complicated by the fact that each niche has two aspects that can be varied independently. In many situations, however, it is natural to assume that the requirement and impact niches of a given organism are related, so that knowledge of one immediately determines the other. In this paper, we show that this assumption leads to a powerful new principle for finding the ecologically stable states (ESS) of ecosystems at any level of environmental complexity. As show in Table I, our derivation is applicable to a wide variety of ecological models including those with biotic and abiotic resource dynamics, substitutable or essential resources, Type I and Type II response functions, environmental feedbacks, and predation and/or resource limited ecosystems.

RESULTS
The ZNGI's, supply vectors and impact vectors all play a role in the construction of the optimization principle, as illustrated in Figure 1(b). For a given regional species pool, the ZNGI's fix the boundaries of the 'ESS region' Ω, which is the region of resource space where ecologically stable states are possible. Outside of this region, the growth rate of at least one species from the regional pool is positive, guaranteeing that the system is susceptible to invasions. Only when all the growth rates g i are zero or negative is it possible for the populations sizes and resource abundances to stably remain constant, with the negative growth rate species going extinct.
Any point in the ESS region can be made into the actual steady state R * of the ecological dynamics with a suitable choice of the supply vector field h(R). In the absence of any consumer species, the environment will relax to its own steady state R 0 where the supply vector h(R 0 ) vanishes. If this point happens to lie in the interior of Ω, then all species in the regional pool will go extinct, and the environment always returns to its unper-  have zero net growth rates, all extinct species have negative net growth rates, and the combined impact of all consumers on the environment perfectly balances its intrinsic dynamics. Contemporary niche theory provides graphical and mathematical statements of these criteria in terms of the impact vectors and zero net-growth isoclines (ZNGI's) of each consumer species, and a supply vector that encodes the intrinsic dynamics of the environment. Two standard choices for the supply vector are logistic growth of self-renewing resources, and a linear model of externally supplied resource fluxes. (b) The Minimum Environmental Perturbation Principle (MEPP) provides an alternative procedure for finding an ESS, by minimizing a dissimilarity measure d(R 0 , R) of the size of the environmental perturbation caused by the consumers. The minimization is performed over resource concentration vectors R within the shaded "ESS region" Ω where all the per-capita growth rates gi(R) are negative or zero. Different supply vector fields generate different measures d(R 0 , R), with self-renewing and externally supplied resource dynamics giving rise to the Euclidean distance and Kullback-Leibler divergence respectively. (c) MEPP readily generalizes to larger numbers of resources. turbed state. If R 0 is outside the ESS region, however, the ESS will lie somewhere on the boundary, and at least one consumer species will survive. The surviving species perturb the environmental state, and we can quantify the magnitude of this perturbation in terms of Euclidean distance, Kullback-Leibler divergence, or any number of other dissimilarity measures d(R 0 , R).
To determine which point on the boundary is chosen by the ecological dynamics as the true steady state, we also need to know the impact vectors q i of the surviving species. In many applications of niche theory, one is primarily interested in the conditions for coexistence of a set of similar species, such as the different species of warblers in MacArthur's original field work [5], or the paramecia of Gause's competitive exclusion experiments [6], with all other organisms treated as components of the environment. The functional similarity of the competitors justifies the assumption of a universal relationship between the ZNGI's and the impact vectors, such that for any new species, the direction of the impact vector can be inferred from the ZNGI. This assumption significantly simplifies the analysis even of two-resource ecosystems, and is present in much of the early work on consumer-resource models [1,2,7]. Later studies showed that strong violations of this postulate lead to limit cycles, heteroclinic cycles and chaos when the number of resources is sufficiently large [8][9][10]. The recent microbiome-related resurgence of interest in consumer-resource models has almost exclusively focused on cases where a strict ZNGI/impact relationship is maintained [11][12][13][14][15].
Our main result is a mathematical proof that this same assumption converts niche theory into an optimization problem. With the impact vectors no longer constituting an independent set of parameters, the ESS can be found from the ZNGI's and the supply vectors alone. The Minimum Environmental Perturbation Principle (MEPP) then says that any ESS locally minimizes the magnitude of the environmental perturbation within the ESS region. The intrinsic environmental dynamics encoded in the supply vector h determine the form of the dissimilarity measure d(R 0 , R). As summarized in Figure 1(b), MacArthur's original assumption of logistic growth for the resources makes the Euclidean distance the appropriate choice for d(R 0 , R), while a linear model of chemostat-like dynamics leads to the Kullback-Leibler divergence [16].  Figure 2 illustrates the MEPP for both of these paradigmatic cases, with contour lines of the relevant perturbation measure d(R 0 , R). In these two-dimensional examples, the ESS can be visually identified as the point lying on the contour line closest to R 0 . Numerical simulations of the ecological dynamics with different initial conditions always end up at this point, even if the transient behavior is complicated. Since R 0 is outside the ESS region in both of these examples, the ESS must lie on the boundary of Ω, so we can also visualize the minimization by plotting the variations of d(R 0 , R) along the one-dimensional boundary. This curve (or surface in higher dimension) is analogous to the free energy landscapes studied in statistical physics, and familiar to biologists from the study of protein folding [17].
In these examples, the environmental perturbation has a unique local minimum along the ESS boundary, guaranteeing a unique ESS. If there are multiple local minima, however, the MEPP implies that each of these minima is a possible ESS. Figure 3(a) shows the ESS region of a classic bistable ecosystem of two species competing for two essential resources, while panel (c) shows how d(R 0 , R) varies along the boundary. The latter plot reveals two local minima, one at the intersection of the two ZNGI's, where the two species coexist, and one at another point, where only one of the species survives. Numerical simulations in panel (a) with two different sets of initial population sizes confirm that each of the two minima can be reached under a suitable choice of initial conditions. So far we have only considered cases of competition between pure consumers. But real ecosystems often exhibit more complex environmental feedbacks, sometimes referred to as 'bioengineering' interactions [3]: large mammals can destroy plant species by trampling or wallowing, beavers can create entirely new habitats for aquatic species, fish excrement can influence algal growth. In general, these effects break the required relationship between growth rates and impact vectors [3]. In some important cases, however, an extended version of the principle still applies. One concrete example is cross-feeding among microbial species. Microbes generically release metabolic byproducts into their environment, serving as both consumers and producers of resources [18,19]. Since core carbon metabolism is largely conserved across many microbial species, these metabolic conversions can plausibly be approximately encoded into a universal stoichiometry matrix [20,21].
Under this "universal chemistry" assumption, as illustrated in Figure 3(b) and (c), byproduct secretion shifts the unperturbed environmental state from R 0 (where the supply vector vanishes) to an "effective" locationR 0 (see Methods for details). ComputingR 0 requires prior knowledge of the steady-state resource concentrations R * , and so the steady state cannot be found by straightforward optimization algorithm in this case. However, one can construct an iterative algorithm that starts with an arbitrary guess R * 0 for the steady-state environment vector, as sketched in Figure 3(f ). The loop terminates when the self-consistency condition for the steady state is achieved, that is, when R * t minimizes d(R 0 t , R). This local version of the MEPP principle is intimately related to Expectation-Maximization algorithms used in statistical inference and machine learning [16,22].

DISCUSSION
The Minimum Environmental Perturbation Principle provides a new perspective on community ecology, by recasting the coexistence conditions of niche theory as the solution to an optimization problem. The minimized quantity has a natural physical interpretation as the size of the environmental perturbation induced by the organisms. Different models for the intrinsic environmental dynamics give rise to different perturbation measures: selfrenewing resources such as algae or plants minimize the Euclidean distance, while externally supplied resources like sugars in a microbial culture minimize the Kullback-Leibler divergence.
The MEPP provides a powerful new way of extending the graphical methods of ecosystem analysis popularized by Tilman to more complex situations. The 'R * ' criterion remains a necessary one for species coexistence [3,23]. Independently of any other assumptions, coexistence requires the ZNGI's to intersect, which means each species can survive at a lower concentration of a least one resource than any of its competitors. But the criterion for resource supply in the graphical coexistence analysis requires that the supply vector field is linear (so that the supply vectors point towards the supply point), and is easiest to implement if the impact vectors are also fixed. These assumptions are valid for self-renewing resources, but are violated in even the simplest models of externally supplied resources, as illustrated in Supplementary Figure S2. MEPP makes it possible to perform a geometric analysis of a much broader class of systems, including all the examples discussed above, and many others that are listed in Table 1. Furthermore, as Chase and Leibold point out, the relationships among the impact vectors are increasingly difficult to visualize in higher dimensions [3]. But high-dimensional constrained optimization is a long-established tool in many engineering disciplines, which have built up a deep pool of theorems, algorithms and intuition that can now be appropriated for ecological reasoning [24].
The validity of MEPP requires that impact vectors and ZNGI's be related in a species-independent way. The correlation between these two niche components has long been known as a condition for ecosystem stability [3], and recent theoretical work has demonstrated that perfect correlation is indeed sufficient for stability in MacArthur's Consumer Resource Model [11]. The fact that this assumption also leads to an optimization principle is further evidence of its importance for understanding ecological dynamics. Non-trophic interactions mediated by natural 'bioengineers' can easily break this postulate, underlining the ecological significance of environmental feedbacks from leaf litter to excrement decomposition [3,[25][26][27].
If different resources are supplied at different rates or have different nutritional value, these factors re-weight the resource axes as described in the Methods. The natural emergence of these dissimilarity measures as optimized quantities suggests that they may be particularly informative ways of quantifying environmental change in field and laboratory experiments. In particular, an immediate consequence of MEPP is that successful invasion by mutation or immigration necessarily increases the appropriate measure of environmental perturbation. This prediction can be directly tested, using only the observables accessible to K. Rotthaupt in his classic zooplankton experiments of the 1980's [28].

A. General Derivation
To derive the MEPP, we begin from the full set of conditions for an ESS to exist at environmental state R * : where the first condition must hold for all α = 1, . . . M and the last three conditions must hold for all i = 1, . . . , S. These equations are very similar in form to the Karush-Kuhn-Tucker (KKT) conditions that must be satisfied at any local minimum R * of a function f (R) subject to inequality constraints g i (R) ≤ 0 [29][30][31]. The KKT conditions are: where the λ i are generalized Lagrange multipliers (known in the literature as KKT multipliers) that enforce the inequality constraints. The KKT conditions have a straightforward and intuitive explanation. At the optimum R * , either g i (R * ) = 0 and the constraint is active λ i ≥ 0, or g i (R * ) ≤ 0 and the constraint is inactive λ i = 0.
The ESS conditions become identical to the KKT conditions when for some function f (R). We are free to add a constant to f (R) while still satisfying the conditions, and can therefore always make f (R 0 ) = 0 at its unconstrained minimum R 0 , which is the location of the unperturbed environmental state. Since f increases as we move away from the minimum, it can be used (at least locally, but usually globally) as a measure of the size of the perturbation away from R 0 . The identity of the ESS conditions with the KKT conditions then guarantees that the perturbation d(R 0 , R) ≡ f (R) is locally minimized in any ESS.
To apply MEPP to a given ecological model, one must transform the steady-state equation for the environment so that q i = −∇g i . The transformed supply vector h then gives rise to the natural dissimilarity measure d.
The central assumption stated in the main text, that the direction of q i can be inferred from the ZNGI, guarantees that this is possible. For two resources, this assumption reduces to the simple statement that the angle between the ZNGI and the impact vector is species-independent (but can depend on the location in resource space). More generally, it means that the dependence of the impact vectors on R can be written as The arbitrary positive-valued functions a i (R) change the magnitude of the impact vectors, but not the direction, while another set of arbitrary positive-valued functions b α (R) provide a universal, species-independent mapping from the normal vector ∇g i of the ZNGI to the impact vector. In terms of the original supply vector field h(R), the optimized function d(R 0 , R) obeys Such a function only exists if the vector h α /b α is curlfree. It can be easily constructed when h α (R) = h α (R α ) and b α (R) = b α (R α ), which is the case in all examples considered here, giving:

B. Self-renewing resources
In MacArthur's Consumer Resource Model [2,7,32], the environmental state is fully described by the abundances R α of M substitutable resources, which are consumed by organisms of species i at a rate −q iα = c iα R α . The per-capita growth rate g i of each species is proportional to a weighted sum of per-capita consumption rates, minus a "maintenance cost" m i that determines the consumption required to maintain a constant population size. The weights w α in this sum measure the "quality" of each resource (e.g., energy density or carbon content). The resources are "self-renewing," with intrinsic dynamics that produce exponential growth up to a finite carrying capacity R 0 α . These definitions and assumptions lead to the following set of steady state conditions: To achieve the correct relationship between the impact vectors and the growth rates, we divide the resource equation by R α /w α , as illustrated in Supplementary Figure S1: The supply vector field can now be written as the negative gradient of a weighted Euclidean distance: In the main text, we considered for simplicity the case where w α = 1 and r α = 1 for all α, so that this becomes the ordinary Euclidean distance. When these factors differ from one another, the natural dissimilarity measure gives more weight to resources with higher energy density or faster resupply rates. Dividing by R α requires some care, however, since it is possible for some R α to vanish in the steady state. Examining the dynamics, we see that the steady state is only stable against invasion by extinct resources if whenever R α = 0. Since the gradient of the objective function points in the positive direction along this axis, all positive values of R α are further away from the optimum. We can therefore accommodate this possibility within the constrained optimization framework, by simply adding a new set of constraints:

C. Externally supplied resources
A straightforward extension to MacArthur's Consumer Resource Model replaces the logistic dynamics of the resources with linear dynamics, which better capture the behavior of externally supplied abiotic nutrients: Performing the same transformation as before to make q i = −∇g i , we obtain The supply vector field can now be written as the negative gradient of a weighted Kullback-Leibler divergence: In the main text, we considered for simplicity the case where w α = 1 and τ α = 1 for all α, making d into the usual KL divergence. As before, when these factors differ from one another, the natural dissimilarity measure gives more weight to resources with higher energy density or faster resupply rates.

D. Essential resources
In Figure 3, we consider the sequential incorporation of non-substitutable resources described in [33]: where γ is a constant with units of inverse time, and the c iα are constants with the same units as R α , based on the ratios of resources required for building biomass.
For the example presented in Figure 3(a), we set the consumption rate per unit resource equal to the marginal benefit accrued from a slight increase of resource α: which becomes small if any resource is in short supply, because all resource types are required for growth. This resulted in q i = −∇g i after dividing the resource equation by R α , so that the KL divergence is minimized under linear ("externally supplied") resource dynamics. Another reasonable choice for the impacts is to make them equal the growth rate times the biomass stoichiometry: Comparing this to the expression for the gradient of g i , we find that in the general form of Equation (3) above. This gives: which is the negative gradient of

E. Environmental Feedbacks
The environmental feedbacks presented in Figure 3(b) are modeled by making each species return a fraction l α of the "energy" acquired from each resource back to the environment, after transforming this energy into other resource types as specified by a universal chemical matrix D αβ [19,21]: The resource equation can be rewritten in matrix form as: We can now put the impact vectors into the required form, by operating with the matrix inverse Q −1 of Q, whose elements satisfy β (Q −1 ) αβ Q βγ = δ αγ , and then dividing through by R α : The transformed supply vector field h(R) = now has a nonvanishing curl in general, as shown in Supplementary Figure S3, and cannot be written as the gradient of an objective function. But we can locally recover MEPP in the vicinity of a given point in resource space by approximating the supply vectors with a function of the same form that arises in the absence of environmental feedbacks: where the weighting factors arẽ w α = (Q −1 ) αα (28) and the effective unperturbed environmental state is given bỹ This state is a function of R, but it needs to be a constant for h(R) to reduce to the desired form. We can make it constant by evaluating it at a fixed value of R, but the point that minimizes the resulting perturbation measure d(R 0 , R) is not necessarily the steady state. If we now re-evaluateR 0 at the new value of R, we will in general find a new optimum. WhenR 0 (R) is evaluated at the true steady state location R = R * , however, then the approximate resource equation with constant R 0 (R) also balances at R * . This means that the steady state does minimize d(R 0 , R) for this choice ofR 0 , and an algorithm that iteratively optimizes d(R 0 , R) and up-datesR 0 will halt when it finds this point. We note that this algorithm is intimately related to Expectation Maximization algorithms in statistical inference and machine learning [16,22]. The pseudocode for this algorithm is given in Figure 3(f ) of the main text.

ACKNOWLEDGMENTS
This work was supported by NIH NIGMS grant 1R35GM119461, Simons Investigator in the Mathematical Modeling of Living Systems (MMLS) to PM. We would also like to thank Ching-Hao Wang, Jacob Ferguson and William Ludington for useful discussions.   (a) Impact vectors and ZNGI's for two species with environmental feedbacks, such that consumption of one resource is associated with production of the other, along with supply vectors corresponding to externally supplied resource dynamics. (b) Transformed impact and supply vectors, with the impact vectors now given by the gradient of the per-capita growth rate, and therefore orthogonal to the ZNGI's. The transformed supply vector field has a nonvanishing curl, and can therefore not be described as the gradient of any perturbation measure. It can be approximated by a curl-free field shown in gray, which becomes exact at the location of the steady state. The unperturbed environmental state corresponding to this approximate supply field is indicated with the gray 'x.'