Pulse Generation in the Quorum Machinery of Pseudomonas aeruginosa

Pseudomonas aeruginosa is a Gram-negative bacterium that is responsible for a wide range of infections in humans. Colonies employ quorum sensing (QS) to coordinate gene expression, including for virulence factors, swarming motility and complex social traits. The QS signalling system of P. aeruginosa is known to involve multiple control components, notably the las, rhl and pqs systems. In this paper, we examine the las system and, in particular, the repressive interaction of rsaL, an embedded small regulative protein, employing recent biochemical information to aid model construction. Using analytic methods, we show how this feature can give rise to excitable pulse generation in this subsystem with important downstream consequences for rhamnolipid production. We adopt a symmetric competitive inhibition to capture the binding in the lasI–rsaL intergenic region and show our results are not dependent on the exact choice of this functional form. Furthermore, we examine the coupling of lasR to the rhl system, the impact of the predicted capacity for pulse generation and the biophysical consequences of this behaviour. We hypothesize that the interaction between the las and rhl systems may provide a quorum memory to enable cells to trigger rhamnolipid production only when they are at the edge of an established aggregation.


Introduction
Pseudomonas aeruginosa is a common Gram-negative bacterium responsible for a wide range of infections, including those of the urinary and gastrointestinal tract, the skin, and, most prominently, the respiratory system in immunocompromised hosts and sufferers of cystic fibrosis (CF). P. aeruginosa is a well-studied opportunistic pathogen in many contexts; it is well known for its ability to form biofilms (O'Loughlin et al. 2013;Singh et al. 2015), its swarming behaviour (Daniels et al. 2004;Shrout et al. 2006), its rapid acquisition of resistance to antibiotics (Shih and Huang 2002) and its quorum sensing (QS) behaviour (Fuqua et al. 2001). QS in P. aeruginosa is of particular interest because the mechanism is more complex than the originally discovered, prototypical Lux homolog positive-feedback loop (e.g. James et al. 2000;Shadel and Baldwin 1991) and the number of genes regulated by QS is large (Sitnikov et al. 1995), especially those associated with virulence (O'Loughlin et al. 2013). Mathematical models of QS in Pseudomonas aeruginosa have received a lot of attention. They provide the formalism to summarize current understanding as well as the means to explore mechanisms and evaluate emergent solution behaviour. Here, we develop a model description, employing recent genomic information and bioinformatic techniques, and explore mechanisms for the generation of pulses and memory effects for downstream rhamnolipid production.
In P. aeruginosa quorum sensing is governed by a hierarchical Luxl/LuxR system, which consists of two homolog pairs: LasI/LasR and RhlI/RhlR (Miller and Bassler 2001). Under this process, formation of the HSL autoinducers N-(3-oxododecanoyl)-HSL and N-(butyryl)-HSL is synthesised by LasI and RhlI, respectively (see Fig. 1). It should be noted, however, that signalling systems of las and rhl are specific in their activation of autoinducers, i.e. N-(3-oxododecanoyl)-HSL is unable to activate RhlR and, similarly, LasR cannot be activated by N-(butyryl)-HSL (Latifi et al. 1995;Pearson et al. 1997). Although biochemically independent, the las system is able to exert control of the rhl system through the transcriptional promotion of the RhlR gene by LasR/N-(3-oxododecanoyl)-HSL Latifi et al. 1996). As well as gene regulation effects, the rhl system has an important function of modulating rhamnolipid production via rhlAB. Rhamnolipids are particularly important in swarming motility where they are postulated to lower surface tension and allow expansion of the colony through their surfactant and wetting properties, driving the bacteria to swarm on surfaces (Glick et al. 2010;Kohler et al. 2000). In addition, a quinolone system (Dubern and Diggle 2008) may also modulate these interconnecting feedback loops (for simplicity, we do not model this aspect of QS here).
The first models of QS in P. aeruginosa were of the Lux (James et al. 2000) and the Las systems (Dockery and Keener 2001). Both descriptions, and subsequent models, highlight the existence of a fold bifurcation structure for the concentration of the response regulator in response to bulk cell concentration. The seminal paper by Dockery and Keener (2001) provides the foundation for the emergence of QS based on formal mass action arguments. However, there have been significant increases in biochemical knowledge of this system in the last 15 years. Fagerlind et al. (2003) constructed a large mass action model of the coupled Las and Rhl systems, including the effects of both RsaL and Vfr. This work also emphasizes the existence of  Fig. 1 The quorum sensing signalling system in Pseudomonas aeruginosa is composed of las and rhl systems. Arrows and barred arrows indicate activating (positive) and inhibiting (negative) regulatory interactions, respectively. Shapes on the diagram depict autoregulation terminology. Letters associated with each arrow reflect the associated time scale (ms = millisecond, s = second, and min = minute). Symbols associated with each shape are detailed in Table 1. Adopted from Van Delden and Iglewski (1998) (Color figure online) the classic twofold bifurcation diagram for the activation level of the system (typically for the levels of liganded LasR) with respect to the external concentration of HSL. Subsequent work Fagerlind et al. (2005) then explored how anti-virulence drugs  Skindersoe et al. 2008) are able to quorum quench this system. The qualitative model of Viretta and Fussenegger (2004) does not include the effect of the rsaL negative loop and lacks the capacity to deal with the nonlinear effects predicted here. The production of rhamnolipid was modelled and tested empirically by Chen et al. (2004). The exhaustive rule-based approach of Schaadt et al. (2013) includes the effect of rsaL but does not include the kinetic possibilities that emerge from nonlinear interactions. There are many other studies on this system but, increasingly, these are based on the perspective of computational rather than mathematical modelling (e.g. Dockery and Keener 2001;Fagerlind et al. 2005;Schaadt et al. 2013), typically with a non-mechanistic emphasis.

S L
In this article, we focus in detail on the Las system and its internal regulation and modulation for individual cell. The las system is composed of LasI, autoinducer N-(3oxododecanoyl)-HSL, LasR and RsaL . Importantly, biochemical evidence has now firmly established that LasR exists as a dimer in solution, with each monomer liganded by a single HSL (N-(3-oxododecanoyl)-HSL) molecule with additional evidence supporting higher multimers upon DNA binding (Schuster and Greenberg 2006). This is in contrast to assumptions of Dockery and Keener (2001) who assume a simpler mechanism for binding. Mathematically, the biochemical evidence is consistent with a Hill number of at least two and possibly much higher. In contrast, the RsaL transcriptional repressor is a helix-turn-helix protein that binds the promoter of lasI (De Kievit et al. 2002) and exists as a monomer in the cell (Rampioni et al. 2007a), leading to a Hill number of one. The transcription of both genes is promoted and regulated via binding of the two proteins to the same intergenic region between the lasI and rsaL operons, so except for rates the functional form for the transcription is likely to be identical (but we discuss variations of this in "Appendix 1"). This improvement in biochemical knowledge offers clear guidance for expected mathematical forms in the equations.
We begin with an investigation of the behaviour of gene regulation by constructing a mathematical model of a single cell. In order to focus on a biologically plausible region of parameter space, we conduct an extensive search of published data, noting deviations from modelling choices in the literature. The dynamical system exhibits a range of interesting solution behaviour. We find that bistable solutions and oscillations are possible and explore the bifurcation structure. Of significant interest is the potential for excitability in the system, in the presence of both single and multiple steady states, such that a modest perturbation from a low-level steady state triggers a large-amplitude excursion around phase space that eventually returns the system to low steady state. The extracellular concentration of HSL can have a significant impact on the dynamics. As signal molecules can accumulate and diffuse in the external environment, we demonstrate how pulse generation of the las system can be propagated in space.
This paper is organized as follows. In Sect. 2, we describe in greater detail the biological system as well as the mathematical approach. In Sect. 3, we conduct a numerical exploration of our model, highlighting the key parameters in determining LasI and RsaL equilibria. We then employ extracellular concentration as a driving factor in the dynamical behaviour of the system. We conclude by discussing our findings, identify challenges and suggest future work in Sect. 4.

Methods
We construct the governing equations using mass action kinetics, except where noted, guided by the literature (e.g. Fagerlind et al. 2003). First, consider the LasR regulator (R L ), its binding activator, the autoinducer 3O-C12-HSL (H L ), and the complex LasR/3O-C12-HSL (R L H ). If LasR associates at rate k + L , dissociates at rate k − L , is produced at a rate β 0 and degrades or dilutes at rate γ L , then we may write In a similar manner, we formulate an equation for the complex, with degradation or dilution rate γ RL , such that where β 0 , γ L , γ RL are positive constants. New biochemical data have indicated approximately steady cellular levels (Ishihama et al. 2014) of many transcription factors and, therefore, we hypothesize it is the change in activation rather than production per se that dictates the dominant dynamics; unlike previous studies that consider the regulated production of LasR. Therefore, we shall assume the association and dissociation processes occur over a sufficiently fast time scale that the production and loss terms can be safely neglected: we disregard terms multiplied by β 0 , γ L and γ RL .
The autoinducer 3O-C12-HSL (H L ) is created in the system via the activity of the LasI synthase (I L ), which we take to be at rate β H L , and is naturally lost from the system at rate γ H L . The most significant loss of the autoinducer from the cell is via diffusion through the cell membrane, a process we account for separately. Taking a simplified description of diffusion we can express the diffusive term as being proportional to the concentration difference across the membrane of H L , where we take the extracellular concentration to be H ex . Therefore, D H L represents an additional loss rate, which is multiplied by the concentration difference, The enzyme LasI (I L ) is produced by the lasI gene through a transcription and translation process of lasI mRNA (Î L ) at rate α L and degrades at rate γ I L , such that In a similar fashion, the inhibitor RsaL (S L ) is produced by rsaL genes through transcription and translation of rsaL mRNA (Ŝ L ) at rate α S and degrades at rate γ S , providing Transcription at the lasI promoter site (Î L ) is activated by the LasR/3O-C12-HSL complex (R L H ). The production process is assumed to follow a Hill form with a Hill number p. Recent biochemical evidence strongly suggests a Hill number > 1 in contrast to the arguments of Dockery and Keener (2001); the activated form of lasR is at least dimeric (Schuster and Greenberg 2006), and it is possible that it forms a tetramer on the DNA. For mathematical simplicity, we adopt the smaller potential value, p = 2, in the analysis that follows. In addition, the RsaL transcriptional repressor has been shown to bind to the lasI-rsaL intergenic region but, unlike LasR, RsaL is a helix-turnhelix protein that exists as a monomer in the cell (Rampioni et al. 2007a), and thus we will adopt Hill number q = 1 for the inhibition factor for production of RsaL (Ŝ L ), mathematically equivalent to a Michaelis-Menten equation. The transcription of both genes is promoted and regulated via binding of the two proteins to the same intergenic region between the lasI and rsaL operons, so the functional form for the transcription is identical, with the exception of the numerical values of the transcription and loss rates. Note that this implies a negative feedback relation between the RsaL protein and its own production, which has not been historically represented on graphical depictions of the las system, but strongly implied by analysis of Rampioni et al. (2007a). Biochemically there is insufficient evidence to determine whether the binding in the intergenic region results in competitive, uncompetitive or non-competitive binding system, or indeed whether there is a symmetry in the expression rates in each direction with all configurations of binding at the intergenic region. In the analysis below, we adopt a symmetric competitive binding form (we choose K L = K S in the analysis below), but show in "Appendix 1" that qualitatively the results are not dependent on this functional choice. With basal expression of β L0 and a loss rate of γ mL , this leads to the following expression for lasI: Similarly, with basal expression of β S0 and a loss rate of γ S for rsaL we obtain The numerical values of the temporal processes are such that we are able to make a number of assumptions regarding the timescales to simplify the system. We assume that the dominant, slowest processes are protein production from mRNA via translation and folding. Therefore, other processes, namely the liganding of regulators, DNA binding, synthetase operation, and the transcription of DNA, are much faster and we can assume that four of our differential equations are at a quasi-steady state, such that From Eqs. (1) and (2), we have a Moiety conservation equation Initially, we make the simplifying assumption that HSL diffusion is rapid and, therefore, that the equation for HSL, H L can also be written in a quasi-steady state, providing Later, we relax this constraint. With the simplifications above, the system of equations for the las system, constructed for the LasI and RsaL loops, becomes just two differential equations. By assuming there is negligible basal production of lasI and rsal genes (β L0 = 0 and β S0 = 0), the governing equations become and with H L as in Eq. (11). We nondimensionalize this model by writing so that (12) and (13) become where Here, a 1 , a 2 , a 3 , a 4 , and a 5 are positive constants. The biological interpretation of this model is that ξ inhibits the expression of both η and ξ , which describes negative feedback from competitive inhibition by RsaL to the expression of both lasI and rsaL genes. The binding of RsaL to the bidirectional lasI-rsaL counters positive feedback, thereby balancing levels of HSL. In addition, both η and ξ degrade exponentially. The model thus constructed is a type of incoherent feed forward motif (Alon 2006) as can be seen from Fig. 1.

Results
The non-dimensional set of differential Eqs. (15) and (16) have been investigated analytically with the assistance of Mathematica (10; Wolfram) and solved numerically using MATLAB (R2016a; MathWorks). We note that the fixed points of this system lie upon the straight line a 3 a 4 η = a 5 ζ but substitution of this expression into either nullcline leads to a quintic equation for the fixed points that does not provide sufficient further simplification to yield general expressions for stability bounds. Table 2 lists parameters that have either been adopted from the literature based on experimental evidence or estimated, as stated. Moreover, some parameters are chosen for the following reasons. The basal production rate of genes-mRNA can be considered as similar to the basal transcription rate of a protein. This is because the transcriptional regulator protein activates genes-mRNA in a very fast process before encoding the protein. We take typical values of total concentration of LasR and LasR/3O-C12-HSL to be 200 nM, as for the concentration of QseB in E. coli (Ishihama et al. 2014). As described in the previous section, we simplify the governing equations by assuming quasi-steady states for the fast reactions. This results in just 5 nondimensional parameters (see Table 3).

Parameter Ranges in the System
The dimensionless variable a 1 is proportional to H ex , the extracellular HSL concentration, an important factor in controlling both the intracellular HSL production and cell-cell communication. We shall see that a 1 influences the location of the key bifurcation.
Parameter a 2 depends on K L and represents the relative binding strength of the LasR dimer, and a 4 assigns binding control of RsaL to LasI through parameter K SL ; consequently, it gives the relative production of RsaL to LasI. Both parameters a 2 and a 4 have a wide potential range, 25 × 10 −6 − 25 and 0 − 25 × 10 3 , respectively. This is because the affinity constant between transcriptional regulator protein to genes-mRNA is largely unmeasured and depends on the chain structure of the signal molecule (short or long chain; 1 − 1000 nM) (Alon 2006).
Parameter a 3 is inversely proportional to β H L , the HSL production rate. It is related to a 5 , which describes the degradation of RsaL relative to HSL production. It is important to note that we set a 5 < a 3 due to the transcription factor being more stable than the synthase. We shall see that the parameters a 3 and a 1 allow for fold and Hopf bifurcations in the system, leading to bistability and oscillations, respectively, as well as excitable dynamics, and provide suitable control parameters that can be varied in experiments.

Excitable Dynamics in the LasI and RsaL Phase Plane
The simplified system (15) and (16) involves just two variables, which we investigate with phase plane analysis. By varying a 3 , four phase portraits of interest can be identified, as depicted in Fig. 2a. We plot the two nullclinesη = 0 andξ = 0, where η and ξ represent the dimensionless of LasI and RsaL concentration, respectively. Fixed points for the system lie at their intersection.
A qualitative description of the model behaviour is now possible, with reference to classical excitable systems such as the Fitzhugh-Nagumo model (Murray 1993). Cases  Welch et al. (2000) α L Rate at which LasI produced by lasI mRNA 0.5 min −1 0.5 2 min to translate protein, Alon (2006) Fig. 2a), which correspond to small and large values of the bifurcation parameter a 3 = 0.1 and 0.6, respectively, possess a single stable fixed point with the nullclines positioned so that there is no possibility of excitable behaviour (see below). For case II (Fig. 2a), there are two further intersections, which yield the central zone of the well-known S-shaped bifurcation diagram with three fixed points in a stable-unstable-stable configuration. This phenomenon has been reported in many experiments on autoregulation of genes (e.g. Alves and Dilo 2005;Angeli et al. 2004;Poignard 2014); more detail is presented in Fig. 3. The feedback loop that consists of both positive and negative autoregulation creates two possible LasI production states, "on" and "off" (at the large and small stable steady states, respectively), affecting signal molecule concentration. This is a familiar pattern that represents quorum sensing; the cells have a low steady state ( not reversible and the system must trace a hysteretic loop to return to the low steady state via changes in the external parameters. In Case III, there is a single fixed point and the nullclineη = 0 has a positive local maximum (which distinguishes it from cases II and IV); see Fig. 2a. This case allows for excitability in the LasI and RsaL phase plane caused by theη nullcline having cubiclike shape with a positive local maximum (a 3 < 0.55). By considering the behaviour of a perturbation to the steady-state solution (see Fig. 2b), if the production of the (N-(3-oxododecanoyl)-HSL) molecule is sufficiently small, then global stability of the single fixed point ensures that the trajectory returns to the steady state. The variable ξ will engage, rapidly deplete and return the system to the small fixed point with a short loop. However, interesting behaviour occurs when there is a greater perturbation of (N-(3-oxododecanoyl)-HSL). For then the variable η pushes the system to the right in the phase plane, so that there is an excitable pulse of (N-(3-oxododecanoyl)-  Table 3 (Color figure online) HSL). Consequently, there is no short route back to the stable fixed point and the trajectory undergoes a large excursion before returning. A small region of parameter space in Case II can also result in excitability (0.27 ≤ a 3 ≤ 0.28). Here there are three intersections between η and ξ nullclines, which are low stable fixed point, an intermediate saddle and a high unstable fixed points. In this region of parameter space, the saddle-unstable fixed point region is not accessible to trajectories that initiate in a zone close to the stable fixed point (see "Appendix 2", particularly Fig. 8).
The initial condition η(0) = η 0 and concentration of extracellular HSL have a large effect on the amplitude and duration of the excited pulse, as can be seen in Fig. 2. The results in Fig. 2b are plotted as a function of time in Fig. 2c, d, where the purple dashed line and red dashed line represent either excitable or non-excitable trajectories for large and small perturbations, respectively (for a 1 = 0.01 or, equivalently, H ex = 10nM). By increasing external HSL concentration in the system via the variable a 1 (equivalent to increasing H ex from 10 nM to 300 nM), we find solutions tending to the single high steady state irrespective of the initial perturbation (shown in Fig. 2c, d by blue solid and black dotted lines; discussed later).

Bifurcation Structure
The QS circuity of P. aeruginosa is complex, with the las system itself consisting of both positive and negative feedback loops. The lasI gene is activated by LasR/3O-C12-HSL leading to increased (N-(3-oxododecanoyl)-HSL) production, positive feedback. However, LasR/3O-C12-HSL also activates the rsaL gene and induces expression of RsaL. RsaL inhibits expression of lasI and rsaL genes, and this process constitutes negative feedback. In general, this negative feedback loop is employed to maintain a homeostatic balance in the system. However, the interaction of both positive and negative feedback can yield bistability (and associated hysteresis; Pfeuty and Kaneko 2009). There are many other studies on the dynamical behaviour of gene regulatory networks involving positive and negative feedback loops. Song et al. (2007) demonstrated that interlocked positive and negative feedback loops play essential roles in bistability and oscillations. More complex bifurcations of co-dimension one or two have also been explored (Hat et al. 2016). Varying a 1 > 0, the effect of extracellular signal molecule alters the qualitative dynamical behaviour of LasI and RsaL. Thus, initially we investigated the dynamics of our system with a one-parameter bifurcation analysis of LasI versus a 1 . We used continuation methods to track the evolution of solutions for η versus a 1 . Figure 3a-e depicts five qualitatively different dynamical behaviours, each at a different fixed value of a 3 , the degradation of LasI relative to signal molecule production. As shown in Fig. 3a-d, there are at least four different bifurcation diagrams. Bifurcations of steady states do not appear in Fig. 3e. Figure 3a, for a 3 = 0.26, depicts the potential for irreversible bistability. The left-hand saddle-node bifurcation can cross the η-axis for sufficiently small a 3 , preventing a drop from the upper to the lower branch. In addition, in Fig. 3b-d, reversible bistability of LasI is evident in the range of a 1 bounded by the two limit points, SN 1 and SN 2 , which denote saddle-node bifurcations. In the low steady-state branch, the concentration of LasI is low until a 1 level exceeds the critical value SN 1 (a 1 = 0.029, 0.032, and 0.0329, for Fig. 3b-d, respectively), at which point the concentration of LasI increases abruptly to a high value. In similar manner, starting with a 1 high, the concentration of LasI does not drop significantly until a 1 reduces below a critical value, either SN 2 (a 1 = 0.023, 0.0318, and 0.0328, for Fig. 3b-d, respectively) or an earlier bifurcation such as a subcritical  (a 1 , a 3 ). The bifurcation lines divide the parameter domain into six regions R 1 ,…, R 6 . Each of these regions is explained in the main text. The blue lines depict the saddle-node lines, defined by the locus of saddle-node bifurcation points (or limit points), with subscripts 1 or 2 as in Fig. (3). The Hopf lines are constructed by subcritical Hopf and supercritical Hopf bifurcation points, which are presented by brown lines. The bistable region, R 3 , consists of reversible (above of grey-dash lines) and irreversible (below of grey-dash lines). The Hopf bifurcation. With a 3 = 0.3 (Fig. 3b), subcritical Hopf (sub-Hopf) bifurcation points exist for the high fixed point in the bistability region, at a 1 = 0.024. The upper steady-state solution is unstable between SN 2 and sub-Hopf point and stable beyond the sub-Hopf point where a 1 > 0.24. For a 3 = 0.32 (Fig. 3c), the two saddle-nodes (SN 1 and SN 2 ) move close together. Consequently, the bistability region is narrow.
Here, the Hopf bifurcation becomes supercritical (sup-Hopf; at a 1 = 0.032). For a 3 = 0.323 (Fig. 3d), the sup-Hopf point disappears and the saddle-node bifurcations collide in Fig. 3e for a 3 ≥ 0.325; the bistable regime ceases to exist. Figure 4 provides a two-dimensional bifurcation diagram in the (a 1 , a 3 )-plane. In general, for a value of a 1 smaller than at the cusp point, decreasing a 3 from a large value we move through a region with one relatively small stable fixed point for η to one small stable fixed point with excitable dynamics. Decreasing a 3 further provides a bistable regime before eventually reaching a region with one large stable fixed point for sufficiently small a 3 . The bistability region collapses through the collision of the two saddle-node lines at the cusp (a 1 = 0.033, a 3 = 0.324), a possibility that might be inferred from Fig. 2a. For a 1 > 0.033, the bistable regime ceases to exist.
Three co-dimension-2 singular points, a cusp point, a Bogdanov-Takens bifurcation point (BT), and a generalized-Hopf bifurcation point (GH), are identified in Fig. 4. The  Bogdanov-Takens bifurcation, at (a 1 , a 3 ) = (0.032, 0.321), is caused by coalescence of a saddle-node point (SN 2 ) and a Hopf bifurcation point (sup-Hopf). A sup-Hopf line extends to the lower right of the saddle-node lines close to the cusp in the bifurcation diagram. The generalized-Hopf bifurcation is at (a 1 , a 3 ) = (0.317, 0.002), a transition point from sup-Hopf to sub-Hopf. It is clear that the bifurcation structure is complex, especially close to the cusp point, although the solution space is dominated by the fold bifurcation and excitability.
The response diagram (Fig. 4) illustrates that the bifurcation lines divide the (a 1 , a 3 ) parameter plane into six regions, as depicted by region R 1 ,…, R 6 . Detailed phase diagrams associated with each regions are given in "Appendix 2". In summary, R 1 has either a monostable fixed point (one stable steady state) with low concentrations of LasI and RsaL for a 3 ≥ 0.55 (Fig. 8a), or excitable solutions, as discussed in the previous section (Fig. 8b). In the narrow region R 2 (see Fig. 4b for detail figure), bounded by SN 2 and sub-Hopf, three steady states arise, but only the lower steady state is stable (Fig. 8c, d). Excitable solutions are also possible. Region R 3 is divided into four domains with different behaviour. Three steady states arise in which the upper and lower steady states are stable, but the middle state is unstable. In the first domain, the upper state in R 3 is surrounded by an unstable limit cycle (Fig. 8e). A homoclinic bifurcation arises to give the second domain (Fig. 8f). Upon decreasing a 3 further in this region, a high stable spiral exists with a large domain of attraction (Fig. 8g, h). Region R 4 is monostable with high concentrations of LasI and RsaL (Fig. 8i). In the area approaching the cusp point R 5 and R 6 , two small distinct regions arise (see Fig. 4b).
Region R 5 provides stable oscillations and R 6 bistability. In summary, the bifurcation lines divide the (a 1 , a 3 )-plane into regions of distinct solution behaviour: monostability (R 1 and R 4 ), bistability or excitation (R 2 , R 3 , and R 6 ), and oscillatory R 5 .

Travelling Wave of a Pulse in a Linear Chain of Cells
We demonstrated the potential for the las system to act as a pulse generator for a single cell in the previous section. To illustrate how this effect may translate into cell-cell communication within a colony, we consider a simplified linear chain of cells. The goal is to provide proof of principle that it is possible to set up a pulse train when the individual cells are coupled to the dynamics of their neighbours. For simplicity, we assume that diffusion across each cell membrane is such that the intracellular concentration of HSL is at a kinetic equilibrium (i.e. dH L dt = 0, as before) and, additionally, assume that there is a neighbourhood surrounding each cell where the extracellular concentration can be modelled simply by where Δ is the isotropic loss of H ex from this loosely defined neighbourhood. Writing in a dimensionless form using Eqs. (11) and (14) we obtain where Additionally, we construct a chain of otherwise identical systems, (η i , ζ i ), for i = 1, 2, . . ., coupled via the external concentrations of HSL in the local neighbourhood of each cell.
where ζ is the rescaled H ex variable and d is now the fraction of the isotropic single loss that is retained in the single cell-cell interchange. In this simplified system, we can investigate whether it is possible for a single cell i = 1 to propagate a signal. To do this, we trigger an arbitrary gain in the HSL, sufficient to trigger an excitable response from system i = 1, and examine the downstream impact. To do this, we increase the levels of HSL by two orders of magnitude from the low steady state (QS off) for the equivalent lifetime of a single cell (corresponding to one dimensionless time unit) after the cells reach steady state. Three types of solution are presented in Fig. 5, illustrating a range of behaviour from no propagation of the signal for small coupling, propagation of a pulse for intermediate coupling, to propagation from a mutually supported high steady state (QS on) for large coupling.

Discussion
We have described a model of the QS system in P. aeruginosa by considering nonlinear positive and negative feedback loops associated with the production of the synthase LasI and the regulator RsaL (Fagerlind et al. 2005;Pearson et al. 1997;Rampioni et al. 2007b). These nonlinear effects create the possibility of novel dynamical behaviours in the model. We have created a dimensionless set of equations to describe these behaviours and explored how the five dimensionless parameters affect the results whilst maintaining biological plausibility. Where possible we have taken parameters from the biological literature and this has led to significant deviation in our parameter choices from existing mathematical biology manuscripts, notably the work of Fagerlind et al. (2005);Fagerlind (2008), which needs commenting upon. When the parameters in these articles, which subsequently have been adopted in other later texts, were compared to biological estimates (presented for example in Alon (2006)), they were found to be significantly different. For example, the parameters used in  (Fagerlind et al. 2005) suggest that the typical lifetime of a transcription factor is of the order of seconds, when biological estimates typically describe transcription factors as stable proteins with lifetimes of the same order as the cellular turnover time, i.e. hours for Pseudomonas aeruginosa. Furthermore, it suggests mRNA molecules are more stable than the proteins they are translated to, which is clearly at odds with biological knowledge about the bursty nature of protein production (Xie et al. 2008).
We are confident that our dimensionless parameters thus lie within a biologically plausible region, though we acknowledge that the effects we report here are sensitive to the values of these parameters. An additional complication is the possibility of other functional forms for the rate of transcription of mRNA from the lasI and rsaL genes as a function of the active forms of both LasR and RsaL. The form we present is one from a family of different choices that result in the same qualitative behaviour: our analysis, presented in the appendix, suggests that competitive binding between RsaL and LasR is a requirement for excitability but symmetrical effects, resulting in the negative action of RsaL on its own production, are not required (as is shown on many diagrams of the las system Dockery and Keener 2001). In the absence of specific biochemical data of this relationship at the intergenic region, our analysis represents a significant step forward in incorporating biochemical knowledge in mathematical models of quorum sensing. There are biological instances where the binding in the intergenic region decouples and the rsaL and lasI genes are promoted and repressed independently, notably in the genetic isolates obtained from cystic fibrosis patients (Rampioni et al. 2007b). In this case, the mathematical modelling of the system becomes trivial as the two systems are decoupled, and there is no possibility of excitation.
A survey of a realistic region of parameter space revealed a range of interesting numerical solution behaviour, such as limit cycles. The full equations, without approximations to reduce the number of equations, were also solved numerically and revealed similar solution behaviour (not shown). Continuation methods were employed for the reduced system in the a 1 -a 3 plane to track bifurcations of the system of co-dimension one and two. The parameters a 1 and a 3 represent the information outside and inside the cell, respectively. We found fold and Hopf bifurcations, both of co-dimension one. Furthermore, there are Bogdanov-Takens and generalized-Hopf bifurcations, each of co-dimension two. For example, at a Bogdanov-Takens bifurcation, Hopf and fold curves in the plane intersect. However, most of the complex bifurcation structure is confined to a relatively small region of parameter space; the dominant behaviour is that of the fold bifurcation outside which are regions of excitable solutions.
We have demonstrated the potential for the las system to act as a pulse generator. As there is no explicit feedback from rhl to las, we may consider the las system as a blackbox controller for the rhl system. We have demonstrated that this can lead to the generation of a pulse train when the individual cells are coupled by diffusion of the HSL signal molecule. These observations allow for a novel stigmergy (Gloag et al. 2015) in the coupled las/rhl system; that of a quorum memory. The las system is coupled to the rhl system in two distinct ways; the active LasR transcription factor promotes the transcription of its counterpart in the rhl system, RhlR, and due to competitive binding on the RhlR molecule of the two homoserine lactones the presence of 3-oxo-C12-HSL acts to prevent the activation of RhlR by C4-HSL. This creates an interesting unreported effect in the combined system: what one might call a 'handbraked acceleration'. Assuming the rhl system follows the same Lux-like dynamics of the las system, the additional production of RhlR increases the likelihood of the system switching to its higher steady state, but this is tempered, or handbraked, by the presence of significant amounts of 3oxo-C12-HSL preventing the activation of the rhlR system and its downstream consequences, notably rhamnolipid production. However, the effect of the pulses (the 'revs' of the 'acceleration') is still to increase the amount of RhlR, sustained due to its relative stability compared to the diffusing homoserine lactones, so that when the handbrake of 3-oxo-C12-HSL is removed the system has an increased likelihood of activating the rhl system. In effect the cells have been primed by the pulses and thus have a memory of experiencing higher densities of cells (or low diffusion regions Redfield 2002). Cells will lose memory only when they no longer experience a local environment rich in the 3-oxo-C12-HSL from the las system. Therefore, the competition between the two signal molecules, one produced by the las system and the other produced by the rhl system, may enable cells to trigger rhamnolipid production only when they are at the edge of an established aggregation. This suggests a previously unreported reason for the coupled nature of the las and rhl systems and a mechanism for how they act in tandem to create a sophisticated control system for sociality and virulence in this important pathogenic organism.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1
RsaL is thought to work as an inhibitor in the las system on P. aeruginosa and, therefore, the Michaelis-Menten equation can be applied. Whilst there are a considerable number of studies about the QS system for P. aeruginosa, to the best knowledge of the authors, none of these studies provide sufficient biochemical evidence to determine whether the binding of RsaL in the intergenic region of the las system is competitive, uncompetitive, non-competitive, or mixed-competitive inhibition (using the standard biochemical descriptors). Here, we restrict attention to competitive, uncompetitive and non-competitive inhibition (mixed-competitive inhibition represents a combination of both competitive and uncompetitive inhibition, such that the RsaL inhibitor binds LasR/3O-C12-HSL or LasR/3O-C12-HSL+lasI genes with different affinities; mixed-competitive is the same as non-competitive inhibition if genes are bound with equal affinity). Furthermore, there is insufficient biochemical evidence whether the expression rates for lasI and rsaL genes are symmetric or not (Fig. 6).

Competitive Inhibition
Competitive binding occurs when both the LasR/3O-C12-HSL enzyme and RsaL inhibitor compete for binding to the active sites of the lasI gene (see Fig. 6a). Thus, the governing equations result from a blockade of either positive or negative lasI gene interaction by either LasR/3O-C12-HSL or RsaL, respectively. Furthermore, this competitive inhibition can be classified into non-symmetric and symmetric binding based on expression rates for LasR/3O-C12-HSL to lasI and rsaL genes. Due to the overlapping of binding sites, LasI is completely inactive when the RsaL inhibitor binds lasI genes. For the non-symmetric case, the affinity constant between LasR/3O-C12-HSL and rsaL is relatively small. Meanwhile, in the symmetric case the affinity constant between LasR/3O-C12-HSL and rsaL gene is identical to the affinity constant between LasR/3O-C12-HSL and the lasI gene.

Symmetric Competitive Inhibition
In this case, we find that yielding the non-dimensional model

Uncompetitive Inhibition
Uncompetitive binding occurs when the RsaL inhibitor binds to a site which only becomes available after the LasR/3O-C12-HSL enzyme has bound to the active site of the lasI gene (see Fig. 6b). It can also be classified into non-symmetric and symmetric binding types, as indicated below.

Non-symmetric Uncompetitive Inhibition
Here, so that the non-dimensional model becomes where ξ is a rescaled ξ .
Here, K SL represents the dissociation constant for the RsaL inhibition of binding of lasI genes to the LasR/3O-C12-HSL enzyme. It is different to K SL , which represents the dissociation constant of inhibitor RsaL on lasI genes.

Non-competitive Inhibition
Non-competitive inhibition occurs when the RsaL inhibitor can bind to the free lasI gene or the lasI-bound LasR/3O-C12-HSL enzyme (see Fig. 6c). As explained above, the non-competitive inhibition form is a special case of mixed-competitive inhibition, where it binds to the target with two equal inhibition constants.
Using the simplification steps highlighted in the main text, we determine appropriate nullclines for each binding type (see Fig. 7). The equations resulting from non-symmetric and symmetric competitive behaviour are presented in Fig. 7a, b, revealing the same qualitative behaviour as for the system in the main text. This suggests that competitive binding between RsaL and LasR/3O-C12-HSL with or without symmetry is sufficient for excitable behaviour. Therefore, although (Rampioni et al. 2007b) suggest there may be negative feedback from RsaL to its a η (LasI) before returning to the stable fixed point (see Fig. 9a, b). Other solutions are also possible analogous to those in the reduced system. These results provide confidence that the behaviour of the full system has been adequately captured by the reduced model.