Ecological feedback in quorum-sensing microbial populations can induce heterogeneous production of autoinducers

Autoinducers are small signaling molecules that mediate intercellular communication in microbial populations and trigger coordinated gene expression via ‘quorum sensing’. Elucidating the mechanisms that control autoinducer production is, thus, pertinent to understanding collective microbial behavior, such as virulence and bioluminescence. Recent experiments have shown a heterogeneous promoter activity of autoinducer synthase genes, suggesting that some of the isogenic cells in a population might produce autoinducers, whereas others might not. However, the mechanism underlying this phenotypic heterogeneity in quorum-sensing microbial populations has remained elusive. In our theoretical model, cells synthesize and secrete autoinducers into the environment, up-regulate their production in this self-shaped environment, and non-producers replicate faster than producers. We show that the coupling between ecological and population dynamics through quorum sensing can induce phenotypic heterogeneity in microbial populations, suggesting an alternative mechanism to stochastic gene expression in bistable gene regulatory circuits. DOI: http://dx.doi.org/10.7554/eLife.25773.001


Introduction
Autoinducers are small molecules that are produced by microbes, secreted into the environment, and sensed by the cells in the population (Keller and Surette, 2006;Hense and Schuster, 2015). Autoinducers can trigger a collective behavior of all cells in a population, which is called quorum sensing. For example, quorum sensing regulates the transcription of virulence genes in the Grampositive bacterium Listeria monocytogenes (Gray et al., 2006;Garmyn et al., 2011;da Silva and De Martinis, 2013) and the transcription of bioluminescence genes in the Gram-negative bacterium Vibrio harveyi (Xavier and Bassler, 2003;Anetzberger et al., 2009), and it may also autoregulate the transcription of autoinducer synthase genes (Fuqua and Greenberg, 2002;Waters and Bassler, 2005). When the concentration of autoinducers reaches a threshold value, a coordinated and homogeneous expression of target genes may be initiated in all cells of the population (Waters and Bassler, 2005;Hense and Schuster, 2015;Papenfort and Bassler, 2016), or a heterogeneous gene expression in the population may be triggered at low concentrations (Anetzberger et al., 2009;Williams et al., 2008;Boedicker et al., 2009;Garmyn et al., 2011;Pérez and Hagen, 2010;Ackermann, 2015;Grote et al., 2014, Grote et al., 2015Papenfort and Bassler, 2016;Pradhan and Chatterjee, 2014). To implement all of these functions and behaviors, a microbial population needs to dynamically self-regulate the average autoinducer production.
Within a given population, the promoter activity of autoinducer synthase genes may vary between genetically identical cells (Garmyn et al., 2011;Grote et al., 2014;Anetzberger et al., 2012;Plener et al., 2015;Cárcamo-Oyarce et al., 2015;Grote et al., 2015). For example, during the growth of L. monocytogenes under well-mixed conditions two subpopulations were observed, one of which expressed autoinducer synthase genes, while the other did not (Garmyn et al., 2011). Such a phenotypic heterogeneity was associated with biofilm formation (Garmyn et al., 2011;da Silva and De Martinis, 2013;Hense and Schuster, 2015;Cárcamo-Oyarce et al., 2015). The stable coexistence of different phenotypes in one population may serve the division of labor or act as a bet-hedging strategy and, thus, may be beneficial for the survival and resilience of a microbial species on long time scales (Ackermann, 2015).
The mechanism by which a heterogeneous expression of autoinducer synthase genes is established when their expression is autoregulated by quorum sensing has remained elusive. For example, expression of the above mentioned autoinducer synthase genes in L. monocytogenes is up-regulated through quorum-sensing in single cells (Garmyn et al., 2011, Garmyn et al., 2009Waters and Bassler, 2005). From an experimental point of view it is often not known, however, whether autoinducer synthesis is up-regulated for all autoinducer levels or only above a threshold level. To explain phenotypic heterogeneity of autoinducer production, currently favored threshold models of quorum sensing typically assume a bistable gene regulation function (Fujimoto and Sawai, 2013;Pérez-Velázquez et al., 2016;Goryachev et al., 2005;Dockery and Keener, 2001). For bistable regulation, cellular autoinducer synthesis is up-regulated above a threshold value of the autoinducer concentration in the population, whereas it is down-regulated below the threshold ('allor-none' expression); see Figure 1B. Stochastic gene expression at the cellular level then explains the coexistence of different phenotypes in one population. If, however, cellular autoinducer synthesis is up-regulated for all autoinducer concentrations (monostable up-regulation), the mechanism by which phenotypic heterogeneity can arise and is controlled has not been explained.
Here we show that the coupling between ecological and population dynamics through quorum sensing can control a heterogeneous production of autoinducers in quorum-sensing microbial eLife digest Bacteria and other microbes can communicate with each other using chemical languages. They release small signaling molecules called autoinducers into their surroundings and sense the levels of the autoinducers in the environment. The response to these autoinducersknown as quorum sensing -can regulate how whole communities of microbes grow and behave; for example, autoinducers can alter the ability of microbes to infect humans or enable the microbes to collectively switch on light production.
Recent experiments suggest that, in a population of genetically identical microbes, some individuals may produce autoinducers while others do not. The coexistence of these different "phenotypes" in one population may enable different individuals to perform different roles, or act as a "bet-hedging" strategy that helps the population to survive if it is later exposed to a stressful situation.
It is not clear how microbes regulate autoinducer production so that only some individuals produce these molecules. Bauer, Knebel et al. developed a theoretical model to address this question. In the model, the microbes shape their environment by producing autoinducers and can respond to this self-shaped environment by changing their level of autoinducer production. Bauer, Knebel et al. found that this establishes a feedback loop that can result in autoinducers being produced by some individuals but not others.
The next step following on from this work is to carry out experiments to test the assumptions and predictions made by the theoretical model. These findings may help to understand how the coexistence of different phenotypes affects collective behaviors, and vice versa, in populations of microbes that use quorum-sensing.
populations. At the same time, the overall autoinducer level in the environment is robustly self-regulated, so that further quorum-sensing functions such as virulence or bioluminescence can be triggered. We studied the collective behavior of a stochastic many-particle model of quorum sensing, in which cells produce autoinducers to different degrees and secrete them into the well-mixed environment. Production of large autoinducer molecules (for example oligopeptides) and accompanied gene expression are assumed to reduce fitness such that non-producers reproduce faster than producing cells. Moreover, it is assumed that quorum sensing enables up-regulation of autoinducer production, that is, individuals can increase their production in response to the sensed average production level in the population ( Figure 1). As a central result, we found that the population may split into two subpopulations: one with a low, and a second with a high production rate of autoinducers. This phenotypic heterogeneity in the autoinducer production is stable for many generations and the autoinducer concentration in the population is tightly controlled by how production is upregulated. If cellular response to the environment is absent or too frequent, phase transitions occur from heterogeneous to homogeneous populations in which all individuals produce autoinducers to the same degree. To capture these emergent dynamics, we derived the macroscopic mean-field equation (1) from the microscopic stochastic many-particle process in the spirit of the kinetic theory in statistical physics, which we refer to as the autoinducer equation. The analysis of the autoinducer equation explains both phenotypic heterogeneity through quorum sensing and the phase transitions to homogeneity.
The key aspect of our work is how the composition of a population changes in time when its constituents respond to an environment that is being shaped by their own activities (see Box 1). This ecological feedback is mediated by quorum sensing and creates an effective global coupling between the individuals in the population. Such a global coupling is reminiscent of long-range interactions in models of statistical mechanics, such as in the classical XY spin model with infinite range interactions (Antoni and Ruffo, 1995;Yamaguchi et al., 2004;Barré et al., 2002;Choi and Choi, 2003;de Buyl et al., 2010;Campa et al., 2009;Pakter and Levin, 2013). Our analysis suggests that quorum sensing in microbial populations can induce and control phenotypic heterogeneity as a collective behavior through such a global coupling and, notably, does not rely on a bistable gene regulatory circuit (see Discussion). or or reproduction before after response no response

Response function
Population average bistable monostable Figure 1. The quorum-sensing model for the production of autoinducers in microbial populations. (A) Sketch of a typical update step. Individuals are depicted as disks and the degree of autoinducer production (p i 2 ½0; 1) is indicated by the size of the green fraction. Non-producers (orange disks) reproduce fastest, full producers (green disks) slowest. Individual i with p i ¼ 1=6 divides into two offspring individuals, one of which replaces another individual j. Both offspring individuals sense the average production level in the population (hpi ¼ 1=3), and may either respond to this environment, with probability l, by adopting the value RðhpiÞ of the response function (¼ 2=3 here, see (B)) or, with probability 1 À l, retain the production degree from the ancestor (¼ 1=6). Here, offspring individual i responds to the environment while j does not (denoted by gray shading). (B) Quorum sensing is characterized by the response function. Perception of the average production level in the population (hpi) enables individuals to change their production degree to the value RðhpiÞ 2 ½0; 1. Sketched are a monostable response function (stable fixed point at 1, unstable fixed point at 0), and a bistable response function (stable fixed points at 0 and 1, unstable fixed point at an intermediate threshold value). Stable fixed points of the response function are depicted as black circles while unstable fixed points are colored in white. For the sketched bistable response function, autoinducer production is down-regulated with respect to the sensed production level in the population below the threshold value, and up-regulated above this threshold. For the monostable response function, autoinducer production is up-regulated at all sensed production levels. DOI: 10.7554/eLife.25773.003 Set-up of the quorum-sensing model We now introduce the quorum-sensing model for a well-mixed population of N individuals ( Figure 1). The phenotype of each individual i ¼ 1; . . . ; N is characterized by its production degree p i 2 ½0; 1, that is, the extent to which it produces and secretes autoinducers. In an experiment with microbes, the promoter activity of autoinducer synthase genes or their enzymatic activity could be a proxy for the production degree. The limiting case p i ¼ 0 denotes a non-producer, and p i ¼ 1 denotes a full producer. The state of the population p ¼ ðp 1 ; . . . ; p N Þ changes stochastically ( Figure 1A): An individual i reproduces with rate f i , which we refer to as the individual's fitness. We assume that fitness decreases with incurring metabolic costs of induction and synthesis of autoinducers, and with other metabolic burdens in the cell's phenotypic state (Ruparell et al., 2016;Diggle et al., 2007;He et al., 2003). For simplicity, we choose f i ¼ fðp i Þ ¼ 1 À sp i . The selection strength 0 s < 1 scales the fitness difference with respect to the non-producing phenotype (fð0Þ ¼ 1). Thus, the larger an individual's production, the smaller its reproduction rate. This assumption is discussed in detail further below (see Discussion).
Whenever an individual divides into two offspring individuals in the stochastic process, another individual from the population is selected at random to die such that the population size N remains constant. Qualitative results of our model remain valid if only the average population size is constant, which may be assumed, for example, for the stationary phase of microbial growth in batch culture. One recovers the mathematical set-up of frequency-dependent Moran models for Darwinian selection (Moran, 1958;Ewens, 2004;Blythe and McKane, 2007;Nowak et al., 2004) if one restricts the production degrees to a discrete set, for example, to full producers or non-producers only, p i 2 f0; 1g. The mathematical set-up of the well-known Prisoner's dilemma in evolutionary game theory is recovered if, in addition, the secreted molecules would confer a fitness benefit on the population (Nowak et al., 2004;Traulsen et al., 2005;Melbinger et al., 2010;Assaf et al., 2013). Since we are interested in the mechanism by which heterogeneous production of autoinducers might be induced and do not study the context under which it might have evolved, we do not include any fitness benefits through signaling, for example at the population level, into the modeling here (see Discussion).
A central feature of our model is the fact that individuals may adjust their production degree via a sense-and-response mechanism through quorum sensing, which is implemented as follows. After reproduction, both offspring individuals sense the average production level of autoinducers hpi ¼ 1=N P i p i in the well-mixed population. With probability l, they independently adopt the value RðhpiÞ 2 ½0; 1 as their production degree in response to the sensed environmental cue hpi, whereas they retain the ancestor's production degree with probability 1 À l through non-genetic inheritance. In an experimental setting, the response probability l relates to the rate with which cells respond to the environment (Kussell and Leibler, 2005;Acar et al., 2008;Axelrod et al., 2015) and regulate their production through quorum sensing. We refer to the function RðhpiÞ as the response function, which is the same for all individuals. The response function encapsulates all biochemical steps involved in the autoinducer production between perception of the average production level hpi and adjustment of the individual production degree to RðhpiÞ in response (He et al., 2003;Williams et al., 2008;Drees et al., 2014;Hense and Schuster, 2015;Maire and Youk, 2015); see Figure 1B. For example, it may be a bistable step or bistable Hill function, which is often effectively assumed in threshold models of phenotypic heterogeneity (Fujimoto and Sawai, 2013;Pérez-Velázquez et al., 2016;Goryachev et al., 2005;Dockery and Keener, 2001). For a bistable response function, cellular production is up-regulated above a threshold value of hpi, whereas it is down-regulated below the threshold. For the bistable response function sketched in Figure 1B, both values hpi ¼ 0 and hpi ¼ 1 are stable fixed points. In this work, however, we particularly focus on monostable response functions RðhpiÞ to model microbial quorum-sensing systems in which autoinducer synthesis is up-regulated at all autoinducer production levels in the population (Garmyn et al., 2009;Waters and Bassler, 2005). In other words, cellular production always increases with respect to the sensed production level in the population (stable fixed point at hpi ¼ 1 and unstable fixed point at hpi ¼ 0). The sense-and-response mechanism is further discussed in the Discussion section.
From a mathematical point of view, the introduced sense-and-response mechanism through quorum sensing constitutes a source of innovation in the space of production degrees because an individual may adopt a production degree that was not previously present in the population. Thus, a continuous production space with p i 2 ½0; 1 as opposed to a discrete production space is a technical necessity for the implementation of the quorum-sensing model. The coupling of ecological dynamics (given by the average production level of autoinducers hpi) with population dynamics (determined by fitness differences between the phenotypes) through quorum sensing results in interesting collective behavior, as we show next. We emphasize that, as long as this coupling is present, the effects of the quorum-sensing model that we found and report next are qualitatively robust against noise at all steps; see below.
Box 1. An ecological feedback can control phenotypic heterogeneity in quorum-sensing microbial populations.
Our work demonstrates that the coupling of ecological and population dynamics through quorum sensing cannot only lead to homogeneously producing populations, but can also control a heterogeneous production of autoinducers in microbial populations. Phenotypic heterogeneity becomes manifest in the quorum-sensing model as long-lived, bimodal states of the population that are dynamically stable; see sketch below and Equation (2). In the quorum-sensing model, ecological dynamics are determined by the average production level of autoinducers, while population dynamical changes are determined by fitness differences between non-producers and producers of autoinducers. Because individuals sense and respond to autoinducers in the environment, the ecological dynamics are coupled with the population dynamics. In other words, an ecological feedback loop is established when cells respond to an environment that is being shaped by their own activities. When fitness differences between nonproducers and producers of autoinducers balance with cellular response to autoinducers in the environment, separated production degrees stably coexist in one population. Therefore, we expect that a heterogeneous production of autoinducers may be induced and controlled by such an ecological feedback in real microbial populations, suggesting an alternative mechanism to stochastic gene expression in bistable gene-regulatory circuits to control phenotypic heterogeneity (see Discussion).
Box 1-figure 1. Effective picture of robust phenotypic heterogeneity through an ecological feedback. The coupling of fitness differences between non-producers and producers (selection strength s) and sense-andresponse to the self-shaped environment through quorum sensing (response probability l and up-regulation of production with response function RðhpiÞ) ensures the stable coexistence of the two subpopulations at the phenotypic states p low and p high ; see Equation (2). The value b ¼ 2l=s quantifies this coexistence. In one subpopulation (fraction y ¼ 1 À b=RðbÞ of the total population), individuals do not produce (p low ¼ 0), while in the other (fraction 1 À y) individuals produce autoinducers to the degree p high ¼ RðbÞ. The average production level in the population is robustly adjusted to the value hpi ¼ b. States of phenotypic heterogeneity arise for a broad range of initial distributions and are robust against noisy inheritance, noisy perception, and noisy response (see Results of mathematical analysis and Appendix 1-figures 1 and 2). DOI: 10.7554/eLife.25773.005

Results of numerical simulations
The quorum-sensing model was numerically simulated by employing Gillespie's stochastic simulation algorithm (Gillespie, 1976(Gillespie, , 1977 for a population size of N ¼ 10 4 individuals and an exemplary selection strength s ¼ 0:2, such that sN ) 1. In this regime, demographic fluctuations are subordinate (Nowak et al., 2004;Wild and Traulsen, 2007;Blythe and McKane, 2007). Within the scope of our quorum-sensing model, the precise value of the selection strength s that scales the fitness differences is not important for the reported mechanism by which phenotypic heterogeneity can be induced, see below. We tracked the state of the population p over time, and depict the histogram of production degrees and the population average in Figure 2.
First, we studied the stochastic many-particle process without sense-and-response (l ¼ 0); see Figure 2A,D and Video 1. In this case, non-producers always proliferate because they reproduce at the highest rate in the population, which is well-studied in evolutionary game theory (Taylor and  Homogeneous and heterogeneous production of autoinducers in the quorum-sensing model. Temporal evolution of autoinducer production in the quorum-sensing model depicted as histograms of production degrees (normalized values), (A-C); and average production level of autoinducers in the population (D-F); see also Videos 1-3. (A) In the absence of sense-and-response (l ¼ 0), only non-producers proliferate. The approach to stationarity is asymptotically algebraically slow for a quasi-continuous initial distribution of production degrees (D). The black line hpi~t À1 serves as a guide for the eye. (B) Sense-and-response through quorum sensing (l ¼ 0:2 here) promotes autoinducer production, and the population becomes homogeneous (ultimately, fixation at a single production degree, data not shown). The response function used here, RðhpiÞ ¼ hpi þ 0:2 Á sin ðphpiÞ, was chosen such that an individual's production degree is always up-regulated through quorum sensing (see Figure 1B). Approach to stationarity is exponentially fast (E), but timescales may diverge at bifurcations of the response function (see Appendix 1- figure 3). The dashed line in (E) shows fit to an exponential decay. (C) When l is small (l ¼ 0:05 here), the population becomes heterogeneous: quasi-stationary states arise in which the population splits into two subpopulations, one of which does not produce autoinducers, while the other does. The same monostable response function was chosen as in (B). Therefore, heterogeneity may arise without bistable response. For very long times, one of the two absorbing states (A, B) is reached, data not shown (see Figure 3A). Heterogeneous, quasi-stationary states arise for a broad class of initial distributions (see Appendix 1-figure 1 and our mathematical analysis). At the same time, the average production level of autoinducers in the population is adjusted by the response probability l if s is fixed (F) or vice versa (data not shown). Bimodal, quasi-stationary states also arise when noisy inheritance, noisy perception, and noisy response are included in the model set-up (see Appendix 1-figure 2). Mean-field theory agrees with all observations (autoinducer equation (1)). The time unit Dt ¼ 1 means that in a population consisting solely of non-producers, each individual will have reproduced once on average. Ensemble size M ¼ 100, s ¼ 0:2, N ¼ 10 4 . DOI: 10.7554/eLife.25773.006 The following source data is available for figure 2: Source data 1. Source data accompanying  , 1978;Maynard Smith, 1982;Hofbauer and Sigmund, 1998). Thus, the initially uniform distribution in the population shifts to a peaked distribution at low production degrees. Ultimately, a homogeneous (unimodal) stationary state is reached in which all individuals produce autoinducers to the same low degree p low ' 0. Such a stationary state is absorbing (Hinrichsen, 2000), that is, the stochastic process offers no possibility of escape from this state of the population.

Jonker
With quorum sensing (l > 0), absorbing states are reached if, again, all individuals produce to the same degree p Ã and, in addition, the value of this production degree is a fixed point of the response function (Rðp Ã Þ ¼ p Ã ); see Figure 2B,E and Video 2. In such a homogeneous absorbing state with hpi ¥ ¼ p Ã , an offspring individual can no longer alter its production degree. It either takes over the production degree p Ã from its ancestor or it adopts that same degree Rðhpi ¥ Þ ¼ hpi ¥ ¼ p Ã through sense-and-response. Thus, all individuals continue to produce with degree p Ã and the state of the population remains homogeneous (unimodal).
Surprisingly, for small response probabilities l, we found that the population may get trapped in heterogeneous (bimodal) states for long times before a homogeneous absorbing state is reached. The temporal evolution of such a heterogeneous state is shown in Figure 2C,F and Video 3 for l ¼ 0:05. A monostable response function was chosen with RðhpiÞ > hpi for all hpi 2 ð0; 1Þ (unstable fixed point at 0, and stable fixed point at 1) such that the production degree is always up-regulated through quorum sensing; see sketch in Figure 1B. After some time has elapsed, the population is composed of two subpopulations: one in which individuals produce autoinducers to a low degree p low , and a second in which individuals produce to a higher degree p high that is separated from p low by a gap in the space of production degrees. Only through strong demographic fluctuations can the population reach one of the homogeneous absorbing states (hpi ¥ ¼ 0 or 1 for the response function chosen above). The time taken to reach a homogeneous absorbing state grows exponentially with N ( Figure 3A). Therefore, states of phenotypic heterogeneity are quasi-stationary and long-lived. These heterogeneous states arise for a broad class of response Video 1. Video accompanying Figure 2A -Homogeneous production of autoinducers in the population if sense-and-response is absent in the quorum-sensing model (l ¼ 0). DOI: 10.7554/eLife.25773.008 Video 2. Video accompanying Figure 2B -Homogeneous production of autoinducers in the population if sense-and-response is frequent in the quorum-sensing model (l ¼ 0:2 here). DOI: 10.7554/eLife.25773.009 Video 3. Video accompanying Figure 2C -Heterogeneous production of autoinducers in the population if sense-and-response is rare in the quorumsensing model (l ¼ 0:05 here). DOI: 10.7554/eLife.25773.010 functions and initial distributions (Appendix 1-figure 1), and they are robust against demographic noise that is always present in populations of finite size ( Figure 3A); see our mathematical analysis below. We demonstrated that states of phenotypic heterogeneity are also robust against changes of the model set-up, which might account for more biological details (see, for example, Papenfort and Bassler, 2016 and references therein). Upon including, for example, noisy inheritance of the production degree, noisy perception of the environment, and noisy response to the environment into the quorum-sensing model, heterogeneous states still arise; see Appendix 1-figure 2. Furthermore, the average production in the heterogeneous state is finely adjusted by the interplay between the response probability l and the selection strength s ( Figure 2F).
The establishment of long-lived, heterogeneous states induced by quorum sensing is one central finding of our study. We interpret this phenotypic heterogeneity as the result of the robust balance between population and ecological dynamics coupled through quorum sensing (see Box 1). On the one hand, fitness differences due to costly production favor non-producers. On the other hand, sensing the population average and accordingly up-regulating individual production enables producers to persist. Remarkably, fitness differences and sense-and-response balance such that separated production degrees may stably coexist in one population; the population does not become homogeneous at an intermediate production degree as one might naively expect. Heterogeneity of the autoinducer production is a robust outcome of the dynamics (and not a fine-tuned effect), and the average production level in the population is adjusted by the interplay of the response probability l and the selection strength s. Phenotypic heterogeneity does not rely on a bistable response function, but arises due to the global intercellular coupling of ecological and population dynamics through quorum sensing, as we show next. The relevance of quorum sensing for phenotypic heterogeneity in microbial populations is further explored below (see Discussion).

A B C
Variance at long times Amplitude of upregulation Population size Time to absorbing state Time to absorbing state Figure 3. Characterization of phenotypic heterogeneity in the quorum-sensing model. (A) For small response probability l, populations get stuck in heterogeneous quasi-stationary states. The time taken to reach a homogeneous absorbing state, T abs , increases exponentially with the population size N (filled circles denote the mean, gray bars denote the range within which 95% of the data points lie closest to the mean; dashed lines show fit to T abs~e gN ). (B) Heterogeneous states are long-lived only if l is small and the response function is nonlinear (in particular, up-regulation is required for some average production level such that RðhpiÞ > hpi). Here, the monostable response function RðhpiÞ ¼ hpi þ k sinðphpiÞ was chosen such that k 2 ½0; 1=p scales the magnitude of up-regulation. As k increases, the gap between the low-productive and high-productive peaks of the heterogeneous state becomes larger such that it takes longer to reach the absorbing state. Mean-field theory (1) predicts the existence and local stability of heterogeneous stationary distributions for 0 < l < l up ¼ s=2 (regime below the black line). Deviations between the stochastic process and mean-field theory are due to demographic fluctuations that vanish as N ! ¥. (C) The variance of production degrees in the population reveals whether the population is in a homogeneous (VarðpÞ ¼ 0) or heterogeneous state (VarðpÞ > 0). The variance was averaged over long times in the quasi-stationary state. Mean-field theory (1) (black line) agrees with our numerical observations (red filled circles); see Methods and materials. Ensemble size M ¼ 100, The following source data is available for figure 3: Source data 2. Source data accompanying Figure 3. DOI: 10.7554/eLife.25773.012

Results of mathematical analysis
In the following, the observed long-lived states of phenotypic heterogeneity in the quorum-sensing model are explained. First, we derived the macroscopic mean-field equation (the autoinducer equation (1)) from the microscopic dynamics of the quorum-sensing model. Second, we analyzed this mean-field equation and characterized phenotypic heterogeneity of autoinducer production.
The microscopic dynamics of the quorum-sensing model are captured by a memoryless stochastic birth-death process as sketched in Figure 1. Starting from the microscopic many-particle stochastic process, we derived a mean-field equation for the probability distribution of finding any individual at a specified production degree p at time t in the spirit of the kinetic theory in statistical physics (Kadar, 2007). We call this one-particle probability distribution the production distribution ; Figure 2 shows the corresponding histogram numerically obtained from the stochastic many-particle process. The mean-field equation for , which we refer to as the autoinducer equation, is obtained as: where Á t denotes averaging with respect to at time t. The details of the derivation of the autoinducer equation from the microscopic dynamics are given in the Methods and materials section and in Appendix 2.
The autoinducer equation (1) involves two contributions: the sense-and-response term with prefactor 2l, and the replicator term with prefactor 1 À 2l. Through the replicator term, probability weight at production degree p changes if the fitness fðpÞ is different from the mean fitness in the population f t (here fðpÞ À f t ¼ Àsðp À p t Þ). Without quorum sensing (l ¼ 0), Equation (1) reduces to the well-known replicator equation of the continuous Prisoner's dilemma (Bomze, 1990;Oechssler and Riedel, 2001;Hofbauer and Sigmund, 2003;Cressman, 2005;McGill and Brown, 2007). The sense-and-response term, on the other hand, encodes the global feedback by which individuals adopt the production degree Rðp t Þ upon sensing the average p t through quorum sensing at rate 2l. The difference between the current state and the state in which all individuals have this production degree Rðp t Þ determines the change in at every production degree. Through the replicator term and the sense-and-response term, the ecological dynamics (average production level p t ) are coupled with the dynamics of .
We now present our results for the long-time behavior of the autoinducer equation (1). First, the autoinducer equation (1) admits homogeneous stationary distributions. Without quorum sensing (l ¼ 0), the initially lowest production degree in the population, p low , constitutes the homogeneous stationary distribution ¥ ðpÞ ¼ dðp À p low Þ, which is attractive for generic initial conditions. With quorum sensing (l > 0), fixed points of the response function p Ã ¼ Rðp Ã Þ yield homogeneous stationary distributions as ¥ ðpÞ ¼ dðp À p Ã Þ, which are attractors of the quorum-sensing dynamics (1) for all initial distributions if l > s=2; see analysis below. These homogeneous stationary distributions confirm our observations of homogeneous absorbing states in the quorum-sensing model, in which all individuals produce to the same degree; see Figure 2A,B. Time scales at which stationarity is approached are discussed in the Methods and materials section.
Second, to analytically characterize long-lived heterogeneous states of the population, we decomposed into a distribution at low production degrees and a remainder distribution at higher degrees. We found that such a decomposition yields the bimodal, heterogeneous, stationary distribution of the autoinducer equation (1): with p high ¼ RðbÞ and y ¼ 1 À b=RðbÞ ; ( if the conditions 0 < p high 1 and 0 < y < 1 are fulfilled; see Box 1 for an illustration and Appendix 3 for the derivation. The parameter b ¼ 2l=s quantifies the balance between fitness differences and sense-and-response mechanism through quorum sensing. Heterogeneous stationary distributions (2) are constituted of a probability mass y at the low-producing degree p low ¼ 0 and a coexisting dpeak with stationary value 1 À y at a high-producing degree p high separated from p low by a gap. Such heterogeneous stationary distributions have mean p ¥ ¼ b and variance VarðpÞ ¥ ¼ bðRðbÞ À bÞ. Therefore, the interplay between selection strength s and response probability l adjusts the average production of autoinducers in the population ( Figure 2F). For simplicity, we assumed in Equation (2) that the initially lowest production degree in the population is p low ¼ 0; generalized bimodal distributions for arbitrary initial distributions 0 are given in Appendix 3. From the conditions on p high and y below Equation (2), one can derive the following conditions on the response function and the value of the response probability l (for given selection strength s) for the existence of heterogeneous stationary distributions: (i) The response function needs to be nonlinear with Rðp ¥ Þ ¼ p high > p ¥ ; that is, quorum sensing needs to up-regulate the cellular production in some regime of the average production level. Therefore, both monostable and bistable response functions depicted in Figure 1B may induce heterogeneous stationary distributions through the ecological feedback. (ii) The response probability needs to be small with l < l up ¼ s=2; that is, to induce phenotypic heterogeneity, cells must respond only rarely to the environmental cue p. This estimate of an upper bound on l is confirmed by our numerical results of the stochastic process ( Figure 3A-C). Vice versa, for a given response probability, the selection strength needs to be big enough to induce heterogeneous stationary distributions. As we show in the Methods and materials section, phase transitions in the space of stationary probability distributions govern the longtime dynamics of the autoinducer equation (1) from heterogeneity to homogeneity as the response probability changes (l ! 0 and l ! l up ); see Figure 3C.
For small l, the coexistence of the low-producing and the high-producing peaks in solution (2) is stable due to the balance of fitness differences and sense-and-response through quorum sensing. In Appendix 3 we show that the heterogeneous stationary distributions (2) are stable up to linear order in perturbations around stationarity. As our numerical simulations show, these bimodal distributions are the attractor of the mean-field dynamics (1) for a broad range of initial distributions when l is small; see Appendix 1- figure 1 for some examples. They are also robust against noisy inheritance, noisy perception, and noisy response as demonstrated in Appendix 1-figure 2. We interpret the stability of the bimodal stationary distributions (2) as follows (see also Box 1). Fitness differences quantified by the selection strength s increase probability mass at production degree p low , whereas nonlinear response to the environment with probability l pushes probability mass towards the upregulated production degree p high ¼ Rðp ¥ Þ. The gap p high À p low > 0 ensures that the exponential time scales of selection and sense-and-response stably balance the coexistence of both peaks; see Methods and materials. Because heterogeneous stationary distributions (2) are attractive and stable, heterogeneous states of the stochastic many-particle process arise and are quasi-stationary. Consequently, the time to reach a homogeneous absorbing state in the stochastic process through demographic fluctuations scales exponentially with the population size N (Elgart and Kamenev, 2004;Kessler and Shnerb, 2007;Assaf and Meerson, 2010;Frey, 2010;Hanggi, 1986); see Figure 3A. Thus, phenotypic heterogeneity is long-lived.
In summary, our mathematical analysis explains how phenotypic heterogeneity in the autoinducer production arises when quorum sensing up-regulates the autoinducer production in microbial populations (Box 1). As an emergent phenomenon, the population may split into two subpopulations: one in which cells do not produce autoinducers ('off' state, p low ¼ 0) and a second in which cells produce autoinducers ('on' state, p high ¼ Rð2l=sÞ>0), but grow slower. The fraction of individuals in the 'off' state is given by the value of y in Equation (2). If quorum sensing is absent (l ¼ 0), the whole population is in the 'off' state (y ¼ 1), whereas all individuals are in the 'on' state (y ¼ 0) if quorum sensing is frequent (l ! l up ). Only when response to the environment is rare (0 < l < l up ) can the two phenotypic states, p low and p high , coexist in the population (0 < y < 1). The transitions from heterogeneous to homogeneous populations are governed by nonequilibrium phase transitions when the response probability changes (l ! 0 and l ! l up ). Our mathematical analysis shows that phenotypic heterogeneity arises dynamically, is robust against perturbations of the autoinducer production in the population, and is robust against noise at the level of inheritance, sense, and response.

Discussion
Summary: Phenotypic heterogeneity in the quorum-sensing model as a collective phenomenon through an ecological feedback In this work, we studied a conceptual model for the heterogeneous production of autoinducers in quorum-sensing microbial populations. The two key assumptions of our quorum-sensing model are as follows. First, production of large autoinducer molecules and accompanied gene expression in the cell's phenotypic state are negatively correlated with fitness such that non-producers reproduce faster than producers. Second, cells sense the average production level of autoinducers in the population and may accordingly up-regulate their production through quorum sensing. As a result, not only does the interplay between fitness differences and sense-and-response give rise to homogeneously producing populations, but it can also induce a heterogeneous production of autoinducers in the population as a stable collective phenomenon. In these heterogeneous states, the average production level of autoinducers in the population is adjusted within narrow limits by the balance between fitness differences (selection strength s in the model), and the rate with which cells respond to the environment and up-regulate their production through quorum sensing (response probability l and response function RðhpiÞ in the model). Due to this robust adjustment of the production level in the population, the expression of other genes (for example, bioluminescence and virulence genes) can be regulated by quorum sensing even when the production of autoinducers is heterogeneous in the population.
In the following, we discuss the assumptions of our model in the light of the empirical reality for both quorum sensing and phenotypic heterogeneity. Furthermore, we indicate possible directions to experimentally test the ecological feedback that is suggested by the results of our theoretical work.

Does autoinducer production reduce individual growth rate?
In our quorum-sensing model, it is assumed that the individual's production degree of autoinducers is negatively correlated with its growth rate (f i ¼ 1 À sp i ). Is this assumption of growth impairment for producing phenotypes justified (Parsek and Greenberg, 2005)? This would be the case if cellular production of autoinducers directly causes a reduction of the cell's growth rate. For example, in L. monocytogenes populations, heterogeneous production was observed for an autoinducer oligopeptide that is synthesized via the agr operon (Garmyn et al., 2011, Garmyn et al., 2009. This signaling oligopeptide incurs high metabolic costs through the generation of a larger pre-protein. For the oligopeptide signal synthesized via the agr operon in Staphylococcus aureus, the metabolic costs were conservatively estimated by Keller and Surette to be 184 ATP per molecule (metabolic costs for precursors were disregarded in this estimate); see Keller and Surette, 2006 for details. In contrast, basically no costs (0-1 ATP) incur for the different signaling molecule Autoinducer-2 (AI-2) that is considered as a metabolic by-product. As to what extent the production of oligopeptides for signaling reduces an individual's growth rate has, to our knowledge, not been studied quantitatively.
For quorum-sensing systems that involve N-acyl homoserine lactones (AHLs) as signaling molecules, however, a reduced fitness of producers has been reported for microbial growth in batch culture (Ruparell et al., 2016;Diggle et al., 2007;He et al., 2003). Even though metabolic costs for the synthesis of C 4 -HSL (one of the simplest AHL signaling molecules that is synthesized via the rhl operon) were conservatively estimated with only 8 ATP per molecule (Keller and Surette, 2006), a growth impairment was experimentally reported only recently for a C 4 -HSL-producing strain (Ruparell et al., 2016). Furthermore, a strain producing a long-chain AHL (OC 12 -HSL, synthesized via the las operon) showed a reduced fitness in both mono and mixed culture compared with a non-producing strain. The reduced fitness of AHL-producers was attributed to (i) metabolic costs of autoinducer production, in particular also to metabolic costs of precursors that were disregarded in the estimates by Keller and Surette, 2006, and (ii) accumulation of toxic side products accompanying the synthesis of autoinducers (Ruparell et al., 2016). As another example, the strain Sinorhizobium fredii NGR234 synthesizes AHLs via both the ngr and the tra operon (Schmeisser et al., 2009), and it was shown that gene expression related to autoinducer production reduces the strain's growth rate in mono culture (He et al., 2003). On the other hand, a heterogeneous expression of the corresponding autoinducer synthase genes was observed during growth of NGR234 only recently (Grote et al., 2014). As to what extent the production of AHLs reduces fitness of NGR234 in mixed culture and, thus, whether the phenotypic heterogeneity observed in Grote et al., 2014 could be explained through the ecological feedback proposed by our quorum-sensing model, remains to be explored experimentally.
In the quorum-sensing model, even small growth rate differences between producer and non-producer, which are quantified by the ratio (growth rate of producer) / (growth rate of non-producer) ¼ 1 À s, may give rise to a bimodal production of autoinducers in the population. Furthermore, it would be interesting to track the expression level of autoinducer synthase genes of a microbial strain during growth for which growth differences between the producing and the non-producing phenotype are known such as in the study of Ruparell et al., 2016. We emphasize that it would be desirable to report the full distribution of expression levels in the population in order to detect whether a population splits into several subpopulations; note that variance or percentiles are not suitable measures to characterize and compare the bimodality of distributions. A bimodal expression of autoinducer synthase genes in the population together with a tightly controlled average expression level could be a signature of the feedback between ecological and population dynamics underlying the observation of phenotypic heterogeneity as suggested by our results.
A question of spatio-temporal scales: How stable and how dispersed are autoinducers in the environment?
Autoinducers are secreted into the environment where they get dispersed and are degraded. For simplicity and to facilitate our mathematical analysis, we assumed in the quorum-sensing model that individuals respond to the current average production level of autoinducers in the whole population. Temporal availability and spatial dispersal of autoinducers determine whether this assumption is valid or not. On the one hand, temporal availability of autoinducers in the environment for signaling depends on many factors. For example, pH and temperature influence the stability of autoinducers (Yates et al., 2002;Byers et al., 2002;Decho et al., 2009;Grandclément et al., 2016;Hmelo, 2017). Biochemical mechanisms that inhibit or disrupt the functioning of signaling molecules (commonly referred to as 'quorum quenching') further determine the time scales at which autoinducers are degraded in the environment (LaSarre and Federle, 2013;Grandclément et al., 2016;Hmelo, 2017). On the other hand, spatial dispersal of autoinducers in the population depends, for example, upon cellular mechanisms that import and export autoinducers into the cell from the environment and vice versa, and upon the spatial structure of the microbial population (Platt and Fuqua, 2010;Hense and Schuster, 2015). The degree of dispersal determines whether autoinducers remain spatially privatized to a single cell, diffuse to neighboring cells, or are spread evenly between all cells of the population. Consequently, the spatio-temporal organization of the microbial population determines as to what extent microbes sense rather the current average production level or a time-integrated production of autoinducers, and to what extent they sense rather the global or a local average production level. Our quorum-sensing model assumes that autoinducers are uniformly degraded in a well-mixed environment. These assumptions do not hold true for a spatially structured microbial biofilm, but should be fulfilled during the stationary phase of microbial growth in a wellmixed batch culture (Yates et al., 2002;Byers et al., 2002).
How is production of autoinducers up-regulated at the single-cell level? Monostable or bistable up-regulation of autoinducer synthesis at the singlecell level Our theoretical results also relate to the question of how cells regulate the production of autoinducers upon sensing the level of autoinducers in the environment. In this work, we showed that positive feedback loops and, thus, up-regulation of cellular autoinducer production may give rise to phenotypic heterogeneity. Positive feedback loops are mathematically introduced in our model as a stable fixed point at the producing phenotype of the response function (up-regulation to the stable 'on' state at p ¼ 1; see Figure 1B). Such a positive feedback is not present in all autoinducer synthase systems, but was reported for the strains L. monocytogenes and S. fredii NGR234 (Waters and Bassler, 2005;Garmyn et al., 2009;He et al., 2003;González and Marketon, 2003) that showed a heterogeneous synthesis of autoinducers at the population level (Garmyn et al., 2011;Grote et al., 2014). From an experimental point of view it is often not known, however, whether autoinducer synthesis is up-regulated for all autoinducer levels or only above a threshold level. Up-regulation at all production levels in the population corresponds to a monostable response function with an unstable fixed point at the 'off' state at p ¼ 0, whereas up-regulation only above a threshold level corresponds to a bistable response function with a stable fixed point at the 'off' state at p ¼ 0 and an additional unstable fixed point at the threshold value (see Figure 1B). Most models of quorum-sensing microbial populations explicitly or implicitly assume a bistable gene regulation for positive feedback loops without experimental verification; see (Hense and Schuster, 2015) for further discussion. Why might it be relevant to distinguish between bistable (for example, a Hill function with Hill coefficient > 1) and monostable (for example, a Hill function with Hill coefficient 1) regulation of autoinducer synthesis -apart from the insight on how regulation proceeds at the molecular level? As the results of our quorum-sensing model show, the qualitative form of the regulation could discriminate between different mechanisms that control phenotypic heterogeneity of the autoinducer production at the population level as we describe in the following.

Heterogeneity through stochastic gene expression only for bistable gene regulation
In recent years, a deeper mechanistic understanding of phenotypic heterogeneity has been achieved by exploring how the presence of different phenotypes in a population of genetically identical cells depends upon molecular mechanisms and stochasticity at the cellular level (Ackermann, 2015). For example, a bistable gene regulation function enables cells to switch between an 'on' and an 'off' state with respect to the expression of a certain gene or operon. Depending on environmental cues, cells are either in the stable 'on' or in the stable 'off' state. A noisy expression at intermediate concentrations of an environmental cue may then cause some cells to be in the 'on' state while others are still in the 'off' state. Thus, stochastic gene expression explains the coexistence of different phenotypic states in one population in many experimental situations (Novick and Weiner, 1957;Ozbudak et al., 2004;Kaern et al., 2005;Dubnau and Losick, 2006;Smits et al., 2006;Raj and van Oudenaarden, 2008;Eldar and Elowitz, 2010). In the context of quorum sensing, the level of autoinducers in the population is the environmental cue that triggers the stochastic switch between 'on' and 'off' state explaining heterogeneous autoinducer production when the response function is bistable (Fujimoto and Sawai, 2013;Pérez-Velázquez et al., 2016;Goryachev et al., 2005;Dockery and Keener, 2001). In other words, bistable regulation together with stochastic gene expression can explain a bimodal autoinducer synthesis in the population. If, however, regulation of autoinducer synthesis is monostable, an explanation of phenotypic heterogeneity in the autoinducer production in terms of stochastic gene expression appears questionable to us.

Heterogeneity through an ecological feedback for monostable and for bistable gene regulation
The analysis of our quorum-sensing model suggests that an alternative mechanism could explain a heterogeneous production of autoinducers in quorum-sensing microbial populations. Our results show that phenotypic heterogeneity may also arise dynamically as a collective phenomenon for monostable regulation of autoinducer production when quorum sensing creates an ecological feedback by coupling ecological with population dynamics. Cells need to up-regulate their expression with respect to the sensed production level in the population. A threshold-like, bistable response function does not need to be assumed in the quorum-sensing model, but would work as well, to establish a bimodal production of autoinducers in the population. Therefore, if phenotypic heterogeneity of autoinducer synthesis is observed in a microbial population and if cellular growth rate is correlated with the cell's production degree of autoinducers, then it would be worth testing experimentally whether regulation of autoinducer synthesis is monostable or bistable. Monostable regulation would be an indicator that heterogeneity on the population level is not caused by stochastic gene expression, but actually is caused by a different mechanism such as the ecological feedback proposed here.

On which timescales do microbes respond to autoinducers in the environment?
Furthermore, in our implementation of the quorum-sensing model, individuals respond to the environment with response probability l upon reproduction. The rule that offspring individuals can only respond at reproduction events represents a coarse-grained view in time to facilitate the mathematical analysis and to identify the ecological feedback. The response probability can actually be interpreted as the rate with which individuals respond to autoinducers in the environment. This cellular response rate is then effectively measured in units of the cell's reproduction rate (f i ) in the quorumsensing model. Phenotypic heterogeneity of autoinducer production arises in the quorum-sensing model if the time scale at which cells respond to autoinducers in the environment is of similar order as or larger than the time scale at which growth rate differences affect the population dynamics. This can be inferred from the prefactors of the sense-and-response term and the replicator term in the autoinducer equation (1): Effective changes of the distribution of autoinducer production in the population occur (i) through cellular response to autoinducers in the environment at rate~2l and (ii) through growth rate differences at rate~s. Both contributions need to balance each other such that a bimodal production in the population is established (quantified in our model by the ratio b ¼ 2l=s; see also Box 1 for an illustration). This balance is robust against several kinds of perturbations and noise as discussed above; see Appendix 1-figures 1 and 2. To understand how bacteria respond to changes of autoinducer levels in the environment and to quantify response rates, experiments at the single-cell level seem most promising to us at present.

Single-cell experiments
Some of the questions raised above may be addressed most effectively with single-cell experiments. For example, it would be desirable to simultaneously monitor, at the single-cell level, the correlations between autoinducer levels in the environment, the expression of autoinducer synthase genes, and the transcriptional regulators that mediate response to quorum sensing. Upon adjusting the level of autoinducers in a controlled manner, for example in a microfluidic device, one could characterize how cells respond to autoinducers in the environment. This way, it might be possible to answer questions of (i) how the cellular production of autoinducers is regulated (monostable or bistable regulation, or a different form of regulation), (ii) whether response times to environmental changes are stochastic and whether response rates can be identified, (iii) as to what extent cellular response in the production of autoinducers depends on both the level of autoinducers in the environment and on the cell's present production degree, and (iv) how production of autoinducers is correlated with single-cell growth rate. In the context of the quorum-sensing model, the results of such single-cell experiments would help to identify the form of the fitness function f and the response function R, to quantify the selection strength s and response probability l, and to refine the model set-up.
Different mechanisms at the cellular (microscopic) level may yield the same behavior at the population (macroscopic) level. Therefore, observations at the population level might not discriminate between different mechanisms at the cellular level. Is phenotypic heterogeneity in the production of autoinducers an example of such a case? In this work, we discussed that phenotypic heterogeneity in the autoinducer production could be the result of stochastic gene expression in bistable gene regulation or, as suggested by our model, the result of the feedback between ecological and population dynamics. We believe that the above-mentioned single-cell experiments could elucidate the mechanisms that allow for phenotypic heterogeneity in quorum-sensing microbial populations, and help to understand how population dynamics and ecological dynamics influence each other.
What is the function of phenotypic heterogeneity in autoinducer production?
The purpose of the quorum-sensing model presented here is to explain how phenotypic heterogeneity in the autoinducer production arises and how it is controlled in quorum-sensing microbial populations. With the current model set-up, however, we did not address its function. Why might this phenotypic heterogeneity in the autoinducer production be beneficial for a microbial species on long times? From an experimental point of view, the evolutionary contexts and ecological scenarios under which this phenotypic heterogeneity may have arisen are still under investigation (Garmyn et al., 2011;Grote et al., 2014, Grote et al., 2015. From a modeling perspective, one could extend, for example, our chosen fitness function with a term that explicitly accounts for the benefit of signaling either at the cellular or population level, and study suitable evolutionary contexts and possible ecological scenarios (Pollak et al., 2016;Dandekar et al., 2012;Czárán and Hoekstra, 2009;Carnes et al., 2010;Hense and Schuster, 2015). Such theoretical models together with further experiments might help to clarify whether heterogeneous production of autoinducers can be regarded as a bet-hedging strategy of the population or rather serves the division of labor in the population (Ackermann, 2015).

Conclusion
Overall, our analyses suggest that feedbacks between ecological and population dynamics through signaling might generate phenotypic heterogeneity in the production of signaling molecules itself, providing an alternative mechanism to stochastic gene expression in bistable gene-regulatory circuits. Spatio-temporal scales are important for the identified ecological feedback to be of relevance for microbial population dynamics: growth rate differences between producers and non-producers need to balance the rate at which cells respond to the environment, degradation of signaling molecules should be faster than time scales at which growth rate differences affect the population composition significantly, and signaling molecules should get dispersed in the whole population faster than they are degraded. In total, if microbes sense and respond to their self-shaped environment under these conditions, the population may not only respond as a homogeneous collective as is typically associated with quorum sensing, but may also become a robustly controlled heterogeneous collective. Further experimental and theoretical studies are needed to clarify the relevance of the different mechanisms that might control phenotypic heterogeneity, in particular for quorum-sensing microbial populations.

Materials and methods
Derivation of the autoinducer equation (1) The microscopic dynamics are captured by a memoryless stochastic birth-death process (a continuous-time Markov process) as sketched in Figure 1. The state of the population p is updated by nongenetic inheritance and sense-and-response through quorum sensing such that at most two individuals i and j 6 ¼ i change their production degree at one time. The temporal evolution of the corresponding joint N-particle probability distribution Pðp; tÞ is governed by a master equation for the stochastic many-particle process (Gardiner, 2009;Van Kampen, 2007;Weber and Frey, 2017), whose explicit form is derived from Figure 1 and given in Appendix 2. This master equation tracks the correlated microscopic dynamics of the production degrees of all N individuals. To make analytical progress, we focused on the reduced one-particle probability distribution ð1Þ ðp; tÞ ¼ 1=Nh P i dðp À p i Þi P in the spirit of a kinetic theory (Kadar, 2007) starting from the microscopic stochastic dynamics. ð1Þ denotes the probability distribution of finding any individual at a specified production degree p at time t; the numerically obtained histogram of ð1Þ was plotted in Figure 2. The temporal evolution of ð1Þ is derived from the master equation, and couples to the reduced two-particle probability distribution and to the full probability distribution P through quorum sensing. By assuming that correlations are negligible, one may approximate ð1Þ by the mean-field distribution , which we refer to as the production distribution. The mean-field equation (1) for is derived in Appendix 2 and referred to as the autoinducer equation. Note that Equation (1) conserves normalization of , that is, R 1 0 dp q t ðp; tÞ ¼ 0. We also proved that ð1Þ converges in probability to as N ! ¥ for any finite time if initial correlations are not too strong. In other words, the autoinducer equation (1) captures exactly the collective dynamics of the stochastic many-particle process for large N. To show this convergence, we introduced the bounded Lipschitz distance d between and ð1Þ , applied Grö nwall's inequality to the temporal evolution of d, and used the law of large numbers; see  for details. Similar distance measures and estimates have been used, for example, to prove that the Vlasov equation governs the macroscopic dynamics of the above-mentioned classical XY spin model with infinite range interactions (Braun and Hepp, 1977;Dobrushin, 1979;Spohn, 1991;Yamaguchi et al., 2004).

Analysis of homogeneous stationary distributions of the autoinducer equation (1)
Without quorum sensing (l ¼ 0), one finds the analytical solution for by applying the method of characteristics to Equation (1) in the space of moment and cumulant generating functions as: ðp; tÞ ¼ 0 ðpÞe Àstp = R 1 0 dp e Àstp 0 ðpÞ; see Appendix 3 for details. Thus, the initially lowest production degree in the population, p low , constitutes the homogeneous stationary distribution ¥ ðpÞ ¼ dðp À p low Þ, which is attractive for generic initial conditions. Only d-peaks at production degrees greater than p low are stationary as well, but they are neither attractive nor stable. The temporal approach to the homogeneous stationary distribution is algebraically slow for continuous initial distributions 0 , and exponentially fast if p low is separated from all greater degrees by a gap in production space; see Appendix 3 and Figure 2D.
With quorum sensing (l > 0), fixed points of the response function p Ã ¼ Rðp Ã Þ yield homogeneous stationary distributions of the autoinducer equation (1) as ¥ ðpÞ ¼ dðp À p Ã Þ. In particular, stable fixed points of the response function (R 0 ðp Ã Þ < 1) constitute homogeneous stationary distributions that are stable up to linear order in perturbations around stationarity. For l > s=2, these distributions are also attractors of the mean-field dynamics (1) for all initial distributions; see Appendix 3. The temporal approach towards homogeneous stationary distributions with quorum sensing is generically exponentially fast ( Figure 2E). This exponentially fast approach is illustrated for the special case of a linear response function and l ¼ 1=2, for which one finds the analytical solution as: ðp; tÞ ¼ yðtÞ 0 ðpÞ þ ð1 À yðtÞÞdðp À p 0 Þ with yðtÞ ¼ e Àf 0 t . However, time scales at which stationarity is approached may diverge at bifurcations of the response function. Such can be seen, for example, if one chooses a supercritical pitchfork bifurcation of a polynomial response function and l ¼ 1=2; see Appendix 1- figure 3 and Appendix 3.
Phase transitions from heterogeneity to homogeneity in the autoinducer equation (1) Here we discuss how the long-time behavior of the quorum-sensing model changes from heterogeneous to homogeneous populations as the response probability l vanishes or reaches the upper threshold l up while the selection strength s is kept fixed. For small response probabilities, 0 < l < l up , the heterogeneous stationary distributions of the autoinducer equation (1) explain the long-lived, heterogeneous states of the stochastic quorum-sensing process. The coexisting d-peaks at the lowproducing and high-producing degree in the heterogeneous stationary distribution are separated by a gap in production space, which gives rise to the non-vanishing variance VarðpÞ ¥ in the phase of heterogeneity ( Figure 3C). As l ! l up , the gap closes, p high ! Rðp high Þ, and y ! 0, such that a homogeneous stationary distribution with VarðpÞ ¥ ¼ 0 is recovered in a continuous transition. This nonequilibrium phase transition from heterogeneity to homogeneity proceeds without any critical behavior. As l ! 0, and under the assumption that 0 is an unstable fixed point of the response function (Rð0Þ ¼ 0 and 1 < R 0 ð0Þ; we further assume R 0 ð0Þ < ¥), the gap between the low-producing and the high-producing peak closes as well because p high ! 0. However, y does not approach 1, but the value 1 À 1=R 0 ð0Þ < 1. The probability weight at the low-producing mode jumps by the value 1=R 0 ð0Þ and the homogeneous stationary distribution with VarðpÞ ¥ ¼ 0 is recovered in a discontinuous transition. Therefore, a discontinuous phase transition in the space of stationary probability distributions governs the long-time dynamics of the autoinducer equation (1) from heterogeneity to homogeneity as the response probability l vanishes (for fixed selection strength s). The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. Phenotypic heterogeneity in the quorum-sensing model is robust against noisy inheritance, noisy perception, and noisy response. Upon including either noisy inheritance of the production degree (A-C), or noisy perception of the average production level and noisy response to it (D-F), or both percx="">(G-I) into the model set-up, bimodal quasistationary states still arise in the relevant parameter regimes (see Figure 2C). Depicted are representative single realizations of the modified stochastic process (histogram over normalized values of production degrees to make the comparison with Figure 2 possible).

Author contributions
(A-C) Noisy inheritance is implemented at reproduction events. Production degree p i is passed on to an offspring as p i 7 ! p i þ h p with noise h p~N ð0; s p Þ sampled from a Normal distribution (and are cut off such that p i þ h p 2 ½0; 1), emulating noisy inheritance of the phenotype. s p ! 0 characterizes the strength of the noise (s p ¼ 0 recovers noiseless inheritance). As s p increases, bimodal quasi-stationary states still arise, but the two peaks become broader than in the noiseless case. (D-F) Noise in the sensing apparatus is implemented as noisy perception of the average production level hpi 7 ! hpi þ h hpi with Gaussian noise h hpi~N ð0; s hpi Þ, and noise in the response is implemented at the level of the response function as RðhpiÞ 7 ! RðhpiÞ þ h R with Gaussian noise h R~N ð0; s R Þ. Therefore, the production degree of an individual is updated through sense-and-response to the environment as p i ¼ RðhpiÞ 7 ! Rðhpi þ h hpi Þ þ h R in the quorum-sensing model. Again, as the strength of both sense and response noise increase, bimodal quasi-stationary states still arise, but the two peaks become broadened compared with the noiseless case. We emphasize that s hpi ¼ s R ¼ 0:1 corresponds to very strong noise on the interval ½0; 1. (G-I) Combined effect of noisy inheritance and noisy sense-and-response. Representative trajectories demonstrate that bimodal quasi-stationary states also arise in the presence of noise at all update steps. Thus, phenotypic heterogeneity in the quorum-sensing model is qualitatively robust against noise at all steps. Initial distribution: p i~U niformð0; 1Þ, independent and identically distributed; Parameters: selection strength s ¼ 0:2, response probability l ¼ 0:05, response function RðhpiÞ ¼ hpi þ 0:2 Á sinðphpiÞ, and population size N ¼ 10 4 .

C D
Appendix 1-figure 3. Time scales at which stationarity is approached may diverge. The response probability was set to l ¼ 1=2, and the nonlinear response function RðhpiÞ ¼ hpi þ 40 Á hpiðhpi À ð0:5 À ÞÞðhpi À 0:5Þðhpi À ð0:5 þ ÞÞðhpi À 1Þ with bifurcation parameter was chosen, see Equation (48); controls a supercritical pitchfork bifurcation of the response function at the fixed point p cr ¼ 0:5 (Rðp cr Þ ¼ p cr ): For > 0, the fixed point at p cr is unstable and non-degenerate (sketch in (B)), and becomes stable and threefold degenerate (z ¼ 3) as ¼ 0 (sketch in (A)). (D) Away from the bifurcation of the response function ( > 0), the approach of an absorbing state in the stochastic many-particle system is exponentially fast (see inset of (D) for an exemplary measurement of hpiðtÞ À hpi ¥ for ¼ 0:1, dashed line denotes fit to exponential decay). The exponentially fast approach of stationarity is confirmed by mean-field theory (p t À p ¥~e Àt=t ), see main text and Equation (51). Mean-field theory also predicts that the time scale of this exponentially fast relaxation diverges as t~ À2 as the bifurcation is approached ( ! 0), indicated by the black line in (D). This prediction agrees with the numerical simulations of the stochastic quorum-sensing model, see (D) (blue crosses denote values of the decay constants obtained from the exponential fits and black dashed line indicates fit to t~1= g with g ¼ 1:95). The divergence of time scales reflects critical slowing down as ! 0. (C) At the bifurcation of the response function ( ¼ 0), the approach of an absorbing state is algebraically slow, p t À p ¥~t À1=n with critical exponent n ¼ z À 1 ¼ 2 obtained from mean-field theory (black line), see Equation (53). This prediction agrees with our numerical simulations of the stochastic quorum-sensing model (black dashed line in (C) indicates fit to hpiðtÞ À hpi ¥~t a with a ¼ À0:50). Initial distribution: unimodal p i~B etað1; 10Þ, independent and identically distributed; Parameters: Ensemble size M ¼ 100, selection strength s ¼ 0:1, population size N ¼ 10 4 . DOI: 10.7554/eLife.25773.015 with reproduction rate of individual i (fitness) given by (and selection strength 0 s <1): and death rate of individual j given by (random death): The transition probabilities A i and A j account for the ensuing changes of at most two production degrees in the population due to non-genetic inheritance and sense-andresponse (see detailed description below Equation (6) and Equation (7)). The initial condition to the master equation (3) is given as Pðp; t ¼ 0Þ ¼ p 0 ðpÞ.
The master equation (3) involves two contributions: gain terms yielding an increase and loss terms yielding a decrease of the probability weight in state p at time t. Loss terms occur when the population is in state p and an individual reproduces. The probability of finding the population in this state is given by Pðp; tÞ. Individual i is selected for reproduction at rate f i ðpÞ and splits into two offspring individuals, and a different individual j 6 ¼ i is removed with probability 1=ðN À 1Þ at the same time (random death). Gain terms involve all events that take the population from an arbitrary state e p i;j to state p, and involve again reproduction for individual i and neutral death for individual j. The transition probabilities A i and A j account for these changes due to non-genetic inheritance and sense-and-response through quorum sensing, and are given as: We abbreviate he p i;j i ¼ 1=N P k ðe p i;j Þ k as the average production degree before the update step. Both transition probabilities A i and A j quantify the probability of attaining the production degrees p i and p j , respectively, for the two offspring individuals of ancestor i. The first summand in both A i and A j captures the response to the perceived average production (p i=j attains the value Rðhe p i;j iÞ) with probability l as the updated production degree, and the second summand accounts for the non-genetic inheritance of the production degree from the ancestor i (p i=j attains the value e p i ) with probability 1 À l. Note that for the transition probability A j , also p j attains the value e p i due to our convention that individual i is labeled as the reproducing individual and individual j is chosen for the death event, see Figure 1 of the main text.
In our prescription of the master equation (3), the introduced gain and loss terms also involve terms that actually do not change the state of the population. Such is the case, for example, when the two individuals i and j have the same production degree (e p i ¼ e p j ) and both offspring individuals retain the production degree from their ancestor i (p i ¼ p j ¼ e p i , that is, both offspring individuals do not update their production through sense-andresponse). Such events do not change the state of the population (e p i;j ¼ p), but are included in the master equation (3). However, these terms always occur both in the gain and loss terms. Therefore, they cancel each other and the master equation can be written in form of Equation (3).
Coarse-grained description: Reduced one-particle probability distribution The reduced one-particle probability distribution ð1Þ is defined as: ¼ Z ½0;1 NÀ1 dp 2 dp 3 . . . dp N Pðp; tÞ ¼ P ð1Þ ðp; tÞ ; and agrees with the marginal probability distribution for the production degree of the first individual P ð1Þ . The equality between the normalized reduced one-particle distribution ð1Þ and the one-particle probability distribution P ð1Þ follows from the symmetry of P with respect to permutation of identical (that is indistinguishable) individuals (Kadar, 2007).
Towards the macroscopic dynamics: Temporal evolution of the reduced one-particle probability distribution In the following we show that the temporal evolution equation of the reduced one-particle probability distribution is obtained from the master equation (3) as: ðp; tÞ ¼ 2l Z ½0;1 N dp 1 dp 2 . . . dp N ðNÞ ðp; tÞ 1 À sp ð Þdðp À RðhpiÞÞ (12) À 2l ð1Þ ðp; tÞ À s Z 1 0 dp 2 ð2Þ ðp; p 2 ; tÞ p 2 þ ð1 À 2lÞs Z 1 0 dp 2 ð2Þ ðp; p 2 ; tÞ p 2 À p ð1Þ ðp; tÞ : To derive the temporal evolution equation for ð1Þ , we specify the production degree of one particular individual (here p 1 ), and integrate out the production degrees of the other N À 1 individuals in the master equation (3): de p i de p j Pðe p i;j ; tÞ fð e p i Þ N À 1 A i ðe p i;j ; iÞA j ðe p i;j ; iÞ (13) À Z ½0;1 NÀ1 dp 2 . . . dp N Pðp; tÞ dp 2 . . . dp N gain À loss ð Þ ¼: I gain À I loss : For the loss term, we split the sum P i * ðiÞ into two contributions: and deal with both contributions separately to obtain: I loss ¼ NP ð1Þ ðp; tÞ À spP ð1Þ ðp; tÞ À sðN À 1Þ Z 1 0 dp 2 P ð2Þ ðp; p 2 ; tÞ p 2 : For the gain term, we split up the sum P i P j6 ¼i * ði;jÞ that occurs in the master equation (3) into three terms as follows: We also introduce the notation de p i;j;k :¼ dp 1 dp 2 . . . de p i . . . de p j . . . dp k . . . dp N in which variables in the superscript are labeled with a tilde in the product (indices i and j in the example), and variables with a hat in the superscript are missing in the product (that is, they are not integrated over; index k in the example). This way, the integral measure in the gain term can be decomposed as follows: j6 ¼i dp1 ;i;j dp i dp j : Upon plugging in the specific form of the transition probabilities and decomposing the integral measure into the three contributions, the gain term can be written as follows (note the asymmetry between the first summand (i ¼ 1 term) and the second summand (j ¼ 1 term); integration over suitable d-functions of the transition probabilities was carried out as well, for example, R 1 0 dp j A j ðe p i;j ; iÞ ¼ 1): de p 1;j Pðe p 1;j ; tÞfðe p 1 Þ ldðp 1 À Rðhe p 1;j iÞÞ þ ð1 À lÞdðp 1 À e p 1 Þ À Á Making use of the fact that P is symmetric with respect to permutation of individuals (individuals are identical), carrying out possible integrals over d-functions, plugging in the explicit form of the fitness function (4), and relabeling variables, one obtains for the gain term: ðp; tÞ þ ðN À 2ÞP ð1Þ ðp; tÞ À sðN À 2Þ Z 1 0 dp 2 P ð2Þ ðp; p 2 ; tÞ p 2 : Combining loss terms I loss and gain terms I gain leads to the result for the equation of motion of the reduced one-particle probability distribution ð1Þ that is given in Equation (12).
Heuristic derivation of the macroscopic dynamics: Mean-field approximation Upon assuming that correlations are negligible, one may approximate ð1Þ by its mean-field approximation , which we refer to as the production distribution. As described in the main text, the temporal evolution equation for ð1Þ serves as a suitable starting point to guess the mean-field equation for , which is the mean-field approximation of ð1Þ . Thus, we naively approximate ð1Þ » and ðNÞ » Q N . From the temporal evolution of ð1Þ in Equation (12), the mean-field equation for is suggested as: q t ðp; tÞ » 2l Z ½0;1 N dp 1 dp 2 . . . dp N Y N i¼1 ðp i Þ 1 À sp t ð Þdðp À Rðp t ÞÞ À 2l ðp; tÞ À s Z 1 0 dp 2 ðp; tÞðp 2 ; tÞp 2 þ ð1 À 2lÞs Z 1 0 dp 2 ðp; tÞðp 2 ; tÞp 2 À pðp; tÞ ; where Á t denotes averaging with respect to at time t. Further collection of terms yields the mean-field equation (1) in the main text: with initial condition ðp; t ¼ 0Þ ¼ 0 ðpÞ, fðpÞ ¼ 1 À sp, f t ¼ 1 À sp t , and p t ¼ R 1 0 dp pðp; tÞ. Alternatively, this mean-field equation can also be written as: We emphasize that the mean-field equation (1) is to be understood in distributional sense, that is, it needs to be integrated over observables (for example, suitable test functions g : ½0; 1 ! R, g smooth) and is interpreted as a linear functional on the space of these observables. This way, can be a continuous probability density function or a discrete probability mass function, or a probability distribution with both density parts and mass parts. To keep notation accessible for a broad readership, we avoid a measure-theoretic notation in this manuscript.
The proof that ð1Þ converges in probability to as N ! ¥ for any finite time if initial correlations are not too strong will be presented in a forthcoming publication .
Without sense-and-response (l ¼ 0): Analytical solution and approach of the homogeneous stationary distribution of non-producers For the case without sense-and-response through quorum sensing, l ¼ 0, it is readily seen from Equation (1) that stationary production distributions are given by d-peaks as ¥ ðpÞ :¼ ðp; t ! ¥Þ ¼ dðp À p low Þ for all p low 2 ½0; 1. However, the distribution with solely non-producers, p low ¼ 0, is the only asymptotically stable solution of the mean-field equation (1); see below.
When sense-and-response is absent, the analytical solution of the mean-field equation (1) for can be obtained by applying the method of characteristics to Equation (27) as outlined above. The implicit solution is given by: Cðu; tÞ ¼ C 0 ðu À stÞ þ sthpi t ; with hpi t :¼ 1=t as the temporal average of the mean production p t . Back-transformation and exploiting normalization of yields: For example, if the initial production distribution 0 is a uniform distribution on ½0; 1, evolves in time as ðp; tÞ ¼ st=ð1 À e Àst Þe Àstp , which is plotted in Figure 2A of the main text (black, solid lines). Every production degree that is different from p ¼ 0 decays exponentially fast and the time scale of the decay is set by the inverse of the value of that production degree. As p ! 0, this time scale diverges and, hence, the stationary distribution, ¥ ðpÞ ¼ dðpÞ ; is approached algebraically slowly; see Figure 2D of the main text.
To quantify the dependence of the time scales to approach stationarity on the initial distribution in more generality, we analyzed the temporal solution of the mean p t , which is obtained from the solution for the cumulant generating function as: Therefore, the temporal evolution of the mean production depends only on the initial distribution 0 via its Laplace transform L½ 0 . For the asymptotic behavior of Laplace transforms it is known that if 0 ðpÞ~p as p ! 0 with > À 1, then L½ 0 ðvÞ~1=v ðþ1Þ for v ) 1 (Doetsch, 1976). Therefore, it follows that the mean evolves in time as p t~1 =t for t ) 1 if the initial production distribution is a continuous probability density with non-vanishing weight at p low ¼ 0 (chosen for simplicity as the lowest production degree). The condition that the exponent satisfies > À 1 is always fulfilled for a continuous probability distribution to ensure integrability at zero. In the same manner, the decay of the variance is shown to evolve in time algebraically as VarðpÞðtÞ~1=t 2 for t ) 1.
In contrast, if the lowest production degree is separated from all other degrees in the population by a gap D > 0 in production space, mean and variance approach their stationary value exponentially fast at a time scale set by D. To see this qualitative difference in the approach of stationarity, we consider an initial probability distribution with probability mass y 0 > 0 at degree p low ¼ 0 (chosen again for simplicity) and a remainder probability distribution e 0 with support on ½D; 1: 0 ðpÞ ¼ y 0 dðpÞ þ ð1 À y 0 Þe 0 ðpÞI ½D;1 ðpÞ (here I ½D;1 denotes the indicator function, which takes value 1 on the interval ½D; 1 and 0 otherwise, and highlights the support of e 0 on ½D; 1). Using this form for 0 and plugging in its Laplace transform into the solution for the mean in Equation (34), one estimates p t < ð1 þ DÞe ÀsDÁt for t ) 1. This result generalizes the exponentially fast approach of stationarity that is known, for example, from the discrete Prisoner's dilemma in evolutionary game theory (Nowak et al., 2004;Traulsen et al., 2005;Melbinger et al., 2010;Assaf et al., 2013).
In total, p t vanishes exponentially fast if and only if the production degree at the smallest production degree is separated by a gap D from all other production degrees that are present in the population. On the other hand, if the lowest production degree is part of an interval with continuously distributed production degrees (that is, D ¼ 0), p t decreases algebraically slowly.
With sense-and-response (l > 0): Homogeneous stationary distributions For the case with sense-and-response through quorum sensing, l > 0, one obtains from Equation (1) or from the cumulant equations (30) that stationary production distributions are given by d-peaks as: In other words, fixed points of the response function give rise to homogeneous stationary distributions. Whether these stationary distributions are stable against small perturbations around stationarity depends on the stability of the fixed points (see linear stability analysis of homogeneous stationary distributions below). Whether they are approached for long times does not only depend on the stability of the fixed points, but also on the initial distribution, the response function, and the value of l (see heterogeneous stationary distributions).

Linear stability analysis of homogeneous stationary distributions
Here, we supplement the statements from the main text on the stability of homogeneous stationary distributions in the linear approximation around stationarity if sense-and-response is present (l > 0). For the sake of simplicity and feasibility, we carry out the stability analysis in the space of cumulants. To this end, we define the vector: Mðu; tÞ ¼ M 0 ðuÞe Àf 0 t þ e up 0 1 À e Àf 0 t ; which yields after back-transformation: ðp; tÞ ¼ yðtÞ 0 ðpÞ þ ð1 À yðtÞÞdðp À p 0 Þ ; with yðtÞ ¼ expðÀf 0 tÞ : The initial production distribution 0 decays exponentially fast on a time scale that is set by the average initial fitness in the population f 0 , whereas a singular probability mass at the initial mean production degree p 0 builds up concomitantly due to sense-and-response through quorum sensing. The population approaches the stationary distribution ¥ ðpÞ ¼ dðp À p 0 Þ exponentially fast.
With sense-and-response (l ¼ 1=2) and polynomial response function: Divergence of time scales at bifurcations of parameters of the response function For l > s=2, the approach of stationarity is typically exponentially fast. However, upon fine-tuning parameters of the response function one observes an algebraically slow approach of stationarity. We exemplify this qualitative change in the temporal evolution by setting the response probability to l ¼ 1=2 and by considering the following nonlinear response function, see Appendix 1-figure 3 (for the sake of readability, we label the argument of R by p instead of hpi): RðpÞ ¼ p þ A Á pðp À ðp cr À ÞÞðp À p cr Þðp À ðp cr þ ÞÞðp À 1Þ ; with some real constant A > 0. The chosen response function (48) is a polynomial of fifth order with Rð0Þ ¼ 0 and Rð1Þ ¼ 1, and parameter 0 < p cr < 1, which is set to p cr ¼ 1=2 in Appendix 1figure 3. The bifurcation parameter 0 minðp cr ; 1 À p cr Þ controls a supercritical pitchfork bifurcation of the response function (48) at p Ã ¼ p cr : Whereas p Ã ¼ 0 and p Ã ¼ 1 are unstable fixed points for all , the fixed points at p Ã ¼ p cr AE are stable for > 0 and merge with p Ã ¼ p cr for ¼ 0. The fixed point p Ã ¼ p cr is unstable for > 0 and is a three-fold degenerate, stable fixed point for ¼ 0, see Appendix 1-figure 3A,B.
For l ¼ 1=2 and upon plugging in the explicit form of the response function (48), the temporal evolution equation of the mean (30) is given by the ODE: q t C 1 ¼ Að1 À sC 1 ÞC 1 ðC 1 À ðp cr À ÞÞðC 1 À p cr ÞðC 1 À ðp cr þ ÞÞðC 1 À 1Þ ; with initial condition C 1 ðt ¼ 0Þ ¼ p 0 . From integrating this temporal evolution equation, one obtains the implicit solution for the mean p ¼ C 1 as: The sum is performed over all non-degenerate fixed points of the right hand side of the equation for the mean (49), that is over the roots p Ã 2 f0; p cr À ; p cr ; p cr þ ; 1; 1=sg of both the response function (48) and the mean fitness f t ¼ 1 À sp t . The coefficients a p Ã arise from the partial fraction decomposition with a pcr;pcrAE~O ð1= 2 Þ and a 0;1;1=s~O ð 0 Þ. Therefore, one concludes that: jp t À p ¥ j~e Àt=a ; for > 0 ; for large times and with a decay constant a that diverges as the bifurcation is approached as a~1= 2 . In other words, stationarity is approached exponentially fast when all fixed points of the response function (48) are non-degenerate, see Appendix 1- figure 3D inset. Which of the two stable fixed points p Ã ¼ p cr AE constitutes the stationary distribution ¥ ðpÞ ¼ dðp À p Ã Þ depends on the initial distribution (and demographic fluctuations of the initial dynamics in the stochastic process). The prediction that the decay constant t diverges as the bifurcation of the response function is approached ( ! 0) is in good agreement with numerical simulations of the stochastic process, see Appendix 1- figure 3D.
In contrast to the exponentially fast approach away from the bifurcation, stationarity is approached algebraically slowly at the bifurcation of the nonlinear response function, that is, for ¼ 0. Since the stable fixed point p Ã ¼ p cr is three-fold degenerate, one finds by integration of Equation (50) the implicit solution for the mean as: In addition to the sum over the non-degenerate fixed points (p Ã 6 ¼ p cr ), a second sum accounts for the degeneracy z ¼ 3 of the fixed point p cr , which is reflected by the singularities in the integrand up to order z. Consequently, the mean production approaches its stationary value as: jp t À p ¥ j~t À1=n ; for ¼ 0 ; for large times with critical exponent n ¼ z À 1 ¼ 2, that is À1=n ¼ À1=2. Appendix 1- figure  3C shows the excellent agreement of our theoretical predictions with numerical simulations of the stochastic process for the algebraically slow approach of stationarity at the bifurcation.
The temporal evolution equation (56) for low has the form of the continuous replicator equation (see Equation (1) with l ¼ 0) with renormalized selection strength sð1 À 2lÞ.
Following the analysis that resulted in Equation (32), the solution for low is given by: low ðp; tÞ ¼ low;0 ðpÞe Àsð1À2lÞtp =L½ low;0 ðsð1 À 2lÞtÞ ; with low;0 ðpÞ ¼ low ðp; t ¼ 0Þ ; if l 1=2. As shown in the main text, the condition l 1=2 is consistent with the condition for the upper threshold of the response probability l s=2 < 1=2, above which heterogeneous stationary distributions cannot occur. For the mean p low;t , one obtains: p low;t ¼ Àq v ln L½ low;0 ðvÞj v¼sð1À2lÞt : In other words, low approaches a stationary d-distribution: low ðp; t ! ¥Þ ¼ low;¥ ðpÞ ¼ dðp À p low Þ ; with p low ¼ p low;¥ ¼ minðsuppð low;0 ÞÞ ¼ minðsuppð 0 ÞÞ : The temporal evolution equation (57) for high has a similar form as the original mean-field equation (1): it involves the sense-and-response term with prefactor 2l, and the replicator term with prefactor 1 À 2l. The sense-and-response term, however, couples to the full production distribution through the argument Rðp t Þ in the d-function and the prefactor ð1 À sp t Þ=ð1 À yðtÞÞ, whereas the replicator term does not couple to low or y. Equation (57) is most suitably analyzed in the space of moment and cumulant generating functions with: M high ðu; tÞ :¼ Z 1 0 dp e up high ðp; tÞ ; and C high ðu; tÞ :¼ ln M high ðu; tÞ À Á ; u 2 ðÀ¥; ¥Þ : The moments and cumulants of high are obtained as M high;k ðtÞ :¼ q k u M high ðu; tÞj u¼0 and C high;k ðtÞ : ¼ q k u C high ðu; tÞj u¼0 for k ! 1. With this notation, it is p high;t ¼ M high;1 ðtÞ ¼ C high;1 ðtÞ. By applying these transformations to the temporal evolution equation (57) of high , one obtains: q t M high ðu; tÞ ¼ Àð1 À 2lÞs q u M high ðu; tÞ À M high;1 ðtÞM high ðu; tÞ À Á þ 2l 1 À sp t 1 À yðtÞ e uRðp t Þ À M high ðu; tÞ ; q t C high ðu; tÞ ¼ Àð1 À 2lÞs q u C high ðu; tÞ À C high;1 ðtÞ À Á þ 2l 1 À sp t 1 À yðtÞ e uRðp t Þ e ÀChighðu;tÞ À 1 ; in which the coupling of high to low and y is apparent explicitly through the occurrence of the factor 1 À yðtÞ and implicitly through the occurrence of p t ¼ yðtÞp low;t þ ð1 À yðtÞÞp high;t . The corresponding equations of motion for the first three cumulants are, thus, obtained as: