Noise regulation by quorum sensing in low mRNA copy number systems

Abstract

Background

Cells must face the ubiquitous presence of noise at the level of signaling molecules. The latter constitutes a major challenge for the regulation of cellular functions including communication processes. In the context of prokaryotic communication, the so-called quorum sensing (QS) mechanism relies on small diffusive molecules that are produced and detected by cells. This poses the intriguing question of how bacteria cope with the fluctuations for setting up a reliable information exchange.

Results

We present a stochastic model of gene expression that accounts for the main biochemical processes that describe the QS mechanism close to its activation threshold. Within that framework we study, both numerically and analytically, the role that diffusion plays in the regulation of the dynamics and the fluctuations of signaling molecules. In addition, we unveil the contribution of different sources of noise, intrinsic and transcriptional, in the QS mechanism.

Conclusions

The interplay between noisy sources and the communication process produces a repertoire of dynamics that depends on the diffusion rate. Importantly, the total noise shows a non-monotonic behavior as a function of the diffusion rate. QS systems seems to avoid values of the diffusion that maximize the total noise. These results point towards the direction that bacteria have adapted their communication mechanisms in order to improve the signal-to-noise ratio.

Background

Gene regulation at the transcriptional level is one of the corner stones of molecular and cellular biology [1]. Recent studies in prokaryotes have revealed the existence of antisense and alternative transcripts and multiple regulators per gene that imply a highly dynamic transcriptome more similar to that of eukaryotes than first thought [2]. Still, prokaryotic gene regulation mainly relies on the binding of regulatory proteins that attach to DNA for either stimulating or repressing transcription. These binding/unbinding events are intrinsically probabilistic because of the significance of thermal fluctuations at that scale and the low number of molecules involved in the process. In this regard, over the past years a growing number of experiments have indeed characterized not only the levels of randomness in cellular biochemical processes but also their functionality [3–8].

Technical advances such as the use of fluorescent tags in single-cell experiments have allowed for quantitative measurements of the noise in protein concentration and have shed light on the mechanisms of gene expression that lead to cell-to-cell variability [9–11]. Moreover, the advent of experimental approaches that permit to count individual mRNA and protein molecules in single cells has further evidenced the role played by fluctuations and their characteristics [12, 13]. Thus, in E. coli, the direct measurement of integer-valued numbers of mRNA as a function of time has revealed transcriptional bursts with Poissonian statistics [14]. The latter is in agreement with the two-state gene expression model where switching between the active and inactive transcriptional regimes occurs with constant probability [15]. It is worth noticing that these noisy sources, far for being a nuisance, have been recognized to play a constructive role in many gene regulatory processes. Examples in this direction include the efficiency of the phage lambda switch [16] or the differentiation into the competence state in bacteria [17]. All in all, it is now accepted that stochastic and non-linear approaches are required for understanding the randomness in the dynamics of biochemical reactions and the effects of fluctuations in gene regulatory networks [7, 8].

While a lot of modeling studies have focused on the single cell level [18, 19], few of them have addressed the role played by noise at the colony level [20–23]. Our recent contributions within this topic, that illustrate the constructive role of stochasticity, include the noise-induced coherence resonance phenomenon in multicellular circadian clocks [24] and the interplay between the stochasticity of the cell cycle duration and the protein expression noise in bacterial colonies [25]. One relevant question within this context is how cellular populations deal with stochasticity in communication processes. In particular, we focus on the simplest cell-to-cell communication mechanism in prokaryotes: the so-called quorum sensing (QS). The term QS generically refers to the mechanism that allow bacteria to count their number, i.e. the colony size, by producing, exporting/importing into/from the environment, and detecting a diffusive signaling molecule, namely, the autoinducer [26, 27]. As the population of bacteria grows, the autoinducer accumulates both in the extracellular medium and inside the cells. When the autoinducer concentration surpasses a threshold the expression of QS-controlled genes starts. This mechanism ultimately results in a response of the colony in a cell density dependent manner. Importantly, QS has opened the door to the design of gene circuits using synthetic biology approaches that control cell populations at the collective level [28–31].

Recent studies have shown that diffusion in QS reduces the noise at the level of the autoinducer [23]. However, to the best of our knowledge, the role played by different sources of stochasticity and their contribution to the dynamics of the signaling molecule has not been characterized yet. Moreover, while in eukaryotes the diffusion seems to contribute for enhancing the precision of regulatory processes, similar effects have not been reported in the context of QS [32]. We point out that a deep understanding of these issues is key in order to design robust synthetic circuits based on such bacterial communication system. Herein, we address these problems by studying the interplay between the QS communication and the transcriptional noise in bacterial populations. First, we aim at understanding how that interaction determines the dynamics of the autoinducer. Second, we aim at shedding light on the mechanisms that confer robustness to noise in QS communication. Transcriptional noise is expected to be particularly relevant when transcription events are short and rare. Under these conditions two main sources of stochasticity naturally arise: the dichotomous fluctuating dynamics of mRNA and the intrinsic noise due to low copy number of species. Thus, we restrict ourselves to the study of the aforementioned problems near the QS activation threshold where we can assume that the transcription events produce basal constitutive levels of mRNA as low as one molecule per cell at a time. We note that such mRNA production processes have been experimentally validated in prokaryotes revealing that the statistics of proteins bursts originates from the translation of a single mRNA molecule [33].

Our main findings are threefold. First, we show how the diffusion process leads to a repertoire of dynamics in regards of the signaling molecule. Second, we demonstrate that, for a large range of diffusion rate values, the main contribution to the total noise of the autoinducer concentration is the mRNA fluctuations. Finally, we show that the total noise exhibits a non-monotomic behavior as a function of the diffusion rate in contrast to previous results [23].

The paper is organized as follows. In the Methods section we introduce our modeling approach, the analytical calculations, and the parameter values used in our in silico experiments. In the next section, Results, we present our findings and compare the results from stochastic simulations with those from analytical calculations. Finally, we discuss further implications of our results in the Discussion, and summarize our findings in the Conclusion.

Methods

Modeling Approach

A large class of gram-negative bacteria use acyl homoserine lactones (AHLs) as signaling molecules [27, 34]. These autoinducer molecules are typically synthesized by enzymes of the LuxI family and can freely diffuse across the cell membrane, i.e. by means of passive diffusion. When the concentration of autoinducer surpasses a critical threshold, it binds to its receptor (a cytoplasmic protein of the LuxR family) which then activates the expression of target genes, e.g. in Vibrio fisheri species the luxICDABE operon, that is responsible for the expression of luxI and luciferase. In contrast, when the autoinducer concentration is below the activation threshold, the transcription of the luxI gene occurs at a low basal rate, thus producing low levels of the enzyme. In this regime, the feedback regulation of the luxI gene leading to autoinduction can be disregarded. As a matter of fact, a number of QS systems lack autoinduction [35]. Consequently, to neglect the feedback in those cases is a valid approximation even above the activation threshold. Herein, for the sake of simplicity, we focus on this situation and describe the dynamics of the autoinducer near, but below, the activation threshold when feedback loops can be neglected and the downstream QS genes are not activated. We stress that, generally speaking, our results are not applicable above the activation threshold since the dynamics of the autoinducer is obviously conditioned by feedback effects.

Following previous approaches we consider a two-stage model for gene expression/regulation [15, 36–38]. Thus, we assume that during the transcription events a single mRNA molecule, e.g. a luxI transcript, is produced and its dynamics can be then described by means of a Markovian dichotomous process [39],

$$M_0^i\underset{\alpha }{\overset{\beta }{\rightleftarrows }}{M}_1^i$$

(1)

where $M_{0,1}^i=0,1$ stands for the number of mRNA molecules at cell i and α and β for the transition rates between these states; i.e. α and β account for the probabilities per unit of time of mRNA degradation and transcription frequency respectively. Notice that the stochastic alternation of the mRNA between the values 0 (no mRNA) and 1 (a mRNA molecule) is not memoryless, i.e. white. Once a mRNA molecule is produced, and until it becomes degraded, the cell keeps noticing its presence and keep producing the autoinducer. That is, the transcriptional noise is a colored noise, and its autocorrelation decays exponentially with a characteristic time scale τ _c = (α + β)^-1[39].

Once a mRNA molecule is produced the translational, and post-translational processes (if any), leads to the appearance of functional LuxI synthetases. Yet, our interest here focuses on the dynamics of the signaling molecule. It has been shown that the amount of the synthetase substrate is not a limiting factor for the production of the autoinducer [40, 41]. As a consequence, the levels of the signaling molecule depends directly on the expression levels of the synthetase. Ignoring intermediate biochemical steps in the autoinducer synthesis reduces the number of noise sources and may even change, under some circumstances, the observed dynamics [42]. Still, it is a valid approximation in many situations, and here we assume that the translation of the synthetase and the subsequent synthesis of the autoinducer, A, can be effectively described by a single chemical step with rate k₊. In addition, we consider that the autoinducer becomes degraded at a rate k_-, that is,

$$M_1^i\stackrel{k_{+}}{\to }{M}_1^i+A^i$$

(2)

$$A^i\stackrel{k_{-}}{\to }$$

(3)

Passive diffusion of the autoinducer can be implemented by considering a new species, A_ext, that accounts for the number of signaling molecules in the extracellular medium such that,

$$A^i\underset{rD}{\overset{D}{\rightleftarrows }}{A}_{\text{ext}}$$

(4)

where D stands for the diffusion rate and r = V/V_ext represents the ratio of the volume of a cell to the total extracellular volume. We consider all cells to have the same value of the diffusion rate. In addition, we assume a well-stirred system where spatial effects can be neglected.

As the bacterial population grows the autoinducer accumulates in the media. In experiments, in order to keep the concentration of the autoinducer below the activation threshold, such growth is compensated by means of a dilution protocol. As detailed below (see Parameters) the latter constitutes the main source of effective degradation of the signaling molecule. Thus, hereinafter we assume the degradation rate of the signaling molecule to be the same inside and outside the cell,

$$A_{\text{ext}}\stackrel{k_{-}}{\to }$$

(5)

Figure 1 schematically represents the biochemical processes considered in our approach. The set of reactions (1)-(5) characterizes the stochastic dynamics of the autoinducer and that of the mRNA. Their probabilistic description is given by the corresponding master equation that is exactly sampled by means of the Gillespie algorithm in a N-cells system [43].

Figure 1

Scheme of a simplified biochemical network of QS systems near the activation threshold. Schematic representation of the biochemical processes considered in our approach for describing the dynamics of the signaling molecule, A, in cell i. The mRNA dynamics satisfies a dichotomous process characterized by the states M_0,1 corresponding to zero and one molecules respectively. Once the autoinducer has been produced, it can diffuse in and out the cell leading to cell communication (see text).

Analytical Calculations: Null Intrinsic Noise Approximation

Further insight into the dynamics of the signaling molecule can be obtained by analytical means as follows. Two stochastic contributions drive the dynamics of A: the mRNA fluctuations due to the random switching (mRNA present or not) and the molecular, i.e. intrinsic, noise due to low copy number of the autoinducer. As for the latter, it can be neglected if over the course of time Aⁱ /(Aⁱ + 1) ≃ 1 (large number of autoinducer molecules). While in our system this approximation is not totally justified (see parameters values below), it is useful to implement it in order to discriminate between the effects caused by different stochastic contributions and to obtain analytical expressions. In this case, it is straightforward to demonstrate that the dynamics of the autoinducer, Eqs. (1)-(5), can be described by the following coupled stochastic equations,

$${\dot{c}}_{A^i}=k_{+}{c}_{M_1^i}(t)-k_{-}{c}_{A^i}+D(c_{A_{\text{ext}}}-c_{A^i})$$

(6)

$$\begin{array}{c}{\dot{c}}_{A_{\text{ext}}}=-k_{-}{c}_{A_{\text{ext}}}+rD{\displaystyle \sum _{i=1}^N(}{c}_{A^i}-c_{A_{\text{ext}}})\\ =-k_{-}{c}_{A_{\text{ext}}}+rDN(〈c_A〉-c_{A_{\text{ext}}}),\end{array}$$

(7)

where $c_{A^i}=A^i/V$, $c_{M_1^i}\left(t\right)=M_1^i/V,$ and $c_{A_{\text{ext}}}=A_{\text{ext}}/V_{\text{ext}}$ stand for the concentration of species A and $M_1^i$ at cell i and for species A_ext at the extracellular medium respectively, N is the colony size (number of cells), and 〈·〉 represents the population average. In Eq. (6) the term $c_{M_1^i}\left(t\right)$ accounts for a dichotomous stochastic process characterized by the rates and states (α, β) and (0, V⁻¹) respectively, and describes the fluctuating mRNA dynamics. We point out that in case that D = 0, Eq. (6) has been proposed to study graded and binary responses in stochastic gene expression. Interestingly, it has been shown that despite its simplicity it can actually reproduce some gene expression phenomena [37, 44].

We can further proceed with the analytical calculations by implementing, as in previous studies e.g. [45], a quasisteady approximation for the dynamics of the external autoinducer, i.e. ${\dot{c}}_{A_{\text{ext}}}=0,$ so that,

$$c_{A_{\text{ext}}}=〈c_A〉\frac{1}{1+\frac{k_{-}}{NDr}}.$$

(8)

By substituting Eq. (8) into Eq. (6) we obtain an equation for the concentration of the signaling molecule inside a given cell that depends on the average 〈c _A 〉(the index i has been dropped),

$$\begin{array}{c}{\dot{c}}_A=k_{+}{c}_{M_1}\left(t\right)-D\left(1+\frac{k_{-}}{D}\right)c_A\\ +〈c_A〉\frac{D}{1+\frac{k_{-}}{NDr}}\end{array}$$

(9)

In the absence of diffusion, Eq. (9) reveals that the concentration of the signaling molecule reaches a maximum value of $c_A^{+}=k_{+}/\left(k_{-}V\right)$ when $c_{M_1}(t)=V^{-1}$. In terms of $c_A^{+}$ and the time scale t _c = 1/k_-, the typical lifetime of a signaling molecule, the dimensionless version of Eq. (9) reads

$${\dot{\tilde{c}}}_A={\widehat{c}}_{M_1}\left(\tilde{t}\right)+k_{+}^{\text{eff}}\left(〈{\tilde{c}}_A〉\right)-k_{-}^{\text{eff}}{\tilde{c}}_A,$$

(10)

where $\tilde{D}=D/k_{-},k_{\_}^{\text{eff}}=1+\tilde{D},$ and $k_{+}^{\text{eff}}(〈{\tilde{c}}_A〉)=〈{\tilde{c}}_A〉\tilde{D}/\left(1+\left(1/N\tilde{D}r\right)\right);{\widehat{c}}_{M_1}(\tilde{t})$ being a Markovian dichotomous noise characterized by the states $\left\{{\widehat{c}}_{M_1}\right\}=0,1$ and the rates $\tilde{\alpha }=\alpha /k_{-}$ and $\tilde{\beta }=\beta /k_{-}$ Equation (10) can be formally closed by invoking the following self-consistency condition:

$$〈{\tilde{c}}_A〉={\displaystyle {\int }_{\tilde{\Omega }}{\tilde{c}}_{A\rho }}({\tilde{c}}_A;〈{\tilde{c}}_A)d{\tilde{c}}_A,$$

(11)

$\rho ({\tilde{c}}_A;〈{\tilde{c}}_A〉)$ being the probability density solving Eq. (10) and $\tilde{\Omega }$ its support [39]:

$$\begin{array}{c}\rho ({\tilde{c}}_A;〈{\tilde{c}}_A〉)=\mathcal{N}{(k_{-}^{\text{eff}}{\tilde{c}}_A-k_{+}^{\text{eff}}(〈{\tilde{c}}_A〉))}^{\frac{\tilde{\beta }}{k_{-}^{\text{eff}}}-1}\\ .{(1+k_{+}^{\text{eff}}(〈{\tilde{c}}_A〉)-k_{-}^{\text{eff}}{\tilde{c}}_A)}^{\frac{\tilde{\alpha }}{k_{-}^{\text{eff}}}-1}\end{array}$$

(12)

with

$$\mathcal{N}=\frac{(1+\tilde{D})\Gamma \left[\frac{\tilde{\alpha }+\tilde{\beta }}{1+\tilde{D}}\right]}{\Gamma \left[\frac{\tilde{\alpha }}{1+\tilde{D}}\right]\Gamma \left[\frac{\tilde{\beta }}{1+\tilde{D}}\right]}$$

(13)

being the normalization constant. The condition (11) can be exactly solved and leads to the following value for the average concentration:

$$\begin{array}{c}〈{\tilde{c}}_A〉=\frac{1+\tilde{D}Nr}{1+\tilde{D}Nr+\tilde{D}}\frac{\tilde{\beta }}{\tilde{\alpha }+\tilde{\beta }}\\ =\frac{1+\tilde{D}Nr}{1+\tilde{D}Nr+\tilde{D}}〈{\tilde{c}}_A〉{|}_{\tilde{D}=0}\end{array}$$

(14)

where $〈{\tilde{c}}_A〉{|}_{\tilde{D}=0}=\tilde{\beta }/\left(\tilde{\alpha }+\tilde{\beta }\right)$ is the average concentration of signaling molecules in the absence of diffusion. Thus, as expected, $〈{\tilde{c}}_A〉<〈{\tilde{c}}_A〉{|}_{\tilde{D}=0}$. For the sake of concision, on what follows we drop the argument term $〈{\tilde{c}}_A〉$ from the notation of $\rho ({\tilde{c}}_A;〈{\tilde{c}}_A〉)$. Note that $\rho ({\tilde{c}}_A)$ has two states (barriers) that define its support. That is, the minimum and maximum values that the concentration of the autoinducer can reach as a function of the diffusion are:

$$c_A^{-}=\frac{{\tilde{D}}^{2}Nr}{(1+\tilde{D})(1+\tilde{D}+\tilde{D}Nr)}\frac{\tilde{\beta }}{\tilde{\alpha }+\tilde{\beta }}$$

(15)

$${\tilde{c}}_A^{+}={\tilde{c}}_A^{-}+\frac{1}{1+\tilde{D}}$$

(16)

It is easy to prove that the probability density $\rho ({\tilde{c}}_A)$ shows a single extremum if,

$$\tilde{\alpha },\tilde{\beta }\lessgtr k_{-}^{\text{eff}},$$

(17)

Where the extremum is a maximum if $\tilde{\alpha },\tilde{\beta }>k_{-}^{\text{eff}}$ and a minimum if $\tilde{\alpha },\tilde{\beta }<k_{-}^{\text{eff}}$. In the other cases the probability density does not display any extrema. Therefore, as a function of $\tilde{\alpha }$ and $\tilde{\beta }$, the probability density $\rho ({\tilde{c}}_A)$ may show four different behaviours depending on the value of the diffusion coefficient as schematically represented in Figure 2A. However, a constraint in our modeling restricts the regions, i.e. behaviors, accessible to the autoinducer dynamics. We have assumed a low constitutive expression such that only a single mRNA molecule can be transcribed at a time. The latter implies that $\tilde{\beta }<\tilde{\alpha }$

Figure 2

Probability densities of the signaling molecule and parameter space. Panel A: Schematic representation of the different probability densities of the autoinducer concentration depending on the value of $\tilde{\alpha }$ and $\tilde{\beta }$ with respect to that of $\tilde{D}$. Given a set of values ($\tilde{\alpha }$, $\tilde{\beta }$) the dynamics of the autoinducer shows different behaviors depending on the value of the diffusion parameter since the transitions lines are located at $\tilde{\alpha }$, $\tilde{\beta }$ = 1 + $\tilde{D}$. The constraints of our modelling in terms of the parameters values make the region on the top-left corner non-accessible (see text). Panel B: Parameter space diagram ($\tilde{\alpha }$, $\tilde{\beta }$) indicating the sets of parameters used in the simulations (solid squares): γ₁ = (8, 2), γ₂ = (15, 5), γ₃ = (8, 0.5), γ₄ = (15, 0.5). The experimental values reported for the degradation rate of the mRNA leads to a biological meaningful range for $\tilde{\alpha }$ (blue region). The low constitutive expression assumption is prescribed by the constraint $\tilde{\alpha }>2\tilde{\beta }$(red colored region).

(the degradation rate of the mRNA is larger than the transcription rate) in order to assure that a maximum of one mRNA molecule is present in a cell at a given time. As a consequence, and independently of the diffusion value, the dynamics leading to the probability density shown at the top-left region of Figure 2A cannot be considered as physical in the context of our modeling approach. We stress that this constraint is not a fundamental ingredient for obtaining our results (see Discussion).

Finally, the noise of the autoinducer concentration reads,

$$\begin{array}{c}{\eta }_{{\tilde{c}}_A}^2=\frac{{\sigma }_{{\tilde{c}}_A}^2}{{〈{\tilde{c}}_A〉}^2}\\ =\frac{\tilde{\alpha }{\left(1+\tilde{D}+\tilde{D}Nr\right)}^2}{\beta \left(1+\tilde{D}\right){\left(1+\tilde{D}Nr\right)}^2\left(1+\tilde{D}+\tilde{\alpha }+\tilde{\beta }\right)}\end{array}$$

(18)

where ${\sigma }_{{\tilde{c}}_A}^2=〈{\tilde{c}}_A^2〉-{〈{\tilde{c}}_A〉}^2$.

Nonetheless we have disregarded the intrinsic noise in the analytical calculations, Eq. (18) will allow us to elucidate the contributions of different sources of noise as follows. By means of the numerical simulations (Gillespie) of the set of reactions (1)-(5) in a N-cell system, we can evaluate the total (intrinsic+transcriptional) noise. Hence, by substracting from that quantity the contribution of the transcriptional noise, i.e. Eq. (18), we obtain the levels of intrinsic noise (see Results).

Parameters

We are interested in the role played by the fluctuations of the signaling molecule, A, when its concentration is close to the activation of the QS switch. Therefore, we fix the mean concentration of the autoinducer and modulate the rest of the parameters in order to keep constant this value. Pai and You [46] have recently studied the core architecture of the QS mechanism for a comprehensive set of systems. These authors have estimated that the critical concentration of autoinducer needed for the activation of the QS genes ranges from 10 to 50 nM for most of the bacterial species. In our model, we set the average concentration of A to a typical value of $C_A^0$ = 25 nM. Yet, our results do not depend on the particular value we choose within that range. As shown below, see Results, this value fixes the level of intrinsic noise of the system. However, the interplay between diffusion and transcriptional noise does not depend on that. Moreover, by defining the so-called sensing potential, ν = (rN)^-1 Pai and You estimated the range of critical cell densities for the QS activation. They concluded that its characteristic value is ν ~ 10³ - 10⁴. In our simulations we set this parameter to ν - 10³. In the experimental setups the cells are typically present in a volume of a few milliliters and the total number of cells is of the order of 10⁸ - 10⁹. Therefore, the concentration of autoinducer in the medium is determined by the exchange of signaling molecules coming from many cells. In contrast, the behavior of QS systems with a very low number of cells can be significantly different, as shown by microfluidic confinement of cells in picoliter droplets [47–49]. In our study, in order to discard small system size effects, we choose a sufficiently large number of cells in the numerical simulations, N = 10². Since the typical volume of an E. coli cell is V = 1.5 μ m³ then V_ext = 10⁵V (i.e. r = 10⁻⁵). We point out that keeping ν to a constant value necessarily requires an external dilution protocol for maintaining constant the cell density. In experiments, the control of the dilution rate is usually achieved by the use of chemostats or microfluidic devices [50]. The rate of dilution should compensate for the cell growth, ~ 2 · 10⁻² min⁻¹ (i.e. cell cycle duration ~50 min). In our modeling, by keeping constant the number of cells and the average concentration of the autoinducer, we tacitly assume a dilution protocol too. Importantly, the dilution rate effectively modifies the degradation rate of the signaling molecule. In this regard, while some bacteria species have hydrolytic enzymes that degrade AHLs, generally speaking, bacteria that synthetize AHLs do not degrade them enzymatically. In fact, AHLs are chemically stable species in aqueous solutions [51]. The degradation rate of the homoserine lactone 3-Oxo-C6-AHL has been measured in vitro revealing that this autoinducer is rather stable: ~ 3 · 10⁻⁴ min⁻¹[51]. Measurements of the degradation rate of other AHL autoinducers show similar results. Based on experimental data and mathematical modeling, the degradation rate of the signaling molecule in vivo has been also estimated [46]. Depending on the pH of the medium, the latter ranges from ~ 5 · 10⁻³ min⁻¹ to ~ 2 · 10⁻² min⁻¹. Consequently, the dilution process constitutes the main source of effective degradation of A, both inside and outside the cell, and here we set k₋ = 2 · 10⁻² min⁻¹.

As for the value of the diffusion rate, the coefficient of passive diffusion has been estimated for the 3-Oxo-C6-AHL autoinducer based on the measure of the diffusion of glucose and lactose through the outer membrane of E. coli[23]. For a typical cell volume of 1.5 μ m³ the estimated coefficient of diffusion is ~ 10³ min⁻¹. Under these conditions the typical value for the dimensionless parameter $\tilde{D}$ is of the order of 10⁴. Yet, active transport mechanisms for the autoinducer leads to much smaller effective diffusion values (see Discussion) and we explore the role played by this parameter.

In regards to the mRNA dynamics, α, the degradation rate, depends on the cell degradative machinery. To this respect, the half-lives of all mRNAs of Staphylococcus aureus have been recently measured during the mid-exponential phase. Most of the transcripts, 90%, have half-lifes shorter than 5 minutes [52, 53].

According to these studies we restrict the mRNA degradation rate to the range ln(2)/5 min⁻¹ < α < ln(2)/2 min⁻¹. Consequently, $\tilde{\alpha }$ > 1. As for the frequency of the transcription events, β is determined by particular characteristics of the gene regulatory process under consideration, e.g. the affinity of the regulatory proteins to the operator site and the initiation rate of transcription. Due to the assumption of a low constitutive expression, we choose values of parameter β satisfying the relation α > β. In particular, we implement the more restrictive condition α > 2β. Figure 2B recapitulate these constraints and show the different sets of $\tilde{\alpha }$ and $\tilde{\beta }$ values that we have used in our simulations and analytical calculations.

Summarizing, N, r, and k₋ are kept fixed in our simulations and analytical calculations and we explore the parameter space α, β, and D within the ranges, and satisfying the constraints, mentioned above. In every particular situation, once a set of those parameters is prescribed, we set the value of k₊ by using Eq. (14) in order to keep the average value of c _A near its critical concentration value (25 nM).

Results

Comprehensive Study of the Autoinducer Dynamics as a Function of the Diffusion Rate

The distribution of c _A at the steady-state is computed for the different parameter sets according to the ranges and constraints described above (section Methods). In order to explore the role of the diffusion in the dynamics of the signaling molecule we first study the case $\tilde{D}=0$. According to the analytical calculations, see Eq. (17), in this case two possible distributions for the concentration of c _A can be observed depending on the value of $\tilde{\beta }$. Since $\tilde{\alpha }$ > 1 we can expect a maximum only if $\tilde{\beta }$ > 1 (note that $k_{-}^{\text{eff}}=1$ if $\tilde{D}=0$), otherwise extrema are not expected. The results of the numerical simulations (Gillespie), Figure 3A, reveal that scenario. Note that in all cases the histogram obtained from the simulations fits fairly well to the expression (12) except for deviations due to the intrinsic noise that are not taken into account by the analytical approach. The differences among dynamics are evidenced by the trajectories, Figure 3B. Thus, for $\tilde{\beta }$ < 1 the dynamics of the autoinducer shows a burst-like behavior. If $\tilde{\beta }$ > 1 the frequency of bursts is high enough to maintain the concentration of signaling molecules near the average and a single-peak distribution develops.

Figure 3

Distributions and dynamics of the signaling molecule in a diffusionless system. Panel A: Distributions of c _A at steady-state for different sets of parameters ($\tilde{\alpha }$, $\tilde{\beta }$) as indicated in Figure 2B. In all cases $\tilde{D}$ = 0. The histogram obtained in simulations (blue bars) compares well with the distribution from the analytical calculations (blue line). Yet, deviations are observed due to intrinsic noise (see text). Panel B: The dynamics of the autoinducer show different behaviors depending on the region of the parameters phase space (see Figure 2). Two typical trajectories are shown with a grey-shaded background indicating the presence of a mRNA molecule in the cell.

If $\tilde{D}$ > 0 we expect a more fruitful phenomenology since the transition lines between behaviors in the parameter space ($\tilde{\alpha }$, $\tilde{\beta }$) shift as a function of the diffusion (see Figure 2A). According to the analytical calculations we can anticipate that, for a given parameter set and as $\tilde{D}$ increases, the system explores different dynamical regimes. By taking as a reference the case γ₂, that is ($\tilde{\alpha }$, $\tilde{\beta }$) = (15, 5), Figure 4 shows the effect of the diffusion on the distribution (left column) and dynamics (center column) of c _A in a given cell. The system initially displays a single-peak distribution for $\tilde{D}$ = 1. By increasing the diffusion coefficient we observe transitions to the other phases (monotonically decreasing and double-peak distributions). The corresponding dynamics of c _A (right panels) show how the diffusion, acting as an additional effective degradation on A, first increases the sharpness of the bursts of production. For $\tilde{D}$ = 10, the diffusion is large enough to remove signaling molecules between consecutive burst events, thus leading to a monotonically decreasing distribution. Increasing the diffusion rate to $\tilde{D}$ = 100 leads to the situation where both $\tilde{\alpha }$ and $\tilde{\beta }$ becomes smaller than 1 + $\tilde{D}$ and a bistable dynamics develops. Under these circumstances the concentration of autoinducer alternates between two states that correspond to a low concentration, when there is no mRNA production, and a high concentration, following the mRNA synthesis. As the diffusion further increases, e.g. $\tilde{D}$ = 2 · 10³, the autoinducer molecules diffusing from the external medium into the cell set a constitutive level of this species. The latter explains the presence of A molecules in the cell even if no mRNA is produced. Finally, at very large values of $\tilde{D}$, e.g. $\tilde{D}$ = 5 · 10⁴, the low constitutive concentration of the autoinducer increases due to the influx of molecules when no mRNA is present whereas the concentration of A that is internally produced decreases due to the efflux of molecules. In this case, the whole N-cells system can be considered as a single volume with no diffusive barriers between cells. Thus, the burst events average out and, as a consequence, a single effective peak appears.

Figure 4

Distributions and dynamics of the signaling molecule in a system with diffusion. Distributions (left column) and dynamics (center column) of c _A at steady-state for different values of $\tilde{D}$. The right most column stands for a density plot of the distribution of $c_{A_1}$ as a function of $c_{A_2}$ for discerning a putative increase in the molecular noise (see text). In all cases the parameters set ($\tilde{\alpha }$, $\tilde{\beta }$) is γ₂ (see Figure 2B). The production rate ${\tilde{k}}_{+}$ is modulated as a function of ($\tilde{\alpha }$, $\tilde{\beta }$, $\tilde{D}$) in order to maintain constant the average 〈c _A 〉 = 25 nM. The histograms obtained in the stochastic simulations (blue bars, left column) are in qualitative agreement with the probability densities from the analytical calculations (blue line, left column). When increasing the diffusion coefficient, the system explores different dynamics as revealed by the trajectories shown in the center column. The grey-shaded background shown in the trajectories of c _A indicates the presence of a mRNA molecule in the cell. The density plots (right column) reveals that the diffusion does not contribute to an increase of the intrinsic noise since the spreading of the distributions in a direction perpendicular to the diagonal does not grow.

Diffusion and Intrinsic Noise

Figure 4 shows that the theoretical distribution captures the essential features of the dynamics obtained in numerical simulations (Gillespie). The noticeable deviations are due to the intrinsic noise that are not considered in the theoretical analysis. Notice that as the diffusion increases those deviations seem to be larger. We stress that in our simulations we keep constant the average concentration of the autoinducer by modulating the effective production rate (see Eq. (14)). That is, as $\tilde{D}$ becomes larger, we increase the production rate k₊ so that the average number of autoinducer molecules per cell remains the same. Consequently, the deviations between simulations and the theoretical analysis cannot be ascribed to a putative intracellular decrease of the number of A molecules (i.e. to an increase of intrinsic noise). Moreover, the deviations cannot be attributed either to a failure of the quasisteady approximation introduced in Eq. (7) because the larger the diffusion, the more accurate that approximation is. As indicated by equation (18), the transcriptional noise behaves as ${\tilde{D}}^{-1}$. Thus, for large values of the diffusion rate, the transcriptional noise level decreases. Therefore, we must conclude that the deviations between the theoretical and the numerical approaches become accentuated as the diffusion increases because there is a drop of the fluctuations related to the mRNA dynamics. This also indicates that for large enough diffusion rate, the intrinsic noise constitutes the main source of stochasticity.

In order to ensure that the intrinsic fluctuations are not actually increasing due to diffusion we perform the following in silico experiment. We consider a modification of our system such that a single mRNA molecule transcript leads to two autoinducer molecules that are considered to be distinguishable. The latter can be experimentally achieved by placing two consecutive copies of the encoding sequence of the autoinducer synthetase labeled with different fluorescent tags in the operon. Thus, we double the set of equations (2)-(5) in order to account for $A_1^i$ and $A_2^i$ molecules synthesis at cell i due to a single mRNA transcript $M_1^i$. Following [9], by plotting the distribution of $c_{A_1}$ as a function of $c_{A_2}$ a putative increase of the intrinsic fluctuations can be discerned by these means. Right column of Figure 4 displays the results in this regard. The width of the distribution in a direction perpendicular to the diagonal is a measure of the intrinsic fluctuations (see [9] for details). As shown, as the diffusion increases there is no amplification of this quantity.

Diffusion and Total Noise

It is interesting to place the previous result in the context of the total noise of the autoinducer concentration. Figure 5 reveals that ${\eta }_{c_A}^2$ shows a non-monotonic behavior. As a function of $\tilde{D}$ the total noise first increases, reaches a maximum, and then decreases as the diffusion becomes larger. While Figure 5 represents data for the γ₂ parameter set, this behavior applies for all the values of $\tilde{\alpha }$ and $\tilde{D}$ explored in our simulations (data not shown). We point out that if ${\eta }_{c_A}^2$ > 1 then the dispersion is larger than the mean. Under this circumstance the fluctuations can lead to catastrophic events. E.g., fluctuations are able to remove all the signaling molecules within the cell. Note that the analytical calculations, that just account for the transcriptional noise, are in agreement with the numerical simulations, that account for both the transcriptional and the intrinsic noise, for a large range of $\tilde{D}$ values. This indicates that the main contribution to the total fluctuations within a large range of diffusion values is the transcriptional noise. Yet, as mentioned above, the latter diminishes as the diffusion increases while the intrinsic fluctuations remain constant. Consequently, the contribution of the intrinsic noise must become more relevant than the mRNA stochasticity beyond some value of $\tilde{D}$.

Figure 5

Noise of the signaling molecule as a function of the diffusion coecient. Noise ${\eta }_{c_A}^2$ as a function of diffusion coefficient $\tilde{D}$ for the set of parameters γ₂ (see Figure 2B): stochastic simulations (circles) and analytical expression. Eq. (18), (solid line). By using the decomposition ${\eta }_{c_A}^2={\eta }_{c_A,\mathrm{int}}^2+{\eta }_{c_A,\text{tran}}^2$ the differences between the computational and the theoretical distributions quantifies the amount of intrinsic noise (squares). As evidenced by the linear regression (blue short-dashed line) the later remains constant and is the main contribution to the total noise only for large diffusion values, $\tilde{D}$ > 10⁴ (see text).

We address this point quantitatively by calculating the relative importance of these noisy sources. To this end, we make use of the decomposition ${\eta }_{c_A}^2={\eta }_{c_A,\mathrm{int}}^2+{\eta }_{c_A,\text{tran},}^2$ where ${\eta }_{c_A,\mathrm{int}}^2$ and ${\eta }_{c_A,\text{tran}}^2$, stand respectively for the intrinsic and the transcriptional contributions to the total noise [54]. Thus, by subtracting the analytical expression given by Eq. (18) to the total noise obtained in the numerical simulations, we are able to compute the intrinsic noise as a function of the diffusion. By performing a linear regression of the points that corresponds to the intrinsic noise we obtain that the slope of the curve is indeed zero in practical terms, ~ 2 · 10⁻⁷. Therefore, in agreement with the qualitative results obtained in Figure 4 (right column), the intrinsic noise remains constant as the diffusion increases, ${\eta }_{c_A,\mathrm{int}}^2$ = 0.054 ± 0.003, and is the main stochastic component if $\tilde{D}$ > 10⁴.

Discussion

As a matter of discussion, our modeling consider passive diffusion as the mechanism for the transport of the signaling molecule. This is indeed true in many QS systems. However, in other cases the autoinducer is actively transported in and out of the cell. For example, in the bacterial species Pseudomonas aeruginosa, C4-HSL can freely diffuse but C12-HSL, a larger signaling molecule, is subjected to active influx and efflux at rates of ~ 10⁻² min⁻¹ and ~ 10⁻¹ min⁻¹ respectively [55]. Other example corresponds to the AI-2 signaling molecule. The latter is present in many Gram-positive and Gram-negative species and it is believed to allow for interspecies communication [56]. In E. coli and Salmonella enterica extracellular AI-2 accumulates during the exponential phase, but then decreases drastically upon entry into the stationary phase. This reduction is due to the import and processing of AI-2 by the Lsr transporter [56, 57]. Moreover, excretion from the cell of this autoinducer also appears to be an active process involving the putative transport protein YdgG (or alternatively named TqsA) [58]. In the case of E. coli these rates have been estimated by computational and experimental means: D _out ≃ 10^-1 min⁻¹ and D _in ≃ 10^-3 - 10^-2 min⁻¹[59]. All in all, the diffusion rates when driven by active processes are four orders of magnitude smaller than the diffusion rates of small molecules through the membrane. Hence, the diffusion rates in QS systems can be categorized into two main, well separated, classes: small diffusion rates due to active transport, and large diffusion rates due to passive mechanisms.

In principle our model does not account for active diffusion processes, but transport driven by concentration differences. Still, it can be demonstrated that if the following condition holds, r ≃ D _in /D _out , then our simple model fairly describes the active diffusion with D = D _out . However, that regime is not accessible in our simulations near the critical sensing potential since that would imply a non-physical situation (N < 1). If in any case we consider the rates of active transport of these QS systems, they would fit in our model with a normalized effective diffusion coefficient in the range $\tilde{D}$ ∈ [10^-1, 10]. Note that, either driven by passive diffusion or by active transport, the region that maximizes the total noise, ${\eta }_{c_A}^2$, is not accessible: $\tilde{D}$ ~ 5 · 10¹ − 10². Notice also that this range of diffusion values corresponds to the dynamics that produce separated peaks in the distribution of the autoinducer (see Figure 4). While our modeling is certainly very simple and the derived consequences should be carefully taken, this observation suggest that bacteria have developed mechanisms for coping with the noise and keep their functional QS regime away from the region where ${\eta }_{c_A}^2$ > 1. In this regard, let us point out that in every form of information exchange the precision is key. If the precision of the information is fuzzy then the related biological function lacks robustness. Thus, our results point towards the direction that bacteria have adapted their communication mechanisms in order to improve the signal to noise ratio.

Diffusion can lessen the effects of fluctuations, both in eukaryotes [32] and prokaryotes [23, 60]. This behavior is certainly obtained in our model: if the diffusion rate is larger than a given amount, then the noise decreases as the diffusion increases. However, to the best of our knowledge, the reverse effect, i.e. the noise increases as the diffusion increases, had not been reported. These opposed behaviors are responsible of the non-monotonic comportment of the noise described herein. One may wonder the reason underlying this phenomenology. By manipulating Eq. (18), it is easy to demonstrate that a) the slope of the noise at $\tilde{D}$ = 0 is always positive and b) in the white noise limit, τ _c → 0, the slope becomes null and the transcriptional noise decreases monotonously as the diffusion increases. Consequently, the observed behavior is due to the colored character of the transcriptional noise, i.e. due to a competition of temporal scales: those of the diffusion and the noise correlation time. This result opens the possibility of finding a similar phenomenology in other systems subjected to relevant levels of transcriptional noise and diffusive signals. In addition, it points up the relevance of non-memoryless noisy sources in biological systems [11]. Moreover, it shows that our results do not depend on details as the precise number of transcripts (one in our case). As long as there are distinct transcriptional phases inside cells with characteristic time scales, that is, presence versus absence of transcripts, the observed phenomenology is, qualitatively, the same. Yet, those phases are prone to appear when the autoinducer is produced at constitutive levels and the number of transcripts is small. Herein we restrict ourselves to the case of one transcript simply because we are able to obtain analytical results in that situation.

Conclusions

Herein we have explored the role played by cell-cell communication and transcriptional noise in QS systems near the activation threshold. Within this context, we have shown that the dynamics of the signaling molecule exhibits different behaviors depending on the diffusion coefficient. When increasing the rate of diffusion, the probability distribution of the autoinducer changes from single-peak distribution (sustained bursts dynamics), to monotonically decreasing distribution (bursts dynamics), to double-peak distribution (bistable dynamics), and finally to narrow single-peak distribution (diffusion-averaged dynamics).

In addition, we have shown that the mRNA dynamics plays a crucial role for regulating the total amount of molecular noise of the signaling molecules. Transcriptional noise is the main contribution to the total noise for a large range of diffusion values, $\tilde{D}$ < 10⁴. Only for very large values of the diffusion the intrinsic noise is the major source of stochasticity. Due to a competition of temporal scales, the total noise shows a non-monotonic behavior as a function of the diffusion rate. For large values of the diffusion coefficient, the total noise decreases as the diffusion rate increases. In this regard, our results are to be compared to previously reported noise reduction mechanisms, as for example in the case of bistable genetic switches coupled by QS communication [60]. On the other hand, when the diffusion rate is small enough compared to the characteristic rate of transcription events, the total noise increases as the diffusion becomes larger. The values of the diffusion rates in QS systems fall into two distinctive categories: either large values corresponding to passive transport mechanism or small values when an active transport mechanism applies. Surprisingly, these two QS classes avoid diffusion rates that maximize the total noise. According to these results, we conjecture that bacteria have engineered the communication mechanism for reducing the signal-to-noise ratio and produce a more reliable information exchange.

Our final comment refers to the possibility of considering other sources of stochasticity. Cell-to-cell variability and extrinsic noise have been proved to act as an important contribution in many cell processes [16, 25, 61, 62]. In the context of the problem studied herein, we can envision that variability, either at the level of the mRNA dynamics or at the level of the diffusion rate, can effectively lead to significant changes in the reported phenomenology. In addition, by considering additional steps in the synthesis of the autoinducer, the levels of intrinsic noise would increase. Whether or not this extra level of fluctuations, coupled with feedback regulation, may generate new effects in the framework of QS is not known. Work in those directions is in progress.