Analysis of Lauffenburger-Kennedy bacterial infection model for tissue inflammation dynamics

ABSTRACT In this paper, we analyze a mathematical model for an inflammatory response to bacterial infection of homogeneous tissues. Specifically, we provide a detailed analysis of the Lauffenburger-Kennedy bacterial infection model and show that the model exhibits three possible equilibria corresponding to a bacteria-free and two endemic compromised steady states. Asymptotic results of the steady states along with the existences of saddle-node connection Hopf bifurcations are shown under certain conditions of the parameters. Within the biological ranges of the parameter values, we observe that the system can exhibit both forward and backward bifurcation. In addition, in both cases, the larger compromise bacterial infection steady state can either approach an equilibrium or can oscillate around it via Hopf bifurcation depending on the value of the ratio of leukocyte mortality to phagocytosis rates. Numerical results are used to provide illustrative examples of these different dynamical patterns observed in the model.


Introduction
In recent decades, mathematical models investigating dynamical patterns of inflammatory response to bacterial invasion of tissue have been considered [2,[4][5][6][7][8][14][15][16]18,23,[25][26][27]29]. We consider the Lauffenburger-Kennedy model describing bacterial (B) and leukocyte (L) densities in a localized tissue region [16]: where h 0 represents the normal emigration coefficient, h 1 denotes the increase per unit local bacterial density-i.e. enhanced emigration coefficient, μ is the mortality rate of leukocyte, and c b is the localized leukocyte density circulating in the vasculature. The surface area of venule wall is denoted by A and the tissue volume is represented by V. The bacterial growth rate is denoted by K g , K i is the inhibition constant, K d is maximum per capita death rate of leukocyte, and the saturation rate is given by K b . The first term of leukocyte density in Equation (1) represents the combined recruitment or emigration rate, and the second expression is the mortality rate of leukocyte. In this expression, the death rate is described as an exponential decay and emigration assumes a random process where production of both vasoactive substances (e.g. histamine) and chemical factors enhancing adherence of circulating leukocyte to blood vessel wall increase leukocyte emigration.
In the bacteria population density equation, the first term denotes the growth rate, which increases exponentially for small bacterial density and saturates as density increases. The second term represents the bacterial death due to leukocyte phagocytosis (i.e. phagocytosis rate). Model (1) describes the tissue inflammation dynamics.
In this paper, we analytically investigate the following non-dimensionalized bacterial infection model for tissue inflammation dynamics: where u = L/L 0 and v = B/K i are normalized leukocyte and bacterial density. Here L 0 = h 0 c b A/μV and ρ = μK i /L 0 K d is the ratio of leukocyte death rate to maximum phagocytosis rate, τ = h 1 K i /h 0 is a dimensionless ratio of enhanced phagocyte emigration rate to normal emigration rate and measures the sensitivity of leukocyte tissue infiltration rate during inflammation to the local bacterial density, β = K g K i /K d L 0 is the ratio of maximum bacterial growth rate to maximum leukocyte phagocytosis rate, and k = K b /K i represents the ratio of inhibitory effect of increasing bacteria density on bacterial growth to inhibitory effect of phagocytosis rate, respectively. Heret = K d L 0 t/K i , where in Equation (2) we have dropped the˜int. System (2) describes the inflammatory response to a bacterial invasion of tissue. Table 1 summarizes the parameters and state variables, where the parameter values were obtained from [16] and refs. therein. The dimensionless parameters are then calculated using the relationships between the parameters as outlined above. The values of these parameters play an important role in the analysis section and are used in the numerical simulation to illustrate different dynamics behaviours, where parameter values are selected within these ranges to show these different regimes.
There have been numerous works studying the above system to model bacterial infection. In particular, in [15,18], the authors used a simple mathematical model to describe the dynamical of local tissue inflammatory response to bacterial invasion, where spatially distributed conservation equations for tissue bacterial and leukocyte cell density were considered as a Metchnikoff-type predator-prey system [21]. Similarly, in [2], the authors developed a system of two coupled nonlinear ordinary differential equations, and use scaling arguments and singular perturbation techniques to study the dependency of the dynamic features of the inflammatory response on certain key parameters. Furthermore, Jing [11] showed the existence of bifurcation of saddle connection and the creation of limit cycles under certain conditions of the parameters. In [17], the author extended the model in [16] by incorporating diffusion and chemotaxis. Using a linear stability analysis, the author showed that when a phagocyte chemotactic response is smaller than a critical value the system can lead to a non-uniform state. However, if the chemotactic response is greater than this critical value, the system stabilizes to the uniform state. In [5], the author further extended the model constructed in [17] by incorporating bacterial taxis and provided similar steady states and stability analyses. In addition, the author discussed the possibility of having a non-uniform steady state, due to the combination of both chemotaxis and diffusion. Numerical simulations were used to check the effect of bacterial chemotaxis on the system's ability to exhibit Turing regimes. Assuming that the diffusion rate and the chemotactic rate are both very small compared to the growth rate, the authors in [22] derived a new equation to describe the time-evolution of the aggregating region and showed the conditions for the existence and stability of radially symmetric equilibrium solutions of the equation. Here, our focus is to revisit the model in Equation (2) to provide detailed analytic results underlying the observations in [16].
In the proceeding sections, we revisit the system described by Equation (2) with particular emphasis on complete and rigorous mathematical results of the stability of the steady states and on the bifurcation results and orbital stability of periodic solutions. In Section 2, we consider positivity and boundedness of the solutions and construct the feasible region where the model is biologically relevant. In Section 3, we establish the conditions for the existence of multiple endemic equilibria. In addition, the asymptotic stability conditions are provided for these equilibria in Section 4, and furthermore, bifurcation results and the existence of at least one limit cycles are provided. Section 5 provides illustrations of the various numerical simulations of the dynamic patterns predicted by our analytical results. Section 5 summarizes these results and provides a further research direction.

Positive invariant feasible region
The model (2) describes dynamics of bacterial infection for tissue inflammation, particularly the interaction between bacteria and phagocyte; and it is therefore critical to show that the model behaviours lie in a feasible region for which the model is biologically relevant. We need to show that the solutions to Equation (2) with natural initial non-negative conditions (u(0), v(0)) = (u 0 , v 0 ) will be confined in a biologically meaningful region (i.e. the solution will remain non-negative) for all time t ∈ R + , where R + := [0, ∞). We have the following results: Theorem 2.1: Solutions to system (2) with non-negative biologically relevant initial conditions u(0) ≥ 0 and v(0) ≥ 0 are non-negative for all time t ∈ R + . That is, solutions u(t) and v(t) remain in the set Furthermore, the following holds: then for all t ∈ [0, T 0 ], the expression for v(t) becomes Therefore, v(t) ≥ v(0) e At . Then evaluating v(t) at time t = T 0 , we observe that v(T 0 ) ≥ v(0) e AT 0 > 0, which is a contradiction. This implies v(t) > 0 for all time t ∈ R + . Similarly, since v(t) ≥ 0 for all t ∈ R + , we observe that, Hence and v(t) are non-negative for all t ∈ R + . That is, u(t) and v(t) remain in the region This completes the proof.

Steady states analysis
In this section, we consider the steady state solutions of the model described by Equation (2) augmented with the following initial conditions: We observe that the system given by Equation (2) with the above initial data admits a semitrivial equilibrium given by E 0 = (u, v) = (1, 0); and we call this steady state the infectionfree (elimination) equilibrium, where bacteria is absent or eradicated from the tissue. It should be noted that the infection-free equilibrium exists. In addition, we observe that Equation (2) can also allow non-trivial steady states or endemic (compromise) equilibria, depending on the model parameters. Following [5,16] and solving we can obtain a quadratic expression, which has two solutions where the system (2) has a positive steady state (u, v) if and only if and u = 1 + τ v > 0. Clearly, we can directly calculate the exact expression using quadratic formula. However, because our aim herein is to obtain the necessary conditions for the various dynamical patterns, we instead use Equation (10) to obtain the necessary conditions for the existence of endemic equilibria in order to simplify the calculations. In particular, we observe that h(v) has the following properties: and h(v) attains its minimal positive value h(λ 0 ) at v = λ 0 .
Noting that h(0) = 1/k and from the aforementioned properties of h(v), we state the following results on the existence of a positive equilibrium of the system (2). (10) and λ 0 be defined by (11), then the following conclusions hold:
The proof of Lemma 3.1 is trivial and will be omitted herein. From Lemma 3.1, the steady states are determined by only k, β and τ . In addition, we observe that the system (2) can undergo both forward (see Lemma 3.1(1)) and backward (see Lemma 3.1(2)) steady state bifurcations. In the proceeding discussion, we will consider and establish the necessary conditions for the stability of these equilibria. Throughout this paper, we denote the infection-free and the two endemic equilibria by E 0 , E λ and E μ , respectively.

Infection-free equilibrium, E 0
Examining the characters of the eigenvalues of the Jacobian matrix obtained from the linearization of Equation (2) at the infection-free equilibrium, we state the following conclusions: Proof: The linearized operator evaluated at (1, 0) is given by Therefore, the eigenvalues are given by −ρ and −1/k + β. Hence, if R 0 = βk < 1, then the infection-free equilibrium E 0 is locally asymptotically stable; and it is otherwise unstable if R 0 = βk > 1. In addition, at R 0 = βk = 1, the eigenvalues of system (2) linearized at the infection-free equilibrium are −ρ and 0.
To establish results under the condition of R 0 = βk = 1, we use the center manifold methods ( [9]), and reformulate the problem using the abstract setting. In particular, we first translate the infection-free equilibrium E 0 to the origin, whereũ = u − 1,ṽ = v. For convenience, we drop the˜inũ andṽ, and use u and v, respectively. Hence, the local system (2) becomes From the center manifolds theory of [9], there exists a center manifold for system (13) which can locally be represented as follows for δ > 0 sufficiently small. To compute W c (0), we assume that h(v) has the form Then, the equation for the center manifold is given by which can be rewritten as follows, Hence, Because the coefficient of v in the left-hand side of Equation (18) is Then the system (2) restricted to the center manifold is, for v sufficiently small, given by the following expression It is easy to check that when k ≥ 1, then the zero solution of (19) is stable on one hand. If 0 < k < 1, then the zero solution of (19) is unstable on the other hand. Therefore, from the center manifold theorem of [9], it follows that if 1/β = k ≥ 1, then E 0 is stable in the original system (2), while unstable when 1/β = k ∈ (0, 1). This completes the proof.
It follows from Lemma 3.1 and Theorem 4.1 immediately that the following result holds.

Endemic equilibria, E μ and E λ
In the previous section, we establish that, when R 0 < 1, the infection-free equilibrium is locally asymptotically stable. Moreover, when R 0 > 1, the infection-free equilibrium is unstable, and there is at least one non-trivial solution. In this section, we consider the stability of endemic equilibria.

Proof:
To study the existence of Hopf bifurcation occurring at the unique equilibrium point E λ = (u * , v * ) = (u λ , λ) = (1 + τ λ, λ), we choose λ as the bifurcation parameter, and letλ be defined asλ Then when λ =λ, the Jacobian matrix L(λ) has a pair of imaginary eigenvalues σ = ±i h (λ)ρλ/(1 +λ). Let σ = ψ(λ) ± iω(λ) be the roots of the characteristic polynomial P λ (σ ) = σ 2 − σ T + D = 0. Then, and By the application of Poincaré-Andronv-Hopf bifurcation theorem [9], we note that the system (2) undergoes a Hopf bifurcation at E λ when λ =λ. However, the detailed nature of the Hopf bifurcation needs further analysis of the normal form of the system. To that end, we translate the equilibrium For the sake of convenience, we drop˜inũ andṽ and use u and v, respectively. After some calculations, the local system (2) becomes Rewriting Equation (43), we have ⎛ where and In addition, let us define the matrix P as follows, with the following inverse where and When λ =λ, we havē (51) Using the following transformation Where L(λ) is defined as follows Besides, and The coefficients in F 2 (x, y, λ) are given by 3 , Rewriting Equation (53) in polar coordinates, we have the following expressionṡ r = ψ(λ)r + a(λ)r 3 + higher-order-terms, θ = ω(λ) + c(λ)r 2 + higher-order-terms.
Then, using Taylor expansion of Equation (58) at λ =λ, we havė Using [9,28], we determine the stability of the periodic solution. In particular, we determine the sign of a(λ), where a(λ) is given by We note that all the partial derivatives of F are evaluated at (x, y, λ) = (0, 0,λ). It is easy to Hence, This completes the proof.
The proof of the above corollary is easy and has been ignored. The Corollary 4.6(1) follows from the use of inequalities and comparison technique, where we can use the fact that k ≥ 1 and R 0 ≤ 1, which imply that β ≤ 1. Corollary 4.6(2) is essential the same as Theorem 4.5.
In addition, for completeness, we also discuss the boundedness of the solution of system (2) and the existence of limit cycles in , and obtain the following theorems. Proof: We consider the following auxiliary system By Theorem 4.5, it follows that the solution (ũ(t),ṽ(t)) of system (70) converges to (u λ , λ) globally. Because 0 < k < 1, we have Hence, by comparison principle, we have This together with the global stability of (u λ , λ) indicates that (u(t), v(t)) of system (2) is bounded. This completes the proof. Proof: To show that the system does not have periodic orbits or limit cycles, we observe that, provided that β − ρ ≤ 0. Hence, there is no closed orbit in . This completes the proof.

Discussion
In this paper, we have considered a mathematical model describing the bacterial infection of a homogeneous tissue. Our analysis of the model has shown the existence of a bacteriafree and two endemic compromised steady states. By using steady state analysis and Hopf bifurcation theory, we show the existence of saddle-node connections and limit cycles. In this section, we discuss our analytic results and its biological implication.
In the case where the ratio of inhibition effects of increased bacterial density on its growth to that of phagocytosis is greater or equal to one (i.e. k ≥ 1), the model exhibit two possible situations. When R 0 = βk ≤ 1, the bacteria-free equilibrium is globally asymptotically stable. That is, whenever R 0 ≤ 1, the elimination of bacterial infection is always sustained after initial transient proliferation of infection. The infection-free steady state is unstable when R 0 > 1, and one of the endemic compromise steady state (i.e. the endemic compromise steady state E λ ) emerges and becomes globally asymptotically stable. The bacteria proliferates gradually before stabilizing to a positive steady state. Figure 1 shows the two possible dynamical features present in the model for k ≥ 1. This case corresponds to the existence of forward bifurcation. In particular, when R 0 ≤ 1, the system approaches the bacteria-free equilibrium which is globally asymptotically stable (see Figure 1(a)); when R 0 > 1, the bacteria-free equilibrium is unstable and model converges to the unique endemic equilibrium (see Figure 1(b)).
For k < 1, the dynamical features are more involved. In particular, we observed that when k ≥ 1/(1 + τ ), the model exhibits forward bifurcation with unique endemic equilibrium and when k < 1/(1 + τ ) the system can have backward bifurcation where there are two endemic equilibria. Unlike their counterparts in standard epidemic models, these  steady state bifurcations structures showed quite interesting dynamical behaviours. For example, in the case of forward, the unique endemic equilibrium can lose stability or can be globally asymptotically stable depending on other secondary conditions on the parameters; and similar results are observed in the case of backward bifurcation.
For k < 1/(1 + τ ) < 1, we have two endemic steady states namely, E λ and E μ . In addition, the endemic equilibrium E μ is always unstable, and the stability of E λ depends on the parameters. In particular, E λ is locally asymptotically stable for λ ∈ (λ * * , ∞) if ρ > A(λ * * ). For R 0 < R c , the bacteria-free equilibrium is globally stable; and the endemic equilibrium E λ is globally asymptotically stable and the bacteria-free E 0 is unstable for R 0 > 1. For R c < R 0 < 1, the situation is different. That is, both E λ and E 0 are locally asymptotically stable. The solution of the model approaches different steady state depending on the initial conditions, and there is no oscillatory behaviours. The dynamical features are similar to the dynamics seen in Figures 1, 2 and 3 (a, b).
For ρ < A(λ * * ), there are three scenarios for the dynamics of the steady states. For the first case, when R 0 < R c < 1, only the bacteria-free equilibrium exists and it is globally stable; and for R c < R 0 < 1 (and R c < R 2 < R 1 < 1 and λ 0 > λ 2 > λ 1 ), both endemic E λ and the bacteria-free E 0 equilibria are stable, and the initial conditions determine which equilibrium is selected. When R 0 > 1, E λ is globally asymptotically stable, and E 0 is unstable; and there is no Hopf bifurcation (see similar behaviours are shown in Figures 1-2). For the second case, we have that R c < R 1 < 1 < R 2 and λ 0 < λ 1 < λ 2 . When R 0 < R c < 1, E 0 is again globally asymptotically stable, and unstable when R 0 > 1. For R c < R 0 < R 1 < 1, both E 0 and E λ are locally asymptotically stable and the initial conditions again determine which steady-state is selected. When R c < R 1 < R 0 < 1 and R c < R 1 < 1 < R 0 < R 2 , both E 0 and E λ are unstable; and E 0 remains unstable for R 0 > 1. For R 0 > R 2 > 1, endemic E λ is globally asymptotically stable. At R 0 = R 1 and R 0 = R 2 , the endemic equilibrium E λ undergoes Hopf bifurcation (see similar dynamical behaviours are shown in Figure 3(c)). The results are equivalent to condition in Theorem 4.3 when λ 0 < λ 1 < λ 2 where at λ = λ 1 and λ = λ 2 the system (2) undergoes Hopf bifurcations (see similar dynamical behaviours are shown in Figure 3(c)). For the third and final scenario, we note that λ 1 < λ 0 < λ 2 and R c < R 1 < 1 < R 2 . When R c < R 1 < R 0 < 1, the bacteria-free is again globally asymptotically stable and unstable if R 0 > 1; and the system (2) undergoes Hopf bifurcations only at λ = λ 2 and R 0 = R 2 (see similar dynamical behaviours are shown in Figure 3(c)). For R 0 > R 2 > 1, the endemic equilibrium E λ is globally asymptotically stable.
In previous discussion, we have focused mostly on the dynamic behaviours of the model under different conditions, which mostly depend on the values of four key dimensionless parameters ρ, β, τ and k, where ρ is the ratio of leukocyte death rate to maximum phagocytosis rate, τ represents ratio of enhanced phagocyte emigration rate to normal emigration rate, β denotes the ratio of maximum bacterial growth rate to maximum leukocyte phagocytosis rate, and k is ratio of inhibitory effect of increasing bacteria density on bacterial growth to inhibitory effect of phagocytosis rate. On one hand, suppose that scale measuring the inhibitory effects of bacterial density on its growth relative to phagocytosis rate is greater or equal to one (k ≥ 1). Then, we observe that only two possible outcomes can be seen, in which the bacterial infection is eliminated or is permanent. For instance, when the maximum bacterial growth rate is smaller relative to the maximum leukocyte phagocytosis rate (e.g. β = 0.2) such that R 0 < 1 as shown in Figure 1, the bacterial infection cannot persist, and leukocytes lead to elimination of these infections. However, when the maximum bacterial growth rate is greater than the maximum leukocyte phagocytosis rate (e.g. β = 2), the leukocytes are unable to effective eliminate the bacterial infection, and therefore, there is a persistence of endemic compromised disease state and the reproductive success of the infection is greater than one (i.e. R 0 > 1).
On the other hand, when the inhibitory effects of bacterial density on its growth relative to phagocytosis rate is less than one (k < 1), the system exhibits two types of steady-state behaviours, whose dynamic features are dependent on whether 1/(1 + τ ) < k < 1 or k < 1/(1 + τ ) < 1. In particular, when 1/(1 + τ ) < k < 1, and the ratio of leukocyte death rate to phagocytosis rate (ρ) is greater than a threshold condition (i.e. ρ > A(λ * * ) := ρ * ), we observe similar dynamic features as in the case of Figure 1 (see Figure 2), where the bacterial infection is always eliminated or persists, depending on the reproduction success R 0 . For R 0 < 1, the bacterial infection is always eliminated, and for R 0 > 1, the compromise state is permanent. However, when the ratio of leukocyte death rate to phagocytosis rate (ρ) is less than a threshold condition (e.g. ρ < ρ * ), bacterial infection is always eliminated whenever R 0 < 1. But, when R 0 > 1, the disease is always permanent. In addition, there exists regions, depending on the initial bacterial density, where the infection always approaches a steady state or oscillates around the endemic compromise solution.
For k < 1/(1 + τ ) < 1, the system exhibits backward bifurcation where the smaller endemic compromise steady state is always unstable. The stability of the larger endemic equilibrium depends on various parameter thresholds. In particular, when the ratio of leukocyte death rate to phagocytosis rate (ρ) is greater than a threshold condition (i.e. ρ > A(λ * * ) := ρ * ), we observe standard backward bifurcation where the larger endemic compromise steady state is stable where the bacterial infection persists for all initial condition whenever R 0 > 1. For R 0 < 1, the bacterial can either be eliminated or can persist, where reducing R 0 below one is no longer effective. When ρ < ρ * , the larger compromise endemic steady state can approach to an equilibrium or oscillate around the equilibrium for some parameter values (see Figure 3). All these behaviours are seen within the biological ranges of the parameters obtained in the literature.

Conclusion
In this paper, we provide a rigorous mathematical analysis of a model describing the bacterial infection of a homogeneous tissue. In particular, we presented steady-state analysis of all the equilibria, namely bacterial infection-free and two endemic compromised steady states. We showed a threshold-like condition for the existence and stability of each of the equilibria, where the bacteria-free equilibrium is globally asymptotically stable whenever k ≥ 1, when R 0 = βk ≤ 1, and is unstable when R 0 > 1. Similar results were established for both of the endemic compromise equilibria, where one of the endemic equilibria is always unstable, and the dynamic property of the second endemic equilibrium is very rich, where different dynamic behaviours can be observed. In particular, we showed that the model exhibits rich dynamic features, including the existence of bubble bifurcation and Hopf bifurcation. Moreover, we showed the existence of saddle-node connections and limit cycles. In this paper, we did not show the uniqueness of the limit cycle, and we leave this analysis as an open problem. Using parameter values obtained from literature, we simulated various dynamic behaviours predicted by our analysis.