Ecological feedback in quorumsensing microbial populations can induce heterogeneous production of autoinducers
 Cited 3
 Views 848
 Annotations
 Article
 Figures and data
 Jump to
Abstract
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 quorumsensing microbial populations has remained elusive. In our theoretical model, cells synthesize and secrete autoinducers into the environment, upregulate their production in this selfshaped environment, and nonproducers 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.
https://doi.org/10.7554/eLife.25773.001eLife 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 autoinducers – known 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 “bethedging” 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 selfshaped 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 quorumsensing.
https://doi.org/10.7554/eLife.25773.002Introduction
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 Gramnegative 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., 2015; Papenfort and Bassler, 2016; Pradhan and Chatterjee, 2014). To implement all of these functions and behaviors, a microbial population needs to dynamically selfregulate 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árcamoOyarce et al., 2015; Grote et al., 2015). For example, during the growth of L. monocytogenes under wellmixed 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árcamoOyarce et al., 2015). The stable coexistence of different phenotypes in one population may serve the division of labor or act as a bethedging 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 upregulated through quorumsensing in single cells (Garmyn et al., 2011, Garmyn et al., 2009; Waters and Bassler, 2005). From an experimental point of view it is often not known, however, whether autoinducer synthesis is upregulated 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érezVelázquez et al., 2016; Goryachev et al., 2005; Dockery and Keener, 2001). For bistable regulation, cellular autoinducer synthesis is upregulated above a threshold value of the autoinducer concentration in the population, whereas it is downregulated below the threshold ('allornone' 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 upregulated for all autoinducer concentrations (monostable upregulation), 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 quorumsensing microbial populations. At the same time, the overall autoinducer level in the environment is robustly selfregulated, so that further quorumsensing functions such as virulence or bioluminescence can be triggered. We studied the collective behavior of a stochastic manyparticle model of quorum sensing, in which cells produce autoinducers to different degrees and secrete them into the wellmixed environment. Production of large autoinducer molecules (for example oligopeptides) and accompanied gene expression are assumed to reduce fitness such that nonproducers reproduce faster than producing cells. Moreover, it is assumed that quorum sensing enables upregulation 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 meanfield equation (1) from the microscopic stochastic manyparticle 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 longrange 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).
Setup of the quorumsensing model
We now introduce the quorumsensing model for a wellmixed population of $N$ individuals (Figure 1). The phenotype of each individual $i=1,\mathrm{\dots},N$ is characterized by its production degree ${p}_{i}\in [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 nonproducer, and ${p}_{i}=1$ denotes a full producer.
The state of the population $\mathbf{p}=({p}_{1},\dots ,{p}_{N})$ changes stochastically (Figure 1A): An individual $i$ reproduces with rate ${\varphi}_{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 $\varphi}_{i}=\varphi ({p}_{i})=1s{p}_{i$. The selection strength $0\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}s\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$ scales the fitness difference with respect to the nonproducing phenotype ($\varphi (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 setup of frequencydependent 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 nonproducers only, $p}_{i}\in \{0,1\$. The mathematical setup of the wellknown 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 senseandresponse mechanism through quorum sensing, which is implemented as follows. After reproduction, both offspring individuals sense the average production level of autoinducers $\u27e8p\u27e9=1/N{\sum}_{i}{p}_{i}$ in the wellmixed population. With probability $\lambda $, they independently adopt the value $R(\u27e8p\u27e9)\in [0,1]$ as their production degree in response to the sensed environmental cue $\u27e8p\u27e9$, whereas they retain the ancestor’s production degree with probability $1\lambda $ through nongenetic inheritance. In an experimental setting, the response probability $\lambda $ 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(\u27e8p\u27e9)$ 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 $\u27e8p\u27e9$ and adjustment of the individual production degree to $R(\u27e8p\u27e9)$ 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érezVelázquez et al., 2016; Goryachev et al., 2005; Dockery and Keener, 2001). For a bistable response function, cellular production is upregulated above a threshold value of $\u27e8p\u27e9$, whereas it is downregulated below the threshold. For the bistable response function sketched in Figure 1B, both values $\u27e8p\u27e9=0$ and $\u27e8p\u27e9=1$ are stable fixed points. In this work, however, we particularly focus on monostable response functions $R(\u27e8p\u27e9)$ to model microbial quorumsensing systems in which autoinducer synthesis is upregulated 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 $\u27e8p\u27e9=1$ and unstable fixed point at $\u27e8p\u27e9=0$). The senseandresponse mechanism is further discussed in the Discussion section.
From a mathematical point of view, the introduced senseandresponse 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}\in [0,1]$ as opposed to a discrete production space is a technical necessity for the implementation of the quorumsensing model. The coupling of ecological dynamics (given by the average production level of autoinducers $\u27e8p\u27e9$) 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 quorumsensing model that we found and report next are qualitatively robust against noise at all steps; see below.
An ecological feedback can control phenotypic heterogeneity in quorumsensing 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 quorumsensing model as longlived, bimodal states of the population that are dynamically stable; see sketch below and Equation (2).
In the quorumsensing model, ecological dynamics are determined by the average production level of autoinducers, while population dynamical changes are determined by fitness differences between nonproducers 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 generegulatory circuits to control phenotypic heterogeneity (see Discussion).
Results of numerical simulations
The quorumsensing model was numerically simulated by employing Gillespie’s stochastic simulation algorithm (Gillespie, 1976, 1977) for a population size of $N={10}^{4}$ individuals and an exemplary selection strength $s=0.2$, such that $sN\gg 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 quorumsensing 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 $\mathbf{\mathbf{p}}$ over time, and depict the histogram of production degrees and the population average in Figure 2.
First, we studied the stochastic manyparticle process without senseandresponse ($\lambda =0$); see Figure 2A,D and Video 1. In this case, nonproducers always proliferate because they reproduce at the highest rate in the population, which is wellstudied in evolutionary game theory (Taylor and Jonker, 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}_{\mathrm{l}\mathrm{o}\mathrm{w}}\simeq 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.
With quorum sensing ($\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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 ${\u27e8p\u27e9}_{\mathrm{\infty}}={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({\u27e8p\u27e9}_{\mathrm{\infty}})={\u27e8p\u27e9}_{\mathrm{\infty}}={p}^{*}$ through senseandresponse. Thus, all individuals continue to produce with degree ${p}^{*}$ and the state of the population remains homogeneous (unimodal).
Surprisingly, for small response probabilities $\lambda $, 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 $\lambda =0.05$. A monostable response function was chosen with $R(\u27e8p\u27e9)\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}\u27e8p\u27e9$ for all $\u27e8p\u27e9\in (0,1)$ (unstable fixed point at $0$, and stable fixed point at $1$) such that the production degree is always upregulated 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}_{\mathrm{l}\mathrm{o}\mathrm{w}$, and a second in which individuals produce to a higher degree $p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ that is separated from $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ by a gap in the space of production degrees. Only through strong demographic fluctuations can the population reach one of the homogeneous absorbing states (${\u27e8p\u27e9}_{\mathrm{\infty}}=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 quasistationary and longlived. These heterogeneous states arise for a broad class of response 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 setup, 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 quorumsensing 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 $\lambda $ and the selection strength $s$ (Figure 2F).
The establishment of longlived, 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 nonproducers. On the other hand, sensing the population average and accordingly upregulating individual production enables producers to persist. Remarkably, fitness differences and senseandresponse 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 finetuned effect), and the average production level in the population is adjusted by the interplay of the response probability $\lambda $ 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).
Results of mathematical analysis
In the following, the observed longlived states of phenotypic heterogeneity in the quorumsensing model are explained. First, we derived the macroscopic meanfield equation (the autoinducer equation (1)) from the microscopic dynamics of the quorumsensing model. Second, we analyzed this meanfield equation and characterized phenotypic heterogeneity of autoinducer production.
The microscopic dynamics of the quorumsensing model are captured by a memoryless stochastic birthdeath process as sketched in Figure 1. Starting from the microscopic manyparticle stochastic process, we derived a meanfield 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 oneparticle probability distribution the production distribution $\rho $; Figure 2 shows the corresponding histogram numerically obtained from the stochastic manyparticle process. The meanfield equation for $\rho $, which we refer to as the autoinducer equation, is obtained as:
where ${\overline{\cdot}}_{t}$ denotes averaging with respect to $\rho $ 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 senseandresponse term with prefactor $2\lambda $, and the replicator term with prefactor $12\lambda $. Through the replicator term, probability weight at production degree $p$ changes if the fitness $\varphi (p)$ is different from the mean fitness in the population ${\overline{\varphi}}_{t}$ (here $\varphi (p){\overline{\varphi}}_{t}=s(p{\overline{p}}_{t})$). Without quorum sensing ($\lambda =0$), Equation (1) reduces to the wellknown 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 senseandresponse term, on the other hand, encodes the global feedback by which individuals adopt the production degree $R({\overline{p}}_{t})$ upon sensing the average ${\overline{p}}_{t}$ through quorum sensing at rate $2\lambda $. The difference between the current state $\rho $ and the state in which all individuals have this production degree $R({\overline{p}}_{t})$ determines the change in $\rho $ at every production degree. Through the replicator term and the senseandresponse term, the ecological dynamics (average production level ${\overline{p}}_{t}$) are coupled with the dynamics of $\rho $.
We now present our results for the longtime behavior of the autoinducer equation (1). First, the autoinducer equation (1) admits homogeneous stationary distributions. Without quorum sensing ($\lambda =0$), the initially lowest production degree in the population, $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$, constitutes the homogeneous stationary distribution ${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}_{\mathrm{l}\mathrm{o}\mathrm{w}})$, which is attractive for generic initial conditions. With quorum sensing ($\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$), fixed points of the response function ${p}^{*}=R({p}^{*})$ yield homogeneous stationary distributions as ${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}^{*})$, which are attractors of the quorumsensing dynamics (1) for all initial distributions if $\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}s/2$; see analysis below. These homogeneous stationary distributions confirm our observations of homogeneous absorbing states in the quorumsensing 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 longlived heterogeneous states of the population, we decomposed $\rho $ 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):
if the conditions $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}1$ and $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}y\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$ are fulfilled; see Box 1 for an illustration and Appendix 3 for the derivation. The parameter $\beta =2\lambda /s$ quantifies the balance between fitness differences and senseandresponse mechanism through quorum sensing. Heterogeneous stationary distributions (2) are constituted of a probability mass $y$ at the lowproducing degree ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$ and a coexisting $\delta $peak with stationary value $1y$ at a highproducing degree $p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ separated from $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ by a gap. Such heterogeneous stationary distributions have mean ${\overline{p}}_{\mathrm{\infty}}=\beta $ and variance $\mathrm{Var}{(p)}_{\mathrm{\infty}}=\beta (R(\beta )\beta )$. Therefore, the interplay between selection strength $s$ and response probability $\lambda $ 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}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$; generalized bimodal distributions for arbitrary initial distributions ${\rho}_{0}$ are given in Appendix 3.
From the conditions on $p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ and $y$ below Equation (2), one can derive the following conditions on the response function and the value of the response probability $\lambda $ (for given selection strength $s$) for the existence of heterogeneous stationary distributions: (i) The response function needs to be nonlinear with $R({\overline{p}}_{\mathrm{\infty}})={p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\overline{p}}_{\mathrm{\infty}}$; that is, quorum sensing needs to upregulate 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 $\lambda \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\lambda}_{\mathrm{u}\mathrm{p}}=s/2$; that is, to induce phenotypic heterogeneity, cells must respond only rarely to the environmental cue $\overline{p}$. This estimate of an upper bound on $\lambda $ 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 ($\lambda \to 0$ and $\lambda \to {\lambda}_{\mathrm{u}\mathrm{p}}$); see Figure 3C.
For small $\lambda $, the coexistence of the lowproducing and the highproducing peaks in solution (2) is stable due to the balance of fitness differences and senseandresponse 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 meanfield dynamics (1) for a broad range of initial distributions when $\lambda $ 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}_{\mathrm{l}\mathrm{o}\mathrm{w}$, whereas nonlinear response to the environment with probability $\lambda $ pushes probability mass towards the upregulated production degree ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=R({\overline{p}}_{\mathrm{\infty}})$. The gap ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}{p}_{\mathrm{l}\mathrm{o}\mathrm{w}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ ensures that the exponential time scales of selection and senseandresponse 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 manyparticle process arise and are quasistationary. 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 longlived.
In summary, our mathematical analysis explains how phenotypic heterogeneity in the autoinducer production arises when quorum sensing upregulates 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}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$) and a second in which cells produce autoinducers (‘on’ state, ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=R(2\lambda /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 ($\lambda =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 ($\lambda \phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}{\lambda}_{\text{up}}$). Only when response to the environment is rare ($0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\lambda \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\lambda}_{\text{up}}$) can the two phenotypic states, $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, coexist in the population ($0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}y\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$). The transitions from heterogeneous to homogeneous populations are governed by nonequilibrium phase transitions when the response probability changes ($\lambda \to 0$ and $\lambda \to {\lambda}_{\mathrm{u}\mathrm{p}}$). 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 quorumsensing model as a collective phenomenon through an ecological feedback
In this work, we studied a conceptual model for the heterogeneous production of autoinducers in quorumsensing microbial populations. The two key assumptions of our quorumsensing 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 nonproducers reproduce faster than producers. Second, cells sense the average production level of autoinducers in the population and may accordingly upregulate their production through quorum sensing. As a result, not only does the interplay between fitness differences and senseandresponse 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 upregulate their production through quorum sensing (response probability $\lambda $ and response function $R(\u27e8p\u27e9)$ 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 quorumsensing model, it is assumed that the individual's production degree of autoinducers is negatively correlated with its growth rate ($\varphi}_{i}=1s{p}_{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 preprotein. 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 Autoinducer2 (AI2) that is considered as a metabolic byproduct. 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 quorumsensing 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}$HSLproducing strain (Ruparell et al., 2016). Furthermore, a strain producing a longchain AHL (OC${}_{12}$HSL, synthesized via the las operon) showed a reduced fitness in both mono and mixed culture compared with a nonproducing strain. The reduced fitness of AHLproducers 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 quorumsensing model, remains to be explored experimentally.
In the quorumsensing model, even small growth rate differences between producer and nonproducer, which are quantified by the ratio (growth rate of producer) / (growth rate of nonproducer) $=1s$, 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 nonproducing 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 spatiotemporal 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 quorumsensing 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 spatiotemporal organization of the microbial population determines as to what extent microbes sense rather the current average production level or a timeintegrated production of autoinducers, and to what extent they sense rather the global or a local average production level. Our quorumsensing model assumes that autoinducers are uniformly degraded in a wellmixed 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 upregulated at the singlecell level?
Monostable or bistable upregulation 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, upregulation 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 (upregulation 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 upregulated for all autoinducer levels or only above a threshold level. Upregulation 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 upregulation 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 quorumsensing 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 $>\phantom{\rule{thinmathspace}{0ex}}1$) and monostable (for example, a Hill function with Hill coefficient $\le \phantom{\rule{thinmathspace}{0ex}}1$) regulation of autoinducer synthesis – apart from the insight on how regulation proceeds at the molecular level? As the results of our quorumsensing 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érezVelá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 quorumsensing model suggests that an alternative mechanism could explain a heterogeneous production of autoinducers in quorumsensing 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 upregulate their expression with respect to the sensed production level in the population. A thresholdlike, bistable response function does not need to be assumed in the quorumsensing 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 quorumsensing model, individuals respond to the environment with response probability $\lambda $ upon reproduction. The rule that offspring individuals can only respond at reproduction events represents a coarsegrained 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 (${\varphi}_{i}$) in the quorumsensing model. Phenotypic heterogeneity of autoinducer production arises in the quorumsensing 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 senseandresponse 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 $\sim 2\lambda $ and (ii) through growth rate differences at rate $\sim 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 $\beta =2\lambda /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 singlecell level seem most promising to us at present.
Singlecell experiments
Some of the questions raised above may be addressed most effectively with singlecell experiments. For example, it would be desirable to simultaneously monitor, at the singlecell 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 singlecell growth rate. In the context of the quorumsensing model, the results of such singlecell experiments would help to identify the form of the fitness function $\varphi $ and the response function $R$, to quantify the selection strength $s$ and response probability $\lambda $, and to refine the model setup.
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 abovementioned singlecell experiments could elucidate the mechanisms that allow for phenotypic heterogeneity in quorumsensing 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 quorumsensing model presented here is to explain how phenotypic heterogeneity in the autoinducer production arises and how it is controlled in quorumsensing microbial populations. With the current model setup, 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 bethedging 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 generegulatory circuits. Spatiotemporal scales are important for the identified ecological feedback to be of relevance for microbial population dynamics: growth rate differences between producers and nonproducers 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 selfshaped 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 quorumsensing microbial populations.
Materials and methods
Derivation of the autoinducer equation (1)
The microscopic dynamics are captured by a memoryless stochastic birthdeath process (a continuoustime Markov process) as sketched in Figure 1. The state of the population $\mathbf{\mathbf{p}}$ is updated by nongenetic inheritance and senseandresponse through quorum sensing such that at most two individuals $i$ and $j\ne i$ change their production degree at one time. The temporal evolution of the corresponding joint $N$particle probability distribution $P(\mathbf{\mathbf{p}},t)$ is governed by a master equation for the stochastic manyparticle 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 oneparticle probability distribution $\rho}^{(1)}(p,t)=1/N\u27e8\sum _{i}\delta (p{p}_{i}){\u27e9}_{P$ in the spirit of a kinetic theory (Kadar, 2007) starting from the microscopic stochastic dynamics. ${\rho}^{(1)}$ denotes the probability distribution of finding any individual at a specified production degree $p$ at time $t$; the numerically obtained histogram of ${\rho}^{(1)}$ was plotted in Figure 2. The temporal evolution of ${\rho}^{(1)}$ is derived from the master equation, and couples to the reduced twoparticle probability distribution and to the full probability distribution $P$ through quorum sensing. By assuming that correlations are negligible, one may approximate ${\rho}^{(1)}$ by the meanfield distribution $\rho $, which we refer to as the production distribution. The meanfield equation (1) for $\rho $ is derived in Appendix 2 and referred to as the autoinducer equation. Note that Equation (1) conserves normalization of $\rho $, that is, ${\int}_{0}^{1}\mathrm{d}p\text{}{\mathrm{\partial}}_{t}\rho (p,t)=0$.
We also proved that ${\rho}^{(1)}$ converges in probability to $\rho $ as $N\to \mathrm{\infty}$ 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 manyparticle process for large $N$. To show this convergence, we introduced the bounded Lipschitz distance $d$ between $\rho $ and ${\rho}^{(1)}$, applied Grönwall’s inequality to the temporal evolution of $d$, and used the law of large numbers; see (Frey et al., 2017) for details. Similar distance measures and estimates have been used, for example, to prove that the Vlasov equation governs the macroscopic dynamics of the abovementioned 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 ($\lambda =0$), one finds the analytical solution for $\rho $ by applying the method of characteristics to Equation (1) in the space of moment and cumulant generating functions as: $\rho (p,t)={\rho}_{0}(p){e}^{stp}/{\int}_{0}^{1}\mathrm{d}p{e}^{stp}{\rho}_{0}(p)$; see Appendix 3 for details. Thus, the initially lowest production degree in the population, $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$, constitutes the homogeneous stationary distribution ${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}_{\mathrm{l}\mathrm{o}\mathrm{w}})$, which is attractive for generic initial conditions. Only $\delta $peaks at production degrees greater than $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ 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 ${\rho}_{0}$, and exponentially fast if $p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is separated from all greater degrees by a gap in production space; see Appendix 3 and Figure 2D.
With quorum sensing ($\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$), fixed points of the response function ${p}^{*}=R({p}^{*})$ yield homogeneous stationary distributions of the autoinducer equation (1) as ${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}^{*})$. In particular, stable fixed points of the response function (${R}^{\mathrm{\prime}}({p}^{\ast})\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$) constitute homogeneous stationary distributions that are stable up to linear order in perturbations around stationarity. For $\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}s/2$, these distributions are also attractors of the meanfield 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 $\lambda =1/2$, for which one finds the analytical solution as: $\rho (p,t)=y(t){\rho}_{0}(p)+(1y(t))\delta (p{\overline{p}}_{0})$ with $y(t)={e}^{{\overline{\varphi}}_{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 $\lambda =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 longtime behavior of the quorumsensing model changes from heterogeneous to homogeneous populations as the response probability $\lambda $ vanishes or reaches the upper threshold $\lambda}_{\text{up}$ while the selection strength $s$ is kept fixed. For small response probabilities, $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\lambda \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\lambda}_{\text{up}}$, the heterogeneous stationary distributions of the autoinducer equation (1) explain the longlived, heterogeneous states of the stochastic quorumsensing process. The coexisting $\delta $peaks at the lowproducing and highproducing degree in the heterogeneous stationary distribution are separated by a gap in production space, which gives rise to the nonvanishing variance $\mathrm{Var}{(p)}_{\mathrm{\infty}}$ in the phase of heterogeneity (Figure 3C). As $\lambda \to {\lambda}_{\text{up}}$, the gap closes, ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}\to R({p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}})$, and $y\to 0$, such that a homogeneous stationary distribution with $\mathrm{Var}{(p)}_{\mathrm{\infty}}=0$ is recovered in a continuous transition. This nonequilibrium phase transition from heterogeneity to homogeneity proceeds without any critical behavior. As $\lambda \to 0$, and under the assumption that $0$ is an unstable fixed point of the response function ($R(0)=0$ and $1\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{R}^{\mathrm{\prime}}(0)$; we further assume $R}^{\mathrm{\prime}}(0)\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\mathrm{\infty$), the gap between the lowproducing and the highproducing peak closes as well because ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}\to 0$. However, $y$ does not approach $1$, but the value $11/{R}^{\mathrm{\prime}}(0)\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$. The probability weight at the lowproducing mode jumps by the value $1/{R}^{\prime}(0)$ and the homogeneous stationary distribution with $\mathrm{Var}{(p)}_{\mathrm{\infty}}=0$ is recovered in a discontinuous transition. Therefore, a discontinuous phase transition in the space of stationary probability distributions governs the longtime dynamics of the autoinducer equation (1) from heterogeneity to homogeneity as the response probability $\lambda $ vanishes (for fixed selection strength $s$).
Appendix 1
Appendix 2
From a microscopic description to a macroscopic description of the quorumsensing model
Description of the microscopic dynamics: Master equation of the stochastic manyparticle process
To describe the temporal evolution of the population, we introduced the joint $N$particle probability distribution $P(\mathbf{\mathbf{p}},t)$. The value $P(\mathbf{p},t)\mathrm{d}{p}_{1}\dots \mathrm{d}{p}_{N}$ denotes the joint probability of finding the first individual with a production degree in the interval $[{p}_{1},{p}_{1}+\mathrm{d}{p}_{1}]$, the second individual with a production degree in the interval $[{p}_{2},{p}_{2}+\mathrm{d}{p}_{2}]$, and so on at time $t$. The stochastic dynamics are captured by a coupled birthdeath process (continuoustime Markov process) as described in the main text and in Figure 1 of the main text. An individual $i$ reproduces randomly after a time that is exponentially distributed with rate ${\varphi}_{i}$, which we refer to as the individual’s fitness in the main text. One update step involves reproduction, senseandresponse through quorum sensing, and nongenetic inheritance such that at most two individuals $i$ and $j\ne i$ change their production degree at one time. We denote the state of the population before the update step as ${\stackrel{~}{\mathbf{p}}}^{i,j}=({p}_{1},\dots ,{p}_{i1},\stackrel{~}{{p}_{i}},{p}_{i+1}\dots ,\stackrel{~}{{p}_{j}},\dots ,{p}_{N})$; the production degrees of individual $i$ and $j$, which might change during the update step, are labeled with a tilde. For the sake of readability, we do not distinguish notationally between a random variable and the value that this random variable attains; both are labeled with the same symbol. The master equation for the joint $N$particle probability distribution $P$ for the individuals’ production degrees $\mathbf{p}=({p}_{1},\dots ,{p}_{N})$ at time $t$ can be written as (Gardiner, 2009; Van Kampen, 2007; Weber and Frey, 2017):
with reproduction rate of individual $i$ (fitness) given by (and selection strength $0\le s\text{}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 nongenetic inheritance and senseandresponse (see detailed description below Equation (6) and Equation (7)). The initial condition to the master equation (3) is given as $P(\mathbf{p},t=0)={p}_{0}(\mathbf{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 $\mathbf{\mathbf{p}}$ at time $t$. Loss terms occur when the population is in state $\mathbf{\mathbf{p}}$ and an individual reproduces. The probability of finding the population in this state is given by $P(\mathbf{\mathbf{p}},t)$. Individual $i$ is selected for reproduction at rate ${\varphi}_{i}(\mathbf{\mathbf{p}})$ and splits into two offspring individuals, and a different individual $j\ne i$ is removed with probability $1/(N1)$ at the same time (random death). Gain terms involve all events that take the population from an arbitrary state ${\stackrel{~}{\mathbf{\mathbf{p}}}}^{i,j}$ to state $\mathbf{\mathbf{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 nongenetic inheritance and senseandresponse through quorum sensing, and are given as:
We abbreviate $\u27e8{\stackrel{~}{p}}^{i,j}\u27e9=1/N{\sum}_{k}{({\stackrel{~}{\mathbf{\mathbf{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(\u27e8{\stackrel{~}{p}}^{i,j}\u27e9)$) with probability $\lambda $ as the updated production degree, and the second summand accounts for the nongenetic inheritance of the production degree from the ancestor $i$ ($p}_{i/j$ attains the value ${\stackrel{~}{p}}_{i}$) with probability $1\lambda $. Note that for the transition probability ${A}_{j}$, also $p}_{j$ attains the value ${\stackrel{~}{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 (${\stackrel{~}{p}}_{i}={\stackrel{~}{p}}_{j}$) and both offspring individuals retain the production degree from their ancestor $i$ ($p}_{i}={p}_{j}={\stackrel{~}{p}}_{i$, that is, both offspring individuals do not update their production through senseandresponse). Such events do not change the state of the population (${\stackrel{~}{\mathbf{\mathbf{p}}}}^{i,j}=\mathbf{\mathbf{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).
The master equation (3) conserves normalization of $P$ because ${\mathrm{\partial}}_{t}{\int}_{[0,1{]}^{N}}\mathrm{d}\mathbf{p}\text{}P(\mathbf{p},t)=0$; see analysis below.
Coarsegrained description: Reduced oneparticle probability distribution
The reduced oneparticle probability distribution ${\rho}^{(1)}$ is defined as:
and agrees with the marginal probability distribution for the production degree of the first individual ${P}^{(1)}$. The equality between the normalized reduced oneparticle distribution ${\rho}^{(1)}$ and the oneparticle probability distribution ${P}^{(1)}$ follows from the symmetry of $P$ with respect to permutation of identical (that is indistinguishable) individuals (Kadar, 2007).
We also define the more general reduced $n$particle probability distribution:
which agrees with the marginal probability distribution for the production degrees of the first $n$ individuals, ${P}^{(n)}$. In particular, one also has $P(\mathbf{\mathbf{p}},t)={P}^{(N)}(\mathbf{\mathbf{p}},t)={\rho}^{(N)}(\mathbf{\mathbf{p}},t)$.
Towards the macroscopic dynamics: Temporal evolution of the reduced oneparticle probability distribution
In the following we show that the temporal evolution equation of the reduced oneparticle probability distribution is obtained from the master equation (3) as:
To derive the temporal evolution equation for ${\rho}^{(1)}$, we specify the production degree of one particular individual (here $p}_{1$), and integrate out the production degrees of the other $N1$ individuals in the master equation (3):
For the loss term, we split the sum ${\sum}_{i}{\ast}_{(i)}$ into two contributions:
and deal with both contributions separately to obtain:
For the gain term, we split up the sum ${\sum}_{i}{\sum}_{j\ne i}{\ast}_{(i,j)}$ that occurs in the master equation (3) into three terms as follows:
We also introduce the notation $\mathrm{d}{\stackrel{~}{\mathbf{p}}}^{i,j,\hat{k}}:=\mathrm{d}{p}_{1}\mathrm{d}{p}_{2}\dots \mathrm{d}{\stackrel{~}{p}}_{i}\dots \mathrm{d}{\stackrel{~}{p}}_{j}\dots \mathrm{d}{\hat{p}}_{k}\dots \mathrm{d}{p}_{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:
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 $\delta $functions of the transition probabilities was carried out as well, for example, ${\int}_{0}^{1}\mathrm{d}{p}_{j}\text{}{A}_{j}({\stackrel{~}{\mathbf{p}}}^{i,j};i)=1$):
Making use of the fact that $P$ is symmetric with respect to permutation of individuals (individuals are identical), carrying out possible integrals over $\delta $functions, plugging in the explicit form of the fitness function (4), and relabeling variables, one obtains for the gain term:
Combining loss terms $I}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}$ and gain terms $I}_{\mathrm{g}\mathrm{a}\mathrm{i}\mathrm{n}$ leads to the result for the equation of motion of the reduced oneparticle probability distribution ${\rho}^{(1)}$ that is given in Equation (12).
Heuristic derivation of the macroscopic dynamics: Meanfield approximation
Upon assuming that correlations are negligible, one may approximate ${\rho}^{(1)}$ by its meanfield approximation $\rho $, which we refer to as the production distribution. As described in the main text, the temporal evolution equation for ${\rho}^{(1)}$ serves as a suitable starting point to guess the meanfield equation for $\rho $, which is the meanfield approximation of ${\rho}^{(1)}$. Thus, we naively approximate ${\rho}^{(1)}\approx \rho $ and ${\rho}^{(N)}\approx {\prod}^{N}\rho $. From the temporal evolution of ${\rho}^{(1)}$ in Equation (12), the meanfield equation for $\rho $ is suggested as:
where ${\overline{\cdot}}_{t}$ denotes averaging with respect to $\rho $ at time $t$. Further collection of terms yields the meanfield equation (1) in the main text:
with initial condition $\rho (p,t=0)={\rho}_{0}(p)$, $\varphi (p)=1sp$, ${\overline{\varphi}}_{t}=1s{\overline{p}}_{t}$, and ${\overline{p}}_{t}={\int}_{0}^{1}\mathrm{d}pp\rho (p,t)$. Alternatively, this meanfield equation can also be written as:
We emphasize that the meanfield 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]\to$ ℝ, $g$ smooth) and $\rho $ is interpreted as a linear functional on the space of these observables. This way, $\rho $ 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 measuretheoretic notation in this manuscript.
The proof that ${\rho}^{(1)}$ converges in probability to $\rho $ as $N\to \mathrm{\infty}$ for any finite time if initial correlations are not too strong will be presented in a forthcoming publication (Frey et al., 2017).
Appendix 3
Analysis of the meanfield equation of the quorumsensing model (autoinducer equation)
Meanfield equation for moment and cumulantgenerating functions
The momentgenerating function $M(u,t)$ for the production degree $p$, which is the random variable of interest, and its corresponding cumulantgenerating function $C(u,t)$ are defined as:
with argument $u\in (\mathrm{\infty},\mathrm{\infty})$ at time $t$. The momentgenerating function $M$ is the (onesided) Laplace transform $\mathcal{L}$ of $\rho $ with negative argument at time $t$. Moments and cumulants of the degree distribution $\rho $ are obtained as:
For the mean production, that is, for the expectation value of the production distribution, it holds that $\overline{p}={M}_{1}={C}_{1}$ and the variance is given by $\mathrm{V}\mathrm{a}\mathrm{r}(p)=\overline{{p}^{2}}{\overline{p}}^{2}={M}_{2}{M}_{1}^{2}={C}_{2}$. By applying transformations (23, 24) to the meanfield equation (1) and plugging in the form of the fitness function in Equation (4), one obtains:
Solution strategy for the moment and cumulantgenerating functions: Method of characteristics
This meanfield equation in moment/cumulant space (26) is more conveniently written as a semilinear partial differential equation (PDE) of first order in $t$ and $u$, for example for $C$:
with $F(C,u,t):=(12\lambda )s{C}_{1}(t)+2\lambda (1s{C}_{1})\left({e}^{uR({C}_{1}(t))}{e}^{C(u,t)}1\right)$ and initial condition $C(u,t=0)={C}_{0}(u)$. This PDE admits the straight lines $r(u,t)=u(12\lambda )st$ as characteristics. Restricted to these characteristic curves, the PDE reduces to a nonlinear ordinary differential equation (ODE) of first order in time for $z(r,t)=C(u(r,t),t)$:
with initial condition $z(r,t=0)=C(u(r,t=0),t=0)={C}_{0}(r)$. The solution for the cumulantgenerating function is then obtained from the solution of the above ODE as $C(u,t)=z(r(u,t),t)=z(u(12\lambda )st,t)$. For the two cases $\lambda =0$ and $\lambda =1/2$ with linear response function, an insightful, analytical solution of the meanfield equation (1) for the production distribution for all times $t$ was found this way; see below.
Moment and cumulant equations
A different approach to characterize the dynamics of the quorumsensing model is to analyze the equations of motions for the moments and cumulants. The moment equations are derived from Equation (26) by applying the definition of the moments (25), which yields for $k\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}1$,
The equations for the first three cumulants are obtained as,
For Figure 2E of the main text, the cumulant equations (30) were numerically integrated after applying a Gaussian approximation, that is a cumulant closure with ${C}_{i}(t)=0$ for $i\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}3$ and all $t$, and plotted for ${\overline{p}}_{t}={C}_{1}(t)$.
Without senseandresponse ($\lambda =0$): Analytical solution and approach of the homogeneous stationary distribution of nonproducers
For the case without senseandresponse through quorum sensing, $\lambda =0$, it is readily seen from Equation (1) that stationary production distributions are given by $\delta $peaks as ${\rho}_{\mathrm{\infty}}(p):=\rho (p,t\to \mathrm{\infty})=\delta (p{p}_{\mathrm{l}\mathrm{o}\mathrm{w}})$ for all ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}\in [0,1]$. However, the distribution with solely nonproducers, ${p}_{\text{low}}=0$, is the only asymptotically stable solution of the meanfield equation (1); see below.
When senseandresponse is absent, the analytical solution of the meanfield equation (1) for $\rho $ can be obtained by applying the method of characteristics to Equation (27) as outlined above. The implicit solution is given by:
as the temporal average of the mean production ${\overline{p}}_{t}$. Backtransformation and exploiting normalization of $\rho $ yields:
For example, if the initial production distribution ${\rho}_{0}$ is a uniform distribution on $[0,1]$, $\rho $ evolves in time as $\rho (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\to 0$, this time scale diverges and, hence, the stationary distribution,
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 ${\overline{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 ${\rho}_{0}$ via its Laplace transform $\mathcal{L}[{\rho}_{0}]$. For the asymptotic behavior of Laplace transforms it is known that if ${\rho}_{0}(p)\sim {p}^{\mu}$ as $p\to 0$ with $\mu \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$, then $\mathcal{L}[{\rho}_{0}](v)\sim 1/{v}^{(\mu +1)}$ for $v\phantom{\rule{thinmathspace}{0ex}}\gg \phantom{\rule{thinmathspace}{0ex}}1$ (Doetsch, 1976). Therefore, it follows that the mean evolves in time as ${\overline{p}}_{t}\sim 1/t$ for $t\phantom{\rule{thinmathspace}{0ex}}\gg \phantom{\rule{thinmathspace}{0ex}}1$ if the initial production distribution is a continuous probability density with nonvanishing weight at ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$ (chosen for simplicity as the lowest production degree). The condition that the exponent satisfies $\mu \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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 $\mathrm{Var}(p)(t)\sim 1/{t}^{2}$ for $t\phantom{\rule{thinmathspace}{0ex}}\gg \phantom{\rule{thinmathspace}{0ex}}1$.
In contrast, if the lowest production degree is separated from all other degrees in the population by a gap $\mathrm{\Delta}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ in production space, mean and variance approach their stationary value exponentially fast at a time scale set by $\mathrm{\Delta}$. To see this qualitative difference in the approach of stationarity, we consider an initial probability distribution with probability mass ${y}_{0}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ at degree ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=0$ (chosen again for simplicity) and a remainder probability distribution $\stackrel{~}{\rho}}_{0$ with support on $[\mathrm{\Delta},1]$: ${\rho}_{0}(p)={y}_{0}\delta (p)+(1{y}_{0}){\stackrel{~}{\rho}}_{0}(p){\mathbb{I}}_{[\mathrm{\Delta},1]}(p)$ (here $\mathbb{I}}_{[\mathrm{\Delta},1]$ denotes the indicator function, which takes value 1 on the interval $[\mathrm{\Delta},1]$ and 0 otherwise, and highlights the support of ${\stackrel{~}{\rho}}_{0}$ on $[\mathrm{\Delta},1]$). Using this form for ${\rho}_{0}$ and plugging in its Laplace transform into the solution for the mean in Equation (34), one estimates $\overline{p}}_{t}\phantom{\rule{thinmathspace}{0ex}}\lesssim \phantom{\rule{thinmathspace}{0ex}}(1+\mathrm{\Delta}){e}^{s\mathrm{\Delta}\cdot t$ for $t\phantom{\rule{thinmathspace}{0ex}}\gg \phantom{\rule{thinmathspace}{0ex}}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, ${\overline{p}}_{t}$ vanishes exponentially fast if and only if the production degree at the smallest production degree is separated by a gap $\mathrm{\Delta}$ 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, $\mathrm{\Delta}=0$), ${\overline{p}}_{t}$ decreases algebraically slowly.
With senseandresponse ($\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$): Homogeneous stationary distributions
For the case with senseandresponse through quorum sensing, $\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, one obtains from Equation (1) or from the cumulant equations (30) that stationary production distributions are given by $\delta $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 $\lambda $ (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 senseandresponse is present ($\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}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:
which is at stationarity (see Equation (35)):
With this notation, the equations of motion for the cumulants of $\rho $ are given as follows:
Here, the functions ${F}_{i}$ for $i\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}0$ are defined by the right hand side of the cumulant equations (30). Upon introducing the distance $\mathrm{\Delta}\mathbf{\mathbf{C}}$ to the stationary vector ${\mathbf{\mathbf{C}}}_{\mathrm{\infty}}$, that is $\mathrm{\Delta}\mathbf{\mathbf{C}}=\mathbf{\mathbf{C}}{\mathbf{\mathbf{C}}}_{\mathrm{\infty}}$, one obtains the temporal behavior of $\mathrm{\Delta}\mathbf{\mathbf{C}}$ as:
with Jacobian ${J}_{ij}({\mathbf{\mathbf{C}}}_{\mathrm{\infty}})={\frac{\partial {F}_{i}(\mathbf{\mathbf{C}})}{\partial {C}_{j}}}_{\mathbf{\mathbf{C}}={\mathbf{\mathbf{C}}}_{\mathrm{\infty}}}$, whose entries are obtained after some algebra as:
The eigenvalues of the upper triangular matrix $J$ determine the stability of the stationary distribution up to linear order in perturbations at the level of cumulants around stationarity. Because of the upper triangular structure of the Jacobian $J$, its eigenvalues are given by the diagonal entries of $J$:
Thus, local stability of homogeneous stationary distributions (${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}^{*})$ with $R({p}^{*})={p}^{*}$) is determined by the stability of the fixed points, that is whether ${R}^{\prime}({p}^{*})$ is less or greater than 1.
In total, homogeneous stationary distributions are unstable up to linear order in perturbations at the level of cumulants around stationarity if ${R}^{\mathrm{\prime}}({p}^{\ast})\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$. In other words, stationary distributions located at a fixed point ${p}^{*}$ are linearly unstable if ${p}^{*}$ is an unstable fixed point of the response function (${R}^{\mathrm{\prime}}({p}^{\ast})\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$). On the other hand, linear stability of the response function at ${p}^{*}$ (${R}^{\mathrm{\prime}}({p}^{\ast})\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}1$) yields linearly stable homogeneous stationary distributions located at ${p}^{*}$.
With senseandresponse ($\lambda =1/2$) and linear response function ($R(p)=p$): Analytical solution and approach of homogeneous stationary distribution
For the choice of linear response function ($R(p)=p$, that is, ${R}^{\prime}(p)=1$ for all $p\in [0,1]$) and $\lambda =1/2$, the mean remains constant in time (see Equation (30)). Furthermore, one obtains the analytical solution of the meanfield equation (1) by applying the method of characteristics (most conveniently in the space of moment generating functions) as:
which yields after backtransformation:
The initial production distribution ${\rho}_{0}$ decays exponentially fast on a time scale that is set by the average initial fitness in the population ${\overline{\varphi}}_{0}$, whereas a singular probability mass at the initial mean production degree ${\overline{p}}_{0}$ builds up concomitantly due to senseandresponse through quorum sensing. The population approaches the stationary distribution ${\rho}_{\mathrm{\infty}}(p)=\delta (p{\overline{p}}_{0})$ exponentially fast.
With senseandresponse ($\lambda =1/2$) and polynomial response function: Divergence of time scales at bifurcations of parameters of the response function
For $\lambda \phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}s/2$, the approach of stationarity is typically exponentially fast. However, upon finetuning 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 $\lambda =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 $\u27e8p\u27e9$):
with some real constant $A\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$. The chosen response function (48) is a polynomial of fifth order with $R(0)=0$ and $R(1)=1$, and parameter $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{p}_{\mathrm{c}\mathrm{r}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$, which is set to ${p}_{\mathrm{c}\mathrm{r}}=1/2$ in Appendix 1—figure 3. The bifurcation parameter $0\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}\u03f5\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}\text{min}({p}_{\mathrm{c}\mathrm{r}},1{p}_{\mathrm{c}\mathrm{r}})$ controls a supercritical pitchfork bifurcation of the response function (48) at $p}^{\ast}={p}_{\mathrm{c}\mathrm{r}$: Whereas ${p}^{*}=0$ and ${p}^{*}=1$ are unstable fixed points for all $\u03f5$, the fixed points at ${p}^{\ast}={p}_{\mathrm{c}\mathrm{r}}\pm \u03f5$ are stable for $\u03f5\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and merge with $p}^{\ast}={p}_{\mathrm{c}\mathrm{r}$ for $\u03f5=0$. The fixed point $p}^{\ast}={p}_{\mathrm{c}\mathrm{r}$ is unstable for $\u03f5\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and is a threefold degenerate, stable fixed point for $\u03f5=0$, see Appendix 1—figure 3A,B.
For $\lambda =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:
with initial condition ${C}_{1}(t=0)={\overline{p}}_{0}$. From integrating this temporal evolution equation, one obtains the implicit solution for the mean $\overline{p}={C}_{1}$ as:
The sum is performed over all nondegenerate fixed points of the right hand side of the equation for the mean (49), that is over the roots $p}^{\ast}\in \{0,{p}_{\mathrm{c}\mathrm{r}}\u03f5,{p}_{\mathrm{c}\mathrm{r}},{p}_{\mathrm{c}\mathrm{r}}+\u03f5,1,1/s\$ of both the response function (48) and the mean fitness ${\overline{\varphi}}_{t}=1s{\overline{p}}_{t}$. The coefficients ${\alpha}_{{p}^{*}}$ arise from the partial fraction decomposition with ${\alpha}_{{p}_{\text{cr}},{p}_{\text{cr}}\pm \u03f5}\sim \mathcal{\mathcal{O}}(1/{\u03f5}^{2})$ and ${\alpha}_{0,1,1/s}\sim \mathcal{\mathcal{O}}({\u03f5}^{0})$. Therefore, one concludes that:
for large times and with a decay constant $\alpha $ that diverges as the bifurcation is approached as $\alpha \sim 1/{\u03f5}^{2}$. In other words, stationarity is approached exponentially fast when all fixed points of the response function (48) are nondegenerate, see Appendix 1—figure 3D inset. Which of the two stable fixed points ${p}^{\ast}={p}_{\mathrm{c}\mathrm{r}}\pm \u03f5$ constitutes the stationary distribution ${\rho}_{\mathrm{\infty}}(p)=\delta (p{p}^{*})$ depends on the initial distribution (and demographic fluctuations of the initial dynamics in the stochastic process). The prediction that the decay constant $\tau $ diverges as the bifurcation of the response function is approached ($\u03f5\to 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 $\u03f5=0$. Since the stable fixed point $p}^{\ast}={p}_{\mathrm{c}\mathrm{r}$ is threefold degenerate, one finds by integration of Equation (50) the implicit solution for the mean as:
In addition to the sum over the nondegenerate fixed points ($p}^{\ast}\ne {p}_{\mathrm{c}\mathrm{r}$), a second sum accounts for the degeneracy $z=3$ of the fixed point $p}_{\mathrm{c}\mathrm{r}$, which is reflected by the singularities in the integrand up to order $z$. Consequently, the mean production approaches its stationary value as:
for large times with critical exponent $\nu =z1=2$, that is $1/\nu =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.
With rare senseandresponse ($0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\lambda \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}s/2$): Heterogeneous stationary distributions
To analyze heterogeneous stationary distributions, we decompose the production distribution as follows:
where $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ denote two probability distributions with support on the interval $[0,1]$. Their respective means are denoted as:
such that $\overline{p}}_{t}=y(t){\overline{p}}_{\mathrm{l}\mathrm{o}\mathrm{w},t}+(1y(t)){\overline{p}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},t$; their stationary values are denoted as $\overline{p}}_{\mathrm{l}\mathrm{o}\mathrm{w},\mathrm{\infty}}=:{p}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\overline{p}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},\mathrm{\infty}}=:{p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, respectively. We decompose the initial distribution ${\rho}_{0}(p)={y}_{0}{\rho}_{\mathrm{l}\mathrm{o}\mathrm{w},0}(p)+(1{y}_{0}){\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},0}(p)$ such that $min(\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}({\rho}_{\mathrm{l}\mathrm{o}\mathrm{w},0}))=min(\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}({\rho}_{0}))$. For a numerical integration of the meanfield equation (1) that not only reproduces the stationary distribution, but also the temporal approach towards stationarity, it turns out suitable to choose the following decomposition: ${\rho}_{\mathrm{l}\mathrm{o}\mathrm{w},0}={\rho}_{0},{\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},0}=\delta (\cdot R({\overline{p}}_{0}))$, and ${y}_{0}=1\u03f5$ with $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\u03f5\phantom{\rule{thinmathspace}{0ex}}\lesssim \phantom{\rule{thinmathspace}{0ex}}0.01$.
With decomposition (54), the meanfield equation (1) for $\rho $ can be rewritten in terms of equations for ${\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}},{\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}},$ and $y$ as follows:
We note that the decomposition (54) of $\rho $ with Equations (56–58) is not unique, but this choice of decomposition enables the characterization of heterogeneous stationary distributions and, thus, phenotypic heterogeneity.
The temporal evolution equation (56) for $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ has the form of the continuous replicator equation (see Equation (1) with $\lambda =0$) with renormalized selection strength $s(12\lambda )$. Following the analysis that resulted in Equation (32), the solution for $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ is given by:
if $\lambda \phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}1/2$. As shown in the main text, the condition $\lambda \phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}1/2$ is consistent with the condition for the upper threshold of the response probability $\lambda \phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}s/2\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1/2$, above which heterogeneous stationary distributions cannot occur. For the mean $\overline{p}}_{\mathrm{l}\mathrm{o}\mathrm{w},t$, one obtains:
In other words, $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ approaches a stationary $\delta $distribution:
The temporal evolution equation (57) for $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ has a similar form as the original meanfield equation (1): it involves the senseandresponse term with prefactor $2\lambda $, and the replicator term with prefactor $12\lambda $. The senseandresponse term, however, couples to the full production distribution $\rho $ through the argument $R({\overline{p}}_{t})$ in the $\delta $function and the prefactor $(1s{\overline{p}}_{t})/(1y(t))$, whereas the replicator term does not couple to $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ or $y$. Equation (57) is most suitably analyzed in the space of moment and cumulant generating functions with:
The moments and cumulants of $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ are obtained as $M}_{\text{high},k}(t):={\mathrm{\partial}}_{u}^{k}{M}_{\text{high}}(u,t){}_{u=0$ and $C}_{\text{high},k}(t):={\mathrm{\partial}}_{u}^{k}{C}_{\text{high}}(u,t){}_{u=0$ for $k\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}1$. With this notation, it is ${\overline{p}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},t}={M}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},1}(t)={C}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},1}(t)$. By applying these transformations to the temporal evolution equation (57) of $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, one obtains:
in which the coupling of $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ to $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $y$ is apparent explicitly through the occurrence of the factor $1y(t)$ and implicitly through the occurrence of $\overline{p}}_{t}=y(t){\overline{p}}_{\mathrm{l}\mathrm{o}\mathrm{w},t}+(1y(t)){\overline{p}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},t$. The corresponding equations of motion for the first three cumulants are, thus, obtained as:
At stationarity, it is ${\partial}_{t}y(t)=0$ and $y(t)\equiv {y}_{\mathrm{\infty}}$ with (see Equation (58); recall also that $\overline{p}}_{\mathrm{\infty}}={y}_{\mathrm{\infty}}{p}_{\mathrm{l}\mathrm{o}\mathrm{w}}+(1{y}_{\mathrm{\infty}}){p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$):
Thus, assuming that a stationary value $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{y}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$ exists, it fulfils the selfconsistency relation:
Note that we denoted ${y}_{\mathrm{\infty}}$ simply as $y$ in the main text.
If $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{y}_{\mathrm{\infty}}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}1$ exists, it follows that the stationary solution for $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ can be obtained via Equation (63) in terms of the stationary moment generating function ${M}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},\mathrm{\infty}}(u)={M}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}(u,t\to \mathrm{\infty})$ with:
where the relation between ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}},{p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}},$ and ${y}_{\mathrm{\infty}}$ in Equation (66) was exploited and the definition $p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}={\overline{p}}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},\mathrm{\infty}$ translates into the boundary condition. In total, one obtains $M}_{\text{high},\mathrm{\infty}}(u)={e}^{u{p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}$ with the selfconsistency relation ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}=R({\overline{p}}_{\mathrm{\infty}})$. In other words, $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$ approaches a stationary $\delta $distribution:
For ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=min(\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}({\rho}_{0}))=0$, one recovers from Equations (61, 67, 69) the heterogeneous stationary distribution (2) that was given in the main text.
For Figure 2F of the main text, equations (58, 60, and 65) were numerically integrated with ${C}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},i}(t)=0$ for $i\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}3$ and for all $t$, and initial conditions ${y}_{0}=0.99$, ${\rho}_{\mathrm{l}\mathrm{o}\mathrm{w},0}\sim \mathrm{U}\mathrm{n}\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}(0,1)$, and ${\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h},0}\sim \delta (\cdot R(0.5))$. The choice of initial conditions, however, is not important for the asymptotic behavior, see Appendix 1—figure 1.
Linear stability analysis of heterogeneous stationary distributions
Here, we supplement the statements from the main text on the stability of heterogeneous stationary distributions (2) in the linear approximation around stationarity. 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:
which is at stationarity:
The cumulants of $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ are obtained in the same way as for $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, that is as $C}_{\text{low},k}(t):={\mathrm{\partial}}_{u}^{k}{C}_{\text{low}}(u,t){}_{u=0$ for $k\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}1$ from ${M}_{\text{low}}(u,t):={\int}_{0}^{1}\mathrm{d}p\text{}{e}^{up}{\rho}_{\text{low}}(p,t)$ and ${C}_{\text{low}}(u,t):=\mathrm{l}\mathrm{n}\left({M}_{\text{low}}(u,t)\right)$ for $u\in (\mathrm{\infty},\mathrm{\infty})$. With this notation, the equations of motion for $y(t)$ in Equation (58) and the cumulants of $\rho}_{\mathrm{l}\mathrm{o}\mathrm{w}$ and $\rho}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}$, respectively, are cast into the compact form:
Upon introducing the distance $\mathrm{\Delta}\mathbf{\mathbf{c}}$ to the stationary vector ${\mathbf{\mathbf{c}}}_{\mathrm{\infty}}$, that is $\mathrm{\Delta}\mathbf{\mathbf{c}}=\mathbf{\mathbf{c}}{\mathbf{\mathbf{c}}}_{\mathrm{\infty}}$, one obtains the temporal behavior of $\mathrm{\Delta}\mathbf{\mathbf{c}}$ as follows:
and with Jacobian ${J}_{ij}({\mathbf{\mathbf{c}}}_{\mathrm{\infty}})={\frac{\partial {F}_{i}(\mathbf{\mathbf{c}})}{\partial {c}_{j}}}_{\mathbf{\mathbf{c}}={\mathbf{\mathbf{c}}}_{\mathrm{\infty}}}$, whose entries are obtained after some algebra as:
and,
The eigenvalues of the matrix $J$ determine the stability of the heterogeneous stationary distribution up to linear order in perturbations at the level of cumulants around stationarity. Its eigenvalues are given by:
the two eigenvalues ${\gamma}_{1,2}$ of the $2\times 2$ matrix,
one eigenvalue 0,
and infinitely many pairs of eigenvalues with values 0 and $s(1y)({p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}{p}_{\mathrm{l}\mathrm{o}\mathrm{w}})\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$ (because ${p}_{\mathrm{h}\mathrm{i}\mathrm{g}\mathrm{h}}{p}_{\mathrm{l}\mathrm{o}\mathrm{w}}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and $1y\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ for the considered bimodal distributions).
For simplicity of the discussion, we assume ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=min(\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}({\rho}_{0}))=0$ in the following, and also introduce the parameter $\beta =2\lambda /s$ as in the main text. The two eigenvalues ${\gamma}_{1,2}$ of $\stackrel{~}{J}$ are given by:
Linear stability for small $\lambda $
Under the assumptions $R(0)=0$ and $1\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{R}^{\mathrm{\prime}}(0)\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\mathrm{\infty}$, one checks that for $0\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\lambda \ll 1$ the eigenvalues of the Jacobian $\stackrel{~}{J}$ in Equation (88) are given by:
and, thus, $\mathrm{R}\mathrm{e}({\gamma}_{1,2})\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$ as $\lambda \searrow 0$.
Therefore, for small response probabilities, the heterogeneous stationary distribution (Equation (2) of the main text) is stable up to linear order in perturbations at the level of cumulants around stationarity (here shown under the assumptions ${p}_{\mathrm{l}\mathrm{o}\mathrm{w}}=min(\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}({\rho}_{0}))=0$, $R(0)=0$, and $1\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{R}^{\mathrm{\prime}}(0)\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\mathrm{\infty}$).
Linear stability for the response function $R(\beta )=\beta +\kappa \cdot \mathrm{sin}(\pi \beta )$
Upon choosing the response function $R(\beta )=\beta +\kappa \cdot \mathrm{sin}(\pi \beta )$ with $\beta \in [0,1]$ (that is $\lambda \in [0,s/2]$) and with $\kappa \in [0,1/\pi ]$, one checks that all eigenvalues of the Jacobian $\stackrel{~}{J}$ in Equation (88) have negative real part.
Therefore, for the special choice of the response function that upregulates the cellular autoinducer production for all sensed average productions in the population, all heterogeneous stationary distributions (Equation (2) of the main text) for choices of the parameters $\lambda \in [0,s/2]$ and $\kappa \in [0,1/\pi ]$ are stable up to linear order in perturbations at the level of cumulants around stationarity.
References

1
Stochastic switching as a survival strategy in fluctuating environmentsNature Genetics 40:471–475.https://doi.org/10.1038/ng.110

2
A functional perspective on phenotypic heterogeneity in microorganismsNature Reviews Microbiology 13:497–508.https://doi.org/10.1038/nrmicro3491

3
Heterogeneity in quorum sensingregulated bioluminescence of Vibrio harveyiMolecular Microbiology 73:267–277.https://doi.org/10.1111/j.13652958.2009.06768.x
 4

5
Clustering and relaxation in Hamiltonian longrange dynamicsPhysical Review E 52:2361–2374.https://doi.org/10.1103/PhysRevE.52.2361

6
Extinction of metastable stochastic populationsPhysical Review E 81:021116.https://doi.org/10.1103/PhysRevE.81.021116

7
Cooperation dilemma in finite populations under fluctuating environmentsPhysical Review Letters 111:238101.https://doi.org/10.1103/PhysRevLett.111.238101
 8
 9

10
Stochastic models of evolution in genetics, ecology and linguisticsJournal of Statistical Mechanics: Theory and Experiment 2007:P07018.https://doi.org/10.1088/17425468/2007/07/P07018

11
Microfluidic confinement of single cells of bacteria in small volumes initiates highdensity behavior of quorum sensing and growth and reveals its variabilityAngewandte Chemie International Edition 48:5908–5911.https://doi.org/10.1002/anie.200901550

12
Dynamical aspects of evolutionary stabilityMonatshefte für Mathematik 110:189–206.https://doi.org/10.1007/BF01301675

13
The Vlasov dynamics and its fluctuations in the 1/N limit of interacting classical particlesCommunications in Mathematical Physics 56:101–113.https://doi.org/10.1007/BF01611497

14
Nonenzymatic turnover of an Erwinia carotovora quorumsensing signaling moleculeJournal of Bacteriology 184:1163–1171.https://doi.org/10.1128/jb.184.4.11631171.2002
 15
 16

17
Confinementinduced quorum sensing of individual Staphylococcus aureus bacteriaNature Chemical Biology 6:41–45.https://doi.org/10.1038/nchembio.264

18
Stability and ensemble inequivalence in a globally coupled systemPhysical Review Letters 91:124101.https://doi.org/10.1103/PhysRevLett.91.124101

19
Stability of the replicator equation with continuous strategy spaceMathematical Social Sciences 50:127–147.https://doi.org/10.1016/j.mathsocsci.2005.03.001
 20

21
Current knowledge and perspectives on biofilm formation: the case of Listeria monocytogenesApplied Microbiology and Biotechnology 97:957–968.https://doi.org/10.1007/s0025301246111
 22

23
Statistical mechanics of collisionless relaxation in a noninteracting systemPhilosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369:439–452.https://doi.org/10.1098/rsta.2010.0251
 24
 25

26
Vlasov equationsFunctional Analysis and Its Applications 13:115–123.https://doi.org/10.1007/BF01077243

27
A mathematical model for quorum sensing in Pseudomonas aeruginosaBulletin of Mathematical Biology 63:95–116.https://doi.org/10.1006/bulm.2000.0205

28
Einführung in Theorie und Anwendung der Laplace Transformation, Vol. 3Stuttgart: Birkhaeuser Verlag.
 29

30
Bistability in bacteriaMolecular Microbiology 61:564–572.https://doi.org/10.1111/j.13652958.2006.05249.x
 31

32
Rare event statistics in reactiondiffusion systemsPhysical Review E 70:041106.https://doi.org/10.1103/PhysRevE.70.041106
 33

34
Meanfield equation of a stochastic manyparticle process for the temporal evolution of quorumsensing microbial populations.

35
Evolutionary game theory: theoretical concepts and applications to microbial communitiesPhysica A: Statistical Mechanics and Its Applications 389:4265–4298.https://doi.org/10.1016/j.physa.2010.02.047

36
A design principle of grouplevel decision making in cell populationsPLoS Computational Biology 9:e1003110.https://doi.org/10.1371/journal.pcbi.1003110

37
Listening in on bacteria: acylhomoserine lactone signallingNature Reviews Molecular Cell Biology 3:685–695.https://doi.org/10.1038/nrm907
 38

39
Evidence of autoinduction heterogeneity via expression of the agr system of Listeria monocytogenes at the singlecell levelApplied and Environmental Microbiology 77:6286–6289.https://doi.org/10.1128/AEM.0289110

40
Communication and autoinduction in the species Listeria monocytogenesCommunicative & Integrative Biology 2:371–374.https://doi.org/10.4161/cib.2.4.8610

41
A general method for numerically simulating the stochastic time evolution of coupled chemical reactionsJournal of Computational Physics 22:403–434.https://doi.org/10.1016/00219991(76)900413

42
Exact stochastic simulation of coupled chemical reactionsThe Journal of Physical Chemistry 81:2340–2361.https://doi.org/10.1021/j100540a008

43
Quorum sensing in nitrogenfixing rhizobiaMicrobiology and Molecular Biology Reviews 67:574–592.https://doi.org/10.1128/MMBR.67.4.574592.2003

44
Transition to quorum sensing in an Agrobacterium population: a stochastic modelPLoS Computational Biology 1:e37.https://doi.org/10.1371/journal.pcbi.0010037

45
Quorum quenching: role in nature and applied developmentsFEMS Microbiology Reviews 40:86.https://doi.org/10.1093/femsre/fuv038
 46

47
Evidence of autoinducerdependent and independent heterogeneous gene expression in Sinorhizobium fredii NGR234Applied and Environmental Microbiology 80:5572–5582.https://doi.org/10.1128/AEM.0168914

48
Phenotypic heterogeneity, a phenomenon that may explain why quorum sensing does not always result in truly homogenous cell behaviorApplied and Environmental Microbiology 81:5280–5289.https://doi.org/10.1128/AEM.0090015

49
Escape from a metastable stateJournal of Statistical Physics 42:105–148.https://doi.org/10.1007/BF01010843
 50

51
Core principles of bacterial autoinducer systemsMicrobiology and Molecular Biology Reviews 79:153–169.https://doi.org/10.1128/MMBR.0002414

52
Nonequilibrium critical phenomena and phase transitions into absorbing statesAdvances in Physics 49:815–958.https://doi.org/10.1080/00018730050198152

53
Quorum sensing in marine microbial environmentsAnnual Review of Marine Science 9:257–281.https://doi.org/10.1146/annurevmarine010816060656
 54

55
Evolutionary game dynamicsBulletin of the American Mathematical Society 40:479–520.https://doi.org/10.1090/S0273097903009881
 56

57
Stochasticity in gene expression: from theories to phenotypesNature Reviews Genetics 6:451–464.https://doi.org/10.1038/nrg1615

58
Communication in bacteria: an ecological and evolutionary perspectiveNature Reviews Microbiology 4:249–258.https://doi.org/10.1038/nrmicro1383

59
Extinction rates for fluctuationinduced metastabilities: a realspace WKB approachJournal of Statistical Physics 127:861–886.https://doi.org/10.1007/s1095500793122
 60

61
Exploiting quorum sensing to confuse bacterial pathogensMicrobiology and Molecular Biology Reviews 77:73–111.https://doi.org/10.1128/MMBR.0004612
 62
 63

64
Evolutionary game theory and adaptive dynamics of continuous traitsAnnual Review of Ecology, Evolution, and Systematics 38:403–435.https://doi.org/10.1146/annurev.ecolsys.36.091704.175517

65
Evolutionary game theory in growing populationsPhysical Review Letters 105:178101.https://doi.org/10.1103/PhysRevLett.105.178101

66
Random processes in geneticsMathematical Proceedings of the Cambridge Philosophical Society 54:60–71.https://doi.org/10.1017/S0305004100033193
 67
 68

69
Evolutionary dynamics on infinite strategy spacesEconomic Theory 17:141–162.https://doi.org/10.1007/PL00004092
 70

71
Nonequilibrium dynamics of an infinite range XY model in an external fieldJournal of Statistical Physics 150:531–539.https://doi.org/10.1007/s1095501205769

72
Quorum sensing signalresponse systems in Gramnegative bacteriaNature Reviews Microbiology 14:576–588.https://doi.org/10.1038/nrmicro.2016.89

73
Sociomicrobiology: the connections between quorum sensing and biofilmsTrends in Microbiology 13:27–33.https://doi.org/10.1016/j.tim.2004.11.007
 74

75
Mathematical modelling of bacterial quorum sensing: a reviewBulletin of Mathematical Biology 78:1585–1639.https://doi.org/10.1007/s1153801601606

76
What's in a name? The semantics of quorum sensingTrends in Microbiology 18:383–387.https://doi.org/10.1016/j.tim.2010.05.003
 77
 78

79
Reversible nongenetic phenotypic heterogeneity in bacterial quorum sensingMolecular Microbiology 92:557–569.https://doi.org/10.1111/mmi.12575
 80

81
The fitness burden imposed by synthesising quorum sensing signalsScientific Reports 6:33101.https://doi.org/10.1038/srep33101

82
Rhizobium sp. strain NGR234 possesses a remarkable number of secretion systemsApplied and Environmental Microbiology 75:4035–4045.https://doi.org/10.1128/AEM.0051509

83
Phenotypic variation in bacteria: the role of feedback regulationNature Reviews Microbiology 4:259–271.https://doi.org/10.1038/nrmicro1381
 84

85
Evolutionary stable strategies and game dynamicsMathematical Biosciences 40:145–156.https://doi.org/10.1016/00255564(78)900779

86
Coevolutionary dynamics: from finite to infinite populationsPhysical Review Letters 95:238701.https://doi.org/10.1103/PhysRevLett.95.238701
 87

88
Quorum sensing: celltocell communication in bacteriaAnnual Review of Cell and Developmental Biology 21:319–346.https://doi.org/10.1146/annurev.cellbio.21.012704.131001

89
Master equations and the theory of stochastic path integralsReports on Progress in Physics 80:046601.https://doi.org/10.1088/13616633/aa5ae2

90
The different limits of weak selection and the evolutionary dynamics of finite populationsJournal of Theoretical Biology 247:382–390.https://doi.org/10.1016/j.jtbi.2007.03.015
 91

92
LuxS quorum sensing: more than just a numbers gameCurrent Opinion in Microbiology 6:191–197.https://doi.org/10.1016/S13695274(03)000286

93
Stability criteria of the Vlasov equation and quasistationary states of the HMF modelPhysica A: Statistical Mechanics and Its Applications 337:36–66.https://doi.org/10.1016/j.physa.2004.01.041
 94
Decision letter

Michael DoebeliReviewing Editor; University of British Columbia, Canada
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Ecoevolutionary dynamics in quorumsensing microbial populations can induce heterogeneous production of autoinducers" for consideration by eLife. Your article has been favorably evaluated by Gisela Storz (Senior Editor) and three reviewers, one of whom, is a member of our Board of Reviewing Editors. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
The paper presents a theoretical model for the emergence of phenotypic heterogeneity in autoinducer production, which is essential for quorum sensing in bacteria. The paper shows that a feedback between the total concentration of autoinducers and individual costly production of autoinducers is sufficient to generate heterogeneity in individual production.
Essential revisions:
The reviewers found the paper very interesting and well written. However, they raised a number of concerns that preclude publication of the paper in its present form.
Notably, some comments point to a lack of integration of the theoretical results with empirical reality. On the one hand, the paper should make more specific suggestions about how the proposed theory can be tested. On the other hand, and more importantly, the authors should consider incorporating more biological realism into their modelling. Reviewer 3 makes a number of useful suggestions in that regard. For example, this reviewer suggests incorporating heterogeneity in signal perception, an issue that was also commented on by reviewers 1 and 2. Further, it is unclear how decay of the signal concentration in the environment would affect the results. You will see that the reviewers have made a number of other suggestions about how to improve the paper.
If you find that you are able to address the concerns raised in the reviewer comments we would welcome a resubmission of your article. In the revision, please try to address all of the reviewers' comments in a constructive way. In particular, it would be good if you could extend the theoretical analysis to cover the salient comments.
Reviewer #1:
This is an interesting paper showing how the emergence of phenotypic heterogeneity in autoinducer production can, in principle, be understood as an emergent collective behaviour of quorum sensing in microbial ecosystems.
The model presented is simple and elegant, and results are supported by both numerical simulations and analytical derivations. The main result of long transients at heterogeneous states in the stochastic model is surprising, and the meanfield approximation serves the purpose of explaining this result very well.
In principle, the analytical work generates novel and testable predictions for autoinducer heterogeneity in quorum sensing (although the paper might benefit from being more specific in this regard).
1) The phenotypic response defined by the response function R is assumed to be deterministic (as e.g. seen by the Delta function appearing in the autoinducer equation (Acar et al., 2008)). I think it is important to investigate what happens if the response is instead probabilistic, i.e., when the offspring phenotype is drawn probabilistically from some distribution with mean R(
) and positive variance.
2) I don't agree with the paper's distinction between "ecological" and "evolutionary" dynamics. The authors call the feedback mechanism for autoinducer production "ecological", and they call differences in growth rates "evolutionary". However, there are no genetic differences between individuals, and hence no evolutionary dynamics occur. Instead, individuals simply have different birth rates, which is a purely ecological difference. Thus, both processes, fitness difference and global feedback, are clearly ecological in this model, and any claim that the former is "evolutionary" is misleading. The coupling is not between "ecological" and "evolutionary" dynamics, but between global production and individual birth rates.
Reviewer #2:
This manuscript shows how phenotypic heterogeneity in autoinducer (AI) production may arise in monostable autoregulation from the interplay between sensing and responding to the environment and fitness differences between producer and nonproducer phenotypes. I really enjoyed the manuscript, and feel that this is an important contribution to the study of collective behavior in microbes, mainly because it provides an alternative mechanism to explain AI production phenotypic heterogeneity using "bistable threshold models". The paper is wellwritten and technically very rigorous. In my view, it would be acceptable for publication in eLife after some minor points are clarified.
According to this model, phenotypic heterogeneity relies on (i) fitness differences between producers and nonproducers linked to the metabolic cost associated with AI production, and (ii) AI production, p, is more likely to be inherited from the parental cell than obtained from the mean level of production,
, in the environment. I am wondering if, due to this fitness difference between producers and nonproducers the authors could devise any population level experimental procedure (in addition to the single cell experiments briefly discussed in the last paragraph of the Discussion) to check the existence of the proposed feedback. The authors envision this possibility, but I think that the paper would greatly improve by establishing tighter connections between theoretical results and their possible empirical confirmation. In my opinion, this would be an important point to reach the broad readership of eLife, although the theoretical results are of enough significance by their own.
I have some doubts about the probabilistic mechanism by which newborn cells adopt. In my opinion, this is the most important ingredient of the model, because it allows for phenotypic heterogeneity, and I think that it would be good to discuss whether there is empirical support for such election (is the production level a trait maintained during the whole life cycle of the cell?) or whether that is a choice of the authors. In the latter case, it could be interesting to discuss other possible mechanisms and maybe outline the robustness of the results presented here against these alternatives.
I think that the importance of the paper could be highlighted, especially for nonspecialist readers, if the Introduction is slightly reorganized. As it is now, one may understand that there are experimental proofs of both bistable and monostable autoinducer synthesis regulation, and that the authors provide a mechanism that could explain phenotypic heterogeneity in the latter case. However, if I understood the Discussion correctly (second paragraph), bistable autoregulation has not been experimentally verified, but there are not empirical studies showing monostable autoregulation either. I think this point should appear earlier in the Introduction, accompanied by more references to existent models that utilize a bistable autoregulation (1,2). My feeling is that if there is not experimental verification of monostable regulation in AI synthesis, then this paper opens a much broader and deeper question because it not only provides a new mechanism by which phenotypic heterogeneity can emerge, but also suggests that AI synthesis regulation could be monostable. In either case, I think that this should be clarified.
1) Goryachev AB, Toh DJ, Wee KB, Lee T, Zhang HB, et al. (2005) Transition to Quorum Sensing in an Agrobacterium Population: A Stochastic Model. PLOS Computational Biology 1(4): e37. doi: 10.1371/journal.pcbi.0010037
2) Dockery, Jack D., and James P. Keener. "A mathematical model for quorum sensing in Pseudomonas aeruginosa." Bulletin of mathematical biology 63.1 (2001): 95.
Reviewer #3:
The present study describes theoretical model to provide a possible explanation for the occurrence of phenotypic heterogeneity in quorum sensing. In general phenotypic heterogeneity is an interesting topic which spans from antibiotics persistence to bacterial competence. Resent study has shown that bacteria exhibit reversible heterogeneity in QS response even in the presence of saturating concentration of QS signal, which indicate that bacteria maintain a stochastic homogenous population as a bethedging strategy to counter fluctuating environmental conditions.
1) The model has been thought in the lines of heterogeneous expression of QS synthase gene by epigenetic mechanism, however in QS, sensor plays an important role in the perception and auto regulation of the QS synthase as well as expression of other QS controlled genes. In this model however, the heterogeneity in terms of QS signal perception is not considered which could play an important role in the process.
2) The model is based on the assumption that there appears to be little change in the extracellular signal concentration once it is produced or the intracellular conc. It has now been shown that signal degradation also takes place in some QS system such as in Pseudomonas syringae, in which bacteria secrete QS degrading enzymes (some are secreted and some are intracellular), which can also influence the heterogeneity in the population. However, it appears that in this model it could have been also considered while simulating the fluctuation along with the assumption of nonproducers grows faster.
3) It has been shown that in certain environmental condition, the QS nonproducers has a growth advantage as they take benefit of the social task performed by the responders and also save on signal production cost. However, it has also has been shown in case of Pseudomonas that QS strains has a big disadvantage under certain environmental condition where QS controlled "private goods" are required (Science. 2012 Oct 12;338(6104):2646.). How will this affect the present model?
4) How the stability of auto inducers in the system will influence the cellular response to fluctuating environment and phase transition to homogeneous response?
5) A factor that could play a role in the QS heterogeneity response is the relative affinity of different QS signal with the receptor. In case the affinity is high, the local concentration of QS signal could be maintained for a prolong period and can affect the distribution. In this model, it is better to include the signal perception component also, as ultimately even for the stochastic expression of signal synthase, the signal perception and henceforth regulation of gene expression also contributes to the overall QS process.
6) What happens when competition experiments are performed with signal blind and signal sensing mutants in this model? Does fluctuating the ratio of these variants also influence the heterogeneous distribution?
https://doi.org/10.7554/eLife.25773.016Author response
Essential revisions:
The reviewers found the paper very interesting and well written. However, they raised a number of concerns that preclude publication of the paper in its present form.
Notably, some comments point to a lack of integration of the theoretical results with empirical reality.
We thank you and reviewers #1, #2, and #3 for the careful reading of our manuscript and for the pertinent comments. Following your suggestions, we have made changes to our manuscript in order to integrate our theoretical results with empirical reality. In particular, we greatly appreciate the suggestion to extend the Discussion of our results in the light of experimental microbiology with respect to model assumptions, model predictions, and their robustness. The following list summarizes the essential changes to the manuscript. Further details are provided in the respective pointbypoint response to the reviewer’s comments; see below.
1) We significantly extended the discussion of our results in the “Discussion” section with focus on integrating our theoretical results with empirical reality. We discuss in detail our model assumptions and added experimental studies to the text that support our assumptions. We also indicate experimental directions to verify or falsify all other model assumptions. Furthermore, we discuss the experimental implications of our results for phenotypic heterogeneity in quorumsensing microbial populations and propose directions for future experiments to test these implications. We included the following points to the manuscript:
1.1) Subsection “Does autoinducer production reduce individual growth rate?”. Discussion of the correlation between the cellular production degree of autoinducers and an individual’s growth rate.
We now describe in detail the experimental background of our model assumption that producers of large autoinducer molecules reproduce slower than nonproducers.
We include the recent experimental study of [Ruparell et al., 2016] (DOI: 10.1038/srep33101) into our discussion, in which metabolic costs were quantified in terms of reproduction rates for a microbial strain that uses AHL molecules as autoinducers for signaling. This study carefully quantifies the consequences of signaling on a strain’s fitness by conducting experiments in both mono and mixed culture with signal blind and signal sensing mutants as envisaged by reviewer #3.
Furthermore, we include the study of [He et al., 2003] (DOI: 10.1128/JB.185.3.809822.2003) into the “Discussion” section. This study showed growth impairments due to signaling for the strain Sinorhizobium fredii NGR234 (the same strain for which phenotypic heterogeneity heterogeneity in the expression of autoinducer synthase genes was observed in [Grote et al., 2014]). Therefore, it would be interesting to explore experimentally whether the phenotypic heterogeneity observed in Sinorhizobium fredii NGR234 could be explained through the ecoevolutionary mechanism proposed by our quorumsensing model.
We also rearranged the model definition in the section “Setup of the quorumsensing model” in that the theoretical estimates of [Keller and Surette, 2006] for the metabolic costs due to signaling are now shifted to the “Discussion” section. This way, the model definition becomes more streamlined as suggested by reviewer #2.
1.2) Subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”. Discussion of the spatiotemporal scales determining the stability and the dispersal of autoinducers, and underlying the assumptions in our quorumsensing model.
We assume in the quorumsensing model that autoinducers are uniformly degraded in a wellmixed 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. In particular, we now refer to the experimental studies of [Yates et al., 2002] (DOI: 10.1128/IAI.70.10.56355646.2002) and [Byers et al., 2002] (DOI: 10.1128/jb.184.4.11631171.2002), in which pH and temperaturedependent degradation of autoinducers in the environment are discussed (for example, the halflife time of AHL molecules can be less than 30 min at pH = 8.5 and temperature = 37 degree Celsius; this time scale is, thus, much shorter than typical microbial doubling times). These studies support our model assumptions of a high degradation rate of autoinducers in the environment under conditions that are fulfilled, for example, during stationary phase of microbial growth.
1.3) Subsection “How is production of autoinducers upregulated at the singlecell level?”. Discussion of the regulation of autoinducer production at the singlecell level.
As proposed by the comments of reviewers #2 and #3, we now discuss in greater detail the regulation of autoinducer production at the cellular level. In particular, we detail on the difference between monostable and bistable upregulation of autoinducer synthesis. Our analysis shows that the qualitative form of the regulation could discriminate between different mechanisms that control phenotypic heterogeneity of the autoinducer production at the population level. We describe how phenotypic heterogeneity may arise through stochastic gene expression for bistable gene regulation (but not for monostable regulation). We then contrast this mechanism with the ecoevolutionary mechanism identified in our study, which induces phenotypic heterogeneity for both monostable and for bistable regulation.
As envisaged by reviewer #2, we now mention that observations at the population level of both phenotypic heterogeneity in autoinducer synthesis and of growth rate differences between producers and nonproducers could be an indicator for monostable regulation of autoinducer synthesis at the cellular level.
1.4) Subsection “How is production of autoinducers upregulated at the singlecell level?”, last paragraph. Discussion of timescales on which microbes respond to autoinducers in the environment.
In the implementation of our quorumsensing model, individuals respond to the environment stochastically at reproduction events. Reviewer #2 asked for a discussion of this implementation. The rule that offspring individuals can only respond at reproduction events represents a coarsegrained view in time to facilitate the mathematical analysis and to identify the ecoevolutionary mechanism by which phenotypic heterogeneity can be induced. 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 in the quorumsensing model. This can also be inferred from the meanfield equation and the prefactor 2\lambda of the senseandresponse term. Even though sense and response to the environment are implemented at reproduction events in the model, effective changes of the population due to senseandresponse occurs at rate 2\lambda in the macroscopic description.
To quantify such response rates in a microbiological setting, experiments at the singlecell level seem most promising to us at present (see point 1.5 below). We hope that our work stimulates such experimental studies.
1.5) Subsection “Singlecell experiments”. Discussion of singlecell experiments.
We now discuss why and how singlecell experiments would help to verify or falsify our model assumptions (for example, to determine the form of the fitness and response function, to quantify the selection strength and response probability, and to refine the model setup).
More explicitly, it would be desirable to simultaneously monitor, at the singlecell 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 singlecell growth rate.
1.6) Subsection “What is the function of phenotypic heterogeneity in autoinducer production?”. Outlook for experiments and theory: function of phenotypic heterogeneity in the autoinducer production.
As an outlook for future experimental and theoretical work, and also in response to reviewer #1 and #3, we now explicitly mention that our work focuses on how phenotypic heterogeneity in the autoinducer production may be controlled, but that our work does not address the function of this phenomenon. The evolutionary contexts and ecological scenarios under which phenotypic heterogeneity in the autoinducer production may have evolved are still to be investigated from an experimental point of view (see [Grote et al., 2014], [Garmyn et al., 2011]). From a modeling perspective, it may be possible to 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.
2) We have extensively tested that the observed phenotypic heterogeneity is robust against incorporating more biological realism into our model setup, as suggested by you and all reviewers. More specifically, we have checked both numerically and mathematically that bimodal quasistationary states still arise in the relevant parameter regimes when noisy inheritance of the production degree, noisy perception of the average production level, and noisy response to it are included into the model setup. Noise was always taken as Gaussian noise with zero mean and positive variance. These results show that the ecoevolutionary mechanism for phenotypic heterogeneity is not a finetuned effect, but may arise as well if more biochemical details of quorum sensing are accounted for from a theoretical point of view. Therefore, we expect the ecoevolutionary mechanism to be relevant for microbial experiments.
Corresponding changes to the manuscript are made as follows:
· Representative trajectories of the quorumsensing model including noise at the abovementioned levels (inheritance, perception, response) are now included in the manuscript as Appendix Figure 2. All details of the numerical implementation are provided in the Figure caption.
· We integrated these new numerical simulations with our mathematical analysis by connecting the numerically observed robustness of the quasistationary states with the linear stability analysis of heterogeneous bimodal distributions of the quorumsensing model [Subsection “Results of mathematical analysis”, seventh paragraph].
· We now explicitly mention the robustness of phenotypic heterogeneity through the ecoevolutionary mechanism throughout the text: at the end of the subsection “Setup of the quorumsensing model”, in the subsection “Results of numerical simulations”, fourth paragraph, in the subsection “Results of mathematical analysis”, seventh and last paragraphs.
3) We improved the wording and corrected typos throughout the manuscript.
On the one hand, the paper should make more specific suggestions about how the proposed theory can be tested.
We are happy to implement this proposal. Point 1 in the response above provides all details to experimentally test our theory and lists all according changes to the manuscript.
On the other hand, and more importantly, the authors should consider incorporating more biological realism into their modelling. Reviewer 3 makes a number of useful suggestions in that regard. For example, this reviewer suggests incorporating heterogeneity in signal perception, an issue that was also commented on by reviewers 1 and 2.
We appreciate this comment and have carried out extensive numerical simulations to test the robustness of the results of our theory against incorporating more biological realism into the setup of the quorumsensing model (noisy inheritance, noisy perception, and noisy response). Point 2 in the response above provides all details and lists the according changes to the manuscript.
Further, it is unclear how decay of the signal concentration in the environment would affect the results.
We are thankful for this comment because it led us to formulate more explicitly the involved time scales underlying our theoretical model. Time scales that determine the stability of autoinducers in the environment are discussed in point 1.2 above; according changes to the manuscript are implemented in a separate paragraph [subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”].
You will see that the reviewers have made a number of other suggestions about how to improve the paper.
We thank the reviewers for their detailed reading and the helpful comments. Our responses to their specific comments and suggestions, and the accompanied changes to the manuscript are attached below.
If you find that you are able to address the concerns raised in the reviewer comments we would welcome a resubmission of your article. In the revision, please try to address all of the reviewers' comments in a constructive way. In particular, it would be good if you could extend the theoretical analysis to cover the salient comments.
Reviewer #1:
[…] In principle, the analytical work generates novel and testable predictions for autoinducer heterogeneity in quorum sensing (although the paper might benefit from being more specific in this regard).
1) The phenotypic response defined by the response function R is assumed to be deterministic (as e.g. seen by the Delta function appearing in the autoinducer equation (Acar et al., 2008)). I think it is important to investigate what happens if the response is instead probabilistic, i.e., when the offspring phenotype is drawn probabilistically from some distribution with mean R(<p>) and positive variance.
We agree with the reviewer that it is an interesting question to ask how noisy response affects phenotypic heterogeneity in our quorumsensing model. We now demonstrate in detail that phenotypic heterogeneity in our quorumsensing model is robust against noise at the level of inheritance of the production degree, perception of the average production level, and response; see our detailed response 2 in the answer to the Reviewing Editor.
2) I don't agree with the paper's distinction between "ecological" and "evolutionary" dynamics. The authors call the feedback mechanism for autoinducer production "ecological", and they call differences in growth rates "evolutionary". However, there are no genetic differences between individuals, and hence no evolutionary dynamics occur. Instead, individuals simply have different birth rates, which is a purely ecological difference. Thus, both processes, fitness difference and global feedback, are clearly ecological in this model, and any claim that the former is "evolutionary" is misleading. The coupling is not between "ecological" and "evolutionary" dynamics, but between global production and individual birth rates.
We agree with the reviewer that the term “evolutionary dynamics” is used to mean different phenomena in different fields of (theoretical) biology. “Evolutionary dynamics” is used both in the field of population genetics (to term changes of genotypes) and in the field of population dynamics (to term changes of phenotypes). From a theoretical point of view, it is ultimately a question of time scales whether one terms the temporal change of a dynamic variable attached to an individual in a population as “evolutionary” or not. Since our manuscript deals with phenotypic heterogeneity (that is, the coexistence of several phenotypes in a population of genetically identical individuals), we use the vocabulary that is common in the field of phenotypic heterogeneity; see, for example, the text books of [MaynardSmith, Evolution and the Theory of Games, 1982], [Hofbauer and Sigmund, Evolutionary Games and Population Dynamics, 1998], and [Nowak, Evolutionary Dynamics, 2006].
In order to make clearer that the context of “evolutionary” is meant in our work at the level of phenotypes in a population of individuals, we changed the wording slightly at the following points:
· “The key aspect of our work is how [new: the composition of] a population evolves in time when its constituents respond to an environment that is being shaped by their own activities.”
· “The coupling of ecological dynamics (given by the average production level of autoinducers
) with evolutionary dynamics (determined by fitness differences [new: between the phenotypes]) through quorum sensing results in interesting collective behavior, as we show next.”
· See also the changes to the manuscript in response to the next comment.
Reviewer #2:
[…] According to this model, phenotypic heterogeneity relies on (i) fitness differences between producers and nonproducers linked to the metabolic cost associated with AI production, and (ii) AI production, p, is more likely to be inherited from the parental cell than obtained from the mean level of production, <p>, in the environment. I am wondering if, due to this fitness difference between producers and nonproducers the authors could devise any population level experimental procedure (in addition to the single cell experiments briefly discussed in the last paragraph of the Discussion) to check the existence of the proposed feedback. The authors envision this possibility, but I think that the paper would greatly improve by establishing tighter connections between theoretical results and their possible empirical confirmation. In my opinion, this would be an important point to reach the broad readership of eLife, although the theoretical results are of enough significance by their own.
We thank the reviewer for the comments and the detailed suggestions. As described in detail in response 1 of the answer to the Reviewing Editor, we now make a couple of suggestions as to how our model assumptions are supported by experimental findings or how they might be verified or falsified from an experimental point of view. We now discuss in the manuscript:
· the correlation between the cellular production degree of autoinducers and an individual’s growth rate [subsection “Does autoinducer production reduce individual growth rate?”],
· the spatiotemporal scales that determine the stability and the dispersal of autoinducers and that are assumed in our model [subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”],
· the regulation of autoinducer production at the singlecell level [subsection “How is production of autoinducers upregulated at the singlecell level?”],
· the timescales on which microbes respond to autoinducers in the environment [subsection “How is production of autoinducers upregulated at the singlecell level?”, last paragraph],
· singlecell experiments [subsection “Singlecell experiments”], and
· an outlook for experiments and theory regarding the function of phenotypic heterogeneity in the autoinducer production [subsection “What is the function of phenotypic heterogeneity in autoinducer production?”].
I have some doubts about the probabilistic mechanism by which newborn cells adopt. In my opinion, this is the most important ingredient of the model, because it allows for phenotypic heterogeneity, and I think that it would be good to discuss whether there is empirical support for such election (is the production level a trait maintained during the whole life cycle of the cell?) or whether that is a choice of the authors. In the latter case, it could be interesting to discuss other possible mechanisms and maybe outline the robustness of the results presented here against these alternatives.
We fully agree with the reviewer that our rule how cells adopt is a simplification of the reality. The rule that offspring individuals can only respond at reproduction events is a coarsegrained view in time to facilitate the mathematical analysis and to identify the ecoevolutionary mechanism. The response probability can actually be interpreted as the rate with which individuals respond to autoinducers in the environment; see response 1.4 of the answer to the Reviewing Editor. We now discuss this assumption in detail in the Discussion section [subsection “How is production of autoinducers upregulated at the singlecell level?”, last paragraph].
The robustness of the quorumsensing model is supported by new numerical results, by our meanfield theory, and by the linear stability analysis of heterogeneous, quasistationary states of the population; please see response 2 of the answer to the Reviewing Editor.
I think that the importance of the paper could be highlighted, especially for nonspecialist readers, if the Introduction is slightly reorganized. As it is now, one may understand that there are experimental proofs of both bistable and monostable autoinducer synthesis regulation, and that the authors provide a mechanism that could explain phenotypic heterogeneity in the latter case. However, if I understood the Discussion correctly (second paragraph), bistable autoregulation has not been experimentally verified, but there are not empirical studies showing monostable autoregulation either. I think this point should appear earlier in the Introduction, accompanied by more references to existent models that utilize a bistable autoregulation (1,2). My feeling is that if there is not experimental verification of monostable regulation in AI synthesis, then this paper opens a much broader and deeper question because it not only provides a new mechanism by which phenotypic heterogeneity can emerge, but also suggests that AI synthesis regulation could be monostable. In either case, I think that this should be clarified.
1) Goryachev AB, Toh DJ, Wee KB, Lee T, Zhang HB, et al. (2005) Transition to Quorum Sensing in an Agrobacterium Population: A Stochastic Model. PLOS Computational Biology 1(4): e37. doi: 10.1371/journal.pcbi.0010037
2) Dockery, Jack D., and James P. Keener. "A mathematical model for quorum sensing in Pseudomonas aeruginosa." Bulletin of mathematical biology 63.1 (2001): 95.
We thank the reviewer for the interesting comment. Accordingly, we made the following changes to the manuscript:
· We now mention explicitly in the Introduction: “From an experimental point of view it is often not known, however, whether autoinducer synthesis is upregulated for all autoinducer levels or only above a threshold level.”
· We added the proposed references (Goryachev et al., 2005; Dockery and Keener 2001) to the manuscript. We thank the reviewer for pointing us to these early mathematical models of quorumsensing.
· We now discuss extensively the regulation of autoinducer production at the singlecell level in the subsection “How is production of autoinducers upregulated at the singlecell level?”, and have taken up many points suggested by the reviewer; see also response 1.3 of the answer to the Reviewing Editor for details.
Reviewer #3:
[…] 1) The model has been thought in the lines of heterogeneous expression of QS synthase gene by epigenetic mechanism, however in QS, sensor plays an important role in the perception and auto regulation of the QS synthase as well as expression of other QS controlled genes. In this model however, the heterogeneity in terms of QS signal perception is not considered which could play an important role in the process.
We agree with the reviewer that it is an interesting question to ask how noisy response affects phenotypic heterogeneity in our quorumsensing model. We now demonstrate in detail that phenotypic heterogeneity in our quorumsensing model is robust against noise at the level of inheritance of the production degree, perception of the average production level, and response; see our detailed response 2 in the answer to the Reviewing Editor.
2) The model is based on the assumption that there appears to be little change in the extracellular signal concentration once it is produced or the intracellular conc. It has now been shown that signal degradation also takes place in some QS system such as in Pseudomonas syringae, in which bacteria secrete QS degrading enzymes (some are secreted and some are intracellular), which can also influence the heterogeneity in the population. However, it appears that in this model it could have been also considered while simulating the fluctuation along with the assumption of nonproducers grows faster.
We are thankful for this comment because it led us to formulate more explicitly the involved time scales underlying our theoretical model. Time scales that determine the stability of autoinducers in the environment are discussed in response 1.2 of the answer to the Reviewing Editor; we discuss signal degradation in the environment now in a separate paragraph [subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”]. In essence, we assume in the quorumsensing model that autoinducers are uniformly degraded at a high rate in a wellmixed environment. A high rate refers here to a time scale that needs to be compared to the strain’s doubling time (for example, the halflife time of AHL molecules can be less than 30 min at pH = 8.5 and temperature = 37 degree Celsius; this time scale is, thus, much shorter than typical microbial doubling times). Our assumptions of high autoinducer degradation in the environment may not hold true for a spatially structured microbial biofilm, but may be fulfilled during the stationary phase of microbial growth in a wellmixed batch culture.
3) It has been shown that in certain environmental condition, the QS nonproducers has a growth advantage as they take benefit of the social task performed by the responders and also save on signal production cost. However, it has also has been shown in case of Pseudomonas that QS strains has a big disadvantage under certain environmental condition where QS controlled "private goods" are required (Science. 2012 Oct 12;338(6104):2646.). How will this affect the present model?
The reviewer raises two interesting points with this comment that we are happy to discuss now in the “Discussion” section.
· How does the spatial dispersal of autoinducers in the environment affect the observed phenotypic heterogeneity? We now discuss the scales that determine the dispersal of autoinducers and that are assumed in our model [subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”]. 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. 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 spatial organization of the microbial population determines as to what extent microbes sense rather the global or a local average production level.
· How do benefits of signaling (either to the cellular level or at the population level) affect the observed phenotypic heterogeneity? From a modeling perspective, this question points in the direction of the function of phenotypic heterogeneity in the autoinducer production, which we briefly discuss now as an outlook for experiments and theory [subsection “What is the function of phenotypic heterogeneity in autoinducer production?”]. It might be possible, for example, to extend 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; also privatization of the signal as a public good may be implemented as the reviewer suggests. Heuristically speaking, it seems straightforward from a theoretical point of view that producing traits may proliferate more easily in (quorumsensing) microbial populations if signal molecules confer an explicit benefit to single cells and if signal molecules remain privatized. We have the impression that modeling the ecological and evolutionary contexts under which phenotypic heterogeneity in the autoinducer production might have evolved would go beyond the scope of the current manuscript. Furthermore, these points are still under investigation from an experimental point of view (see [Grote et al., 2014], [Garmyn et al., 2011]). Therefore, we formulated these thoughts as an outlook in the aforementioned subsection.
4) How the stability of auto inducers in the system will influence the cellular response to fluctuating environment and phase transition to homogeneous response?
Again, this is a very interesting question from an experimental point of view. In our model, we assume that autoinducers are degraded at a rate that is higher than the cellular growth rate. We now discuss ecological scenarios under which this assumption may hold true and conclude that the stationary phase of microbial growth in batch culture might provide such a scenario. Time scales that determine the stability of autoinducers in the environment are further explained in response 1.2 of the answer to the Reviewing Editor; we discuss the stability of autoinducers in the environment now in a separate paragraph of the Discussion section [subsection “A question of spatiotemporal scales: How stable and how dispersed are autoinducers in the environment?”].
5) A factor that could play a role in the QS heterogeneity response is the relative affinity of different QS signal with the receptor. In case the affinity is high, the local concentration of QS signal could be maintained for a prolong period and can affect the distribution. In this model, it is better to include the signal perception component also, as ultimately even for the stochastic expression of signal synthase, the signal perception and henceforth regulation of gene expression also contributes to the overall QS process.
We agree with the reviewer that it would be possible to include further molecular details of quorum sensing into our model. At present, we capture all molecular and biochemical details by the response function on a macroscopic level. The response function encapsulates all steps involved in the autoinducer production between perception of the average production level (“input” of signal) and adjustment of the individual production degree in response (“output” of producing the signal). Because many different quorumsensing architectures are known (see, for example, [Waters and Bassler, 2005]) and relevant to our work, a specific choice of regulation architecture will narrow the scope of the results of our study. Therefore, we decided to capture all steps between signal perception and regulation by one effective function, which we introduced as the response function. Such response functions are actually measured, as for example depicted in Figure 3 of [He et al., 2003] (DOI: 10.1128/JB.185.3.809822.2003; termed doseresponse curve for an acylHSL as autoinducer). This way, our model remains more flexible and adjustable to specific experimental situations if one is interested in fitting a model to experimental data. On the other hand, the response function captures the necessary condition of upregulation to explain the mechanism of phenotypic heterogeneity in autoinducer production through the ecoevolutionary mechanism.
6) What happens when competition experiments are performed with signal blind and signal sensing mutants in this model? Does fluctuating the ratio of these variants also influence the heterogeneous distribution?
We thank the reviewer for this comment.
· First, this question anticipates exactly the experiments that we envision, for example, to determine the correlation between the cellular production degree of autoinducers and an individual’s growth rate. We refer the reviewer to point 1.1 in the response to the Reviewing Editor and according changes to the manuscript [subsection “Does autoinducer production reduce individual growth rate?”].
In brief, we now added the study of [Ruparell et al., 2016] (DOI: 10.1038/srep33101), which showed a reduced fitness of a producing strain (producing a longchain AHL; OC$_{12}$HSL, synthesized via the las operon) in both mono and mixed culture compared with a nonproducing strain. It would be interesting to conduct similar experiments for the strain Sinorhizobium fredii NGR234 for which a heterogeneous synthesis of autoinducers was reported.
· Second, we have checked extensively the robustness of phenotypic heterogeneity in our model. We refer the reviewer to point 2 in the response to the Reviewing Editor for a detailed answer and the according list of changes to the manuscript. In essence, perturbations of the ratio between producers and nonproducers in the heterogeneous state of the population decay in our quorumsensing model; phenotypic heterogeneity through the ecoevolutionary mechanism is a robust outcome.
https://doi.org/10.7554/eLife.25773.017Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (SPP1617 through grant FR850/111, 2)
 Matthias Bauer
 Johannes Knebel
 Matthias Lechner
 Erwin Frey
German Excellence Initiative (Nanosystems Initiative Munich)
 Matthias Bauer
 Johannes Knebel
 Matthias Lechner
 Erwin Frey
Qualcomm European Research Studentship (n/a)
 Matthias Bauer
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank the QBio 2014 summer course ‘Microbial Strategies for Survival and Evolution’ at the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara, from which this work originated. We acknowledge fruitful discussions on phenotypic heterogeneity and quorum sensing with Paul Rainey, Kirsten Jung, Kai Papenfort, Wolfgang Streit, Jessica Grote, Vera Bettenworth, Friedrich Simmel, and Madeleine Opitz. We also thank Mauro Mobilia, Meike Wittmann, Florian Gartner, Markus F. Weber, Karl Wienand, Jonathan Liu, and Alexander Dobrinevski for discussions on the quorumsensing model. MB appreciates funding by a Qualcomm European Research Studentship. This research was supported by the German Excellence Initiative via the program ‘Nanosystems Initiative Munich’ (NIM) and by the Deutsche Forschungsgemeinschaft within the framework SPP1617 (through grant FR 850/111, 2) on phenotypic heterogeneity and sociobiology of bacterial populations. MB, JK, ML, PP, and EF designed research, performed research, and wrote the paper. The authors declare no conflict of interest.
Reviewing Editor
 Michael Doebeli, Reviewing Editor, University of British Columbia, Canada
Publication history
 Received: February 5, 2017
 Accepted: June 15, 2017
 Version of Record published: July 25, 2017 (version 1)
Copyright
© 2017, Bauer et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 848
 Page views

 183
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Computational and Systems Biology

 Cell Biology
 Computational and Systems Biology