Mutualistic cross-feeding in microbial systems generates bistability via an Allee effect.

In microbial ecosystems, species not only compete for common resources but may also display mutualistic interactions as a result from metabolic cross-feeding. Such mutualism can lead to bistability. Depending on the initial population sizes, species will either survive or go extinct. Various phenomenological models have been suggested to describe bistability in mutualistic systems. However, these models do not account for interaction mediators such as nutrients. In contrast, nutrient-explicit models do not provide an intuitive understanding of what causes bistability. Here, we reduce a theoretical nutrient-explicit model of two mutualistic cross-feeders in a chemostat, uncovering an explicit relation to a growth model with an Allee effect. We show that the dilution rate in the chemostat leads to bistability by turning a weak Allee effect into a strong Allee effect. This happens as long as there is more production than consumption of cross-fed nutrients. Thanks to the explicit relationship of the reduced model with the underlying experimental parameters, these results allow to predict the biological conditions that sustain or prevent the survival of mutualistic species.

between two different states of the system, in this case survival and extinction, is related to the concept of bistability. In mutualism, this is generated by an Allee effect [19][20][21][22][23] . An Allee effect is characterized by a decreased fitness at low densities so that the individual growth rate reaches a maximum at an intermediate density due to cooperative behavior. In the case of obligate mutualism, the Allee effect leads to a survival dependency as one species cannot survive without the other 22 . This effect is in contrast with the prediction of the classical logistic growth which predicts that an increased population density limits the growth 24 . A distinction between a weak and a strong Allee effect is made 1 . Whereas a weak Allee effect leads to a single stable state (the species always survives), a strong Allee effect is characterized by bistability, whereby a density threshold for survival is present.
Different models showing how mutualism causes bistability via the interdependence between the species have been proposed [25][26][27][28] . In some of these models the species are competitive as well as mutualistic, e.g. when mutualists become competitors for a resource or for available space at high density [29][30][31] . One study showed that the type of interaction could be modulated by varying the resource concentration 17 . Bistability is generated if the benefit of cooperation is counteracted by a cost, for example via the dispersal of the species. This increases the possibility of a collapse of the mutualistic system 23 . In order to predict when such disruptions occur and, if needed, to intervene to prevent the collapse of the community, a deeper understanding of mutualistic interactions and of the occurrence of thresholds in microbial communities is necessary.
A currently unresolved problem is that phenomenological models where interaction mediators, like nutrients, are neglected can behave differently than models where these are explicitly incorporated 32 . As both types of models describe similar phenomena, it should be possible to reduce the mechanistic model into a phenomenological model 33,34 . Such reductions of mechanistic models already exist for the growth of a single species 8,35,36 or for competitive consumer-resource models 37 . In the case of mutualism, one approach described the occurrence of bistability by the saturation of the mutualistic strength at high densities 38 . Nevertheless, it remains unclear how the occurrence of bistability in a mutualistic cross-feeding community is related to nutrient concentrations and to their consumption and production kinetics. This is essential to quantify the effects of prebiotics or biological parameters on the survival of the species.
In this theoretical work, we use a nutrient-explicit model for the growth of mutualistic species in a chemostat reactor and show how an Allee effect is created. This allows to interpret the effect of biological parameters on the dynamics in order to predict when bistability is created and to estimate the density threshold for survival. Nutrient-explicit models of a single species in the chemostat can be reduced to the logistic growth equation 8,35,36 . Using a similar approach, we reduce chemostat equations of a mutualistic system to an appropriate mechanistic model which only involves the species densities. This allows to relate the obtained equations to a generic growth model with an Allee effect. By establishing this analogy, we show that mutualism causes a weak Allee effect, which can be turned into a strong Allee effect under the influence of the dilution in the chemostat. Critical chemostat parameters such as the dilution rate and the influx of nutrients thus allow to manipulate the strength of the Allee effect and therefore of the survival threshold. As a consequence, it is possible to switch between regions of bistability, monostable survival, or monostable extinction. We also show that the production of cross-feeding nutrients needs to be larger than the consumption for an Allee effect to exist. This explicit relationship between experimental parameters and the Allee effect provides a way to bridge the gap between biological experiments and theoretical models.

Results
Bistability in mutualistic systems creates a survival threshold. We study a theoretical system of two mutualistic species with densities ρ 1 and ρ 2 . Mutualism is mediated by cross-feeding: each species consumes a nutrient, with resp. concentrations P 1 and P 2 , produced by the other species. We also assume that each species requires an additional nutrient which we refer to as the substrate, with resp. concentrations S 1 and S 2 . The substrate is a necessary component of the model to avoid a violation of mass conservation when production of the cross-fed nutrients is higher than the consumption. We choose to model the use of two different substrates as it is the least complex, symmetric case. Other scenarios like the exploitative competition for the same substrate for both species or the dependence of one substrate by one of the two species lead to cumbersome algebraic calculations but result in similar dynamics (Supplementary Material, Section S2).
A chemostat consists of a well-mixed growth vessel with an inflow of nutrients, at concentrations  S i and  P i ( = i 1, 2), and an outflow of the suspension. Both inflow and outflow occur at a dilution rate d. The consumption of nutrients is considered proportional to the growth of the species, which corresponds to the conservation of biomass when = d 0 (Supplementary Material S1). In the same way, we also assume the production of nutrients to be proportional to the growth of the species. The equations that describe the mutualistic chemostat system are then:  We focus on obligate cross-feeding. Therefore, we assume obligatory dependence of the nutrients which means that the growth rate of the species (F S P ( , ) ) needs to satisfy the following condition: . Different assumptions for the type of growth rates can be made. We will show that our results require the possibility to linearize the growth rate. A frequently used growth rate is the Monod function 39 , so that we use the following extension of the Monod function for two nutrients as an example 40 : The mutualistic relationship creates a positive feedback loop which can lead to bistability. This phenomenon becomes apparent if we simulate the behavior of the species when cultured in a chemostat reactor (Fig. 1). As an illustrative example, we consider what happens if the dilution is varied with steps of . 0 08 for consecutive time frames of = t 100. When the equilibrium state is reached, the growth rate of each species equals the dilution rate d, allowing to set the equilibrium point by adjusting the value of . This increase in the dilution rate leads to a decrease of the equilibrium density. However, a survival threshold exists: if the density drops below some critical values the species are washed-out. This is an irreversible collapse of the system, in the sense that reversing the dilution rate does not immediately lead to a return of the surviving state. Instead the species density continuous to drop. Survival is only restored by a sufficient decrease of the dilution rate. This concept is often called hysteresis 41 : the state of the system depends on its history. It is an effect related to bistability, meaning that there are two distinct equilibrium states for the same environmental parameters. In this case there is bistability between survival and extinction of the two species. Even though this simulation provides some intuition, it is not straightforward to predict how the different parameters of the system affect bistability. This is essential to determine when the threshold for survival is crossed. To provide an answer, we show how these parameters affect the dynamics by revealing a close connection to the widely studied growth equation with an Allee effect. Dilution allows to switch between a weak and a strong Allee effect. Mutualistic cross-feeding is a form of cooperative behavior and causes a positive feedback loop between the two species. As a consequence, at small densities an Allee effect is created. The per capita growth rate of a species, a proxy for their fitness, increases with the density at small densities and it reaches a maximum at an intermediate density. This is in contrast with the logistic growth, where the per capita growth rate becomes maximal near zero density. There is an important distinction between a weak Allee effect, associated to monostable dynamics, and a strong Allee effect, associated to bistability. Here, we show that the mutualistic cross-feeding creates a weak Allee effect and that increasing dilution in the chemostat is able to turn this into a strong Allee effect. Dilution can thus promote bistability: the two species coexist via cooperation or both die (the community collapses). The separation between the two types of behavior is determined by a threshold of the population size.
Based on the reduction of nutrient-explicit equations for the growth of one species in the chemostat to the logistic equation 8,35,36 (see Supplementary Material, Section S1), we can reduce the explicit chemostat equations to a two-variable mutualistic system. The calculations are detailed in the Supplementary Material (Section S2). In brief, this reduction relies on the following assumptions ( = i 1, 2 and ≠ i j):  Table 1. Definition of the variables and parameters of nutrient-explicit chemostat equations of the mutualistic cross-feeding system, Eq. (1), for = i 1, 2 and ≠ i j.
1. Conservation of mass for t ≫ 1/d: A Taylor approximation of the growth rate at low nutrient concentrations: The second assumption only holds in the case of obligatory dependence of the nutrients , such as the given example of the adapted Monod growth rate (Eq. (2)). Using a newly defined reduced set of parameters (see Table 2), we find that the two-species system can be described by the following equations: The microbial system we study consists of two species, characterized by their densities ρ 1 and ρ 2 , consumption of substrates (concentrations S 1 and S 2 ) and mutualistic cross-feeding via the production of nutrients (concentrations P 1 and P 2 ). We theoretically investigate the growth of these species in a chemostat reactor: a well-mixed growth vessel with an inflow of nutrients with concentrations  S i ,  P i (i 1, 2 = ) and an outflow of the suspension, occuring at an equal rate: the dilution rate d. (B) The behavior of the system is simulated with Eq. (1) using the growth rate Eq. (2). The dilution is gradually increased from = d 1 to = . d 0 24 with steps of . 0 08 for consecutive time frames of = t 100, using a logarithmic scale for the species density for clarity. The equilibrium densities of the species decrease with the dilution. For = .
d 0 24, a threshold is reached and the species will be washed-out. Decreasing the dilution to its previous rate does not lead to the recovery of the initial abundances of the species. The dilution needs to be further decreased to = .
d 0 08 the make the population growing again. This is a phenomenon called hysteresis: the system has memory of the previous state. It is a consequence of bistability between survival and extinction at intermediate dilution, so that a density threshold for survival exists. Initial These equations can be interpreted as follows. Each species (with density ρ i ) grows with a growth rate r i , which can be increased via mutualistic interactions (b i ). The available nutrients increase with the other species population. As this is the case for both species, this results in positive feedback: the growth of each species is positively influenced by its own density. Such growth is limited by two carrying capacities, one related to the available substrate (C si ) and one determined by the available cross-fed nutrients (C pi ). The exact biological mechanism behind the limitation through these carrying capacities is not critical for the dynamics, which remain qualitatively similar even when species compete for the same substrate (see Supplementary Material, Section S2). Additionally growth is further limited by the dilution rate in the chemostat (d).
In order to analyze the dynamics of these reduced Eqs. (4), we first study the symmetric case where all parameters of both species are the same (e.g.: = = r r r 1 2 ). In this situation the dynamics maps to the subspace where the densities of the species are equal: ρ ρ ρ = = 1 2 . The reduced mutualistic system is then found to be described by the following generic equation including cubic growth and an Allee effect 42 : Here, the growth rate r a , carrying capacity C a , Allee parameter a, and loss term caused by the dilution (parameter d) are defined in Table 3. The Allee effect is introduced by the factor ρ + a ( ) and is a consequence of the positive feedback between the species so that at low densities the per capita growth increases. Without dilution ( = d 1) this equation takes the usual form of a growth equation with an Allee effect 42 . Based on the parameter a, the following distinction is made: A weak Allee effect corresponds to monostable growth towards an equilibrium, while a strong Allee effect is associated to bistability between survival and extinction so that a densitiy threshold for survival is generated (Fig. 2B). As the Allee parameter is determined by the inflow of the cross-feeding nutrients and is therefore strictly positive, the mutualistic cross-feeding causes a weak Allee effect. The equilibria of the system are determined by the points where the growth is zero, which happens at ρ = 1 and ρ = C a when = d 1 (Fig. 2C). The extinction state ρ = 1 is unstable, while the coexisting state ρ = C a is stable. For the chemostat case, the dilution is switched on ( > d 0), the per capita growth is reduced by d and a new steady state may arise (Fig. 2D). The new steady state is unstable forming the threshold for survival between the two stable equilibria (survival and extinction). This way, the dilution rate can turn a weak Allee effect into a strong Allee effect.
Bifurcation analysis reveals regions of bistability. Before analyzing the biological consequences of the Allee effect in mutualistic interactions, it is useful to illustrate how the different parameters in Eq. (5) affect the behavior. This can be done via a bifurcation analysis, showing the stability of the different equilibria as a function of the parameters.
The equilibria of the system correspond to the intersections of the growth term and the dilution term, where ρ = d dt / 0, and are given by:   Here, ρ e is the extinction state, ρ C the coexisting state and ρ s the separating state between the two equilibria in the case of bistability. The equilibria are determined geometrically by calculating how many times the per capita growth intersects with the dilution term (Fig. 3D,E). When there is only one intersection, the only stable    5) shows the effect of bistability for different initial abundances ( t ( 0) ρ = varies from 0 and 1 with steps of 0.05). The species may coexist or become extinct. A threshold for survival separates the two outcomes. (C) The dilution rate is crucial for bistability. Without a dilution term, the system is monostable: there is only one stable equilibrium (ρ C ), corresponding to the survival of the species. Under the influence of the dilution, bistability occurs when there are three intersections of the growth and the dilution term: two stable fixed points corresponding to extinction (ρ e ) and survival (ρ C ) and an unstable steady state (saddle point, ρ s ) that separates both regions. (D) An Allee effect is present when the per capita growth increases with the density, which is the case for small densities and is the result of the positive feedback between the mutualistic species. For a small dilution rate, the per capita growth does not become negative for small densities, corresponding to a weak Allee effect. When the dilution rate is increased, the per capita growth rate becomes negative for small densities. This is a strong Allee effect and creates bistability between survival and extinction. www.nature.com/scientificreports www.nature.com/scientificreports/ equilibrium is ρ C (ρ s is negative and thus not physical). When there are two intersections, i.e. when ρ s and ρ C are both positive, the system is bistable: both the extinction state ρ e and the coexisting state ρ C are stable, and ρ s is a so-called saddle point that separates the two stable equilibria. When the dilution rate is too large, there are no intersections and the net growth is always negative so that the only equilibrium is ρ e , corresponding to extinction of the species. In this manner, it is straightforward to calculate the values of the dilution rate which determines the stability of the solutions. This can be represented in bifurcation diagrams, which show the different equilibria ρ e , ρ s and ρ C and their stability as a function of parameters like the dilution rate d (Fig. 3B) or the Allee parameter a (Fig. 3C). For the dilution rate d, the density ρ C decreases and the survival threshold is created at = d r C a a a , leading to bistability. This critical point is called a transcritical bifurcation 41 the coexistence state ρ C and the saddle point ρ s collapse and disappear, leaving ρ e as the only equilibrium state. This leads to the following condition for bistability: a a a a 2 For lower values of the dilution rate there is only survival and for higher values there is only extinction, thereby explaining the observed hysteresis in Fig. 1. The same analysis shows that the Allee parameter a has a counter-acting effect on the dynamics in comparison to d (Fig. 3A,C,E): when the dilution rate is too large for survival, increasing a causes the per capita growth to intersect with the dilution, thereby allowing survival of the community. The complete behavior is summarized in Fig. 3A: for large values of a and low values of d there is always survival via cooperation (regime I), for intermediate values of a and d there is bistability (regime II) and for low values of a and large values of d there is always extinction of the community (regime III).
In order to check the validity of the used simplification of the growth rates, we simulated the original chemostat equations, Eq. (1), with extended Monod growth rates, Eq. (2) (Fig. 4). The used Taylor approximation of the Monod growth rates in the simplified model corresponds to a linear dependence of the nutrients. Both are monotonically increasing functions, but the Monod growth rate yields a smaller value than its linearization. Therefore, we expect extinction to occur at lower dilution than predicted with the simplified model (Fig. 4A). We quantified the differences in the locations of the critical dilution rate ⁎ d where the saddle-node bifurcations (SN) occurs for both the full model and the reduced model (Fig. 4B). We performed the same calculation for different values of www.nature.com/scientificreports www.nature.com/scientificreports/ the Monod constants, whereby we multiplied the Monod constants by a factor F and the maximal growth rate µ by F 2 so that the Taylor approximation of the growth rate remains the same. The bifurcation diagrams for different values of F shows the qualitative correspondence to the bifurcation diagram related to Eq. (5), whereby the density decreases with the dilution and the regime changes from monostable survival to bistability via a transcritical bifurcation and to extinction via a saddle-node bifurcation (Fig. 2B). For = F 1 the similarity is only qualitative, while a better quantitative correspondence is obtained for increased values of F. The dynamic regimes are quantified for = F 1 (Fig. 4C) as a function of the Allee parameter a, highlighting the qualitative correspondence and a lower critical dilution rate ⁎ d than estimated. This shows how the regimes of Eq. (5) map unto the real dynamic regimes. The resemblance is only qualitative. For = F 1 a large error is observed. We estimate the error of the simplification by considering the critical dilution rate ⁎ d for different values of F. Multiplying the Monod constants by a factor = .
F 3 3 significantly decreases the error. A factor = F 10 yields a negligible error, i.e. a good quantitative match of the bifurcation diagrams was obtained with the explicit model and the simplified model (Fig. 4B). By using the extended Monod function (Eq. 2) as an example, we show that two conditions on the experimental growth rates should be fulfilled to expect Eq. (5) to be a valid approximation. First, it should be possible to linearize the growth rate for low nutrient concentrations (  S K). Second, the function should be monotonically increasing. As long as these conditions are fulfilled, the dynamics is qualitatively described by Eq. (5).
Bistability requires a sufficient production of cross-feeding nutrients. Bistability results from the mutualistic relationship between the two species. Having a dilution rate d in the correct parameter range is, by itself, not sufficient to guarantee bistability. Another necessary condition for bistability is obtained by considering the effect of the mutualistic interaction strength b, defined as the ratio of the production to the consumption of cross-fed nutrients = = . This parameter affects the growth rate r a as well as the Allee parameter a (Table 3), so that its overall effect on the dynamics is not straightforward to predict intuitively. Therefore, we simulated the mutualistic system (Eq. (4)) with symmetric parameters, for different values of b and of dilution rate d (Fig. 5A). This analysis shows that the region of bistability is limited to > b 1. This means that bistability only occurs when the production of cross-fed nutrients is higher than the consumption. When < b 1, the production is used to quantify the deviation from the simplified bifurcation curve of Eq. (5). At = F 1 the error is large as K s and K p are of the order of S and P (Table S5), the error is significantly reduced for = .  Table S5).

Scientific RepoRtS |
(2020) 10:7763 | https://doi.org/10.1038/s41598-020-63772-4 www.nature.com/scientificreports www.nature.com/scientificreports/ of cross-fed nutrients is insufficient. The growth is limited by these cross-fed nutrients, and the dynamics becomes similar as in the case of logistic growth, so that there is only one equilibrium state. The carrying capacity is then determined by the inflow of cross-feeding nutrients, via C p , and not by the inflow of substrate, via C s . Thus, for < b 1 there is no Allee effect.
The difference between the logistic growth at < b 1 and the Allee effect at > b 1 is visually interpreted by representing the physical growth (i.e. feasible) regions in the phase plane for < b 1 (Fig. 5B) and for > b 1 (Fig. 5C). These feasible regions originate from the mass conservation laws (Eq. (3)), as the nutrient concentrations need to be positive. For each species, the limitation by the substrate causes the growth region to be bounded by the following function (i 1, 2 = ): i s Similarly for the cross-feeding nutrients the following limitation functions are obtained ( = i 1, 2 and ≠ i j): The interaction strength b determines the slope of Eq. (9), so that for < b 1, the growth is entirely limited by Eq. (9) and thus by the availability of cross-feeding nutrients (Fig. 5B). The species cannot grow sufficiently to become limited by the substrate. In contrast, for > b 1, the species can grow sufficiently. In this case, the feasible region is limited at high densities by the availability of the substrate, Eq. (8) (Fig. 5C).
The overall influence of the parameters is understood by looking at the nullclines: the functions where ). The intersections of the nullclines define the equilibria of the system. The nullclines are hyperbolic functions, with the limiting functions (Eqs. (8) and (9)) as asymptotes. For bistability, the slopes need to be such that the nullclines can intersect twice (Fig. 5C), thereby forming a stable equilibrium corresponding to coexistence of the species and a saddle point which creates the threshold between coexistence and survival. This only occurs when > b 1, so that in this region we observe the same dynamics as described in last section: survival when the hyperbola intersect once at low dilution, bistability when the nullclines intersect twice for intermediate  4). For < b 1, the growth is limited by the crossfeeding nutrient, leading to monostable dynamics similar to logistic growth. There is survival for small dilution ( < d r C a a a ) and extinction for large dilution. For > b 1, the growth only becomes limited by the substrate, allowing for bistability. There is a weak Allee effect, corresponding to monostable survival, for small dilution, a strong Allee effect for intermediate dilution ( ) and monostable extinction for large dilution. (B) For < b 1, the growth is limited by the cross-feeding nutrients (Eq. (9)), so that the physical growth region is small and does not allow bistability. For small dilution the dynamics is similar to logistic growth as the population monotonically grows towards its equilibrium. (C) For > b 1, high population densities are possible, whereby the growth is limited by the substrate S (Eq. (8)). Increased values of b lead to a larger growth region, allowing for bistability when > b 1. Increasing the dilution rate d causes the nullclines to bend into hyperbola with the linear functions at = d 1 as asymptotes. Bistability is obtained when the hyperbolic nullclines intersect twice parameter values as listed in Table S6).  Fig. S3). Bistability is created in the same way when asymmetric parameters are used. A similar analysis was performed for asymmetric values of the mutualistic interaction strength b (Supplementary Material, Fig. S5 and Section S3), resulting in the general necessary condition for bistability: > b b 1 1 2 . Finally, the choice to model two substrates, as it is the least complex case, does not affect this result. When two mutualistic species compete for a substrate, we found that the qualitative effect of the dilution rate and the mutualistic strength remains the same ( Supplementary  Material, S4)). Returning to the chemostat parameters, the obtained condition for bistability > b b 1 1 2 corresponds to ν ν > a a p p 1 2 1 2 , which means that the overall production of cross-fed nutrients needs to be larger than their consumption.

Discussion
Microbial interaction networks are characterized by multiple positive and negative interactions. Species enter in competition for limited resources, but they can also display mutualistic relationships through cross-feeding. Through mutualistic interactions, both species benefit of each others presence. This may be seen as a stabilizing factor. However, mutualism carries the seed of its own instability: under the influence of dilution bistability may occur, causing a critical density threshold for the species to survive. Once the abundances of the species drop below this threshold, the community eventually collapses and the species become extinct.
How biological parameters affect the survival threshold is often unclear. To provide an understanding of the effect of different parameters, we showed how nutrient-explicit equations for two mutualistic cross-feeding species can be reduced to a set of equations which only involve the densities of the species. These could then be related to a growth equation with an Allee effect, which can be analyzed to obtain a deeper understanding of the impact of the different biological parameters.
Density thresholds for survival have previously been observed in mutualistic systems. For example, a survival threshold was found in a spatially cooperating microbial community 18 and in a cross-feeding system which could switch to other interactions like competition and parasitism, depending on the availability of nutrients 17 . Besides cross-feeding, mutualism can also arise via the protection of another species towards antibiotics. In such a cross-protection system it was found that periodic dilution drives oscillatory dynamics, potentially leading to extinction if the survival threshold was crossed 43 .
Our results showed that the overall production rate of cross-fed nutrients needs to be larger than the overall consumption rate to create an Allee effect. The production and consumption rates can experimentally be altered by making use of synthetic cross-feeding systems 17,44,45 , which can be designed to display bistability, or on the contrary, to prevent bistability and thereby the risk of an abrupt collapse of the community. We obtained quantitative results for the case where both species have symmetric parameter values, but showed that this framework still applies to the case of asymmetric parameter values.
By using time traces of microbial growth experiments it is possible to fit the experimental values of biological parameters 46 . This would allow the validation of these theoretical predictions if the microbes are obligate cross-feeders. If the system exhibits a sufficient production of cross-fed nutrients ( ν ν > a a p p 1 2 1 2 ), then bistability is predicted for dilution rates within the range < < + r C a d r C a ( )/4 a a a a 2 . The dilution rate and the influx of nutrients are experimental parameters which can be tuned so that these allow to manipulate the behavior. As a result, this provides experimental guidelines to study the presence of bistability in a microbial system. Moreover, the influx of nutrients (  S and  P ) allow to study the effect of prebiotics, in order to determine the effect of these growth-promoting nutrients on the survival of the species. Furthermore, if antibiotics are used, then determining whether a survival threshold exists would be essential to predict the survival or extinction of the species.
The presence of bistability can significantly affect microbial behavior. For example, the existence of a survival threshold creates susceptibility to cheaters 28 . A cheater is an individual of the species which does not cooperate in the production of cross-feeding nutrients, thereby increasing the risk of a collapse of the system 47 . Addition of a third species can create global stability of the coexistence state if it is a facultative mutualist, providing a solution for the presence of the survival threshold 48 .
Bistability also affects the behavior in the case of spatial expansions. Here, an Allee effect counteracts genetic drift of a species as it creates a pushed wave rather than a pulled wave corresponding to logistic growth 49,50 . Such behavior has been observed experimentally in a system of two cross-feeding species whereby the mutualistic strength was modulated by the inflow of nutrients 51 .
Mutualistic species are usually part of a larger ecosystem. Therefore, when mutualists are decreased under the survival threshold, for example in response to antibiotics, the entire ecosystem can be destabilized 10 . Such critical effects on ecological interactions are often not well characterized 52 . Generalized Lotka-Volterra models are often used to study the interactions in microbial ecosystems [53][54][55] . However, Lotka-Volterra models for mutualism involve the addition of fitness effects 32 so that nonlinear growth rates are not incorporated 56 .
We studied a theoretical microbial system that is mutualistic via cross-feeding, which is a necessary step to obtain a deeper understanding of complex behavior in microbial communities. The obtained phase plane and the nullclines, describing the dynamics of the species, are observed in different mutualistic models [25][26][27][28] , so that we can state this is a general phenomenon. Our results give insights into necessary conditions to obtain bistability in a model for microbial species which are obligate cross-feeders: there needs to be an Allee effect as well as a limiting function. In particular, we show that nonlinear growth rates significantly affect the dynamics of microbial communities.