Bacterial quorum sensing is often mediated by multiple signaling systems that interact with each other. The quorum-sensing systems of Pseudomonas aeruginosa, for example, are considered hierarchical, with the las system acting as a master regulator. By experimentally controlling the concentration of auto-inducer signals in a signal deficient strain (PAO1ΔlasIΔrhlI), we show that the two primary quorum-sensing systems—las and rhl—act reciprocally rather than hierarchically. Just as the las system’s 3‑oxo‑C₁₂‑HSL can induce increased expression of rhlI, the rhl system’s C₄‑HSL increases the expression level of lasI. We develop a mathematical model to quantify relationships both within and between the las and rhl quorum-sensing systems and the downstream genes they influence. The results show that not only do the systems interact in a reciprocal manner, but they do so asymmetrically, cooperatively, and nonlinearly, with the combination of C₄‑HSL and 3‑oxo‑C₁₂‑HSL increasing expression level far more than the sum of their individual effects. We next extend our parameterized mathematical model to generate quantitative predictions on how a QS-controlled effector gene (lasB) responds to changes in wildtype bacterial stationary phase density and find close quantitative agreement with an independent dataset. Finally, we use our parameterized model to assess how changes in multi-signal interactions modulate functional responses to variation in social (population density) and physical (mass transfer) environment and demonstrate that a reciprocal architecture is more responsive to density and more robust to mass transfer than a strict hierarchy.

Fig 1. The Pseudomonas aeruginosa QS regulatory network is typically viewed as a hierarchy with the las system on top.

In general, the regulatory network of two QS systems may be organized in three architectures: (A) Independent, in which the signal of each has no influence on the expression of synthase or receptor in the other. (B) Hierarchical, where one system’s ssignal influences expression of the other’s components but without reciprocation. And (C) reciprocal, where both systems’ signals influence the others’ components. The consensus of the review literature for P. aeruginosa (D) as summarized in 17 review papers published since 1996 (Tables A and B in S1 Text) suggests a hierarchical architecture for the las and rhl systems.

https://doi.org/10.1371/journal.pbio.3003316.g001

The las and rhl systems interact reciprocally

We first evaluate QS behavior under the influence of a single signal. We establish a baseline expression level by measuring reporter luminescence with no signal present. We then observe the increase in luminescence as exogenously controlled signal concentration increases. We chose signal concentrations to reflect experimentally observed concentrations typically ranging from 0.1 to 5 μM [17]. The ratio of luminescence with signal to luminescence with no signal represents the fold-change in expression induced by the defined signal concentration. Fig 2A and 2D shows the results for 3‑oxo‑C₁₂‑HSL. As expected, expression of both lasI and rhlI genes increases as signal concentration increases. The availability of the las signal molecule influences the expression of rhlI as well as lasI, and, therefore, the las system affects the rhl system.

While we find no surprises with 3‑oxo‑C₁₂‑HSL, our experiments with C₄‑HSL challenge the conventional hierarchical view. Fig 2B and 2E shows those results: expression of lasI and rhlI increases with higher C₄‑HSL concentration. The response of lasI (Fig 2B) does not correspond to a strict hierarchy with las as the master. Here we find that the rhl system affects the las system.

To quantify the impact of each signal alone, we model gene expression using Michaelis–Menten kinetics under quasi-steady state assumptions (Methods), allowing us to estimate both the maximal fold-change under signal activation and sensitivity to signal (Table C in S1 Text). Our model fits (Fig 2A, 2B, 2D, 2E and Fig D in S1 Text) indicate that while the las and rhl systems have reciprocal impacts, those impacts are not symmetrical. The las signal 3‑oxo‑C₁₂‑HSL has a substantially greater influence on gene expression than C₄‑HSL. In both cases the potential fold-change from 3‑oxo‑C₁₂‑HSL is approximately six times greater than the potential fold-change from C₄‑HSL. Both lasI and rhlI are also more sensitive to 3‑oxo‑C₁₂‑HSL than to the C₄‑HSL as the concentrations required to reach half of maximal expression are roughly 4 times and 30 times higher for the latter.

Fig 2A, 2B, 2D, and 2E considers the effects of each signal in isolation, but wildtype cells with functioning synthase genes can produce both signals. To understand environments where both signals are present, we use controlled concentrations of both signals in combination. Fig 2C and 2F presents those results in the form of heat maps. The qualitative responses of both genes are similar: raising the concentration of either signal increases expression regardless of the concentration of the other signal. As with our observations of C₄‑HSL alone, these results demonstrate again that the rhl system (via C₄‑HSL) affects the las system (lasI expression).

Quantifying the las and rhl interactions reveals non-linear and synergistic effects

Having established a simple Michaelis–Menten model for each signal in isolation (Fig 2A, 2B, 2D, and 2E), we next consider whether that model is sufficient to explain the effect of the signals in combination (Fig 2C and 2F). Can we estimate total expression as the sum of expression induced by each signal alone? Such a response could result from two independent binding sites in the promoter regions [18], one site for LasR/3‑oxo‑C₁₂‑HSL and a separate site for RhlR/C₄‑HSL. Fig 3A and 3B clearly shows that we cannot. The maximum expression observed, shown as a “ceiling” in that figure’s panels, far exceeds the sum of the signals’ individual influence.

To account for the synergy between the signals, we incorporate a cooperativity term in the gene expression model. Note that the cooperativity term is a multiplication of signals, and it alone cannot explain the full response, as the product is necessarily zero when any signal is absent. This term accounts for any non-additive interaction, for example the ability of one bound transcription factor to recruit the binding of a second transcription factor [19]. Eq 1 (annotated in Fig 4) shows the result. Each gene has a basal expression level, amplification from each signal alone, and additional amplification from each pair-wise combination of signals. The interaction from these pair-wise combinations captures the cooperative enhancement from the combined signals.

(1)

For both lasI and rhlI we again minimize the sum of squared errors to estimate parameter values (non-linear regression using the Gauss-Newton algorithm). The resulting multi-signal models in Table D in S1 Text have R² values of 0.82 and 0.77. Fig 3C and 3D compares the model estimates with observations. For both genes, the model captures the effect of each signal in isolation and also both signals in combination.

The parameter estimates in Table D in S1 Text quantify the relative effect of individual and combined signals. The maximum expression induced by both signals nearly doubles compared to the maximum expression induced by any signal alone. Fig 3E summarizes the model parameters graphically and revises the consensus view of Fig 1D—the rhl system does influence the las system—and it shows the relative magnitudes of the effects.

las and rhl synergy also shapes control of quorum-sensing effector genes

Our results, summarized in Fig 3, establish that both AHL signals influence the expression levels of both synthase genes in a synergistic manner. We note that the methods outlined in Figs 2–3 for the lasI and rhlI synthase genes can be applied to any QS-controlled gene of interest and can quantify dual (or more) signal control over expression levels. Here we look at lasB, a classic QS effector gene that codes for the secreted digestive enzyme LasB and is widely used as a model of QS-mediated virulence [20,21] and cooperation [22–24].

The prior literature provides a clear expectation that lasB expression will be positively influenced by both AHL signal molecules [3,25], but does not provide a quantitative partitioning of the relative importance of each signal’s contribution—alone and via synergistic effects. Using the approach outlined above, we measure luminescence of a lasB reporter in an AHL signal null strain exposed to defined, exogenous concentrations of both AHL signals, revealing maximal activation under dual-signal exposure (Fig 5A and 5B). As illustrated for lasI and rhlI (Fig 3), fitting the baseline single signal model (Eq A in S1 Text) cannot account for the maximal expression of lasB (Fig 5A), while the multi-signal model (Eq 1) can fit the data (Fig 5B). Fitting Eq 1 to the lasB data (Table E in S1 Text) supports previous conclusions that 3‑oxo‑C₁₂‑HSL has a stronger effect than C₄‑HSL on lasB expression [26,27]. Our results add that lasB expression is dominated by the synergistic combination of both signals rather than either signal in isolation. We further note that lasB expression is more sensitive to C₄‑HSL alone than to 3‑oxo‑C₁₂‑HSL alone. (Fig 5C); the maximum increase in expression in response to C₄‑HSL is much less than to 3‑oxo‑C₁₂‑HSL, but C₄‑HSL requires a significantly lower concentration to achieve that maximum increase (Table E in S1 Text).

Mathematical models incorporating multi-signal interactions predict wildtype quorum-sensing response to environmental variation

Our parameterized gene expression model (Eq 1) predicts gene i’s expression E_i as a function of the AHL signal environment S (Figs 3C, 3D, 5B and Fig G in S1 Text), which leaves open the question of how the signal environment S relates to underlying dimensions of environmental variation—the hypothesized sensing targets of QS (population density [28], diffusion [29], efficiency [30], containment [31], genotypes [23,32], combinations [33], etc.)

To connect our gene expression model to critical environmental dimensions of bacterial population density and mass transfer (e.g., diffusion or advective flow), we build on previous models of extracellular signal dynamics [33–37]. While our earlier experimental and modeling analyses (Figs 2, 3, and 5 and Eq 1) treated signal concentration as an experimentally fixed exogenous factor, we now treat signal concentration as a dynamical variable that emerges from the interplay of regulatory feedbacks (both within- and among-signal modules) and environmental conditions.

We assume that signal concentration increases in proportion to the corresponding synthase’s expression level, multiplied by the number of cells expressing synthase, and decreases due to a constant rate of decay and removal via mass transfer. These assumptions lead to the differential equation model of Eq 2 (annotated in Fig 6), where S_i is the concentration of signal i, E_i (S) is the expression level of the synthase for signal i (as a function of both signal concentrations, S, see Eq 1) and c_i a proportionality constant. N is the population density; δ_i is the decay rate of signal i, and m is the rate of mass transfer.

(2)

This equation models the dynamics of extracellular signal concentrations in response to environmental conditions defined by density N and mass transfer rate m. While it is possible to derive analytical equilibrium solutions to Eq 2 for either independent (Fig 2A, 2B, 2D, 2E and Eq A in S1 Text) or synergistic (Eq 1) signal activation E_i(S), the resulting solutions are cumbersome and do not yield clear insights into system behavior. (See “Analytic Solutions for Equilibrium” in S1 Text.)

To numerically solve for equilibrium signal values in Eq 2 that result from given values of population density and mass transfer, we derive parameter estimates for synergistic E_i(S) from our experimental data (Eq 1 and Table D in S1 Text). The remaining parameters, c_i and δ_i, we estimate from published literature as detailed in Table E and Fig H in S1 Text. We can connect these equilibrium environmental signal values to QS-controlled effector behavior using our parameterized model of lasB expression (Fig 5 and Eq 1 as parameterized by Table E in S1 Text). Integrating these steps allows us to predict how lasB expression varies in wildtype bacteria (with intact signal synthases) as the social (N) and physical (m) environment changes.

In Fig 7A and 7B, we illustrate our integrated model prediction of wildtype lasB expression as a function of change in stationary phase population density, N (solid line, Fig 7B), capturing the causal chain (Fig 7A) from environment N to signal levels S (Eq 2) and signals S to lasB expression (Eq 1 and Fig 5 and Table E in S1 Text). To test this prediction using an independent dataset, we use previously published data on wild-type lasB expression as a function of stationary phase density [38]. Our model prediction is in strong agreement with the independently derived data (R² = 0.91), demonstrating not only the utility of our approach but also providing theoretical support for the conclusions of the prior experimental study [38]: QS controlled behaviors in P. aeruginosa are not governed by a critical cell density or ‘quorum’ threshold separating off/on states, but are better described as producing graded or rheostatic responses to changes in stationary phase density.

Multi-signal architectures govern functional responses to environmental variation

Following an initial validation of our quantitative predictions linking environmental properties to wildtype QS-controlled behaviors (Fig 7B), we now turn to a counter-factual approach to assess how changes in model architecture translate to changes in the ability to sense and respond to differences in environmental conditions. Specifically, we ask: how do reciprocal versus hierarchical versus independent architectures (Fig 1) influence QS signal and effector responses to variation in bacterial density and mass transfer? To answer that question, we recognize that Eq 1 is general enough to model hierarchical and independent architectures in addition to reciprocal architectures. Those alternatives emerge when specific ɑ parameter values are equal to zero.

Starting with our full reciprocal activation model (Table D in S1 Text), we can derive a las-controlled hierarchical model by setting ɑ_1,2 and ɑ_1,1,2 to zero (removing influences of the rhl system on the las system). To derive an independent model, we additionally set ɑ_2,1 and ɑ_2,1,2 to zero (removing influences of the las system on the rhl system). These adjustments (parameters summarized in Table F in S1 Text) allow the multi-signal architecture to represent the alternate models Fig 1 depicts, but they also diminish the total weight of QS activation processes. To compensate for this diminished activation, we also examine rescaled hierarchical and rescaled independent parameter sets, where the maximal activation weightings are held constant across models (parameters in Table G in S1 Text).

Using our alternate “counterfactual” models, we examine in Fig 7C–7H how equilibrium lasB expression changes over a range of population densities N and mass transfer rates m, for all three architectures. Note that in Fig 7C (reciprocal architecture) with m = 0, we recover the graded response to increasing stationary density N predicted theoretically and confirmed empirically in Fig 7B. Looking across alternate models without rescaling (Fig 7C–7E), we see that removing activating processes from rhl to las (hierarchical model) and also from las to rhl (independent model) leads unsurprisingly to a general weakening of QS maximal response. In the re-scaled models (Fig 7F–7H), we normalize maximal expression, so we can contrast functional responses to changes in density and mass transfer. Comparing these panels, we can see that changes in architecture change the functional response to environmental variables. In comparison to the hierarchical or independent models, the reciprocal architecture broadens the environmental parameter space where lasB expression is elevated (areas above 5% black line, and above 50% white line) and also reduces the slope of activation contours (angles of black and white lines). The change in contour slope indicates a greater robustness to increasing mass transfer, given a reciprocal activation architecture. In Fig I in S1 Text, we assess the temporal activation dynamics of lasB under the three alternate architectures and highlight that even under the rescaled model, the reciprocal architecture provides the most rapid activation process (given high density and low mass transfer).

Literature review

The PubMed database of the US National Institutes of Health was queried on 20 July 2021 using the query PubMed Search (“review”[Title/Abstract] OR “review”[Publication Type]) AND “quorum sensing”[Title] AND “pseudomonas aeruginosa”[Title/Abstract], resulting in 76 results with publication dates from 1996 to 2021. Papers that included a diagram of the gene transcription networks for the las and rhl QS systems were further analyzed to interpret the interactions present in those diagrams. Tables A and B in S1 Text show the results. Of the papers analyzed, all show the las system positively activating the rhl system, and none show the rhl system positively activating the las system.

Data collection

We used three strains for the experimental observations: lasB:luxCDABE genomic reporter fusion in NPAO1∆lasI∆rhlI, lasI:luxCDABE genomic reporter fusion in NPAO1∆lasI∆rhlI, and rhlI:luxCDABE genomic reporter fusion in NPAO1∆lasI∆rhlI. The transcriptional reporters are based on the mini-CTX-lux reporters [53]. Promoter regions of the genes of interest were cloned upstream of the lux operon, creating a transcriptional fusion reporter. The CTX system ensures a single chromosomally integrated reporter. We streaked out all strains in Luria–Bertani (LB) agar at 37°C for 24 h and then subcultured a single colony in 10 mL LB, incubated at 37°C under shaking conditions (180 rpm) for 24 h.

We prepared 3‑oxo‑C₁₂‑HSL and C₄‑HSL in methanol at 7 different concentrations: 0.1, 0.5, 1, 2, 3, 4 and 5 µM, each diluted from 100 mM stock. We centrifuged all cultures and washed each three times using PBS. We then re-suspended in LB and diluted to an OD (600) of 0.05. We then transferred 200 µl of each culture to a black 96-well plate with a clear bottom and inoculated with signals at the indicated concentrations. We repeated each experiment to generate five replicates. Methanol with no signal was used as a control. The plates were incubated in BioSpa at 37°C for 18 h. Measurements of OD (600) and RLU (Relative Luminescence Units) were collected every hour.

Data analysis

We estimated parameter values in Tables C, D, E, and F in S1 Text with non-linear regression by least squares using the Gauss–Newton algorithm [54]. Observations were limited to time ranges with peak expression values. (See Figs A, B, and C in S1 Text for detailed time course analysis.) Comparisons of model predictions and observed values are shown in Figs E, F, and G in S1 Text. Equilibrium values shown in Fig 7 were computed using a Trust-Region-Dogleg Algorithm [55]. Analyses performed and data visualizations created with MATLAB: Version 9.13 (R2022b) from The MathWorks, Natick, MA and Stata Statistical Software: Release 17 from StataCorp LLC, College Station, TX. All original code is available on GitHub at https://github.com/GaTechBrownLab.

Additional third-party modules:

1. Jann B. “PALETTES: Stata module to provide color palettes, symbol palettes, and line pattern palettes,” Statistical Software Components S458444, Boston College Department of Economics, 2017, revised 27 May 2020.
2. Jann B. “COLRSPACE: Stata module providing a class-based color management system in Mata,” Statistical Software Components S458597, Boston College Department of Economics, 2019, revised 06 Jun 2020.
3. Jann B. “HEATPLOT: Stata module to create heat plots and hexagon plots,” Statistical Software Components S458598, Boston College Department of Economics, 2019, revised 13 Oct 2020.

Custom color schemes adapted from seaboarn:

1. Watson ML. “seaborn: statistical data visualization.” Journal of Open Source Software 2021 Apr 06;6(60): https://doi.org/10.21105/joss.03021.

S1 Text. Supporting information.

Table A. Activation of QS genes by LasR/3‑oxo‑C₁₂‑HSL in review of published literature. Solid dots indicate positive activation in the paper’s diagram of gene transcription, while hollow dots indicate that the diagram shows no effect. No diagrams indicated repression. Note that some papers made no attempt to indicate particular interactions; several, for example, concentrated strictly on the QS genes themselves and did not show the effect on downstream genes such as those for elastase. Table B. Activation of QS genes by RhlR/C₄‑HSL in review of published literature. Same notation as previous table. Table C. Single signal parameter estimates. Estimated fold-change, derived from raw parameters of Eq A as (ɑ + ɑ₀)/ɑ₀, and half-concentration, K, values for gene expression as a function of a single signal in isolation. Values shown with 95% confidence intervals. Table D. Multi-signal parameter estimates. Model parameters for gene expression as a function of multiple signal concentrations. Parameter definitions are the same as in Table C with addition of cooperative fold-change, again derived from raw parameters as (ɑ + ɑ₀)/ɑ₀, and cooperative half-concentration K_Q. Values shown with 95% confidence intervals. Table E. Multi-signal parameter estimates for lasB. Model parameters for lasB expression as a function of multiple signal concentrations. Parameter definitions are the same as in Table D. Values shown with 95% confidence intervals. Half-concentration estimates less than 0.001 μM are below the limits of precision of the experimental data. Table F. Estimated proportionality constants that relate synthase expression levels to per-capita signal production rates. Final column shows adjusted R² of non-linear least squares estimate. Table G. Hierarchical and independent architectures are special cases of the reciprocal architecture. The multi-signal model of Eq 1 (main text) can represent hypothetical, alternative QS architectures by setting appropriate ɑ values to zero. Zero ɑ values result in a corresponding maximum fold-change of 1. For a hierarchical architecture, this setting nullifies the effect of C₄‑HSL on lasI. For an independent architecture, this setting additionally nullifies the effect of 3‑oxo‑C₁₂‑HSL on rhlI. Table H. Models of hierarchical and independent architectures can be normalized to ensure that maximum synthase expression is the same for all architectures. Parameters are the same as those in Table G but with increased values where appropriate. Table I. Model parameters for hypothetical three-signal architectures. Different parameter values result in the different responses of the third QS system’s synthase expression level as population density increases. Fig A. Expression level of lasI over time course of experiment. Shaded regions highlight peak expression and indicate two-hour period used in analysis. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig B. Expression level of rhlI over time course of experiment. Shaded regions highlight peak expression and indicate two-hour period used in analysis. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig C. Expression level of lasB over time course of experiment. Shaded regions highlight peak expression and indicate two-hour period used in analysis. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig D. Effect of each signal in isolation on the expression level of lasI Plotted points are observations and dashed lines show model (Eq A) predictions when parameterized per Table C. Dark blue points are additional observations collected using low signal concentrations. (A single data point identified as a faulty outlier is indicated in light blue and excluded from the analysis.) These data points are not used in estimating model parameters yet still show strong agreement with the model. Coefficient of determination R² between additional observations and initial model predictions is 0.82. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig E. Multi-signal model for lasI expression. Panels compare model predictions to observations for all combinations of signal concentrations. Horizontal bars indicate model predictions, while plotted points show observed values. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig F. Multi-signal model for rhlI expression. Panels compare model predictions to observations for all combinations of signal concentrations. Horizontal bars indicate model predictions, while plotted points show observed values. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig G. Multi-signal model for lasB expression. Panels compare model predictions to observations for all combinations of signal concentrations. Horizontal bars indicate model predictions, while plotted points show observed values. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig H. Equilibrium signal concentration predicted using proportionality constants. Individual data points show experimental observations and dashed lines indicate model predictions. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig I. Time response of lasB expression for reciprocal, rescaled hierarchical, and rescaled independent architectures. Dynamics are those of Eq 2 (main text) with parameters from Table H. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig J. Extracellular signal concentration as a function of density and mass transfer varies based on the quorum-sensing architecture. Heat maps of equilibrium 3-oxo-C₁₂-HSL (A–C) and C₄-HSL (D–F) concentration for three quorum-sensing architectures. Both population density and mass transfer rate are varied over the same ranges for all heatmaps. The lines on each heat map indicate density and mass transfer values for which equilibrium concentration is constant, either 50% of its maximum value (white) or 5% of its maximum value (black). Equilibrium concentrations calculated from Eq 2 model with parameters from Table D and architectural parameters normalized according to Table H. These results follow from the same model parameterization presented in Figs 7F–7H, which showcased the predicted outcome behavior of lasB expression. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig K. Ratio of signal concentrations as a function of density and mass transfer varies based on the quorum-sensing architecture. The figure shows heat maps of the ratio of equilibrium 3-oxo-C₁₂-HSL to C₄-HSL concentration for the reciprocal architecture. Equilibrium concentrations calculated from Eq 2 model with parameters from Table D. (The data underlying this figure and the code used to analyze it can be found in https://doi.org/10.5281/zenodo.15808353.) Fig L. Interaction strength for both induction and repression determines population behavior. The figure considers hypothetical architectures for a quorum-sensing network with three QS systems. The first two, mimicking the architecture of las and rhl, are mutually reinforcing. The third system is both induced and repressed by the other two, matching the reported interactions of pqs. The panels show all three synthase expression levels as a function of population density. As the four panels show, even within the constraints of a particular architecture, a wide variety of responses are possible. (A) Baseline case with weak activation of system 3 by system 1 (low α_3,1). (B) Strong activation (high α_3,1). (C) Limited activation with weak repression of system 3 by system 2 (moderately negative α_3,2). (D) Damped activation (strongly negative α_3,2). (The data underlying this figure and the code used to generate it can be found in https://doi.org/10.5281/zenodo.15808353.) Equation A. Expression level as a function of a single signal’s concentration.

https://doi.org/10.1371/journal.pbio.3003316.s001

(DOCX)