What models are clone-consistent?

In this section, we first define impact functions, which are fundamental ingredients of ecosystem models that allow us to mathematically encode our consistency criteria. We then expose the consequences of our consistency criteria – first for impact functions and then for entire models. In particular, we present recipes for checking and building models.

Defining and constraining impact functions

Impact functions describe the impact of a community on a species, on a resource, or on any other relevant feature of the ecosystem that is captured in a model. Features and phenomena described by impact functions include:

the effective growth rate of a given species,

the remaining size of a niche,

the rate of predation,

the availability of a resource or, if the resource is a dynamical variable, its consumption and production,

reproductive services, e.g. pollination,

the amount of crowding, and

general interaction terms, e.g. the sum in the generalized Lotka–Volterra model [1].

The arguments of impact functions are the abundances of all populations in the ecosystem x = (x₁, x₂, …, x_n) and parameters a = (a₁, a₂, …, a_n) that quantify the impact of the populations. Often the impact of a population i is described by a single number. Yet our results also hold for the more general case that m parameters per population are required, i.e., $$ a_i=(a_{i1},a_{i2},\dots ,a_{im})\in {\mathbb{R}}^m$$.

A prominent example of an impact function is:

where $$ \mathcal{Y}\mapsto \mathcal{Z}$$ denotes the function that maps $$ \mathcal{Y}$$ to $$ \mathcal{Z}$$. This is employed for the interaction term in the generalized Lotka–Volterra model [1] among others. In general, impact functions can take many forms.

We require impact functions to fulfill the following basic criteria (illustrated in Fig 2; see Methods, The functional algebra of impact functions for mathematical formulations):

I1 Commutativity: The properties and idiosyncrasies of a given population are exclusively captured by its associated parameters, as opposed to dedicated mathematical terms in the function. This is equivalent to pairs of abundances and parameters ((x_i, a_i)) being interchangeable as arguments of the impact function.

I2 When a population is absent, its associated parameters have no effect on the value of the impact function.

I3 When each parameter associated with a given population is zero, that population’s abundance has no effect on the value of the impact function.

I4 Clone consistency: If two (or more) populations have identical parameters, the value of the impact function must only depend on their summed abundance and not on its distribution among the two populations.

Note that this criterion is a special case of clone consistency as described in the introduction. Using impact functions complying with this criterion is therefore necessary for a model to be clone-consistent, but not sufficient (we will address clone consistency on the scope of the entire model later). Due to the narrower scope of this criterion, the populations do not need to be identical in all respects, but only in the parameter(s) used by the respective impact function.

Fig 2

Criteria for impact functions exemplified for the grazing impact of animal populations.

Each island represents a community. Each row and color represents one population in the model, with animal heads representing individuals. Numbers on the left represent parameters governing the respective population (grazing rate in the example), and head shapes indicate whether populations have identical properties as per these parameters. The similar sign (∼) indicates that two communities are equivalent as arguments of an impact function, i.e., they should yield the same result (total amount of grazing in the example).

Note that it is often reasonable to choose a considerable portion of parameters to be zero. For example, if our impact function describes predation loss of a given focal species, we would choose a_i = 0 for all populations i that do not prey upon the focal species. In the following we call a function impact function if and only if it satisfies these criteria.

The form of impact functions

While the criteria I1–I4 for impact functions are conceptually simple, it is not straightforward to test directly if a given function complies with them or to devise a new function that does. To address this issue, we investigated the functional form of impact functions and found a set Ω of basic impact functions, which can serve as building blocks for ecosystem models.

These basic impact functions are linear combinations of abundances and (potentially transformed) parameters, i.e., all functions of the form:

where $$ \zeta :\mathbb{R}\to \mathbb{R}$$ is an arbitrary parameter transformation with ζ(0) = 0. Often, ζ is the identity function (id), simplifying Eq 5 to Eq 4. Using methods from functional analysis, we mathematically proved that any impact function can be built from these basic impact functions via addition, multiplication, function composition, and similar operations (see Methods, The functional algebra of impact functions). Conversely, everything built from those elements or other impact functions is again an impact function. Formally, the general functional form of impact functions is:

where ω₁, …, ω_l ∈ Ω are basic impact functions as per Eq 5 and $$ \chi :{\mathbb{R}}^l\to \mathbb{R}$$ is an arbitrary function combining their results, such as a product, sum, or a more complex function.

To illustrate the composition of impact functions from these building blocks, we consider the case of a single population of flowering plants that may be both pollinated and grazed upon by several insect populations [64] as a toy example (also see Fig 3, bottom). We use a function ρ_graz to describe the rate at which insects (and their larvae) graze on the plants:

where a_i is the grazing rate of insect population i. The function ρ_graz has the form of Eq 5 (with ζ = id) and therefore it is a (basic) impact function. We use a function ρ_poll to capture the rate at which insects fertilize flowers:

where b_i is the contribution to pollination provided by insect population i. A Holling type-II response here ensures that the fertilization rate saturates at 1. Again, we can see that this is an impact function as it is built via composing $$ z\mapsto \frac{1}{1+z^{-1}}$$ and the basic impact function $$ \tau (\mathbf{x},\mathbf{b})≔{\sum }_{i=1}^n{b}_i{x}_i$$. Finally, we combine the two impact functions to a function ϕ that describes the relative change of the plant population due to insect influence:

where g_max is the growth rate of the plant population in the absence of death and with maximum fertilization. As ϕ is built from impact functions, it is an impact function itself. We can also see directly that ϕ complies with the general form of an impact function (Eq 6) by choosing $$ \chi (y,z)≔\frac{g}{1+y^{-1}}-z$$, ω₁ = τ, and ω₂ = ρ_graz. As a result, we can be certain that the resulting model satisfies the fundamental consistency criteria. Similarly, the basic building blocks we identified enable constructing consistent models of any other ecosystem.

Fig 3

Recipes for building and validating models using our framework.

In general x_i denotes the abundance of population i, and a_i, b_i, c_i, and d_i are parameters describing its impact. The example for checking is based upon Eq 21 and tailored for covering relevant cases. The example for building extends the one from The form of impact functions. For both examples, the biological background is discussed in more detail in the main text.

Assessing and building models

Our treatment of impact functions allows us to ensure clone consistency when modeling the effects of a community. To ensure clone consistency of an entire population model, we additionally need to require clone consistency when modeling the effects on the growth of a population. Specifically, we require that the total growth of two populations of identical individuals must be the same as if those two populations were joined in the model.

We consider models of the forms:

where the right-hand side $$ {\mathcal{R}}_j$$ describes the change of abundance of population j due to the internal dynamics of the ecosystem. For simplicity, we omit the dependencies of $$ {\mathcal{R}}_j$$ and ignore immigration and other external effects as well as dynamic variables that are not populations (such as explicitly modelled resources); adding them to such a model is straightforward.

We proved that $$ {\mathcal{R}}_j$$ must have the form (see Methods, Non-impact-function contribution to abundance changes must be proportional):

where ϕ is an impact function. This means that all dependencies of $$ {\mathcal{R}}_j$$ on x must either happen within an impact function or in the form of a single factor x_j. Intuitively, each individual multiplies with a rate that is the result of all impacts it experiences in the ecosystem (ϕ) – these impacts include interactions between individuals of the same population, e.g. due to crowding.

This insight provides an easy way to verify if models comply with our consistency criteria. We simply need to check if they have the form of Eq 11. To verify in turn if ϕ is an impact function, we can look for terms of the form of Eq 5. For instance, a common formulation of the generalized Lotka–Volterra model [1] can be rewritten as:

with a_jj = −1, which reflects that a population maximally competes with itself. With this, the model is clearly built from a linear combination and a factor x_j and we can thus be sure that the model is clone-consistent. We can now have another look at our introductory example (Eq 1): As the interaction term features logarithms of abundances, it does not comply with the form of Eq 5. Thus, the model violates our consistency criteria and the observed clone inconsistency (Fig 1) is inevitable. We summarize the recipe for checking a model and provide a more extensive example in the top of Fig 3.

This recipe can be inverted to build a clone-consistent model. In the example from The form of impact functions, we can directly insert Eq 9 into Eq 11 and obtain a model for the change of a plant population in light of pollination and grazing (see Fig 3 bottom). If we stringently apply our framework here, the only thing that distinguishes plants and insects is that the former have a grazing and pollination rate of zero.

Our framework can also be applied to experiments that do not assess the details of ecological interactions (nutrients, toxins, etc.) but only aggregated, phenomenological interaction observables, such as the carrying capacity of a population in the presence of another. This is typical for high-throughput experiments assessing the pair-wise interactions of microbial communities [15, 16, 19–22]. Combining Eqs 10, 11, 6, and 5, we obtain a general ansatz for such a model:

Here, we refrain from transforming the parameters (corresponding to ζ = id in Eq 5), since this would not affect the final model. First we identify the number of basic impact functions m as the number of interaction observables: Using fewer than m basic impact functions would entail that interaction observables would go unused, i.e. available information on the system would be ignored; more than m basic impact functions would make the model overly complex given our limited knowledge of the system. The choice of the combining function χ depends on the application, but we expect that a product of transformations of individual basic impact functions is often appropriate. In this case our ansatz becomes:

with some transformation functions η_k. With χ chosen, the parameters and functions (a_jik and η_k) can be determined using:

the requirement that the model should reproduce key characteristics of the experimental scenarios, e.g., the carrying capacity of a population,

ecological assumptions and facts about the scenario, e.g., that the predation rate should increase with the abundance of predators, or

assumptions of simplicity (Occam’s razor), e.g., that η_k should not be more complex than necessary to fulfill the experimental and ecological constraints.

In New model, we provide an example for this approach.

Higher-order interactions

Our framework is readily extended to models describing second- or higher-order interactions. For this, one simply has to consider parameters that are associated with more than one population. Then, it is convenient to use other basic building blocks (instead of Eq 5), for example for second-order interactions:

which is an impact function independent of whether a_ij is considered a parameter associated with population i or with population j. Analogous building blocks exist for higher interaction orders. Such building blocks are featured in existing models that capture higher-order interactions [12, 13, 17].

Case study: Semi-empirical models for microbial communities

As an instructive example, we apply the impact-function framework to describe the dynamics of a microbial community for which multiple ecological interaction observables were recently measured experimentally. We first describe the experiment and the general modeling challenge. Then we show that an existing model is clone-inconsistent and discuss the specific reasons and implications of this. Finally, we show how our framework can be used to build a new clone-consistent model.

A recent study (co-authored by one of us) used a high-throughput approach to systematically measure ecological interactions in microbial communities consisting of strains isolated from polymicrobial urinary-tract infections (UTI) [16]. For each strain, the exponential growth rate g_j and the carrying capacity c_j in isolation were measured (Fig 4, left). For convenience, abundances of each strain j were normalized such that c_j = 1. Furthermore, for each strain k, a medium partially conditioned by that strain was produced (it contains a fraction v of the supernatant). In each such partially conditioned medium, the conditioned growth rate g_jk and carrying capacity c_jk of each strain j were measured to quantify how strain k affects strain j (Fig 4, right). Creating a model from this dataset is particularly challenging since it has to include two interaction observables: growth rate and carrying capacity in conditioned medium.

Fig 4

Acquisition of data used in our case study.

Measurement of growth characteristics (left) and pairwise interactions (right) of bacterial strains isolated from urinary-tract infections [16]. Left: Each strain j was cultivated for 48 h in artificial urine. Solid bold letters represent individuals of the respective strain. The exponential growth rate g_j as well as the carrying capacity c_j (named yield in the original study [16]) were experimentally determined via optical densities. Right: For each strain k, a conditioned medium was produced by letting the strain grow for 48 h, mechanically removing the bacteria to obtain a supernatant, and mixing the result with fresh medium in a ratio of v ≔ 0.4. Outline letters (“footprints”) indicate to what extent the respective culture consists of supernatant. 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.

For modeling this system, it is a plausible assumption that the abundance of a population also represents its footprint, i.e. the nutrients, toxins, and other relevant substances produced or depleted by that population. The basis for this simplifying assumption is that populations are declining only due to dilution of the entire system, which has the same effect on the footprint. With this assumption, we can treat the medium partially conditioned by strain k as an ecosystem where the abundance of that strain is fixed to the corresponding fraction of its carrying capacity (x_k = vc_k = v).

The general form of an ordinary differential equation describing the (normalized) abundance x_j of population j in this system is $$ {\dot{x}}_j={\mathcal{R}}_j(\mathbf{x},\mathbf{a})$$. Such models should reproduce the observed growth rates and carrying capacities for all situations that were experimentally investigated. For instance, in the absence of other strains, the initial exponential growth rate of strain j in the model should be equal to its experimentally observed exponential growth rate g_j:

Existing model

The original study [16] proposed a model for communities consisting of such strains based on Verhulst’s logistic model for one population [65]:

where both the growth rate g and the carrying capacity c (normalized to 1 here) are modified by interaction terms incorporating the experimentally obtained parameters:

with ⌈z⌉ ≔ max(0, z), $$ a_{ji}≔\frac{g_{ji}}{g_j}-1$$, $$ b_{ji}≔\{\begin{array}{ll}{c}_{ji}-1\hfill & \text{if}{c}_{ji}\ge 1\hfill \\ \frac{1}{v}\left(c_{ji}-1\right)\hfill & \text{if}{c}_{ji}<1\hfill \end{array}$$, and the remaining symbols are as in Fig 4. Looking at this model through the lens provided by our framework, both the carrying-capacity and the growth term should be impact functions, with a and b being the parameters quantifying these impacts. However, it is clear that neither term is built from linear combinations (with complete sums), matching the form of Eq 5. Thus the model must violate at least one of our criteria. As it clearly satisfies Criteria I1–I4, it therefore must violate I4 and be clone-inconsistent which can indeed be shown explicitly (S1 Appendix). We can complete the sums with appropriate choices of a_jj and b_jj, and we can expand the solitary x_j in the numerator to:

where the new parameter d_ji can be interpreted as quantifying the extent to which population i occupies the niche of population j. However, this expansion implicitly assumes that the strains have non-overlapping niches, which is not justified for this study, as it features many communities containing two strains of the same genus or even species. This interpretation of d_ji would also ignore that the conditioned carrying capacities c_ji capture niche overlap. Following the recipe from Fig 3 (top), Question 5 must be answered with no at this point.

The fixed points of this model are characterized by:

Consequently, they are clone-consistent (also see S1 Appendix). Most conclusions of the original study [16] are based on these fixed points and thus unaffected by the clone inconsistency of the model. However, clone inconsistency affects the transient dynamics (see also S1 Appendix), which is relevant as these communities are subject to frequent dilutions (due to bladder voiding) which can happen long before the system has equilibrated.

New model

We use our framework to construct a new population dynamics model for this scenario. Since two experimental interaction observables are available, we make an ansatz using two basic impact functions (see Assessing and building models). As one of the observables and thus one of the impact functions captures the carrying capacity, we choose to multiply (and not add) the impact functions to ensure that an impact function can single-handedly reduce the growth to zero. Our ansatz is thus Eq 14 with m = 2:

where a_j1 ≔ r_j and a_j2 ≔ s_j.

Inserting this ansatz into Eqs 16–19 and making a few choices that do not affect generality already yields strong constraints on the functions ρ_j and ς_j and on how the parameters r_jk and s_jk relate to these and the experimental parameters g_j, g_jk and c_jk, namely (see Methods, Deriving a new model for UTI strains – the legwork):

By further making simple choices for ρ_j and ς_j within these constraints and accounting for singularities and discontinuities (see Methods, Deriving a new model for UTI strains – the legwork), we arrive at the following model (with ⌈z⌉ ≔ max(0, z)):

Like the existing model (Eq 21), this can be understood as an expansion of Verhulst’s logistic model with the following differences: First and foremost, the carrying-capacity term is simplified to achieve clone consistency – interactions affecting the carrying capacity of population j are now captured in a single sum. Second, the additional parameter q governs how abruptly the saturation effect kicks in. Third, the dilution factor v is included consistently without case distinctions. Fourth, as per the initial assumptions, populations cannot decline anymore (unless dilution is added to the model). The particular form of this model illustrates how challenging it can be to write down clone-consistent models from scratch without using the framework presented here.

We find that this model can explain observed species abundances and ecological stability in a small experimental dataset (S1 Appendix) at least as well as the previous model (Eq 21). Moreover, the fixed points of both models are the same under most conditions (S1 Appendix).

Implications

Specific implications: Clone-inconsistency points to implicit assumptions

While many popular models, including most variants of the generalized Lotka–Volterra model [1], comply with our criteria, others do not [50, 52: Eqs 9 and 10, 53, 54: Eqs 1.28–1.30 and 1.50, 55: Eqs 11 and 12, 56, 57, 58: Eq 5, 59: NFR model, 60, 16, 61, 62: Eq 3, 3: Figs 3b and c, 63: UIM and IIM model]. However, as we will elaborate below, this does not mean that these models should be dismissed outright.

Any model can be made clone-consistent by setting the right parameters to zero (Fig 3 and Eq 22). For instance, suppose that u − x_j quantifies the unoccupied portion of the niche of population j. We can extend this term to an impact function:

with a_ii = 1 and a_ji = 0 for j ≠ i. Here a_ji describes the extent to which population i occupies the niche of population j. Thus, a_ji = 0 for j ≠ i implies that population j exclusively occupies its niche. Whether this assumption of an exclusive niche is justified depends on the specific application. The strength of our framework is to systemically point to such (often implicit) assumptions and prompt an explicit justification or an improvement of the model if the assumptions are not justified.

As an instructive example, consider a recently proposed unique-interactions model (UIM) [63]. This model features three interaction terms. Two describe the influence of mutualistic or exploitative interactions, respectively, on a focal population j as:

where a_ji describes the strength of the respective interaction and h is the half-saturation constant. These terms are not impact functions and thus clone-inconsistent. However, a key assumption of this model is that mutualistic and exploitative interactions are unique for each focal population, which is a legitimate approximation for the purposes of that work. Following this assumption, each summand of Eq 28 represents a separate interaction mechanism, such as a specific resource, service, or mode of predation, and only population i affects the focal population j via this mechanism. In our framework, each summand can be expanded to an impact function, with only one parameter different from zero. By contrast the third term in the UIM model, which describes competitive interactions, is clone-consistent. For this term, unique interactions are neither assumed nor would this be justified.

A similar case of assumptions justifying clone inconsistency arises if some features of populations cannot be feasibly encoded in numerical parameters, which we assumed as given so far. For example, compatibility for sexual reproduction is tedious to capture in a parameter; instead, it is often reasonable to assume that populations do not interbreed (i.e., contain species as defined by Mayr). In this case, splitting a population in two equal parts also halves the availability of partners for sexual reproduction. Thus, this availability should not be described by a clone-consistent impact function. These examples illustrates how clone-inconsistent terms can make sense in a model if properly justified.

Many of the aforementioned studies based on clone-inconsistent models primarily make statements about the effects of model properties on population dynamics. If these findings are based on strong assumptions, their generality, relevance, and applicability are considerably diminished. Moreover, if these assumptions are implicit, this increases the risk that the study is misinterpreted and misapplied by others. The framework presented here will help to avoid these problems by forcing key assumptions of ecological models to be explicitly stated and justified.

General implications

While each instance of unjustified clone inconsistency reflects a shortcoming of the model, general statements about the consequences of clone inconsistency for the solutions of a model are most likely impossible. Clone inconsistency is tightly intertwined with the fabric of the model, and thus we cannot study its effect in isolation. Moreover, there is no reason to expect that clone-inconsistent models have any relevant commonalities, as they are extremely diverse. For example, it makes a difference whether joining identical populations increases or decreases some impact (that would be unchanged in a clone-consistent model). For illustration, the diversity of clone-inconsistent models may be compared to that of all numbers not divisible by seven. Similarly, it is likely not possible to draw general conclusions about the dynamic behavior of the clone-consistent models – such as favoring or suppressing oscillatory dynamics. While they are based on linear combinations, they can certainly be non-linear. Not only can a non-linear function still be applied to the linear combination (e.g., χ in Eq 6), but any existing non-linear model can be made clone-consistent by sufficiently strong assumptions. The central question is whether these assumptions are biologically justified. However, insights can be gained by interpreting each impact function as a mechanism of ecological interaction:

Many pure modeling studies use a model of the general form [12, 59, 62]:

where g_j is the unperturbed growth rate of population j. If all η_ji are linear, this kind of model uses a single basic impact function and requires no further assumptions. However, within our framework, non-linear η_ji can only be justified as follows: Each summand in Eq 29 corresponds to one interaction mechanism by which population i uniquely affects population j. In this case, each summand would be expanded to an impact function, where all parameters except one are zero. This would result in a total of n basic impact functions. A notable implication is that models of the above form can either feature 1 or n basic impact functions (depending on whether η_ji is linear), but they cannot capture the middle ground in between. As each basic impact function can be associated with one interaction mechanism, this limitation is relevant beyond the consistency issues addressed by our framework.

To fill this gap, our framework suggests an alternative form for general ecosystem models, such as:

where m is the number of impact functions (compare to Eq 14). Each factor in Eq 30 corresponds to one interaction mechanism, which can involve multiple populations. To capture that most populations do not participate in most interaction mechanisms, we propose to sample the interaction coefficients (a_ijk) from a distribution containing many zeros. By narrowing down or expanding the space of possible models taken into consideration, the alternative form for general models above can inform studies employing random interaction parameters [9–13], generalized modeling [25], or machine learning.

Conclusion

We introduced a framework for building ecosystem models using impact functions as building blocks. This framework is aimed at ensuring the clone consistency of models and thus constrains the possible choices of models. While at first this may seem like at a burden, we anticipate that it will rather facilitate the modeling process by guiding ecologists when choosing from the (still infinitely many) clone-consistent models. Our framework further prompts relevant questions about the underlying assumptions of models. Specifically, the absence of impact functions in a model exposes that it is clone-inconsistent, which may indicate a fundamental problem. Alternatively, clone-inconsistency can reveal implicit assumptions at the heart of the model, which need to be justified and may limit the model’s generality. Our framework also informs the form of more general models by outlining the space of possible models for ecosystems and enables new research directions in this field. Finally, our approach could be extended to implement criteria for specific ecological scenarios such as predation [27–30, 33]. Overall, the framework presented here provides a systematic way to understand models and can form the backbone for a wide range of ecological modeling studies.

Mathematical notation

lowercase italic letters: numbers or parameter configurations (tuples of numbers);

lowercase Greek letters: functions;

boldface letters: vectors or similar;

uppercase letters: sets of respective contents;

n: the number of populations;

m: the number of parameters per population;

$$ {\mathbb{R}}_{+}$$: the non-negative real numbers;

$$ X={\mathbb{R}}_{+}^n$$: the space of all possible population abundances;

$$ A={\mathbb{R}}^{n\times m}$$: the space of all possible parameter configurations of these populations;

X × A: the domain of impact functions;

x = (x₁, …, x_n) ∈ X: an arbitrary first argument of an impact function (abundances), where x_i is the abundance of population i;

a = (a₁, …, a_n) ∈ A: an arbitrary second argument of an impact function (parameters), where $$ a_i\in {\mathbb{R}}^m$$ are the parameter values that describe population i;

non-italic sans-serif letters: modifications of specific components of arguments of an impact function (similar to named arguments in many programming languages). For example: ϕ(x, a, x₂ = y) denotes ϕ((x₁, y, x₃, …, x_n), (a₁, …, a_n)). Here the arguments of the function φ are x and a except for the abundance of the second population (x₂) being changed to y.

$$ \mathcal{Y}\mapsto \mathcal{Z}$$: the function that maps $$ \mathcal{Y}$$ to $$ \mathcal{Z}$$ (anonymous function).

$$ \mathcal{Y}≔\mathcal{Z}$$: $$ \mathcal{Y}$$ is defined as $$ \mathcal{Z}$$.

Proof: Ξ generates Φ requires additional notation that is mostly introduced in that subsection.

The functional algebra of impact functions

In this subsection, we describe our main mathematical result and connect it to the main text, from encoding our criteria mathematically to translating the results back to application.

Expressed in equations, our criteria for impact functions are:

I1 Commutativity:



I2 When a population is absent, its associated parameters have no effect:



I3 When all parameters associated with a given population are zero, that population has no impact:



Note that the parameter value corresponding to no impact could be readily changed from zero to any other value.

I4 Clone consistency:

Note that through commutativity (I1), the other criteria apply to all populations or pairs of populations of the impact function ϕ, respectively (and not just to populations 1 and 2). Clone consistency (I4) of more than two populations is covered by applying the respective criterion repeatedly.

In the terms of functional analysis, impact functions form a functional algebra Φ. This 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. This algebra is also closed, which means that the limit of uniformly converging sequences of impact functions is again an impact function.

To easily build and detect impact functions, it is crucial to find a (small) set of impact functions from which all impact functions can be build, i.e., a generating set of Φ. Our main mathematical result is that Ξ = Λ ∪ Γ is such a generating set, where Γ is the set of constant functions and taking the limit of a uniformly converging sequence is considered amongst the generating operations. Λ is the set of all linear combinations of powers of parameters and abundances, i.e., functions of the form:

We will formally prove this in the next subsection (Proof: Ξ generates Φ), but the essential idea is this: Our criteria (I1–I4) require an impact function to have the same value on given subsets of its domain. For Ξ to be a generating set of Φ, it must reflect this: First, the functions in Ξ must be constant on each such subset, i.e., fulfill our criteria for impact functions. Otherwise Ξ would also generate functions that are not impact functions. (It is straightforward to show that the elements of Ξ fulfill our criteria, however, in our proof this is a by-product.) Second and more crucially, for each pair of points that are not in the same subset, there must be a function in Ξ that differs between these points (this is called separating points). Otherwise some impact functions could not be generated by Ξ. We show the latter by using our criteria (I1–I4) to systematically transform arguments of impact functions to a canonical form, in which populations are ordered by impact and maximally lumped together.

In application, the fact that Φ is a closed functional algebra, i.e., that limits remain within it, is relevant as it addresses the case of a (non-polynomial, continuous) function being applied to the result of an entire impact function or the parameters. We can rewrite:

Proof: Ξ generates Φ

We here state and prove our main mathematical result, namely:

Theorem 1. Let $$ {\Lambda }_k:=\{\mathbf{x},\mathbf{a}\mapsto {\sum }_{i=1}^n{a}_{ik}^p{x}_i|p\in \{1,\dots \}\}$$ denote the set of linear combinations of powers of values of the k-th parameter and abundances. Denote the set of all such functions as $$ \Lambda ={\bigcup }_{k=1}^m{\Lambda }_k$$. Let Ξ ≔ Λ ∪ Γ, where Γ is the set of constant functions. Let Ψ be the generated set of Ξ, i.e., the smallest closed functional algebra that contains Ξ. Then Ψ = Φ, i.e., Φ contains all impact functions as characterized by criteria I1–I4.

To prove it, we apply Bishop’s Theorem [70, 71]. We here only need the reduction to the special case of real-valued functions (as opposed to complex-valued functions):

Bishop’s Theorem. Let Z be a compact Hausdorff space. Let Ψ be a closed unital subalgebra of $$ C(Z,\mathbb{R})$$. Let $$ \varphi \in C(Z,\mathbb{R})$$. Suppose that ϕ|_S is constant for each subset S ∈ Z such that ψ|_S is constant for all ψ ∈ Ψ. Then ϕ ∈ Ψ.

The requirements of Bishop’s Theorem on Ψ are fulfilled since Z can be any sufficiently large compact subset of X × A and the inclusion of Γ ensures unitality. To show that the functional algebra Ψ contains all impact functions, we therefore need to show that for an arbitrary impact function ϕ for any $$ \mathbf{x},\widehat{\mathbf{x}}\in X$$ and $$ \mathbf{a},\widehat{\mathbf{a}}\in A$$:

We prove that Λ separates points for m = 1 with three lemmas, for which we transform the arguments to a canonical form (Definitions 1–4), in which populations are ordered by impact and maximally lumped together. We first show that if all functions from Λ have the same value for two arguments, these arguments have the same canonical form (Lemma 1). We then employ our criteria to show that no impact function will differ for two arguments that have the same canonical form (Lemma 2). Finally, we combine the first two lemmas to show that if all functions from Λ have the same value for two arguments, so do all impact functions (Lemma 3). Thus Λ separates points where it needs to (as per Eq 38).

Definition 1. Let x ∈ X and a ∈ A. Let $$ \mathfrak{A}=({\mathfrak{a}}_1,\dots ,{\mathfrak{a}}_s)$$ be the ordered sequence of non-zero values of a that correspond to a non-zero abundance, i.e.:

As the $$ \mathfrak{a}$$ are unique and ordered, we will directly use them like indices to avoid additional levels of indexing. One can think of them as equivalence classes of parameters.

Definition 2. Let x ∈ X and a ∈ A. For a given $$ \mathfrak{a}$$, let $$ I_{\mathfrak{a}}$$ be the set of indices where this parameter value is assumed and the corresponding abundance is not zero, i.e., the maximal set I such that $$ a_i=\mathfrak{a}$$ and x_i > 0 for all $$ i\in I_{\mathfrak{a}}$$. Consequentially, $$ I_{\mathfrak{b}}≔\left\{\right\}$$ for $$ \mathfrak{b}\notin \mathfrak{A}$$.

Definition 3. Let x ∈ X and a ∈ A. Denote the sums of abundances for one absolute parameter value as $$ \begin{array}{c}\hfill z_{\mathfrak{a}}≔{\sum }_{i\in I_{\mathfrak{a}}}{x}_i\end{array}$$.

Lemma 1. Suppose $$ \mathbf{x},\widehat{\mathbf{x}}\in X$$ and $$ \mathbf{a},\widehat{\mathbf{a}}\in A$$ are such that:

Then:

We show Eq 41 by induction over $$ \mathfrak{A}\cup \widehat{\mathfrak{A}}$$ in descending order of absolute value. We first note that the lemma trivially holds for all $$ \mathfrak{a}\in \left\{\right\}$$. In the following we show that, if the lemma holds for all $$ \mathfrak{b}$$ with $$ \left|\mathfrak{b}\right|>\mathfrak{a}>0$$, it also holds for $$ \mathfrak{a}$$ and $$ -\mathfrak{a}$$. (If one of $$ I_{\mathfrak{a}}\cup {\widehat{I}}_{\mathfrak{a}}$$ and $$ I_{-\mathfrak{a}}\cup {\widehat{I}}_{-\mathfrak{a}}$$ is empty, this does not affect this part of the proof.) To this end, we first show that the linear combinations must also be equal when only considering coefficients $$ \mathfrak{c}$$ with $$ \left|\mathfrak{c}\right|\le \mathfrak{a}$$ (for all p):

If $$ z_{\mathfrak{a}}+z_{-\mathfrak{a}}\ne 0$$ and $$ {\widehat{z}}_{\mathfrak{a}}+{\widehat{z}}_{-\mathfrak{a}}\ne 0$$ the above equality will be dominated by $$ {\mathfrak{a}}^p$$ for p → ∞, which gives us:

If exactly one of $$ z_{\mathfrak{a}}+z_{-\mathfrak{a}}$$ and $$ {\widehat{z}}_{\mathfrak{a}}+{\widehat{z}}_{-\mathfrak{a}}$$ were zero, the above limit would evaluate as either 0 or ∞ instead of 1; hence this cannot be. If both are zero, Eq 43 holds without further ado. Analogously, we obtain:

Definition 4 Define the canonical form of x ∈ X and a ∈ A as:

Lemma 2 Let x ∈ X and a ∈ A and ϕ be an impact function. Then $$ \varphi (\mathbf{x},\mathbf{a})=\varphi (\tilde{\mathbf{x}},\tilde{\mathbf{a}})$$.

We first transform blocks of arguments to the canonical form (with some zero arguments added if necessary) step by step, and show that the value of an impact function is not affected by these transformations. The first kind of block we consider are blocks of equal non-zero parameters and corresponding non-zero abundances, i.e., $$ I_{\mathfrak{a}}=:\{i_1,\dots ,i_v\}$$ for some $$ \mathfrak{a}$$. Then, for some $$ \check{\mathbf{x}},\check{\mathbf{a}}$$:

If a parameter a_i or abundance x_i, respectively, is zero, we transform the single-index block {i} to zero (for some $$ \check{\mathbf{x}},\check{\mathbf{a}}$$):

Second, after all blocks are transformed, we swap abundances and parameters in parallel to match the order in the canonical form. This does not affect the value of the impact function ϕ as it is commutative (I1).

Lemma 3 Suppose $$ \mathbf{x},\widehat{\mathbf{x}}\in X$$ and $$ \mathbf{a},\widehat{\mathbf{a}}\in A$$ are such that:

Let ϕ be an impact function. Then $$ \varphi (\mathbf{x},\mathbf{a})=\varphi (\widehat{\mathbf{x}},\widehat{\mathbf{a}})$$.

To prove this, we only need to note how the canonical forms $$ \tilde{\mathbf{x}}$$ and $$ \tilde{\mathbf{a}}$$ only depend on the parameters values $$ \mathfrak{a}$$ corresponding to non-zero total abundance $$ z_{\mathfrak{a}}$$ and these abundances. Those in turn are equal per Lemma 1. Thus:

Non-impact-function contribution to abundance changes must be proportional

Here, we show that $$ {\mathcal{R}}_j$$ as defined in Eq 10 must have the form $$ {\mathcal{R}}_j=x_j{\varphi }_j(\mathbf{x},\mathbf{a})$$, where ϕ_j is an impact function. To keep the notation simple, we assume that $$ {\mathcal{R}}_j$$ features no delay, noise, explicit time dependency, or similar, and thus $$ {\mathcal{R}}_j:X\times A\to \mathbb{R}$$.

We require every impact of population other than j to be comprised in an impact function ψ_j. We can write $$ {\mathcal{R}}_j$$ in the form:

Deriving a new model for UTI strains—The legwork

We here elaborate the details of creating a new model for the case described in Case Study.

Inserting our ansatz (Eq 24) into our first requirement (Eq 17), we obtain:

Assuming that the two factors do not “take turns” in being zero for different j, this means that either ρ_j(r_jj) = 0 or ς_j(s_jj) = 0. Without loss of generality, we assume that the latter applies, thus assigning ς_j the role of quantifying the carrying capacity. Furthermore, we choose ς_j(0) = 1. These are normalization choices, as they can be compensated by including a respective factor in ς_j or ρ_j respectively. Using this and expanding Eq 19, we obtain:

Assuming that ς_j is again responsible for the product being zero and it has only one root, namely 1, we arrive at: c_jk + s_jk v = 1, and thus: $$ \begin{array}{c}\hfill s_{jk}=\frac{1-c_{jk}}{v}\end{array}$$. Note that since s_jj = 1, this is consistent with our choice of c_jj = 1 − v (see S1 Appendix).

Using the above, we can expand Eqs 16 and 18:

We choose the arguably simplest function to fulfill the criteria for ρ, namely ρ_j(z) ≔ g_j + z. This has the consequence:

A group of functions fulfilling 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 effect of a occupied niche kicks in. Note that this choice results in terms similar to what was named hyperlogistic [72].

Finally, like the original study [16], we constrain the growth and capacity term to be non-negative to avoid the occasional implausible result. For example, we do not allow negative growth because we equate the abundance of a population with its footprint, which cannot be undone, and we lack the data to capture cell death. Putting everything together, we arrive at the model:

A problem with this model is that for 0 < x_k < 1, we have: $$ \begin{array}{c}\hfill {lim}_{c_{jk}\to 0}{\dot{x}}_j={lim}_{c_{jk}\to 0}{\mathcal{R}}_j(\mathbf{x},\mathbf{a})=\infty \end{array}$$. 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. Thus, one might argue that the actual point of the singularity requires a dedicated case distinction anyway. However, $$ \begin{array}{c}\hfill {lim}_{c_{jk}\to 0}{\dot{x}}_j=\infty .\end{array}$$ also means that $$ {\dot{x}}_j$$ becomes arbitrarily large for small c_jk. A way to address this problem is to consider the case q → ∞, or more specifically:

In our simulations, we therefore make a trade-off between complying with Eq 18 and the numerical benefits of a continuously differentiable model by setting q = 10 and approximating ς(1 − c_jk) ≈ lim_p→∞ ς(1 − c_jk) = 1 in Eq 58, thus arriving at: