Bacteria Floc, but Do They Flock? Insights from Population Interaction Models of Quorum Sensing

Uncoupled model.We introduce a general model that focuses on the following four aspects of each cell: (i) x_i, the position of cell i; (ii) v_i, the velocity of cell i; (iii) p_i, the protein expression of cell i; and (iv) A_i, the concentration of autoinducers surrounding cell i. The model is a system of 4N ordinary differential equations that tracks the dynamics of the four above-named state variables. N represents the total number of cells, as follows:dxidt=vi (1)dvidt=λ1N∑j=1Nk1(xi,xj)(vj−vi)+F0(υoe^s−vi) (2)dpidt=λ2N∑j=1Nk2(xi,xj)(pj−pi)+L0(p∞−pi)ψA(Ai) (3)dAidt=λ3∑j=1Nk3(xi,xj)r−ku(xi)Ai−kdAi (4) with the corresponding weights, k1(xi,xj)=k2(xi,xj)=k3(xi,xj)=1(1+|xj−xi|2)β logistic function, ψA(Ai)=1[1+ec0(c1−cAAi)]a and uptake functionku(xi)=a1exp⁡(−{a[γ1(∑j=1N|xj−xi|2β)]−b}22c)+a2 (5)

This model is referred to as “uncoupled” because the protein expression equation (equation 3) is not coupled with the velocity equation (equation 4) (i.e., protein expression does not influence movement), unlike in the coupled model, which will be discussed later.

The weight functions k₁, k₂, and k₃ are symmetric, which will prove useful in the subsequent asymptotic analysis. These functions decrease as the distance between the cells increases. The exponent in the denominator (β) controls how fast the functions decrease from 1. Qualitatively, β represents the scope of influence that the neighboring cells have on the individual. So, low β values (such as β = 0) would model a system in which the individual is influenced by all cells.

The logistic function Ψ_A(A_i) acts as a smoothed piecewise constant function and takes values ranging from 0 to 1 depending on the surrounding autoinducer concentration (A_i) of cell i. The uptake function is modeled as a shifted Gaussian. k_u(x_i) is constructed in such a way so that the uptake rate increases as the density increases until reaching a maximum. Thereafter, the uptake decreases when the density is high. This is to reflect a decrease in metabolic activity caused by overcrowding (i.e., decreased nutrient availability when cells are crowded). With this Gaussian, we assume that there is a density range in which uptake is highest.

Equation 1 provides the definition of velocity. Equation 2 models the changing dynamics of the velocity of cell i by employing a flocking (population interaction) term and a source term. The flocking term addresses how the velocities of neighboring cells influence the individual's velocity. The weight function k₁ as well as the parameter β determines how much influence that the other cells, j, have on cell i. Next, the source term models chemotaxing behavior in which the cell runs toward the higher concentrations of chemoattractant at a velocity of v_oê_s, where ê_s represents the unit direction vector. The strength of this term is controlled by the constant F₀. The combination of both the flocking term and the source term generalizes the movements of these cells in response to the chemical stimuli. Note that depending on F₀, the source term influences the movement of cells much more than the flocking term. In this way, a random perturbation by one individual, a small collection, or a perturbation that reflects an environmental cue can lead to global directional changes for the population, and these may accurately reflect natural circumstances.

Equation 3 tracks the dynamics in autoinducer-triggered protein expression, which we use as a marker of the collective output of the QS phenotype, such as biofilm formation (40) or agglomeration (32). To motivate the existence of equation 3, for the sake of a concrete example, we use the expression of bioluminescence by Vibrio fischeri (41). In this species, the autoinducer AI-1 (an acyl-homoserine lactone [AHL]) controls the regulatory system (LuxI/LuxR) that governs the production of both the autoinducer and various proteins. LuxI is an AI-1 synthase, and LuxR is an activator regulating transcription of the luxCDABEG operon. The luxCDABEG genes code for the regulation of bioluminescence, including the expression of the luciferase enzymes that are responsible for light emission. The flocking term then serves to average the bioluminescence of the surrounding cells. The additional term, which we call the source term, represents the bioluminescence triggered when the surrounding concentration of autoinducers passes a certain threshold. Above this threshold, the logistic function Ψ_A approaches 1, thereby activating this source term. The strength of this term is controlled by the constant L₀. In our model, flocking plays only a small part in the cell dynamics when the source term is in effect (i.e., the logistic function is not zero). That is, flocking behavior is more pronounced before the source term is activated.

The changes in the surrounding concentration of autoinducers of a particular cell are represented by equation 4. The equation is constructed as follows:dAidt=synthesis−uptakes−sloughing off/degradation of autoinducers=λ3∑j=1Nk3(xi,xj)r−ku(xi)Ai−kdAi and this is phenomenologically similar to several previous models (7, 20).

Uncoupled model: simulations and discussion.The model presented was constructed to analyze both the movement of cells toward a chemoattractant and autoinducer-triggered protein expression. To recap, this system consists of 4N differential equations, in which the four state variables, namely, position, velocity, protein expression, and surrounding autoinducer concentration, are assigned to each cell. N represents the total number of cells. The equation governing the surrounding concentration of autoinducer of cell i consists of a weighted synthesis term, an uptake term, and a degradation term. The uptake function can be customized to reflect the system. We initially used a logistic uptake function that was dependent on time (see property iv in Text S1 in the supplemental material (“Asymptotic Analysis”) to prove that the surrounding concentration of autoinducers converged to a constant value, which differs for each cell. Our assumption was that the cells grow sufficiently dense over time and that the uptake increases concomitantly and remains steady at an elevated level. The purpose of this logistic function was to find a simple representation of this general behavior that was with respect to time, instead of position, so that one could perform the asymptotic analysis of the differential equation.

In Fig. 2, we show how a randomized distribution of cells with randomized initial velocities ultimately converges toward a uniform velocity trajectory. In the plot, panel a reveals the rapid alignment of velocities toward the value v_oê_s as time progresses. The velocities are parsed into x and y directions. Panel b includes trajectories of 100 cells in a 2D plane, the color indicating the level of protein expression for each cell, and this shifts from green to yellow as the QS-mediated activity is actuated. The final position in each trajectory is indicated by a dot, the color of which indicates the extent to which QS-actuated protein was accumulated. This simulation is purely a representative case, with constants indicated in the figure caption; it indicates generally the propensity by which the cells reached an asymptotic velocity and exhibited QS-mediated gene expression.

• Open in new tab
• Download powerpoint

FIG 2

Simulated trajectories of a randomized initial population. (a) Cells exhibit a convergence of velocity in the x direction (V_x) and y direction (V_y), with λ₁equal to 5 and β equal to 0.2 for 100 cells and with initial conditions randomly generated from a uniform distribution on the interval ([0, 20] and [0, 30] for V_x and V_y, respectively). (b) The cells are represented by circles. The trails represent the cell trajectories, and the emergence of yellow from green represents the emergence of QS-mediated gene expression. Note that the velocities of cells remain bounded (shown here for a fixed time [the proof is in Text S1 in the supplemental material]). a.u., arbitrary units.

For these simulations (x_i in ℝ², v_i in ℝ², p_i in ℝ², and A_i in ℝ²), we used a shifted Gaussian for the uptake function in which, at low cell densities, the uptake is low (the lowest being a₂). At higher densities, the uptake increases. At even higher densities, the uptake begins to decrease. The reason for this is to account for the effect of overcrowding (e.g., nutrient depletion) on metabolic processes. Also, this allows for more interesting dynamics in the surrounding concentration of autoinducers, which can lead to a higher variance in protein expression among the population, as seen in Fig. 3, where the emergence of subpopulations of cells can be observed. This is more representative of AI-2-mediated QS, wherein there is a negative-feedback loop governing AI-2-mediated protein expression (1). A monotonically increasing uptake function with decreased distance would not result in this behavior, and this would be more representative of AI-1-mediated QS. Note also, however, that the variation in protein expression here is also attributed to the logistic function Ψ_A(A_i), which, together with the autoinducer concentration, governs the per-cell protein expression in the absence of neighbors.

• Open in new tab
• Download powerpoint

FIG 3

Histogram of QS-mediated protein expression in 100 cells. The logistic function Ψ_A(A_i) and the Gaussian uptake expression in the model enable the bimodal expression variation that is often observed for AI-2-mediated QS. The histograms depicted reveal a simulated bimodal distribution at three arbitrary time points (t₃ > t₂ > t₁). Parameter values for the Gaussian uptake function are as follows: a equals 4, b equals 4, and c equals 27.81, with β equal to 0.62 and γ equal to 1/2,000.

The β parameter in the model can be interpreted as the weight placed on the behaviors of the neighbors. When β is equal to 0, the individual cell places equal weight on all other cells. As β increases, the individual cares less about the behaviors of the cells further away. To illustrate, we varied β and illustrate resultant behaviors in order to qualitatively hone in on obtaining realistic β values. This type of simulation also illustrates the ease with which population-based behavior can be modified by simple adjustment of one parameter, β. In Fig. 4, the β values analyzed were 0.55, 0.6, 0.65, and 0.7. For the lower values of β (i.e., 0.55 and 0.6) (Fig. 4a and b), the cells at a higher density (those closer together and not on the periphery) express more protein (they are more yellow) than the cells at a lower density. That is, the cells at lower density are seen on the outskirts of the group, and these are all green, as are their entire trajectories. This is particularly well visualized in Fig. 4a, where the plot’s center is predominantly yellow, the color that is the visual sum of the colors of the individual trajectories that at times intersect. In Fig. 4b, the QS-mediated yellow is a bit more focused in the middle of the plot. This heuristic graph represents the prototypical QS response wherein the cells that are closer to others and experience the highest levels of autoinducer exhibit the greatest QS genetic response. As noted above, this is easily described by parameter β and the Gaussian shape of the uptake function, which is determined by the constants a, b, c, and a₂, providing the range of distances between cells over which the autoinducer uptake is highest. Note that as β increases to 0.65 and 0.7 (Fig. 4c and d), the simulated result is largely reversed; cells at a lower density are observed to express more protein, which is counterintuitive in QS. This too is due to the Gaussian uptake function. As β is increased further, the cells with a greater difference from other cells (i.e., large |x_i(t) – x_j(t)|) are less influential on each other and, importantly, become more individualistic in their phenotype. Also, numerically, the autoinducer uptake rate approaches 0 for these cells. Since autoinducers are continually synthesized according to the model, QS-mediated proteins are still being expressed. The cells at a higher density continue to take up the autoinducers as the distances between cells become smaller. Correspondingly, they experience lower autoinducer levels and lower metabolic activity, etc., and express less QS-mediated protein. Thus, in this configuration at higher values of β, the autoinducer-triggered protein expression would be lower for cells of higher density than for cells of lower density, which is not observed experimentally. Importantly, we include these simulations to illustrate both the model simplicity and the flexibility in the mathematical formalism ascribed to the simple adjustment of parameter β and the Gaussian uptake function.

• Open in new tab
• Download powerpoint

FIG 4

Trajectories of 100 cells initially randomly distributed but with different dependencies on their nearest neighbors. (a to d) Cell paths for a β of 0.55, 0.6, 0.65, and 0.7, respectively. The circles represent the cells, and the lines represent the paths that the cells have travelled (from time zero). The colors represent levels of QS-mediated protein expression (green is low, and yellow is high). At low β values, cells that are closer together are seen to be more yellow; they express more protein as a function of the autoinducer uptake. As the β value increases, however, there is a reversal in behavior. The cells which are at a higher density (middle of plots) express less protein than the cells at a lower density (on the periphery). This functional flexibility in the simulation is due to the Gaussian uptake term, in which the β value decreases the importance of the sum of the distances between cells. The cases observed in natural settings are more aligned with lower β values (e.g., a β of 0.55 or 0.6).

With this model, which phenomenologically describes the natural phenomena, we have introduced a method of modeling quorum sensing that does not a priori require knowledge of the intrinsic gene and protein interactions that are necessary to communicate with other bacteria. Rather, the methodology depends only the distances between cells and simple parameters that provide the relative importance of those distances. That is, the flocking terms and similar density-dependent weights, as introduced here, serve to generalize the behavior of quorum sensing. We assumed here that QS marker protein expression in this uncoupled model formalism does not influence velocity. This keeps the ordinary differential equations simple enough to perform an asymptotic analysis (see Text S1) while encapsulating the dynamics of bacterial communication in tandem with chemotaxis.

Coupled model.The uncoupled model enabled the first mathematical representation of QS-mediated behavior as a form of flocking and was confirmed by asymptotic analyses (see Text S1). We now add a layer of complexity by creating a coupled model. By coupled, we mean the coupling of the velocity equation and the protein expression equation. Chemotaxis is the movement of bacteria triggered in response to a gradient in a concentration of particular molecular substances which are classified as either attractants or repellents. The changes in the concentrations of these substances activate a cascade of signal transduction mechanisms that affect motility (42).

Not much is known about the direct connections between the quorum sensing and chemotaxis systems, including for E. coli. Connections tend to describe the interactions or effects of one key protein in one system influencing the other system in some way. For example, AI-2 is a known chemoattractant for E. coli (43, 44), and the LsrB protein, which is responsible for AI-2 uptake, is necessary for chemotaxing to AI-2 (45). Thus, LsrB of the quorum-sensing system is integral to the chemotaxis system with regard to the movements toward AI-2. Also, upregulation in the motility and production of flagella was observed in E. coli luxS knockouts (46). That is, the regulatory systems in quorum sensing and chemotaxis consist of many genes, effector molecules, and proteins that contribute to the dynamical behaviors, as well as resultant activities (e.g., phosphorylation, induction, etc.) that are involved in these processes (47). Taking into account the various reactions and regulatory mechanisms of these systems can be complicated, as they are ill-defined, so instead, we aimed to generalize the behaviors of the quorum-sensing and chemotaxis systems again with the help of density-dependent weighted terms. Our revised model is represented in equations 6 to 9.dxidt=vi (6)dvidt=λ1N∑j=1Nk1(xi,xj)(vj−vi)+F0[υoe^s+(υ∞−υo)e^sψp(pi)−vi] (7)dpidt=λ2N∑j=1Nk2(xi,xj)(pj−pi)+L0(p∞−pi)ψA(Ai) (8)dAidt=λ3∑j=1Nk3(xi,xj)r−ku(xi)Ai−kdAi (9)

The velocity equation has been modified to be a function of protein expression. The strength of the source term for velocity is still governed by the constant F₀, and ê_s remains the unit direction vector. We have incorporated a logistic equation, Ψ_p(p_i), that takes values between 0 and 1 depending on the protein expression corresponding to that cell. Ψ_p(p_i) couples the protein expression to velocity, whereas Ψ_A(A_i) couples autoinducer activity to protein expression.ψp(pi)=1[1+ec0(c1−cppi)]α1ψA(Ai)=1[1+ec2(c3−cAAi)]α2

If Ψ_P is 0 for all t’s, the source term reduces to F₀(v_oê_s – v_i). If Ψ_P is 1 for all t’s, the source term becomes F₀(v_∞ê_s – v_i). v_o is the magnitude for the initial velocity. After the protein expression approaches a certain level, the cell's velocity increases, approaching v_∞ if protein expression, p_i, continues to increase. Protein degradation is assumed to be minimal. In the scenario when the source term for the protein expression equation is initially turned on [i.e., Ψ_A(A_i) takes on a nonzero value] but then turns off and remains off {Ψ_A[A_i(t)] is 0 for all t’s above t_threshold, where t_threshold is the time when the autoinducer concentration drops below the threshold value}, we note that protein degradation should be taken into account for a more accurate analysis of long-term cell dynamics.

Coupled model: discussion and application.The 4N system of differential equations can be customized for specific quorum-sensing systems by using different threshold values, which are controlled by changing α₁ and α₂ (or, alternatively, c_p and c_A), modifying the scopes of influence (β₁, β₂, β₃), customizing the uptake function in equation 9, and toggling the strength of each of the flocking and source terms (λ₁, λ₂, λ₃, F₀, L₀). First, we use the same uptake function (5) as in the uncoupled model above.

When the α values were equal (e.g., equal influence of velocity on protein expression and influence of autoinducer on protein expression), we observed an interesting dynamic within our cell population: two subgroups of cells had formed. Because this was not observed in the previous model under any condition, we conclude that this behavior was attributed to the coupling of quorum sensing and chemotaxis systems. This was an interesting finding that emerges only because of the dependence of velocity on the QS protein level. That is, in the uncoupled model, the cells exhibited trajectories in one single group, with only a portion of the group expressing more protein than the rest. Since quorum sensing does not influence velocity in that model, cells with higher expression would not move faster. The lengths of the trajectories in Fig. 2 and 3 were essentially similar. When we couple the model such that protein expression influences velocity, the portion with higher expression was seen to move faster, leaving the rest of the group behind. Not only did we observe that different cells had different velocities, in Fig. 5, we found the formation of two distinct subgroups. That is, this figure presents simulations starting from a randomized initial state (Fig. 5a, or similarly, the initial conditions in Fig. 4) and progresses in time (to Fig. 5f). That is, in Fig. 5b, a subpopulation of cells emerges from the group by moving to the upper right. This emerging subgroup of cells remains tightly formed as if in a flock. It is noteworthy that in this configuration, the QS-mediated protein expression is maintained. In essence, this also provides a means by which a subgroup of cells accelerates away from the others. This type of behavior, in general, is reminiscent of a bird flock.

• Open in new tab
• Download powerpoint

FIG 5

Simulations of a coupled model where velocity is linked to QS-mediated protein expression. (a through f) Population distributions as the population progresses in time. Cell populations were randomly distributed initially, as in Fig. 4. Two groups of cells emerged early on with different velocities. (b) The dense portion of cells at the upper right reaches a higher velocity. As it moves through panels c and d, it splits away from the rest of the group. At later times (e and f), the two groups seemingly remain as distinct groups. These phenomena were observed experimentally.

The organization of cells into two groups, here, is consistent with experimental observations (35, 39). In chemotaxis experiments conducted over 50 years ago using capillary tubes (39), E. coli cells were placed at one end, and the chemoattractant (galactose or oxygen) was placed at the other end. E. coli moved toward galactose, as expected. What was unexpected was the appearance of distinct migrating bands. This is exactly analogous to the simulations in Fig. 5. The author, Adler (39), made an attempt to uncover a mechanism for this partitioning. He attributed the results to the varied consumption rates of nutrients and their placement at the distal end of the capillary tube. In one case, the first band of cells depleted oxygen while consuming part of the galactose, and the second band depleted the remaining galactose in what was remaining as an anaerobic environment. Interestingly, there was also a nonmotile fraction that stayed in place at the proximal end. To prove that this behavior was not attributed to individual cell heterogeneity, Adler scraped the cells from each of the two bands and repeated the experiment with those cells. Two groups of cells formed, yet again. He did not discuss or hypothesize how a subfraction of cells and not the entire population had emerged from the initial population in the first place.

Increased QS-dependent protein expression that is caused by the increased autoinducer locally near a high density of cells might lead to the emergence of a small group of cells transiting out of the larger group, especially if their velocity is linked to QS. This mechanism provides for the formation of two sets of cells during the process of chemotaxis, as described by Adler. Our model, for the first time using flocking formalisms, qualitatively reveals this emergent flocking behavior. Our simulations reinforce the idea that quorum-sensing systems, or analogous density-dependent mechanisms, are connected to the chemotaxis system of E. coli.

We next postulated thought experiments in which we interrogated the model seeking to observe the potential for QS-mediated microbial recognition. That is, what kind of response would we see when α₁ was not equal to α₂, reflecting scenarios where the threshold levels required for protein-coupled motility or autoinducer-coupled protein expression are varied relative to each other. Or, in what kind of experiment would we observe a scenario in which α₁ is not equal to α₂? We first respond to the latter question by describing the experiments conducted in a transwell apparatus by Servinsky et al. (35). There are two sets of cells in these experiments (Fig. 6). One set is “sentinel” cells that secrete AI-2 (Fig. 6a). These are placed in capsules that, in turn, are placed in the lower chamber of the transwell. The other set of cells, called “recruitable” cells, are genetically engineered to produce green fluorescent protein (GFP) when transcription is induced by AI-2 (48). Located in the upper chamber, these cells also constitutively express DsRed as a means of indicating their accumulation in the lower chamber of the transwell. These recruitable cells are engineered not to produce AI-2 (lux mutant) so that the only source of AI-2 was the sentinel cells located in the lower chamber. The transwell experiments therefore enabled us to track the movement of recruitable cells toward AI-2 (chemoattractant) by monitoring red fluorescence and, simultaneously, QS-dependent protein expression (AI-2-mediated expression of GFP). The experimental results (Fig. 6b) revealed cells (∼350) with DsRed, but not GFP, in the lower chamber at the 16-h time point. However, both GFP and DsRed had appeared by the 40-h time point (e.g., ∼15% green). The main conclusions in the work of Servinsky et al. (35) were, indeed, that a subset of cells moved toward a new locale that was, in turn, marked by the presence of autoinducer, AI-2, that had been secreted by distant sentinel cells. Also, because of the sequential nature of motility and QS-mediated gene expression, the work demonstrated further that the motility response was more sensitive than the quorum-sensing-triggered protein expression response with respect to AI-2. In such a scenario, α₁ is less than α₂. Such an observation was also noted in the work of Wu et al. (29), wherein AI-2 was generated on the surfaces of cancer cells so that payload-synthesizing bacteria could swim to those cells, synthesize a drug, and deliver it specifically to cancer cells and not healthy cells (devoid of surface-synthesized AI-2). There too, the AI-2 threshold for motility was lower than lsr-mediated gene expression.

• Open in new tab
• Download powerpoint

FIG 6

Encapsulated sentinel cells secrete AI-2 (orange) and draw recruiter cells through a membrane and into the bottom of a transwell apparatus. Data from the work of Servinsky et al. (35) are in panel b. (a) In the transwell experiments, initially, the recruitable cells are in the upper chamber and the sentinel cells are sequestered in alginate beads and placed in the lower chamber; the upper and lower chambers are separated by a membrane. The recruitable cells (ΔluxS) constitutively express DsRed and express enhanced green fluorescent protein (egfp) in the presence of AI-2 (via the lsr promoter). (b) The bar chart shows the number of cells in the lower chamber at hour 16 and hour 40. The fraction of cells expressing GFP is also indicated (horizontal green bars). The cells moved first, and subsequently, QS-mediated protein expression was actuated. (c) Our simulation results are qualitatively similar: the cell number in the lower chamber increases first, and the fraction of cells expressing protein increases later. Note that time is in arbitrary units (a.u.).

Consequently, we modeled an analogous system in which chemotaxis was more sensitive than the (quorum-sensing-mediated) protein expression. This can be modeled by lowering the threshold of Ψ_p(p_i) (i.e., the value of p_i for which the logistic function increases to 1) and raising the threshold of Ψ_A(A_i). That is, α₁ was set to less than α₂. Here, we use a modified uptake function in which the measure of density is the inverse of the sum of k₃.ku(xi)=a˜1exp⁡(−{a˜[γ2∑j=1Nk3(xi,xj)]−b˜}22c˜)+a˜2 (10)

In this simulation (Fig. 6c), we mimicked the upper and lower chambers of the transwell apparatus in the work of Servinsky et al. (35) by assigning a fixed distance from the origin of each cell trajectory to represent the membrane in the experiment. We recorded the number of cells and their protein expression level when they had swum past this membrane. We classified the cells that expressed greater than 0.7 arbitrary unit (AU) in our model as fully expressing GFP. Figure 6c shows the increasing cell count in the lower chamber of our simulation as well as the increasing fraction of cells fully expressing GFP. More specifically, the number of cells that had migrated to the lower chamber had increased dramatically by the second time point (24 AU), and the GFP-expressing fraction increased appreciably only after the third time point (32 AU). We qualitatively compared these simulations to the results of Servinsky et al. (35), shown in Fig. 6b, and we observe similar results: as time passed, the cell count increased but GFP was not expressed until later. By the incorporation of flocking interactions and by accommodating signal-mediated sensitivity for chemotaxis and QS-mediated gene expression, our flocking model for the first time was able to demonstrate the potential for a subset of bacterial cells to seek out new environments and, upon entrée into the new locale, to exhibit classic QS behavior. We know of no other modeling efforts that have described these phenomena.