Axiomatic Ecology: Rules for Building Consistent Ecosystem Models

Non-mechanistic ecosystem models are employed in many ecological studies ranging from purely theoretical to data-driven ones. With such models in mind, we derive fundamental consistency criteria (axioms) from rst principles. ese particularly cover what we call clone consistency: e outcome does not change if a population is split into two with identical properties. We mathematically prove that, these axioms are fullled if and only if the model is based on linear combinations of powers of parameters and abundances. Using this insight, we formulate a framework that allows to quickly assess the consistency of existing models and to build new models. We demonstrate our approach by invalidating a data-based model proposed for polymicrobial urinarytract infections and developing an alternative. We argue that our framework reveals implicit assumptions and informs the general modelling studies by narrowing the space of possible models or pointing to new forms of models.


I. INTRODUCTION
Many theoretical and semi-empirical modelling studies of ecological communities employ general models [ -], such as the popular Volterra model [ ]. Such models are o en not mechanistic or process-oriented, i.e., they do not explicitly feature agents of interactions, such as nutrients and toxins. As a result the equations governing each population all have the same form, and the species of a population only manifests in the values of the associated control parameters. ese parameters may describe the properties of a single population, the interplay of two populations, or higher-order interactions, i.e., e ects involving three or more populations [ , ].
ey are usually chosen randomly [ -] or determined from experiment [ -].
In particular for microbial ecosystems, recent advances in automatising experiments have made it feasible to determine interaction parameters for richer ecosystems [ , -], to quantify the interaction between two populations with more than one parameter [ ], or to measure higher-order interactions [ , ]. ese new experimental scenarios o en demand new models that can incorporate the respective data, in particular as there is no single answer as to how multiparameter or higher-order interactions should be measured [ , , , , ].
One di culty when building a model is to avoid inconsistencies. As an example, we compare two simulations of a bacterial community using the same model (Box ): e rst simulation is straightforward; we represent each strain of bacteria as one population. For the second simulation, we divide one population into two with identical properties. Allegorically, we paint half the individuals in a di erent colour. Although these two simulations describe an identical situation, their outcomes strongly di er (Fig. , solid and dashed lines), which is why we consider the model underlying these simulations as inconsistent. is di erence cannot be explained by numerical noise and sensitivity to initial conditions, as demonstrated by simulations with perturbed initial conditions (do ed lines) exhibiting a much smaller di erence. is model was taken from a study [ ] that featured communities containing two strains of the same species; therefore the case of populations with very similar properties is relevant here. In Appendix B we provide an explicit arithmetic demonstration of this inconsistency that is independent of some choices we had to make for this example (see Appendix A).
We here formulate criteria (axioms) that exclude such fundamental inconsistencies. Applying methods from functional analysis, we explore the consequences of these axioms and derive a framework that allows to easily decide whether a given model is consistent and to build models that are. We demonstrate the la er by building a new semi-empirical model for polymicrobial urinary-tract infections. Moreover, we argue that our theoretical results are relevant to all kinds of ecosystem models and may inform them with respect to their consistency or implicit assumptions required to maintain this consistency.

Box . Our Case Study: Semi-Empirical Models of Urinary-Tract Infections
Ref.
used a high-throughput approach to systematically measure ecological interactions between strains isolated from polymicrobial urinary-tract infections (UTI). Employing this data, Ref. also proposed a model for ecosystems consisting of such strains. We use this scenario and the model as an example for applying our concepts.
We brie y summarise the measurements of the growth characteristics and pairwise interactions of these bacterial strains: Le : Each strain j was cultivated for 48 h in arti cial urine (solid bold le ers represent individuals of the respective strain). e exponential growth rate g j as well as the carrying capacity c j (named yield in Ref. ) were experimentally determined via optical densities. For convenience, abundances of each strain were normalised such that c j = 1. Right: Moreover, for each strain k, a conditioned medium was produced by le ing the strain grow for 48 h, mechanically removing the bacteria to obtain a supernatant, (represented by outline le ers) and mixing the result with fresh medium in a ratio of v := 0.4. In each such medium, each strain j was cultivated, and the conditioned growth rate g jk and carrying capacity c jk were determined as above.
In the model proposed by Ref. using this data, the (normalised) abundance x j of population j is described by the following di erential equation (with z := max(0, z)): ( ) with a jk := g jk gj − 1 and b jk :=

II. IMPACT FUNCTIONS
We consider functions that describe the impact of an ecosystem consisting of n populations on a species in that ecosystem, a resource, or similar. Examples of phenomena described by such impact functions include: • the e ective growth rate of a given species • the remaining size of a niche, • the availability of some nutrient, • reproductive services, e.g., pollination. e arguments of these functions are the abundances of all populations in that ecosystem and m parameters per population quantifying its impact.
Before we formulate our criteria, we formalise this scenario and introduce some helpful notational conventions: We denote by X = R n + the space of all possible population abundances, where R + denotes the non-negative real numbers. Also, we denote by A = R n×m the space of all possible parameter con gurations of these populations. e domain of impact functions is thus X × A. Furthermore, we denote an arbitrary pair of arguments for impact functions by (x, a), where x := (x 1 , . . . , x n ) ∈ X and a := (a 1 , . . . , a n ) ∈ A with a i ∈ R m being the parameter values that describe population i. In general, lowercase italic le ers denote numbers or parameter con gurations (tuples of numbers); Greek le ers denote functions; boldface le ers denote vectors or similar; uppercase le ers denote sets of respective contents. Finally, we use non-italic sans-serif le ers to identify modi cations of speci c components of arguments of these functions (similar to named arguments in many programming languages) to obtain notational abbreviations like the following: φ(x, a, x 2 = y) := φ((x 1 , y, x 3 , . . . , x n ), (a 1 , . . . , a n )).
Here the arguments of the function φ are x and a except for the abundance of the second population (x 2 ) being changed to y.

A. Ensuring Consistency
As our requirements have a very broad scope of application and to distinguish them from model-speci c criteria, we call them axioms. e axioms we require an impact function φ : X × A → R to ful l are: I : ere are no population-speci c mechanisms in the general form of the model -the properties of a given population are completely captured by its associated parameters. A speci c model can still feature populationspeci c mechanisms, e.g., when a speci c parameter has a value of 0 for all but one populations (also see Axiom I ). is way, the axiom can also be ful lled by a typical mechanistic model.
Following mathematical nomenclature, we refer to this axiom as commutativity. φ(x, a, x 1 = 0, a 1 = a 1 ) = φ(x, a, x 1 = 0, a 1 = b 1 ) ∀ a 1 , b 1 ∈ R m I : Parameters quantify the impact of the associated population and are scaled such that all parameters associated to a given population being zero corresponds to no impact of that population ( Fig. , bo om le ): Note that while we chose zero as the parameter value corresponding to no impact, it is straightforward to translate our results to any other choice of this constant.
I : Suppose two populations have identical parameter values. is means that they feature identical individuals (clones) within the limitations of the model. en the impact of these two populations should only depend on their total abundance, and not on how it is distributed onto the two populations ( Fig. , bo om right): If this requirement is not ful lled, the resulting model can produce di erent outcomes when implementing the same scenario in di erent ways, as exempli ed in Fig. . We refer to this axiom as clone consistency.
Note that through commutativity (I ), the other axioms apply to all arguments or pairs of arguments, respectively. Also, clone consistency (I ) of more than two populations is covered by applying the respective axiom repeatedly. Furthermore note that the above axioms do not explicitly capture (but also do not exclude) the case that a parameter is associated with more than one population, which is relevant for a higher-order interaction.

Box . e Functional Algebra of Impact Functions
In this box we explain the mathematical concepts of the following theorem and provide the main parts of a proof: where Γ is the set of constant functions. Let Φ be the smallest closed functional algebra that contains Ξ. en Φ contains all impact functions.
at the impact functions form a functional algebra Φ means that each product or sum of two impact functions is again an impact function and that each multiple of an impact function is an impact function. is algebra being closed means that the limit of uniformly converging sequences of impact functions is again an impact function.
Φ being the smallest such algebra containing Ξ means that all impact functions can be built from elements of Ξ using addition, multiplication, and taking a limit. In mathematical terms, Φ is the generated set of Ξ, and conversely, Ξ is the generating set of Φ, where taking the limit of a uniformly converging sequence is considered amongst the generating operations.
To prove eorem , we apply Bishop's eorem [ , ], which states, when reduced to algebras of real-valued functions: Bishop's eorem. Let Z be a compact Hausdor space. Let Ψ be a closed unitial subalgebra of As Z can be any su ciently large compact subset of X × A, we can identify Ψ with Φ, i.e., the generated set of Ξ. Unitial means that the functional algebra shall contain the constant functions, which is ful lled by construction for Ψ. us, to show that the functional algebra Ψ contains all impact functions, it only remains to be shown that for an arbitrary impact function φ for any x,x ∈ X and a,â ∈ A: or, in the language of functional analysis, Ψ has to point-separate, except where impact functions do not point-separate either. Since point-separations are una ected by algebraic operations of functions and limits, Ψ is point-separating, if and only if Ξ is. Moreover, since the functions form Γ are constant everywhere (and thus point-separating nowhere), this is equivalent to Λ being point-separating. Finally, since for any i = j, the functions from Λ i are constant wherever the functions from Λ j are not, it su ces to consider one Λ i only. is case in turn is covered by the following lemma, which we prove in Appendix C: Lemma . Suppose x,x ∈ X and a,â ∈ A are such that: Examples of impact functions that will be important in the remainder of the paper are: • Linear combinations of parameter values and abundances: x, a → n i=1 a i x i , where Y → Z denotes the function that maps Y to Z. ese are for example featured in the interaction term of most variations of the Volterra model [ ].
• Constant functions: x, a → h, with h ∈ R. ese cover the case that some quantity is not in uenced by the ecosystem at all.

B. Basic Building Blo s of Ecosystem Models
In Box we state and prove a theorem on the composition of impact functions, which facilitates testing whether models comply with the axioms and building models that do. It has two important implications for modelling: e rst is that everything built from impact functions via addition, multiplication, function composition, etc. is again an impact function, formally: If φ 1 , . . . , φ l are impact functions and χ : R l → R is an arbitrary function combining the results of these impact functions, then: x, a → χ(φ 1 (x, a), . . . , φ l (x, a)), is also an impact function. As an illustrative example consider a population that can diauxically live on two nutrients P and Q with e ciency e P and e Q . On the one hand, we can use two impact functions γ P and γ Q to describe the concentrations of these nutrients in dependence of the other populations in the ecosystem that produce or consume them. On the other hand, we can also use the impact function x, a → e P γ P (x, a) + e Q γ Q (x, a), which describes how much our population can grow on the nutrients provided by the ecosystem. e other implication of the theorem is that all impact functions can be wri en using (possibly in nitely many) elements from a set Ξ of basic building blocks, namely linear combinations of powers of parameter values and abundances and constant functions. Conversely, all models complying with our axioms must be decomposable into these basic building blocks. For application, it is arguably more convenient to use the equivalent set of basic building blocks that have the form: Here, the functions ζ and κ may comprise in nite series of elements from Ξ, thus simplifying the representation. On the other hand, all elements from Ξ can be represented in this form. We refer to these as basic impact functions from now on. A corresponding basic building block featuring general second-order interactions is: which is an impact function no ma er whether a ij is considered a parameter associated to population i or to j. Analogous building blocks exist for even higher interaction orders. Such building blocks are featured in existing models of higher-order interactions [ , , ]. Many ecosystem models feature a change of abundances, ). If we assume for simplicity's sake that there is no delay, noise, or similar, we can write the right-hand side R of the equation in the form: If, similar to Axiom I , we consider the case of two populations j and k with identical properties with abundances y and z, their total growth must be the same as if all individuals were assigned to one population: Using that j and k are identical as well as the properties of impact functions, we can conclude from this that: . erefore β must be proportional in its rst argument. Practically this means that all dependencies of R j on x must either happen within an impact function or in the form of a single factor x j .
We can now easily verify for most models whether they comply with our consistency requirements by checking whether they are built from impact functions. is in turn we can do by looking for terms of the above form (Eq. ). For instance, the Volterra model [ ] (with the carrying capacity for each population normalised to 1) can be rewri en as follows: with a jj = −1. is is clearly built from a linear combination and a factor x j and we can thus be sure that it complies with our axioms.
In another example, we can look at the model given by Eq. from Box within our framework: Both the capacity and growth term correspond to impact functions, and a and b are the parameters quantifying these impacts. However, it is clear that neither is built from linear combinations (with complete sums). From this, we can deduce that the model violates at least one of our axioms. Since it clearly ful ls Axioms I -I , it must fail to be clone-consistent (Axiom I ), which we observed in Fig. . Note that completing the sums with appropriate choices of a jj and b jj would not su ce to address this since it still leaves a solitary x j in the numerator of the capacity term.
Beware that a representation using the above basic impact functions may not be immediately obvious, as in the le -hand sides of the following examples: (a 1 − a 2 ) 2 a 1 a 2 x 1 x 2 = (a 1 x 1 + a 2 x 2 ) a 3 1 x 1 + a 3 2 x 2 − a 2 1 x 1 + a 2 2 x 2 2 .
However, we expect such functions to be a rare occurrence in application.

III. BUILDING MODELS
We can use impact functions as an ansatz for constructing consistent ecosystem models tailored to a given experimental scenario. is is particularly relevant when deducing a model from assays of a huge number of populations, which become increasingly available thanks to high-throughput experiments. Such assays usually do not provide su cient information to include all relevant agents (nutrients, toxins, etc.) into the model, thus rendering the models nonmechanistic. Moreover, in such approaches, one can rarely exclude the case of populations with nearly identical properties. When building a model (that does not feature higher-order interactions), it is o en appropriate to assume that there is one building block of the form of Eq. for each of the m experimentally determined interaction parameters. Fewer building blocks would mean unused parameters, while more building blocks would result in models that are usually overly complex. For illustration of the la er, the arguably simplest impact function that requires two building blocks of the form of Eq. while featuring only one parameter is: x 1 , x 2 , a 1 , a 2 → (a 1 x 1 + a 2 x 2 ) a 2 1 x 1 + a 2 2 x 2 . In many applications, such complexity does not level with our knowledge about the system and would thus violate Occam's razor by trying to model observations in an overly complex manner. If we have only one building block for a given parameter, we can simplify Eq. to: is is appropriate if a and κ are beyond the scope of the model in question, for example because they are not experimentally accessible.

A. Example: Deducing a Model for Urinary-Tract Infections
As a demonstration of our framework, we apply it to construct a new model for the scenario of polymicrobial urinarytract infections from Ref. (see Box ). is model shall use the same data and also employ ordinary di erential equations. We have two experimental interaction parameters and therefore make an ansatz using two basic impact functions for the reasons outlined above: where a j1 := r j and a j2 := s j . Like the model from Ref. (Eq. from Box ), we assume that a population's abundance also represents its footprint, i.e., the nutrients, toxins, and other relevant substances produced or depleted by that population. Hence featuring the death of individuals in our model could lead to implausible outcomes, as it would undo its footprint. is simpli cation is justi ed as we assume the major cause of declining populations to be dilution of the entire system, which also a ects the footprint. Note that this is another limitation of the model from Ref. , as it allows populations to decline without dilution (see, e.g., Fig. ). We furthermore assume that there is no lag phase caused by populations adapting to a new environment since we lack the data to quantify it.
In the situations that were experimentally investigated, this model should reproduce the observed growth rates and capacities. First, in the absence of other strains, the initial exponential growth rate of strain j should be g j : where 0 denotes a length-n vector of zeroes (abundances are zero unless speci ed otherwise). Also, if all other strains are absent and strain j has reached its saturation abundance c j = 1, it should not grow anymore: , we assume that the medium conditioned by strain k is equivalent to an ecosystem where the strain k is xed to an abundance of v, as the footprint of strain k is the same in both situations and makes up for most of the interaction between the strains. For the medium conditioned by strain k, the initial exponential growth rate of strain j should be g jk : Moreover, strain j should stop growing when it reaches c jk : Eqs. -, it is straightforward to derive constraints on the functions ρ j and ς j and to determine the parameters r jk and s jk in dependence of these functions, g j , g jk and c jk (see Appendix D). Making simple choices for ρ j and ς j within these constraints and accounting for singularities and discontinuities (see Appendix D), we arrive at the following model (with z := max(0, z)): ( ) To compare this model with the previous one (Eq. from Box ), we use both to simulate the scenario of a growth experiment with regular dilution that has a known outcome (see Fig. ). We note that quantitatively predicting experimental scenarios such as these without in-depth knowledge about the involved strains is a highly di cult challenge. Moreover, the high-throughput interaction data is restricted; for example, it does not feature higher-order interactions, and the supernatant used to determine interactions will not contain toxins whose production is triggered by products of its target. We therefore do not expect either model to have a high absolute predictive power. We assume that the abundances had converged when they were experimentally observed (see Fig. ), and thus the experimental results are not a snapshot of an oscillatory behaviour (beyond the expected e ects or a regular dilution). is is corroborated by the replicates being in good agreement and oscillations only occurring in one of the simulations. We nd that both models are in equally good or bad agreement with the experiment for six communities ( , , , , , and ), while the predictions of Eq. are better for two communities ( and ). Given the low number of samples, we refrain from further quantifying the agreements.
ese results indicate that models satisfying our criteria are at least equally suitable for describing ecosystem dynamics. (cf. Fig. S B) and simulations with our model (Eq. ) and the model from Ref.
(see Box ). Communities are numbered as in Ref.
. For the in-vitro experiments, each strain was inoculated at a xed optical density (OD600 = 0.001) grown for 4 × 30 h and diluted by a factor of 1 36 in between. Each experiment was performed in triplicate. Final abundances were determined from colony-forming units on Chromagar. e simulations mimicked the experimental procedure. Diamonds at the bo om indicate an abundance of the respective population between 10 −3 and 10 −2 . Community featured two strains of Enterococcus faecium, for which we only report the summed populations since they could not be distinguished in experiment. All simulations converged, except for the model from Ref.
and Community , which exhibited strong chaotic uctuations. See Appendix A for details of the simulation.
is expected to some extent since both models have the same xed points if the growth term is ignored and c jk < 1 ∀j, k and thus can be expected to yield similar nal states (see for example Fig. ).

IV. DISCUSSION
While we motivated and exempli ed our framework with applications to non-mechanistic models based on highthroughput data, our axioms can be required for all models that operate on the population level. We consider exceptions from this to be rare and deserving justi cation. For instance, if we can be sure that populations do not interbreed (i.e., contain species as de ned by Mayr), replacing a population with two copies with half the abundance a ects the availability of partners for sexual reproduction. us, the function modelling this availability must not be clone-consistent. While many popular models such as most variations of the Volterra model [ ] comply with our axioms, others do not (the models from Ref. -, Eqs. , , and in Ref. , the NFR model in Ref. , Eq. in Ref. , Figs. b and c in Ref. , and Eq. from Box ). However, for reasons we elaborate in the following, this does not necessarily mean that they should be dismissed outright.
In case of mechanistic models, applying our framework usually means that each mechanism is covered by its own (typically basic) impact function with many zero parameters, namely whenever the respective population is not involved in the mechanism. Taking this approach to the extreme, it is possible to extend every term to an impact function with enough assumptions. For instance, suppose a population j exclusively occupies its niche (which implies no other pop-ulations with identical properties) and u − x j is the size of the remaining niche available to this population. en we can extend the la er term to an impact function: where a jk describes to what extent population k occupies the niche of population j and per our assumptions has the values a kk = 1 and a jk = 0 for j = k. Here, the value of our framework is to prompt the question: What assumptions need to be made to comply with the axioms and are these assumptions justi ed? Many pure modelling studies use a model of the general form:ẋ where g j is the unperturbed growth rate of population j. If all η jk are linear, this kind of model employs a single basic impact function and requires no further assumption. If, on the other hand, the η jk are non-linear, each of them has to correspond to one basic impact function, which requires justi cation that each of them only depends on one population. One assumption yielding this would be that each η jk re ects one interaction mechanism which is exclusive to the impact of population k on population j. Moreover, models of the above form can either feature 1 or n impact functions (depending on whether η jk is linear), but do not capture the middle ground in between. Since these basic impact functions can be roughly associated with interaction mechanisms, this limitation may be relevant beyond our framework. As an alternative lling this gap, our framework suggests general models of the form: where l is the number of impact functions. In Sec. III A we used such a model with l = 2 as an ansatz (g j was incorporated in the rst factor of the big product). To cover higher-order interactions, the approaches of Eqs. and can be combined. Amongst others, the above general model may inform studies employing random interaction parameters (e.g., Ref. ), general modelling [ ], or machine learning by narrowing down or expanding the space of possible models taken into consideration. We note that for a high l, it may make sense to regard a distribution of parameters that contains a considerable amount of zeroes, corresponding to populations not participating in some mechanism. Finally we emphasise that our approach is general and not limited to the types of models featured in our examples: It is not restricted to models employing ordinary di erential equations, but can also be applied to models with discrete time, noise, or delay. Also, higher-order interactions are covered by our framework. Moreover, while we mainly used impact functions to describe the impact of a set of populations on a population, both the targets and the source may be something else, e.g., a resource or toxin concentration or an aggregated observable such as the pH value. For example, the impact of the ecosystem on substances in the model proposed by Ref. can be described in terms of impact functions. We particularly note the parallels to two pharmacological approaches to describe the cumulative e ect of two drugs that are not synergistic or antagonistic, i.e., are not subject to higher-order interactions [ , ]: • Loewe additivity, which is based on arguments similar to clone consistency (I ) and holds if the two drugs a ect the same component of the cell, • Bliss independence, which violates clone consistency (I ) at rst glance and holds if the two drugs affect di erent components of the cell.
In our framework, drugs that target the same cell component would correspond to using the same interaction mechanism and thus would be captured by the same basic impact function. e e ect of a complex drug cocktail could be captured by several Bliss-independent basic impact functions, each of which comprises a series of Loewe-additive components.

V. CONCLUSION
We introduced impact functions, which are basic building blocks of ecosystem models adhering to axioms capturing fundamental consistency requirements. We mathematically proved that such impact functions must have linear combinations of parameters and abundances at their core. We demonstrated that impact functions can be used to easily build models complying with the axioms to capture a speci c experimental scenario -though one should still beware that the model makes sense in respects not covered by the axioms. Our framework also informs the form of more general models, pointing out potential new directions of research in this area or outlining the space of possible models. Conversely, the absence of impact functions in a model points out that it makes implicit assumptions or has consistency problems. We are therefore con dent that our framework can form a backbone for a wide range of ecological modelling studies. .

Appendix A: Simulation Details
For simulations with both, the model from Ref.
(Eq. from Box ) and our model described by Eq. , we use the data from Ref. as is, with the exception of data describing interactions between two identical strains: First, in the experiment, a strain cannot grow on the portion of the medium that is its own supernatant, but only on the portion that is fresh medium, which makes up 1 − v = 0.6 of the medium. We set c ii = 1 − v to adhere to this ideal. Second, in the medium conditioned by itself, a strain's growth rate should at best slightly lower than in an unconditioned medium and at worst be proportional to the concentration of nutrients, and thus to 1 − v. We therefore restrict g ii to the interval [(1 − v)g i , (1 − )g i ] with = 0.01. Without these adjustments, we would obtain implausible results, e.g., in case of c ii > 1 − v, the respective strain could never stop growing since it e ectively increases the size of its own niche.
We performed all simulations with JiTCODE [ ] using the DoPri method. To obtain continuity as required by the integration method, we approximate z = max(0, z) ≈ 1 2 z + √ z 2 + 2 with = 0.001. We converted abundances in the simulation results to optical densities by undoing the respective normalisation of abundances ( xing c j = 1 in Box ). We then converted the optical densities to displayed abundances (in Figs. and ) by approximating that optical density is proportional to biovolume.

Appendix B: Arithmetic Example for Inconsistencies in a Model
We here provide an arithmetic example for an inconsistent behaviour of the model described by Eq. from Box . is example features no growth interaction and thus does not rely on how we chose identical populations to a ect each other's growth (see Appendix A) -a choice that is not clear without experiment.
We consider the case of three populations {1, 2, 3} =: J with the rst two populations having identical properties. We choose g j = 1 ∀j ∈ J, a j,k = 0 ∀j, k ∈ J ×J, i.e., the growth term is not a ected by interaction and is always 1. Finally, we let the coe cients of the capacity term be: b 12 = −1 and b 21 = −1 re ects that two populations with identical properties deplete each other's niches. Now, consider two states of the ecosystem x = 1 4 , 1 4 , 1 and x = 1 2 , 0, 1 . As the rst two populations are indistinguishable, these states describe an equivalent situation. us, they should also evolve equivalently, i.e., the temporal derivative of the summed populations 1 and 2 should be the same in both cases:ẋ 1 +ẋ 2 =ẋ 1 +ẋ 2 =ẋ 1 . However, If z + p + z − p = 0 andẑ + p +ẑ − p = 0 the above equality will be dominated by a j r for j → ∞, which gives us: In case that exactly one of z + p + z − p andẑ + p +ẑ − p is zero, the above yields a contradiction due to the limit evaluating as either 0 or ∞, and hence this cannot be the case. If both are zero, the concluded equation holds without further ado. Analogously, we obtain: We choose the arguably simplest function to ful l the criteria for ρ, namely ρ j (z) := g j +z. is has the consequence: A group of functions ful lling the criteria for ς is: ς j (z) := 1 − z q with q > 0 and z := max(0, z). Here, the free parameter q controls how early and smoothly the saturation e ect of a depleted niche kicks in. Note that this choice results in terms similar to what Ref. named hyperlogistic.
Finally, like Ref. , we constrain the growth and capacity term to be non-negative. Pu ing everything together, we arrive at the model: A problem with this model is that for 0 < x k < 1, we have: lim c jk →0ẋj = lim c jk →0 R j (x, a) = ∞. Now, c jk = 0 means that there is no growth of strain j in the medium conditioned by strain k and thus we already have a problem with experimentally determining g jk . us, one might argue that the actual point of the singularity requires a dedicated case distinction anyway. However, lim c jk →0ẋj = ∞. also means thatẋ j becomes arbitrarily large for small c jk . A way to address this problem is to consider the case q → ∞, or more speci cally: In this case, the term ς(1 − c jk ) in Eq. D can be assumed to be 1 (otherwise, we would have the aforementioned problem of not being able to experimentally determine g jk ). is eliminates the singularity, but also renders the model not continuously di erentiable.
In our simulations, we therefore make a trade-o between complying with Eq. and the numerical bene ts of a continuously di erentiable model by se ing q = 10 and approximating ς(1 − c jk ) ≈ lim p→∞ ς(1 − c jk ) = 1 in Eq. D , thus arriving at: Biosci. , ( ).