Stability and Sensitive Analysis of a Model with Delay Quorum Sensing

Abstract

This paper formulates a delay model characterizing the competition between bacteria and immune system. The center manifold reduction method and the normal form theory due to Faria and Magalhaes are used to compute the normal form of the model, and the stability of two nonhyperbolic equilibria is discussed. Sensitivity analysis suggests that the growth rate of bacteria is the most sensitive parameter of the threshold parameter and should be targeted in the controlling strategies.

1. Introduction

Quorum sensing is a process that enables bacteria to communicate using secreted signaling molecules called autoinducers [1]. It makes bacteria regulate their gene expression collectively and control their behaviors on community scale. Quorum sensing was initially observed in the marine bacterium Vibrio fischeri about 30 years ago [2, 3]. Now, many other species are observed to exhibit quorum sensing behavior, including major human pathogens such as Staphylococcus aureus and Pseudomonas aeruginosa. Quorum sensing has received more and more attention (see [4–15] and the references therein) and some models are formulated to investigate its effect on the transmission of disease. Braselton and Waltman [4] formulated the dynamically allocated inhibitor production. Dockery and Keener [5] were devoted to developing and studying an ODE and a PDE mathematical models for quorum sensing in Pseudomonas aeruginosa and found 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. Koerber et al. [6] presented a mathematical model for the early stages of the infection process by Pseudomonas aeruginosa in burn wounds which accounts for the quorum sensing and the diffusion of signalling molecules in the burn-wound environment, and the effects of important parameters on the dynamic properties of the model are discussed in detail. They gave some sufficient conditions for the global asymptotic stability of two boundary equilibria which, respectively, correspond to the survival of the allelochemical producer species or the susceptible one. Fergola et al. [7] formulated an allelopathic competition model in which a distributed delay term simulates quorum sensing which regulates the delay production process of allelochemicals, and they proved the unique existence of the positive solution and the stability of biologically meaningful steady-state solutions. Anguige et al. [9] constructed a multiphase mathematical model of quorum sensing in a maturing Pseudomonas aeruginosa biofilm to investigate the effect of antiquorum sensing and antibiotic treatments on the exopolysaccharide concentration, signal level, bacterial numbers, and biofilm growth rate. The above articles leave out the immune response to the bacterial invasion. However, the immune status of the hosts has a significant impact on the transmission of an infection in a population. Literatures [11–15] employed a quadratic function to describe the quorum sensing of bacteria and formulate some models to characterize the competition between bacteria and the immune system, in which the existence of periodical solution, chaotic motion, and subharmonic bifurcation, the properties of Hopf bifurcation, and the stability of equilibrium et al. were investigated. As a novelty of this paper, a cubic function is used to express quorum sensing. This makes the model own more general nonlinearity, which results in much wider set of outcomes including the coexistence of multiple positive equilibria and the existence of critical equilibrium with simple zero singularity.

2. Model Formulation

We denote by the concentration of the uninfected target cells, the concentration of the infected target cells, the concentration of the bacteria, the concentration of the innate cells, and the concentration of the adaptive cells. The dynamic relations among them are as follows: the uninfected target cells have a natural turnover and a half-life and they are infected by bacteria with mass-action term ; the infected target cells are cleared by half-life or adaptive immune cells with mass action term ; both the innate and the adaptive immune cells have a source term and a half-life time; for the innate immunity, the source term includes a wide range of cells involved in the first wave of defense of the host such as natural killer cells, polymorphonuclear cells, macrophages, and dendritic cells, and for the adaptive immunity, the source term represents the memory cells, derived from a previous infection or vaccination, a zero source means the first infection with this pathogen and there are no memory cells; both of the two kinds of cells are increased by the signals captured by the bacteria load; the bacteria population has a net growth term represented by a logistic function and it is cleared by the innate immunity with mass action term . Here, we use a function to formulize the bacteria that compete with the immune cells at time , which receive signal molecules time units ago. is a positive constant, is the growth rate of bacteria, and is the effective carrying capacity of the environment. Consequently, the vital dynamics are governed by

Remark 1. The first equation of system (2) suggests that bacteria are controlled and increased by quorum sensing except for their net growth and they are cleared by the innate immune cells. The second equation of system (2) characterizes the dynamics of the uninfected target cells, and the third one reflects the dynamics of the infected target cells. The uninfected target cells are infected by bacteria in mass action law and they have their own constant input flow, and the infected target cells are killed by the adaptive immune cells. The last two equations of system (2) show that each kind of the immune cells has a special source term and their responses are enhanced by the bacteria load. The target cells and the immune cells have their own half-life terms.

3. The Existence and Stability of Equilibria

We introduce

Theorem 2. System (2) always admits a bacteria-free equilibrium . If , then system (2) admits a unique positive equilibrium . If , then system (2) admits a unique positive equilibrium . If and , then system (2) admits no positive equilibrium. If and , then system (2) admits a unique positive equilibrium . If and , then system (2) has two positive equilibria and .
Specifically,

Theorem 3. If , the bacteria-free equilibrium is asymptotically stable, while if , is unstable.

The proofs for Theorems 2 and 3 are trivial, so omit them.

In the sequel, we study the stability of the positive equilibrium . First of all, we transfer it to the origin and get where

The characteristic equation of (5) at the origin is where Equation (7) has two negative roots , , and the other roots can be obtained by solving the following equation: It can be seen that

Theorem 4. For , both and are unstable when they exist.

Obviously, does not own zero eigenvalue singularity; namely, zero is not a eigenvalue of the Jacobi matrix at . By using the Routh-Hurwitz stability criterion, is locally asymptotically stable if , , and are satisfied together.

Theorem 5. For , is asymptotically stable if , , and are small sufficiently.

Proof. Let and . From Theorem 2, exists if and only if and , which leads to and ; that is, and . By using , it follows that approaches if does so. Because we have Note that . We have Seen from above formulae, , , , and are positive when is small, which results in the locally asymptotical stability of .

Clearly, the left side of (9) is continuous in and has roots with positive real parts if and only if it has purely imaginary roots. We will determine whether (9) has purely imaginary roots or not, from which we then will be able to get conditions for all eigenvalues to have negative real parts.

Denote the eigenvalue of the characteristic equation (5) by , where , continually depend on the delay . Under the same conditions as Theorem 5, we have . Since is continuous in , one still has and remains stable if is sufficiently small. If there exists a positive value satisfying , that is, is a purely imaginary root of (9), then loses its stability and eventually becomes unstable when becomes positive. On the other hand, if such a does not exist is always stable.

Obviously, (9) has a purely imaginary root , , if and only if Separating the real and imaginary parts of (14) and adding up the squares of them lead to where , ,  , and .

It can be verified that and (15) has positive roots. Without loss of generality, one assumes that (15) has three positive roots defined by , , and , respectively. Then, (14) has three positive roots , . And then, we have Thus, if we denote where , , then is a pair of purely imaginary roots of (9). Define . We have the following.

Theorem 6. Under the same conditions as Theorem 5, is asymptotically stable for .

Theorem 7. If , system (2) undergoes Hopf bifurcation at the positive equilibrium when .

4. Normal Forms on the Center Manifold

From the discussions in the above section, it can be seen that the Jacobi matrix at has a uniquely simple zero eigenvalue if and . To determine the dynamic properties of , we have to compute the normal forms on the center manifold. The method used is based on the center manifold reduction and normal form theory due to Faria and Magalhaes; see [16, 17].

By means of , it obtains where It is seen from (19) that if and are small enough, which means there exists a unique positive solving (18). Denote it by and define as where is a small parameter. Obviously, if . Next, we transfer to the origin by Normalizing the delay by , denoting by , and neglecting the higher order terms , we have where

Let be the Banach space of continuous functions from into with supremum norm. Define as , . Equation (22) can be written as the functional differential equation where is a neighborhood of zero in space of real numbers, is a parameterized family of bounded linear operators, and is a function with , for all . They have the following respective forms: where

From the Riesz representation theorem the linear map can be expressed in integral form as follows: where is a bounded variation matrix-valued function on . In fact, we can define as where is the Dirac delta function.

Let be the -dimensional vector space of row vectors, and denote . Then, the adjoint bilinear form on can be defined as where and .

Let be the infinitesimal generator of the flow for the linear system with spectrum . The adjoint operator is defined as the infinitesimal generator for the solution operator of the adjoint equation in :

It is well known that the eigenvalues of with zero real parts play an important role in the bifurcation theory of RFDES. Denote by , and let . We have .

Using the formal adjoint theory for FDEs in [18], the phase space can be decomposed by as , where is the generalized eigenspace associated with the eigenvalues in , , and the dual space is the generalized eigenspace for associated with the eigenvalues in . Assume that and are the respective dual bases of and and satisfy . We might as well choose and as follows: where is the transpose of and .

Let . Then, the following equations hold simultaneously:

As shown in [16, 17], an appropriate phase space for considering normal forms of (24) is the Banach space BC of functions from into , which are uniformly continuous on and with a jump discontinuity at 0. Then, the elements of BC have the form , where , , and so that BC is identified with with the norm .

Let BC denote the projection and then the decomposition yields a decomposition of BC by as the topological direct sum BC with the property , where is an infinite dimensional complementary subspace of and as shown above. Now, we decompose in (24) as , where and , is the subset of consisting of continuously differentiable functions.

Next, we rewrite (24) as follows: And, then, under the composition , (24) can be decomposed as a system of ODEs in as follows: where is the restriction of as an operator from into , and

As for autonomous ODEs in , the normal forms are obtained by a recursive process of changes of variables. At a step , the terms of order are computed from the terms of the same order and from the terms of lower orders already computed in previous steps. Assuming that steps of orders , have already been performed leads to where is the terms of order in after the previous transformations of variables and h.o.t stands for the higher order terms. Following the algorithm of   [16, 17] at step , using a change of variables of the form where , and , are homogeneous polynomials of degree in and , after dropping the hats for simplification of notations, (37) can be put into the normal form where

It can be verified that (24) satisfies nonresonance conditions; see [16, 17]. Then, the locally invariant manifold for (24) tangent to at zero must be and the flow on this manifold is given by -dimensional ODE The nonlinear terms in (41) are in normal form in the classical sense with respect to matrix . In application, usually can be determined by the following procedure.

Theorem 8. For , let denote the operator defined in , with values in the same place, by with domain . Here, we use the natation to denote the space of homogeneous polynomials of degree in 2 variables , with coefficients in a Banach space .

According to [16, 17], we derive and then where is the projection of   on .

Since , it can be checked that . Then, we obtain and further then the normal form of (24) on the invariant local center manifold is given by If , by using the definitions, we can verify that

Theorem 9. If , the positive equilibrium is unstable for any .

For , is also an equilibrium with simple zero singularity. To discuss its stability, we employ the following perturbation form: where and is a small parameter. Then, the normal form of (24) near the bacteria-free equilibrium is as follows:

Theorem 10. If , the bacteria-free equilibrium is unstable for any .

5. Sensitive Analysis

Sensitivity indices allow us to measure the relative change in a variable when a parameter changes. The normalized forward sensitivity index of a variable to a parameter is the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differentiable function of the parameter, the sensitivity index may be alternatively defined using partial derivatives. Here, one adopts the following definition as described by Chitnis et al. [19].

Definition 11. The normalized forward sensitivity index of a variable, , that depends differentiably on a parameter, , is defined as .

To clear bacteria in the body, we must take measures to make hold. For this end, it is important to determine how crucial each parameter is to . Table 1 exhibits the analytical sensitivity indices of to the parameters , , , , , , and .

Table 1 

Sensitivity indices of to its parameters.

The negative sign of the sensitivity index for implies that increase in the relevant parameter leads to the decrease in . Note that . The most sensitive parameter is the growth rate , which has a high impact on and should be targeted by intervention strategies. To reduce to ensure , we may develop an inhibitor to control the reproduction of bacteria or to kill bacteria individual, which agrees with an intuitive expectation.

To understand how affects , we might as well assume , , and , and then gets . Clearly, for , only the innate immune cell competes with bacteria. Figure 1 exhibits the relationship between and with , , and . It can be seen from the plot that as decreases increases faster when than the case of , which reveals that vaccination or other strategies adopted to stimulate immunity of the body are beneficial to the clearance of bacteria.

Figure 1 

Plot of against the the effective reproductive rate for .

Figure 2 shows the solutions of system (2) starting from different initial values chosen arbitrarily together with , ,  , , , , , , , , , , , , and . It is easy to see that the bacteria concentration tends to zero, and the smaller growth rate leads to the faster convergence.

(a)

(b)

(c)

Figure 2 

Simulations of for in each figure.

6. Conclusions

This paper formulates the competition between bacteria and immune system by DDEs. Then, the qualitative properties of the model are analyzed. Specially, by virtue of the center manifold reduction and normal form theory due to Faria and Magalhaes [16, 17], the normal form of system (2) associated with zero eigenvalue is computed, from which one deduces that the bacteria-free equilibrium and the positive equilibrium are unstable under the conditions of and , respectively. Next, sensitivity analysis and numerical simulations indicate that the effective reproductive rate is the most sensitive parameter to . Theorem 3 suggests the strategies target the decrease of the growth rate which can be successful in disease elimination.

On the biological viewpoint, the terms and measure the respective strengths of the innate and adaptive immune system defense against the bacterial challenge, while the factor measures the bacteria’s offensive strength. So with , we can compare the strength of the immune system against the bacterial offensive. Thus, Theorems 2 and 3 have the biological explications: in the domain of attraction of , bacteria will be cleared if ; that is, the strength of the immune system defense against the bacteria challenge is not weaker than the bacteria’s offensive strength; in the domain of attraction of , bacteria coexist with immune cells when and the bacterial challenge is weaker than bacteria’s offensive strength.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the Natural Science Foundation of China (nos. 10971064, 11271125, and 11201277); the Mathematics Tianyuan Funds of NSFC (no. 11226260); the Scientific Research Plan Projects of Shaanxi Education Department (nos. 12JK0851, 2013JK0611); the Key Laboratory of Simulation and Control for Population Ecology (Xinyang Normal University), Xinyang 464000, China (no. 201004); the Program for Innovative Research Team of Science and Technology of University of the Henan Province (no. 2010IRTSTHN006); the Innovation Scientists and Technicians Troop Construction Projects of the Henan Province.