Parameters, architecture and emergent properties of the Pseudomonas aeruginosa LasI/LasR quorum-sensing circuit

2.1. Bacterial strains, plasmids and growth conditions

We used the P. aeruginosa PAO1 wild-type strain (WT, Iglewski lineage) [38], and an isogenic mutant containing a markerless in-frame deletion in lasI [18]. We used the promoter-probe vector pPROBE.AT [39] with lasI’-gfp, lasR’-gfp and rsaL’-gfp transcriptional fusions, respectively. The encoded green fluorescent protein (GFP) variant is the stable and fast-folding gfpmut1 [40]. The lasI’-gfp and rsaL’-gfp fusion constructs were described previously [18], whereas the lasR'-gfp fusion was constructed for the present study, using routine molecular cloning techniques. A 333 bp lasR promoter region was PCR-amplified from the PAO1 genome using forward primer 5′-N₆AAGCTTGCGATGGGCCGACAGTGAAC-3′ and reverse primer 5′-N₆GAATTCTCTTAAACTATTAACCAATCAGCCAAATATGGATT-3′. The PCR product was digested with HindIII and EcoRI (restriction sites underlined) and ligated with equally digested pProbe-AT. The resulting construct was confirmed by sequencing. Plasmids were introduced into chemically competent P. aeruginosa strains according to standard procedures [41]. Strains were grown at 37°C with shaking in Lennox Lysogeny Broth (LB) liquid cultures buffered with 50 mM 3-(N-morpholino)-propanesulfonic acid (MOPS), pH 7.0, or on LB agar plates. The antibiotic carbenicillin was added at a final concentration of 200 µg ml⁻¹ for plasmid selection and maintenance. Growth was measured as optical density at 600 nm (OD₆₀₀) in a photometer (Eppendorf Biophotometer), and as colony-forming units (CFU) per ml by dilution-plating. The conversion factor between OD₆₀₀ and CFU ml⁻¹ was determined by linear regression (electronic supplementary material, figure S1).

2.2. Gene expression and AHL measurements

We conducted gene expression and AHL quantitation of P. aeruginosa strains in LB-MOPS liquid cultures. We employed a pre-culturing scheme as described previously to dilute pre-existing GFP in cells to background levels [18]. In short, we inoculated pre-cultures with freshly streaked P. aeruginosa colonies, and grew them to an OD₆₀₀ of 0.1 to 0.2. Using these pre-cultures, we inoculated experimental cultures to an initial OD₆₀₀ of 0.001. Cultures were grown for 3 h before the first sampling. Subsequent samples were taken every 0.5 to 1 h, for up to 12h of total cultivation time. For gene expression analysis, P. aeruginosa WT and lasI mutant strains contained lasI’-gfp, lasR’-gfp or rsaL’-gfp reporter fusions. At each time point, aliquots were transferred to a black-walled microtiter plate to quantify GFP intensity in a fluorescence plate reader (Tecan M200) at 480 nm excitation and 535 nm emission wavelengths. The fluorescence of LB-MOPS medium was subtracted as a blank. The net fluorescence values were divided by the OD₆₀₀ at the respective time point to attain relative fluorescence intensities.

For AHL analysis, at each time point a 5 ml culture aliquot was extracted twice with acidified ethyl acetate as described in [42,43]. Extracts were concentrated to dryness by evaporation and resuspended in fresh ethyl acetate. The 3OC12-HSL content in these extracts was determined by an established Escherichia coli bioassay containing lasR and a lasB’-lacZ reporter [42,43]. ß-galactosidase activity of bioassay cultures was measured using the Galacto Light Plus Reporter Gene Assay (Thermo Fisher). 3OC12-HSL concentrations in the experimental samples were determined by comparison with a standard curve generated with synthetic 3OC12-HSL (RTI International).

2.3. Mathematical model

Mathematical modelling was based on a previously described quantitative model of QS [28,44]. This model consists of a system of ordinary differential equations (ODEs) that describes the change in signal and protein components, as well as cell growth, over time. Parameters are described in table 1. We considered the cellular signalling dynamics with two different versions of this model: (i) a minimal model with positive feedback on lasI expression only, and (ii) an extended model with positive feedback on both lasI and lasR expression, as well as negative feedback on lasI and rsaL expression via RsaL. The minimal model is described as follows:

extracellular 3OC12-HSL, S_e

intracellular 3OC12-HSL, S_i

and LasI synthase, I

The minimal model is virtually identical to that by Fujimoto and Sawai [28]. It makes several assumptions and simplifications. It models the positive feedback of 3OC12-HSL on lasI expression directly through the Hill function $S_i^m/(S_i^m\hspace{0.17em}+K_S^m)$ and does not explicitly model LasR. However, it considers positive cooperativity (expressed by the Hill coefficient m), which primarily manifests through receptor dimerization [49]. The binding of one signal molecule to a receptor monomer, which then drives dimerization, is assumed to be non-cooperative, as is the binding of LasR-3OC12-HSL to target promoters [53,54]. The model omits the concentrations of mRNA species, as they are assumed to be in quasi-steady state with the respective protein concentrations [27]. It considers extracellular but not intracellular signal degradation, because the contribution of intracellular degradation is negligible, given the extremely small intracellular volume fraction. The model also does not distinguish between contributions of diffusion and active transport to the movement of signal in and out of the cell [16], and does not explicitly consider potential effects of growth rate on the concentrations of intracellular species [55]. These simplifications reduce the complexity of the model such that parameter estimation is computationally feasible. The cellular volume fraction (${\rho }_v$) affects the extracellular 3OC12-HSL concentration (S_e) through population growth, as described below.

In addition to the minimal model, we built the following extended model that incorporates additional protein components and feedback regulation:

extracellular 3OC12-HSL, S_e

intracellular 3OC12-HSL, S_i

LasI synthase, I

LasR receptor, R

and RsaL repressor, L The full model includes LasR and RsaL, assuming separate basal and maximum synthesis rates (V_min and V_max) for each protein. Here, the Hill function models the dual positive feedback of signal–receptor complex formation on both signal synthase and receptor expression, following a previously proposed notation [28,44]. It considers positive cooperativity through receptor dimerization as above, and the same induction threshold (K_S) for all three proteins, which is consistent with our experimental data (see below). The repression of lasI and rsaL by RsaL is modelled by the term $K_L^n/(L^n\hspace{0.17em}+K_L^n)$, with separate repression constants (K_L1 and K_L2), respectively. This distinction in repression strength is again motivated by our experimental data and is consistent with the asymmetric promoter architecture (even though there is only a single RsaL binding site in the lasI-rsaL bidirectional promoter, it is positioned such that disproportionate repression is conceivable) [13]. A Hill coefficient (n) allows for cooperativity in the form of RsaL dimerization [50].

Further, the ODE for extracellular 3OC12-HSL is simplified in both models by assuming a quasi-steady state between intracellular and extracellular signal. Hence, S_e can be expressed through S_i

Each cellular model is connected to a growing population as follows.

population density, P

and volume fraction, ρ_v Population density P is modelled with a logistic growth equation. P is defined as cell number N (in colony-forming units, CFU) per total culture volume V_tot (normalized to 1 l). The volume fraction ${\rho }_v$ links population growth to the single-cell model through a change in S_e. Strictly speaking, the volume fraction ${\rho }_v$ is a volume ratio $N_{\hspace{0.17em}}(V_{\mathrm{cell}})/(V_{\mathrm{tot}}-V_{\mathrm{cell}})$, which can be readily shown by mass conservation law. However, because V_cell is very small compared with V_tot, V_tot ≅ V_tot – V_cell.

Finally, we explicitly model GFP intensities to relate our measurements of lasI’-gfp, lasR’-gfp and rsaL’-gfp reporter expression to the actual protein components in the cell (only lasI’-gfp in the reduced model, and all three reporters in the full model).

and GFP_rsaL, G_L The factor τ relates the rate of GFP synthesis to that of the respective protein. We assume here that promoter activity is proportional to protein concentration. We believe this assumption is appropriate given that QS is a transcriptional regulatory pathway, but of course it neglects possible post-transcriptional effects. The apparent degradation rate of GFP $({\gamma }_G)$ is the same in all equations.

2.4. Model analysis and parameter estimation

The model was built, simulated and analysed with the complex pathway simulator software package (COPASI 4.29) [56,57]. Parameters were estimated with the parameter estimation function (method Levenberg-Marquardt), using the experimental data in figure 1 and the guessed initial parameter values shown in table 1. Parameter fitting was done in iterations, starting with the initial parameter set and then again with the parameter set obtained after the initial step. The fitting algorithm calculates the following objective function, to which each experiment contributes with a weighted sum of squares [57]:

where P is the tested parameter set, x_i.j is a point in the dataset, and y_i.j is the corresponding simulated value. The indices i and j denote rows and columns in the dataset, corresponding to the measured time points and quantities, respectively. The weight for each data column is given by ${\omega }_j$. Weights were calculated using the mean squares method.

The Akaike information criterion (AIC) as an estimator of the prediction error was calculated as

where n is the number of data points in the dataset, k is the number of parameters fitted plus one and SSE is the least-squares error, i.e. the sum of the objective function in table 2 [58].

A steady-state analysis was conducted with the parameter scan function in Copasi. The GFP_lasI level at steady state was determined with cell density as the control variable. To maintain a constant cell density, the growth rate was set to zero, and time courses were run for 1000 h at 1000 different cell density values within a given range. Each cell density value was held constant throughout each time course. Two different initial LasI and GFP_lasI concentrations were set to separately investigate the off-to-on transition (I₀ = 0.01, G_{I 0} = 100) and the on-to-off transition (I₀ = 1, G_I0 = 10 000).

2.5. Feedback analysis

To investigate the significance of feedback in the model, we employed the parametrized extended QS model. In addition to S_i, we evaluated effects on the actual circuitry components I, R and L directly, rather than on G_I, G_R and G_L. To essentially eliminate each one of the three feedback loops individually, we set the equilibrium constants in the respective Hill and repression functions to very large values (1 × 10⁶), namely: K_s in the equation for I to eliminate positive feedback on I, K_s in the equation for R to eliminate positive feedback on R, and K_L1 and K_L2 in the equations for I and L to eliminate negative feedback on I and L via RsaL. Further, we added a hypothetical QS-controlled protein T to the model to evaluate induction kinetics and steady-state behaviour independent of individual feedback manipulations that alter the kinetics of I, R and L. It follows I induction kinetics but does not include negative influence from RsaL, and the value for K_s remains at 0.11 µmol² l⁻² (table 1).

3.1. Measurement of gene expression and signal production

As a first step, we collected the experimental data required to describe the las QS circuit. We focused on species relevant to our mathematical model, and on experimental methods that are feasible and quantitative. We measured the expression of lasI, lasR and rsaL in P. aeruginosa via GFP reporter fusions, and we measured the concentration of extracellular 3OC12-HSL by bioassay (figure 2). We used a standard growth medium, buffered Lysogeny broth, that is most commonly used in the field and allows comparison with other studies. We collected time courses across the full range of QS induction dynamics, from basal expression in the off state to maximal expression in the on state.

• Download figure
• Open in new tab
• Download PowerPoint

Consistent with previous studies, we observed a sharp cell density-dependent increase in GFP_lasI and GFP_rsaL from initially very low, basal levels in the WT strain. The expression of GFP_lasR increased concurrently, although basal expression levels were higher. For comparison, the expression of all GFP reporter fusions was much lower in the lasI deletion mutant, confirming that their induction is QS-dependent. By contrast to GFP_lasI and GFP_rsaL, GFP_lasR showed a modest, partially QS-independent increase throughout growth, the implications of which are discussed below. As expected, the observed concentrations of 3OC12-HSL during culture growth of the WT correlated well with the expression levels of GFP_lasI.

3.2. Mathematical modelling and parameter estimation

Next, we sought to fit our experimental data to a mechanistic mathematical model of QS. To establish a baseline, we first considered a minimal QS model with single feedback, introduced by Fujimoto & Sawai [28]. This ODE model includes four species, cell density, LasI concentration, extracellular 3OC12-HSL concentration and intracellular 3OC12-HSL concentration. It models positive feedback on lasI directly via changes in signal concentration (figure 3 top left; see Methods), with an underlying quasi-steady-state assumption. We amended the model with a GFP equation that directly relates gene expression measured by lasI’-gfp fusion to the concentration of LasI in the cell. Hence, three model variables, cell density, extracellular 3OC12-HSL and GFP_lasI, could be directly fit to experimental data. We asked whether this model is sufficiently complex to describe the experimentally observed induction dynamics. We used the modelling software Copasi [56] to implement the model and estimate parameters using a least-squares method. To guide parameter estimation, start values were taken from our own experimental data when available, or taken from the literature (table 1). In addition to the cell growth data, the minimal model was able to fit the timing of induction well, but failed to accurately capture the subsequent decrease in both 3OC12-HSL and GFP_lasI levels (figure 3 and table 2). We reasoned that the negative feedback on lasI by RsaL could largely account for the discrepancy.

• Download figure
• Open in new tab
• Download PowerPoint

As a second step, we conducted a parameter estimation with the extended QS model, which incorporates repression of both lasI and rsaL by RsaL (figure 3 top right). The repression of both genes works by the binding of an RsaL dimer to a conserved sequence element in the lasI-rsaL bidirectional promoter [13,14,50]. The extended model also incorporates LasR and its regulation by a second, positive feedback loop, which is supported by our experimental data. We performed a parameter estimation with the extended model as described above and obtained a greatly improved fit that captured the decline in 3OC12-HSL and GFP_lasI levels late in growth. Compared with the minimal model, the extended model showed decreased RMSE values for 3OC12-HSL and GFP_lasI (table 2). The RMSE for the newly fit GFP_rsaL was equally low, while that for GFP_lasR was somewhat higher, possibly indicating additional layers of regulation not captured with our model (see Discussion). Many estimated parameters were remarkably close to the start values used to seed the model (table 1).

To estimate the relative quality of each model, we determined the Akaike information criterion (AIC), which weighs the goodness of fit against the number of parameters used to achieve that fit. The AIC is −116 for the minimal model, and −181 for the extended model (electronic supplementary material, table S1). Because a smaller number indicates higher model quality, we conclude that the extended model performs better, even when considering the larger parameter space.

3.3. Bistability and hysteresis

Having parametrized the extended model, we sought to investigate two important and related emergent properties of the QS network. We asked whether the extended model, as fit to experimental data, can function as a bistable switch with two stable states (on and off), and whether this switch exhibits hysteresis. To answer these questions, we modelled signal synthase expression, GFP_lasI, at cellular steady state, as a function of different cell densities (figure 4). In support of bistability and hysteresis, we observed two distinct on and off states without a graded transition in between, and we observed sharp transitions at different cell densities, depending on the direction of the change. The cell density range between the low → high and the high → low transition is sizeable, suggesting that hysteresis could significantly improve the stability of the QS response in a fluctuating environment.

• Download figure
• Open in new tab
• Download PowerPoint

Of note, the state switching modelled here is so-called group-level bistability triggered by the accumulation of extracellular signal in a growing population of cells [28]. A concerted group-level response is supported by recent single-cell gene expression data showing unimodal distribution of GFP_lasI under culture conditions identical to those used here [17,18].

3.4. The role of positive and negative feedback

Next, we employed the parametrized extended model to analyse the significance of feedback architecture on the QS response. There are two positive feedback loops (LasR-3OC12-HSL activating lasI and lasR, respectively), and a negative feedback loop (RsaL activated by LasR-3OC12-HSL repressing lasI and rsaL). We removed one of the respective feedback loops in silico, as explained in Methods, and we simulated time courses for the concentrations of the actual QS circuitry components S_e, I, R and L in the modified models (figure 5). To also evaluate these effects on a QS-controlled component that is not itself involved in positive or negative feedback, we added a separate, hypothetical target protein T to the model, and we simulated its behaviour under dynamic and steady-state conditions. In the P. aeruginosa LasR regulon, such a target protein candidate would be, for example, the Type VI secretion effector PAAR4, the alkaline protease AprA, or the nucleoside hydrolase Nuh [7,14,59,60].

• Download figure
• Open in new tab
• Download PowerPoint

The omission of positive feedback on I had a large effect. Compared with WT conditions, it delayed the timing of induction, greatly decreased the concentrations of I and S_e and more modestly those of the other components, and resulted in rheostat-like rather than bistable behaviour. Of note, basal expression of I still resulted in the gradual accumulation of S_e to concentrations near the K_s (0.1 µM), such that R, L and T were induced to substantial levels.

By contrast, the omission of positive feedback on R had a much smaller effect. It did not change the magnitude or timing of induction (except for R, of course) and still produced bistable behaviour, albeit with a narrower margin for hysteresis. The observation that positive feedback on I is essential for bistability is consistent with previous theoretical work [30], as is the observation that positive feedback on R enhances hysteresis [26].

The omission of negative feedback on I and L via L had the expected effect on the magnitude of I, L and S_e, but it did not affect the timing of induction, nor bistability in general. The increase in the levels of I, L and S_e is in agreement with experimental data on lasI, rsaL and 3OC12-HSL expression in a P. aeruginosa rsaL mutant [14].