This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Bulletin of Mathematical Biology (2001) 63, 95–116
doi:10.1006/bulm.2000.0205
Available online at http://www.idealibrary.com on A Mathematical Model for Quorum Sensing in
Pseudomonas aeruginosa
JACK D. DOCKERY∗
Department of Mathematics,
Montana State University,
Bozeman,
MT 59718, U.S.A.
Email: umsfjdoc@math.montana.edu JAMES P. KEENER
Department of Mathematics,
University of Utah,
Salt Lake City,
UT 84112, U.S.A.
Email: keener@math.utah.edu The bacteria Pseudomonas aeruginosa use the size and density of their colonies to
regulate the production of a large variety of substances, including toxins. This phenomenon, called quorum sensing, apparently enables colonies to grow to sufﬁcient
size undetected by the immune system of the host organism.
In this paper, we present a mathematical model of quorum sensing in P. aeruginosa
that is based on the known biochemistry of regulation of the autoinducer that is crucial to this signalling mechanism. Using this model we show that quorum sensing
works because of a biochemical switch between two stable steady solutions, one
with low levels of autoinducer and one with high levels of autoinducer.
c 2001 Society for Mathematical Biology 1. I NTRODUCTION The bacterium Pseudomonas aeruginosa is an increasingly prevalent human
pathogen, responsible for 12% of hospitalacquired urinary tract infections, 10%
of blood stream infections, and 8% of surgical wound infections. It is also one of
the most common and lethal pathogens responsible for ventilatorassociated pneumonia in intubated patients, with directly attributable deaths reaching 38%. Cystic
ﬁbrosis patients are characteristically susceptible to chronic infection by P. aeruginosa, which is responsible for high rates of illness and death in this population
(Van Delden and Iglewski, 1998).
The ability of P. aeruginosa to invade a potential host relies largely on its ability
to control many of its virulence factors by a mechanism that monitors cell density
∗ Author to whom correspondence should be addressed. 00928240/01/010095 + 22 $35.00/0 c 2001 Society for Mathematical Biology 96 J. D. Dockery and J. P. Keener and allows communication between bacteria. That is, isolated production of extracellular virulence factors by a small number of bacteria would probably lead to an
efﬁcient host response neutralizing these factors. However, the typical behavior is
for virulence factors to be expressed by an entire colony only after the colony has
achieved a certain size or density. In this way, cells excrete extracellular virulence
factors only when they can be produced at high enough levels to overcome host
defenses.
Another feature of persistent P. aeruginosa infections in cystic ﬁbrosis patients is
that when cell densities are high enough, the population as a whole begins to produce copious amounts of exopolysaccharide alginate (Govan and Deretic, 1996).
It is thought that these bacteria grow inside a bioﬁlm [deﬁned as microcolonies
surrounded by an exopolysaccharide, see Characklis and Marshall (1990)] as a
survival strategy since the surrounding polymer matrix protects bacteria. P. aeruginosa growing in an alginate ‘slime matrix’ are resistant to antibiotics and disinfectants (Govan and Deretic, 1996).
Recently, a celltocell signaling system has been shown to be involved in the
differentiation of P. aeruginosa bioﬁlms (Davies et al., 1998). A mutant defective
in the production of a certain signaling molecule formed an abnormal bioﬁlm that
in contrast to the wild type bioﬁlm was sensitive to low concentrations of the detergent sodium dodecyl sulfate (SDS). Furthermore, the addition of the signaling
molecule to the batch culture medium restored production of a slimy bioﬁlm that
was SDS resistant (Davies et al., 1998).
The ability of cells to sense their own cell density, to communicate with each
other and to behave as a population instead of individually is called quorum sensing. Quorum sensing is known to occur in a growing number of examples, including sporulation and fruiting body formation by Myxococcus xanthus and antibiotic
production by several species of Streptomyces. Perhaps the best characterized example is in the autoinduction of luminescence in the symbiotic marine bacterium
Vibrio ﬁscheri, which colonizes the light organs of certain marine ﬁshes and squids
(Fuqua et al., 1996).
In all of these examples the fundamental mechanism of quorum sensing is the
same, although the biochemical details are surprisingly different. These celltocell signaling systems are all composed of a small molecule, called an autoinducer,
which is synthesized by an autoinducer synthase, and a transcriptional activator
protein. The production of the autoinducer is regulated by a speciﬁc Rprotein. The
Rprotein by itself is not active without the corresponding autoinducer, but the Rprotein/autoinducer complex binds to speciﬁc DNA sequences upstream of the target genes, enhancing their transcription. Thus, at low cell density, the autoinducer
is synthesized at basal levels and diffuses into the surrounding medium, where it
is diluted. With increasing cell density, however, the intracellular concentration of
the autoinducer increases until it reaches a threshold concentration beyond which
it is produced autocatalytically, resulting in a dramatic increase of product concentrations. The autoinducer, therefore, allows the bacteria to communicate with each A QuorumSensing Model 97 other, to sense their own density, and together with a transcriptional activator to
express speciﬁc genes as a population rather than individually.
Quorum sensing is an unusual signaling mechanism for several reasons. First,
there are no membrane surface receptors for the signaling molecule, and since the
autoinducer is produced on the interior of the cell, the extracellular concentration
of autoinducer is generally smaller than the intracellular concentration. One might
therefore question how the autoinducer can act as a signal.
The purpose of this paper is to develop and study mathematical models of quorum sensing, in order to obtain a deeper understanding of how and when this mechanism works. We focus on the speciﬁc system of P. aeruginosa although the mathematical principles are the same in other systems. Our emphasis here is on the
biochemical mechanism underlying quorum sensing. In a related study, James et
al. (2000) found a simple mathematical model appropriate to V. ﬁscheri that had
multiple stable steady states and was controlled by extracellular autoinducer. By
contrast, in Ward et al. (2000), a population dynamics approach was taken for V.
ﬁscheri, and switching behavior in which there was cell differentiation was found
as the colony grew in size. The models in both of these studies were systems of
three ﬁrstorder, nonlinear, ordinary differential equations. In what follows, we
examine models involving both ordinary and partial differential equations. 2. T HE M ODEL The quorumsensing system of P. aeruginosa is unusual because it has two somewhat redundant regulatory systems. The ﬁrst system described in P. aeruginosa
was shown to regulate expression of the elastase LasB and was therefore named
the las system. The two enzymes, LasB elastase and LasA elastase, are responsible for elastolytic activity which destroys elastincontaining human lung tissue and
causes pulmonary hemorrhages associated with P. aeruginosa infections. The las
system is composed of lasI, the autoinducer synthase gene responsible for synthesis
of the autoinducer 3oxoC12HSL, and the lasR gene that codes for transcriptional
activator protein. The LasR/3oxoC12HSL dimer, which is the activated form of
LasR, activates a variety of genes, but preferentially promotes lasI activity. The
las system is positively controlled by both GacA and Vfr, which are needed for
transcription of lasR. The transcription of lasI is also repressed by the inhibitor
RsaL.
The second quorumsensing system in P. aeruginosa is named the rhl system
because of its ability to control the production of rhamnolipid. Rhamnolipid has
a detergentlike structure and is responsible for the degradation of lung surfactant
and inhibits the mucociliary transport and ciliary function of human respiratory
epithelium. This system is composed of rhlI, the synthase gene for the autoinducer
C4HSL, and the rhlR gene encoding a transcriptional activator protein. A diagram
depicting these two systems is shown in Fig. 1 (Van Delden and Iglewski, 1998). 98 J. D. Dockery and J. P. Keener Vfr
+
lasR + 3oxoC12HSL
LasR + LasI +
lasI +
LasR + GacA rsaL +
rhlR RsaL C4HSL
RhlI RhlR RhrR + rhlI Figure 1. Schematic diagram showing the gene regulation for the las and rhl systems in
P. aeruginosa. The ﬁrst step in our analysis is to ﬁnd equations describing the kinetics of this
system. For simplicity, we consider only the las system. We introduce variables for
all the concentrations (shown in Table 1). We assume that the dimer P is formed
via the law of mass action at rate k R A and degraded at rate k P ,
dP
= k R A R A − k P P.
dt (1) The enzyme LasR is used in the production of P and naturally degrades at rate k R ,
and is produced by the degradation of P , and by lasR mRNA at rate k1 ,
dR
= −k R A R A + k P P − k R R + k1r.
dt (2) Similarly, the autoinducer is used in the production of P and naturally degrades at
rate k A , and is produced by the degradation of P , and by the lasI enzyme at rate k2 ,
dA
= −k R A R A + k P P + k2 L − k A A .
dt (3) A QuorumSensing Model 99 Table 1. Variables used to identify concentrations.
Variable Concentration R
A
P
L
S
r
l
s LasR
3oxoC12HSL
LasR/3oxoC12HSL
LasI
RsaL
lasR mRNA
lasI mRNA
rsaL mRNA The enzyme LasI is produced by lasI mRNA at rate k3 and degrades at rate k L ,
dL
= k3l − k L L .
dt (4) The inhibitor RsaL is produced by rsaL mRNA at rate k4 and degrades at rate k S ,
dS
= k4 s − k S S .
dt (5) All messenger RNAs are produced by DNA at rates that are Michaelis–Menten in
type. Thus, for example, the inhibitor rsaL mRNA is produced at a Michaelis–
Menten rate depending on P , and degrades at some natural rate ks ,
ds
P
= Vs
− ks s ,
dt
Ks + P (6) and similarly, lasR mRNA is produced at a Michaelis–Menten rate depending on P ,
and degrades at some natural rate kr . We also assume that lasR mRNA is produced
at some basal rate r0 .
dr
P
= Vr
− kr r + r 0 .
(7)
dt
Kr + P
Finally, the production of lasI mRNA is activated by P and inactivated by S , degrades at rate ll , and is produced at some basal rate l0 ,
dl
P
1
= Vl
− kl l + l 0 .
dt
Kl + P K S + S (8) We would like to simplify this system by taking into account that some reactions
are fast compared to others. The difﬁculty here is that the rate constants in the
model are not known. There is evidence that many proteins are more stable than
the mRNA that code for them [see, for example, Anderson et al. (1998), Chalﬁe et
al. (1998) and Ehrenberg and Sverredal (1995)]. If this is the case here, then LasR 100 J. D. Dockery and J. P. Keener mRNA and lasI mRNA are much shorter lived than LasR and LasL, respectively,
so that kr and kl are much larger than k L and k R . With this assumption, we take l
and r to be in quasisteady state, so that
kl l = Vl 1
P
+ l0 ,
Kl + P K S + S kr r = Vr P
+ r0 .
Kr + P (9) The variable L can be understood as a ﬁrstorder linear ﬁlter, so that L tracks l
with some delay. The simplest approximation to this behavior is to ignore the delay
and take
k3l = k L L .
(10)
The quantity S inhibits the production of l , but it seems not to have much effect
on quorumsensing behavior. We ignore this variable by eliminating it from the
production term in equation (8).
With these simpliﬁcations, the governing system of equations becomes
dP
= k R A R A − k P P,
dt
dR
P
= −k R A R A + k P P − k R R + V R
+ R0 ,
dt
KR + P
dA
P
= −k R A R A + k P P + V A
+ A0 − k A A ,
dt
KL + P (11)
(12)
(13) (with redeﬁned parameters, of course). Finally, since the production of A and R
involves transcription of mRNA, it is probably slow compared to the binding and
unbinding of R and A to form complex P . Thus, we assume that P is in quasisteady state so that
kP P = kR A R A
(14)
and the governing equations become
dR
P
= −k R R + V R
+ R0 ,
dt
KR + P (15) dA
P
= VA
+ A0 − k A A ,
dt
KL + P (16) P= kR A
R A.
kP (17) Next, we need to determine how the density of organisms controls the activity
of this network. We assume that autoinducer A diffuses across the cell membrane,
and that the local density (volume fraction) of cells is ρ . Then by assumption, the
local volume fraction of extracellular space is 1 − ρ . A QuorumSensing Model 101 The extracellular autoinducer is assumed to diffuse freely across the cell membrane, with conductance δ , and naturally degrades at rate k E . If we suppose that
the density of cells is uniform and the extracellular space is well mixed, then the
concentration of autoinducer in the extracellular space, denoted E , is governed by
the equation
dE
(1 − ρ)
+ k E E = δ( A − E ).
(18)
dt
Here the factor 1 − ρ must be included to scale for the difference between concentration in the extracellular space and concentration viewed as amount per unit total
volume. Similarly, the governing equation for intracellular autoinducer A must be
modiﬁed to account for diffusion across the cell membrane, yielding
ρ dA
P
− VA
− A0 + k A A = −δ( A − E ).
dt
KL + P (19) There is some evidence that the transport of the autoinducer may involve both
passive diffusion and a cotransport mechanism (Pearson et al., 1999). For this
model we assume that diffusion alone acts to transport the autoinducer [as is apparently correct for the rhl system (Pearson et al., 1999)]. The inclusion here of an
additional cotransport mechanism is possible but seems not to be a crucial ingredient.
It is now fairly easy to see that this system of equations exhibits quorum sensing.
One way to see this is to take E to be in quasisteady state. This is probably not
a good assumption since A and E are the same chemical and so probably degrade
at about the same rates inside and outside the cell. However, this assumption allows us to study this system in the phase plane. The behavior is not changed if
a threevariable system is used, but the analysis of the threevariable systems is
substantially more complicated with little increase in insight.
The resulting system is (with redeﬁned parameters)
P
dR
= VR
− k R R + R0 ,
dt
KR + P (20) dA
P
= VA
+ A0 − d (ρ) A ,
dt
KA + P (21) k
ρ)
δ
A
where P = k Rk PR A and d (ρ) = k A + ρ δ+E (1(−−ρ) .
kE 1
Quorum sensing works because of the ρ dependence of d (ρ). In particular, when
ρ is small, d (ρ), the decay rate for A, is large, while when ρ is close to one, d (ρ)
is small. This has the feature of modifying the A nullcline in an important way.
The nullclines for the system (20), (21) are shown in Fig. 2. Here it is readily seen
that for small values of ρ and for large values of ρ there is a unique steadystate
solution in the positive R – A quadrant. For small values of ρ the steadystate values 102 J. D. Dockery and J. P. Keener 4 A 3 2 1 0 0 1 2 3 4 R
RA
RA
Figure 2. Nullclines d R = −k R R + KV R+ R A + R0 = 0 and d A = KV A+ R A + A0 −d (ρ) A =
dt
dt
R
A
0 with three different values of ρ , shown for parameter values V R = 2.0, V A = 2.0,
K R = 1.0, K A = 1.0, R0 = 0.05, A0 = 0.05, δ = 0.2, k E = 0.1, k R = 0.7, and
k A = 0.02. of R and A are small, while for large values of ρ the steady statevalues of R and
A are large.
For intermediate values of ρ , there are three steady solutions in the positive quadrant. The small and large values are stable steady solutions while the intermediate
solution is unstable, a saddle point.
The parameter ρ provides a switch between the two stable steady solutions. For
small values of ρ , there is a unique stable steady state with small values of R and
A. As ρ increases, two more solutions appear (a saddle node bifurcation), and as ρ
increases yet further, the small solution disappears (another saddle node) leaving a
unique solution, thus initiating a ‘switch’.
It is easy to arrange parameter values that have this switching behavior, and the
switch can be adjusted to occur at any desired density level.
An even simpler model that has the same qualitative behavior (but is not correct
from a biochemical perspective) is found by setting R = A to obtain
A2
dA
= VA
+ A0 − d (ρ) A .
dt
K A + A2 (22) A QuorumSensing Model 103 0.2 dA/dt 0.1 0.0 − 0.1 − 0.2 0 1 2 3 A
2
Figure 3. Plot of d A = V A A 2 + A0 − d (ρ) A as a function of A for several values of
dt
K A+A
ρ with parameter values V A = 1.0, K A = 1.0, A0 = 0.05, δ = 0.5, k E = 0.1, and
k A = 0.02. Here, the curve d A , shown in Fig. 3, is a monotone increasing function of ρ and
dt
for small and large values of ρ there is a unique positive steadystate solution. For
small values of ρ the steadystate value of A is small, while for large values of ρ
the steadystate value of A is large.
For intermediate values of ρ , the curve d A is ‘nshaped’ and there are three steady
dt
positive solutions. The small and large values are stable steady solutions while the
intermediate solution is unstable.
The switching behavior is readily observed by varying ρ . For small values of ρ ,
there is a unique stable steady state with a small value of A. As ρ increases, two
more solutions appear, and as ρ increases yet further, the small solution disappears
leaving a unique solution, thus initiating a ‘switch’.
We can now give a verbal explanation of how quorum sensing works. The quantity A, the autoinducer, is produced by cells at some nominal rate. However, the
cell must dump its production of A or else the autocatalytic reaction would turn on.
As the density of cells increases, the dumping process becomes less effective, and
so the autocatalytic reaction is turned on. In other words, quorum sensing takes
place because the ‘drain’ for A backs up. Note that the external concentration of 104 J. D. Dockery and J. P. Keener A is always smaller than the internal concentration of A. Thus, A is not actually a
signaling chemical, and there is no necessity for surface receptors for A, and there
is no need to develop high concentrations of autoinducer in the extracellular space. 3. A PDE M ODEL While the above model shows the desired behavior in a homogeneous environment, such as a chemostat, a more realistic model would take into account the
possibility of inhomogeneous distributions of the autoinducer in the extracellular
space, while the cells themselves remain relatively motionless.
We suppose there is a uniform layer of cells of thickness L and ﬁxed volume
fraction ρ which are attached to a substratum in an aqueous bath. As above, the
cells produce the autoinducer with concentration A which has extracellular concentration E ,
δ
dA
= F ( A) + ( E − A),
dt
ρ (23) ∂E
∂2 E
δ
=
+
( A − E ) − kE E .
2
∂t
∂x
1−ρ (24) The factor ρ in equation (23) and the factor 1 − ρ in (24) are necessary to account
for the differences in volume fraction between intracellular and extracellular space.
We assume that the cells occupy the onedimensional region 0 < x < L . At the
boundary of the cellular domain, x = L , we assume that there is mass transfer into
the bulk ﬂuid. We model this simply by assuming the Robin boundary condition
Ex (L , t ) + α E (L , t ) = 0 (25) where α is a positive parameter. At the substratum we impose the Neumann boundary condition
E x (0, t ) = 0.
(26)
The goal of this section is to ﬁnd the steadystate solutions of this pde, i.e.,
Ex x + δ
( A − E ) − k E E = 0,
1−ρ with F ( A) + δ
( E − A) = 0 (27)
ρ on 0 < x < L with E x (0) = 0 and E x ( L ) + α E ( L ) = 0.
The analysis of this equation is accomplished in the phase plane. First, we must
determine the nature of the solutions of the equation
F ( A) + δ
( E − A) = 0.
ρ (28) A QuorumSensing Model 105 5
4 E 3
2
1
0
0 1 2 3 4 5 A
Figure 4. Plot of the function E ( A) = A − ρ F ( A) for K A = 1.0, V A = 1.0, A0 = 0.05,
δ
δ = 0.1, ρ = 0.3, and k A = 0.2.
0.20
g+(E) g(E)
0.15
0.10 g (E)
0 0.05
g−(E) 0.00
− 0.5 E mn E mx
0.0 0.5 1.0 1.5 E0 2.0 E
δ
Figure 5. Plot of the ‘function’ g ( E ) = 1−ρ ( A−1 ( E ) − E ) − k E E for K A = 1.0,
V A = 1.0, A0 = 0.05, δ = 0.1, ρ = 0.3, k A = 0.2, and k E = 0.1. We begin by stating general assumptions for the function F ( A), namely that it is
‘nshaped’ and strictly positive on the interval 0 ≤ A < A∗ , and that F ( A∗ ) = 0.
A speciﬁc example of the function F ( A) is
F ( A) = V A A2
− k A A + A0 .
K A + A2 (29) This function is ‘nshaped’ if k A < 3V A K3A and has only one zero if A0 is sufﬁ8
ciently large.
It is easy to determine E as a function of A from (28) (shown in Fig. 4). We
further assume that ρ is sufﬁciently large so that the function E ( A) = A − ρ F ( A) is
δ
δ δ
an ‘nshaped’ function. For the function (29), this occurs whenever ρ < 3V A K3A −
8
k A . Since E ( A) = A − ρ F ( A) is ‘nshaped’, there are values Amx and Amn , with
δ 106 J. D. Dockery and J. P. Keener 0.50 dE/dx 0.25 0.00 − 0.25 − 0.50 0.0 0.5 1.0 1.5 E
Figure 6. Phase portrait of the solution of the boundary value problem on the upper branch
g+ ( E ), with parameter values K A = 1.0, V A = 1.0, A0 = 0.05, δ = 0.1, ρ = 0.3,
k A = 0.2, and k E = α 2 = 0.1. Dashed lines are the curves on which E x = −α E and
E x = 0. Amx < Amn < A0 , which are local maximum and local minimum, respectively, for
the function E ( A), and that E mx = E ( Amx ). The value E mn = E ( Amn ) may be
positive or negative, however, it is required that E mn < E mx .
Now we invert the function E ( A) by reversing the axes, and note that the inverse
function has three branches, a lower, middle and upper branch. For each of these
we ﬁnd the quantity
g( E ) = δ
( A−1 ( E ) − E ) − k E E ,
1−ρ where (30) ρ
F ( A−1 ( E )) = E .
(31)
δ
Of necessity, g ( E ) has three branches, which we denote as g− ( E ), g0 ( E ), and
g+ ( E ), shown in Fig. 5. These are deﬁned on the three subintervals 0 ≤ E ≤ E mx
(for g− ( E )), E mn < E < E mx (for g0 ( E )), and E mx < E < E 0 (for g+ ( E )), and
g+ ( E 0 ) = 0. Where they are comparable, g− ( E ) < g0 ( E ) < g+ ( E ).
The ﬁnal restriction placed on the parameters is that we require k E to be small
enough so that g− ( E ) > 0 on the interval 0 ≤ E ≤ E mx .
A−1 ( E ) − A QuorumSensing Model 107 1.6 E(0) 1.2 0.8 0.4 0.0 0 2 4 6 8 10 L
Figure 7. Steady solution branches with E (0) plotted as a function of L shown for parameter values K A = 1.0, V A = 1.0, A0 = 0.05, δ = 0.1, ρ = 0.3, k A = 0.2, and
k E = α 2 = 0.1. Now we look for solutions of the ordinary differential equation
E x x + g( E ) = 0 (32) on the interval 0 < x < L and subject to boundary conditions E x (0) = 0 and
E x ( L ) = −α E ( L ) where α > 0. The easiest way to demonstrate the existence of
such solutions is in the E − E x phase plane, depicted in Fig. 6. In this plane, the
boundary conditions are depicted as two straight lines E x = 0 and E x = −α E . In
Fig. 6, trajectories are shown for the function g+ ( E ) although the trajectories for
g− ( E ) are qualitatively identical.
We use the function g− ( E ) to construct the ﬁrst solution branch. In fact, for any
E − with E − < E mx , the trajectories with E ≤ E − and
12
E−
2x E− g− ( E )d E = 0 (33) E work. It follows that for each E − with 0 < E − < E mx there is a corresponding
value of L , say L = L − ( E − ), for which this solution exists. Furthermore, since 108 J. D. Dockery and J. P. Keener 0.6
0.5 E 0.4
0.3
0.2
0.1
0.0 0.0 0.2 0.4 0.6 0.8 1.0 x
Figure 8. Steady solution proﬁles E (x ) plotted as a function of x on the interval 0 < x <
L = 1.0 for which E (0) = 0.057 on the lower branch and E (0) = 0.53 on the upper
branch, shown for parameter values K A = 1.0, V A = 1.0, A0 = 0.05, δ = 0.1, ρ = 0.3,
k A = 0.2, and k E = α 2 = 0.1. g− ( E ) is positive and bounded away from 0, L − ( E − ) remains bounded on its range
of existence. In other words, there is a number, say L ∗ , which bounds L − ( E − )
above, L ∗ ≤ L − ( E − ). This solution branch is depicted by the smaller of the two
curves in Fig. 7 with E (0) shown plotted as a function of L . This solution branch
terminates at L ∗ since the trajectories for the phase portrait can only be deﬁned
on the range of deﬁnition of g− ( E ). A typical solution proﬁle is depicted by the
smaller of the two curves in Fig. 8.
The second branch of solutions of interest uses g+ ( E ) on the interval max{ E mn , 0}
< E < E 0 . The phase portrait trajectories are solutions of
12
E−
2x E+ g+ ( E )d E = 0 (34) E provided E ≤ E + ≤ E 0 . Now, however, since g+ ( E 0 ) = 0, the point E =
E 0 , E x = 0 is a saddle point in the phase plane. Thus, the trajectory emanating
from this point corresponds to L = ∞. It follows that there is a second solution
branch which can be represented as L = L + ( E + ) with L 1 < L < ∞. Furthermore,
if E mn < 0, then L 1 = 0. This solution branch is depicted by the larger of the two A QuorumSensing Model 109 0.4 g(E) 0.3 0.2 0.1 0.0
0.5 0.0 0.5 1.0 1.5 2.0 E
δ
Figure 9. Plot of the ‘function’ g ( E ) = 1−ρ ( A−1 ( E ) − E ) − ke E for K A = 1.0, V A =
1.0, A0 = 0.05, δ = 0.1, k A = 0.2, k E = 0.1, and three different values of ρ . curves in Fig. 7 with E (0) shown plotted as a function of L . This solution branch
extends to L = ∞ since the trajectory that contains the point E 0 is the unstable
manifold of a saddle point. A typical solution proﬁle is depicted by the larger of
the two curves in Fig. 8.
Quorum sensing in a spatially distributed system can now be described. If the
density ρ is too small or the size L is small, the small steady solution is attained.
However, if the density is sufﬁciently large, then as the colony size L grows, the
small steady solution ceases to exist, and the solution switches to the large steady
solution. Note that this switch fails if k E is too large, since then the curve g− ( E ) is
not strictly positive, and therefore the small solution branch exists for all positive L .
A second way to describe quorum sensing is as a function of density ρ . We
suppose that there is a colony of a ﬁxed size L that is growing so that the density ρ
is slowly increasing.
It is clear from Fig. 9 that for ρ sufﬁciently small there is only the single solution
branch g− ( E ), while for ρ sufﬁciently large the only branch deﬁned for positive E
is the g+ branch. In other words, if δ is sufﬁciently small, there is a value of ρ < 1
above which E mx < 0. However, if E mx < 0, the lower branch of solutions fails to
exits. Thus, for any ﬁxed L , there is a value of ρ above which the lower solution
fails to exist. At this value of ρ , the solution must switch to the upper branch, if
stable, thereby dramatically increasing the production of autoinducer.
We have numerically computed a branch of steadystate solutions of (23) and (24)
as a function of the density ρ for ﬁxed domain size L . These results are shown in
Fig. 10. We have found that the upper and lower branches are stable while the 110 J. D. Dockery and J. P. Keener
6
5 A(0) 4
3
2
1
0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 10. Plot of the bifurcation diagram for the steady states of (23) as a function of the
cell density ρ with L ﬁxed at L = 2. The parameter values are K A = 1.0, V A = 1.0, A0 =
0.05, δ = 0.1, k A = 0.2, k E = 0.1, ρ = 0.3, α = 1.0, and L = 2.0. middle branch is unstable. As the cell density increases, the production of the autoinducer increases slowly until a critical value of the density is reached. At this
point there is a large increase in autoinducer production. Once on the upper branch,
the population maintains a high rate of autoinducer production even if the cell density decreases. However, if the cell density falls below the value at the knee on the
upper branch, production falls back to a nominal rate.
It has been shown (Davies et al., 1998) that autoinducers can increase polymer
production. With increased polymer the cell density would decrease. Thus we see
that hysteresis assures that the switch to high autoinducer production is not easily
reversed by a decrease in cell density due to an increase in polymer production.
It is also easy to speculate that once the cell density becomes low enough that the
autoinducer is reset to basal levels, the genes responsible for polymer production
are turned off. If there is too much polymer, the diffusion of vital nutrients to the
population could be hindered.
In summary, quorum sensing can be accomplished by increasing the density of a
colony of ﬁxed size, or by increasing the size of a colony of large enough density.
In the next section, stability results are discussed. 4. S TABILITY R ESULTS The stability of the steady state can be determined by numerical simulations of
the system or by linearization (Weinberger, 1983).
The following proposition implies that the large and small amplitude steadystate
solutions to our system are stable provided that on the range of the steady solution A QuorumSensing Model
L=2 2 L=2 0.22 A
E 1.8 A
E 0.2 1.6 111 0.18
0.16 1.4 0.14 1.2 0.12
1
0.1
0.8 0.08 0.6 0.06 0.4 0.04 0.2 0 0.2 0.4 0.6 0.8
x 1 0.02 0 0.2 0.4 0.6 0.8
x 1 Figure 11. Plot of two different solutions of (23), (24) at t = 200. The parameter values
are K A = 1.0, V A = 1.0, A0 = 0.05, δ = 0.1, k A = 0.2, k E = 0.1, ρ = 0.3, α = 1.0,
and L = 2.0. ( A0 , E 0 ), g ( E 0 (x )) < 0. We show this by constructing timeindependent upper
and lower solutions that are arbitrarily close to the steadystate solutions.
The system (23), (24) is a cooperative system (quasimonotone), the righthand
side of (23) is an increasing function of E and the righthand side of (24) is increasing in A. This allows one the following Comparison Principle.
D EFINITION 1. The pair ( A, E ) is an upper solution pair for the system (23), (26)
if and only if
At ≥ F ( A) +
Et ≥ δ
( E − A)
ρ ∂2 E
δ
+
( A − E ) − kE E
2
∂x
1−ρ (35)
(36) E x (L ) + α E (L ) ≥ 0 (37) E x (0) ≤ 0. (38) The pair ( A, E ) is a lower solution pair if it satisﬁes the above with each inequality
reversed.
T HEOREM 1 (C OMPARISON P RINCIPLE ). Let ( A(x ), E (x )) be a tindependent
upper solution pair for the system (23)–(26) and let ( A(x , t ), E (x , t )) be the solution to the system with the initial data ( A(x ), E (x )). Then A and E are nonincreasing functions of t and approach the largest steadystate solution ( A∗ (x ), E ∗ (x )) 112 J. D. Dockery and J. P. Keener 3 L
E(0) 2.5
2
1.5
1
0.5
0
0 20 40 60 80 100 120 140 160 180 200 t
Figure 12. Plot of L (t ) and E (0, t ) starting from constant initial data E (x , 0) = 0.0,
A(x , 0) = 0.0, L (0) = 0.01. The parameter values are K A = 1.0, V A = 1.0, A0 = 0.05,
δ = 0.1, k A = 0.2, k E = 0.1, ρ = 0.3, α = 1.0, and L = 2.0. 3 L
E(0) 2.5
2
1.5
1
0.5
0
0 50 100 150 200 250 300 350 400 t
Figure 13. Plot of L (t ) and E (0, t ) starting from constant initial data E (x , 0) = 0.0,
A(x , 0) = 0.0, L (0) = 0.01. The parameter values are K A = 1.0, V A = 1.0, A0 = 0.05,
δ = 0.1, k A = 0.2, k E = 0.1, ρ = 0.3, α = 1.0, and L = 2.0. A QuorumSensing Model 113 that satisﬁes A∗ (x ) ≤ A(x ) and E ∗ (x ) ≤ E (x ). Similarly, let ( A(x ), E (x )) be a
lower solution pair to the system and ( A(x , t ), E (x , t )) a solution with ( A(x ), E (x ))
as initial data. Then ( A(x , t ), E (x , t )) are nondecreasing in t and tend to the
smallest stationary solution, ( A∗ (x ), E ∗ (x )) that lies above ( A(x ), E (x )).
Let ( A0 (x ), E 0 (x )) be a steadystate solution of (23)–(26) and let
but arbitrary. Consider > 0 be ﬁxed ˆ
E (x ) = E 0 (x ) + (39) ˆ
A(x ) = A−1 ( E 0 (x ) + ). (40) ˆˆ
We show that ( A, E ) is an upper solution pair provided certain conditions are satisﬁed. First, we assume that our steady solution lies either on the upper branch of
A−1 or the lower branch. It is clear from Fig. 4 that A−1 is an increasing function
of E on these two branches, thus
ˆ
A(x )q = A−1 ( E 0 (x ) + ) ≥ A−1 ( E 0 (x )) = A0 (x )
ˆˆ
if > 0 is sufﬁciently small. It follows that ( A, E ) satisﬁes (35). It is also clear
that the boundary inequalities (37), (38) are satisﬁed.
Now consider inequality (36). If we assume that g is nonincreasing on the range
of E 0 , then for small enough we obtain
ˆ
ˆ
E x x + g( E ) = E0 x x + g( E0 + ) (41) ≤ E0 x x + g( E0) = 0 (42) and we see that (36) is satisﬁed. The following result follows from Theorem 1.
P ROPOSITION 1. Let ( A0 , E 0 ) be a steadystate solution which lies on the upper
or lower branch of g, g± . If on the range of E 0 (x ), g is nonincreasing, then for
sufﬁciently small and positive ( A0 (x ) + , E 0 (x ) + ) is an upper solution pair and
( A0 (x ) − , E 0 (x ) − ) is a lower solution pair. Furthermore, in this case ( A0 , E 0 )
is a stable solution of (23)–(26).
While in general we cannot verify that g± are decreasing functions for all parameter values, it is clear from Fig. 9 that this is the case for g+ ( E ) for large enough
values of ρ , and for g− ( E ) for small enough values of L .
There are also steady solutions that correspond to the middle branch of g , g0 . We
can show that these solutions are unstable.
P ROPOSITION 2. Steadystate solutions ( A0 (x ), E 0 (x )) for which E ( A0 (x )) <
0 are linearly unstable. 114 J. D. Dockery and J. P. Keener To prove this, consider the linearization of (23), (24) about a steady solution
( A0 , E 0 )
at = F [ A0 (x )]a +
et = ex x + δ
(e − a ),
ρ δ
(a − e) − k E e.
1−ρ (43)
(44) δ
Let α(x ) = F [ A0 (x )] − ρ , and note that α(x ) = − ρ E [ A0 (x )], which is posδ
itive on the middle branch. Suppose that (a , e) is a solution of (43), (44) with
(a (x , 0), e(x , 0)) = (a0 (x ), e0 (x )) and such that e(x , t ) remains bounded for all
t ≥ 0, say e(x , t ) ≤ M . Then from (43) we see that at ≥ α(x )a − M δ
ρ which implies that
a (x , t ) ≥ exp(α(x )t )a0 (0) + Mδ
(1 − exp(α(x )t ))
ρα(x ) so that a (x , t ) cannot remain bounded for all t > 0 if we start with positive initial
data, which implies instability. If e(x , t ) is unbounded for t > 0 starting from
arbitrarily small initial data, then clearly the steady solution is linearly unstable. 5. N UMERICAL S IMULATIONS In this section we report on numerical simulations which illustrate the stability
results of the previous section. We also indicate how these results may be applicable to bioﬁlms.
We take α = 1 and scale the x interval to have unit length. In Fig. 11 are
displayed the two stable steady states that exist when L = 2. One is an order
of magnitude larger than the other. The small solution was generated from the
constant initial data A(x , 0) = 0.0 and E (x , 0) = 0. The larger solution was
also generated from constant initial data, A(x , 0) = 1.0 and E (x , 0) = 1.0. Both
solutions are shown at time t = 200. At this time the time derivatives of both A
and E are less than 10−10 so we take these to be steady states.
One can model a growing population by adding dynamics for the domain length
L . Here we set
dL
= κ( L f − L )
(45)
dt
where κ = 0.04 and L f = 3. A QuorumSensing Model 115 The time series of a solution of (23), (24) and (45) with constant initial data is
shown in Fig. 12. Here we display the length of the interval L and the value of
the extracellular autoinducer at x = 0, E (0, t ) as functions of time. From other
numerical calculations we know that the small amplitude solution does not exist at
L = 2.2. We see that after a lag the solution of the evolution equation approaches
the large amplitude steady solution that exists at L = 3.
To demonstrate the hysteresis indicated in Fig. 7, we use the same dynamics for
L as above in a second run until t = 180 whereupon we set L = 2. This could be
thought of as a major detachment event whereby the bioﬁlm thickness is suddenly
reset to two thirds of its asymptotic thickness. The time series of such a solution
is shown in Fig. 13. We see that upon reset, the autoinducer concentration remains
high. The hysteresis assures that the switch to high autoinducer production is not
easily reversed. This is very important for a bioﬁlm since these populations are
often subjected to major detachment and slufﬁng events whereby there are large
and sudden decreases in population size.
6. D ISCUSSION We have presented a simple model of quorum sensing using autoinduction that is
based on the known biochemistry of P. aeruginosa. Using this model we demonstrate that quorum sensing works because the rate of elimination of autoinducer
depends on the colony size and density. Thus, the autoinducer production switches
to its high state when the elimination of autoinducer from the extracellular space is
decreased.
This biochemical switch is hysteretic so that autoinducer production switches on
at different (higher) levels than it switches off. This hysteresis is possibly important
in the regulation of the production of exopolysaccharide, which tends to decrease
bacterial density. The hysteresis predicted by this model investigation has not been
veriﬁed experimentally.
ACKNOWLEDGEMENTS
JDD was supported in part by NSF grant DMS9805701. JPK was supported in
part by NSF grant DMS99700876. R EFERENCES
Anderson, J. B., C. Sternberg, L. K. Poulsen, M. Givskov and S. Molin (1998). New unstable variants of green ﬂuorescent protein for studies of transient gene expression in
bacteria. Appl. Environ. Microbiol. 64, 2240–2246.
Chalﬁe, M., Y. Tu, G. Euskirchen, W. W. Ward and D. C. Prasher (1994). Green ﬂuorescent
protein as a marker for gene expression. Science 263, 802–805. 116 J. D. Dockery and J. P. Keener Characklis, W. G. and K. C. Marshall (1990). Bioﬁlms, New York: John Wiley & Sons,
Inc.
Davies, D. G., M. R. Parsek, J. P. Pearson, B. H. Iglewski, J. W. Costerton and E. P.
Greenberg (1998). The involvement of celltocell signals in the development of bacterial bioﬁlm. Science 280, 295–298.
Ehrenberg, M. and A. Sverredal (1995). A model for copy number control of the plasmid
R1. J. Mol. Biol. 246, 472–485.
Fuqua, C., S. C. Winans and E. P. Greenberg (1996). Census and concensus in bacterial
ecosystems: The luxRluxI family of quorumsensing transcriptional regulators. Annu.
Rev. Microbiol. 50, 727–751.
Govan, J. R. and V. Deretic (1996). Microbial pathogenesis in cystic ﬁbrosis: mucoid
pseudomonas aeruginosa and Burkolderia cepacia. Microbiol. Rev 60, 539–574.
James, S., P. Nilsson, G. James, S. Kjelleberg and T. Fagerstr¨ m (2000). Luminescence
o
control in the marine bacterium Vibrio ﬁscheri: An analysis of the dynamics of lux
regulation. J. Mol. Biol. 296, 1127–1137.
Pearson, J. P., C. Van Delden and B. H. Iglewski (1999). Active efﬂux and diffusion are involved in transport of Pseudomonas aeruginosa celltocell signals. J. Bacteriol. 181(4),
1203–1210.
Van Delden, C. and B. H. Iglewski (1998). Celltocell signaling and Pseudomonas aeruginosa infections. Emerging Infect. Dis. 4, 551–560.
Ward, J. P., J. R. King and A. J. Koerber (2000). Mathematical modelling of quorum sensing in bacteria. Preprint.
Weinberger, H. F. (1983). A simple system with a continuum of stable inhomogeneous steady states, in Nonlinear Partial Differential Equations in Applied Science
(Tokyo, 1982), Amsterdam: NorthHolland, pp. 345–359.
Received 16 January 2000 and accepted 3 September 2000 ...
View
Full
Document
 Summer '06
 DeLeenheer

Click to edit the document details