GLOBAL DYNAMICS OF A COUPLED EPIDEMIC MODEL

. In this paper, we propose a novel epidemic model coupling direct and indirect transmission of disease and study the global dynamic of the model system. Despite the nonlinearity and complexity of the system, the basic reproduction number exhibits a nice linear property: it is simply the sum of two basic reproduction numbers for direct and indirect disease transmissions respectively. We further demonstrate that the local and global dynamics of the system are related to the basic reproduction number. The new model has the advantage that it generalizes or connects to various disease models on HIV, Zika virus, avian inﬂuenza, H1N1 and so on.


1.
Introduction. We propose and study the following three-dimensional disease model system: where x and y denote the susceptible and infected population sizes, respectively. The constants µ and γ are death rates of these two groups, and b is a constant birth rate. The disease could transmit directly via a standard incidence function β 1 xy/(x+y) and indirectly via a mass-action infection term β 2 xz, where z accounts for the vector of the indirect transmission. We assume in the third equation that the growth of disease vector is proportional to the number of infected individuals, and its decay rate is a constant δ. The units of parameters β 1 , µ, γ, p, δ are the same as the reciprocal of time unit, but the unit for b is the unit of β 1 multiplied by the population size/density unit, while the unit of β 2 is the unit of β 1 divided by the population size/density unit. If β 1 = 0 (namely, there is no direct transmission), the above system (1-3) reduces to the classical HIV model [17], where x, y, z represent the densities of target cells, infected cells and virus, respectively. It is well-known that the global dynamics of this reduced system is fully determined by the basic reproduction number [5,21] To be more specific, if R i 0 ≤ 1, then the infection-free equilibrium is globally asymptotically stable; if R i 0 > 1, then there exists a unique positive equilibrium which is globally asymptotically stable. Here, the superscript i means indirect disease transmission.
On the other hand, if the indirect transmission is ignored (i.e., β 2 = 0), then the equation (3) can be decoupled from the original system, and the remaining system (1)(2) is the same as the Kermack-McKendrick epidemic model [8] with standard incidence function, while the infected class is assumed to be removed from the social activity after being recovered/quarantined; see [22]. The basic reproduction number for the reduced direct-transmission-only model is given by [5,21] where the superscript d corresponds to the direct disease transmission. A similar dichotomy result holds: if R d 0 ≤ 1, then the disease-free equilibrium is globally asymptotically stable; if R d 0 > 1, then there exists a unique positive equilibrium which is globally asymptotically stable. For global stability analysis on epidemic models with more general direct transmission rates, we refer to [9,10,15] and references therein.
The proposed system (1-3) is motivated by the observations of cell-to-cell transmission of HIV [20] and human-to-human transmission of Zika virus [7]. It is of both biological and mathematical interests to study a disease model that couples direct transmission with indirect transmission. The coupled virus models with both transmissions being assumed to be mass-actions have been studied in [13,14,18], where the basic reproduction number for the coupled system is simply the sum of two basic reproduction numbers for the subsystems with only direct or indirect disease transmissions respectively. Similar results were also obtained for cholera models [19], pathogen models [2,3] and treatment models [24]. The coupled model proposed in [1] incorporated two virions with sensitive and resistant strains, respectively, and the basic reproduction number for the full system is the maximum of two basic reproduction numbers of viral strains. However, the basic reproduction number corresponds to each viral strain is still the sum of two basic reproduction numbers for the direct and indirect transmissions when both transmission functions are chosen as mass-actions. One nature question is whether this linear property of basic reproduction number of nonlinear coupled systems is universal, regardless of the choices of nonlinear incidence functions. In our system (1-3), we assume that the contact probability between a susceptible individual and an infected individual is decreasing as total population of two groups is increasing [4,11], which leads to a standard incidence function for the direct disease transmission. On the other hand, the contact probability between a susceptible individual and a disease vector is assumed to be a constant, so the indirect disease transmission is characterized by a mass-action incidence (or bilinear function). Biologically, we may still define the basic reproduction number of this coupled system as the sum of R i 0 and R d 0 : We intend to explore the relation of local and global dynamics for the model system with the basic reproduction number.
2. Main results. We first establish the nonnegativeness and boundedness of the solutions to the system (1-3) with given nonnegative initial conditions. Proof. Note from (1) that x (t 0 ) > 0 whenever x(t 0 ) = 0 for some t = t 0 . This implies that x(t) > 0 for all t > 0 if x(0) ≥ 0. Now, we claim y(t) and z(t) are nonnegative for all t ≥ 0. If not, then there exists a t 0 ≥ 0 such that y(t) ≥ 0 and z(t) ≥ 0 for all t ≤ t 0 and either y(t) or z(t) becomes negative right after t 0 . On one hand, if y(t) < 0 for t close to the right of t 0 , then y(t 0 ) = 0 and y (t 0 ) < 0. However, we obtain from (2) that y (t 0 ) = β 2 x(t 0 )z(t 0 ) ≥ 0, a contradiction. On the other hand, if z(t) < 0 for t close to the right of t 0 , then z(t 0 ) = 0 and z (t 0 ) < 0. However, we obtain from (3) that z (t 0 ) = py(t 0 ) ≥ 0, a contradiction. Therefore, neither y(t) nor z(t) could be negative for any t ≥ 0. Next, we show that the solutions are ultimately bounded as t → ∞. Adding equations (1) and (2) gives Thus, we have lim sup t→∞ [x(t) + y(t)] ≤ b/ min(µ, γ). Especially, x(t) and y(t) are ultimately bounded as t → ∞. It is then easily seen from (3) that z(t) is also ultimately bounded as t → ∞. Now, we consider the equilibria which satisfy the following algebraic system: It is readily seen that there is a disease-free equilibrium (x 0 , 0, 0) with x 0 = b/µ. Now, we are looking for the endemic equilibrium and assume y = 0. Eliminating the variables y and z from equations (7)-(9) gives z = py/δ, y = (b − µx)/γ and Recall that γ and µ denote the death rates of infected and susceptible individuals, respectively. Biologically, it is nature to assume that γ ≥ µ. In this case, we always have a unique positive solution x = x * to the above equation, where , γ > µ.
We further have y = y * and z = z * , where To make y * > 0 and z * > 0, we need x * < b/µ, which, by a simple calculation, is equivalent with R 0 > 1, where R 0 is defined in (6). The equivalence is obvious when γ = µ. When γ > µ, we rewrite x * < b/µ as Isolating the square root and squaring both sides gives B + 2Ab/µ > 0 and Simplifying the last inequality yields Cµ 2 < Bbµ + Ab 2 . We make use of the expressions of A, B, C and rewrite the inequality as By expanding and canceling, we obtain γδµ < bpβ 2 + β 1 δµ, which is the same as R 0 > 1. Furthermore, coupling R 0 > 1 and γ > µ implies B + 2Ab/µ > 0. Thus, we obtain that the endemic equilibrium exists if and only if R 0 > 1. Next, we will establish the local stability theory of disease-free and endemic equilibria.
then the disease-free equilibrium is locally asymptotically stable. If R 0 > 1, then the disease-free equilibrium is unstable. If further, γ ≥ µ, then there exists a unique endemic equilibrium and it is locally asymptotically stable.
Proof. We calculate the Jacobian matrix of the system (1-3) about the diseasefree equilibrium (x 0 , 0, 0) and endemic equilibrium (x * , y * , z * ), respectively. At the disease-free equilibrium, we have The three eigenvalues of the Jacobian matrix are calculated as −µ and (T r ± √ T r 2 − 4Det)/2, where T r and Det denote the trace and determinant of the submatrix That is, we have T r < 0 and Det > 0, and thus all three eigenvalues have negative real parts, which implies that the disease-free equilibrium is locally asymptotically stable. On the other hand, if R 0 > 1, then Det < 0, which implies that at least one of the three eigenvalues (i.e., (T r − √ T r 2 − 4Det)/2) is real and positive. Consequently, the disease-free equilibrium is unstable.
For the critical case R 0 = 1, the Jacobian matrix J 0 has two negative real eigenvalue −µ and β 1 − γ − δ, and one zero eigenvalue. We introduce the matrix of eigenvectors The Jacobian matrix for the differential equations of (u, v, w) about the zero equilibrium is exactly Λ. To analyze the local asymptotic stability of this zero equilibrium, we need to calculate the restricted dynamical system on the center manifold for u sufficiently small and v = O(u 2 ), w = O(u 2 ); see [23,Ch. 2]. Note that u = (δy + β 2 x 0 z)/(δ 2 + pβ 2 x 0 ). We obtain from equations (2-3) that Since the above restricted system is stable about u = 0, the original system (2-3) is stable about the disease-free equilibrium (x 0 , 0, 0). Finally, we investigate the Jacobian matrix at the endemic equilibrium: The corresponding characteristic equation is λ 3 + c 2 λ 2 + c 1 λ 1 + c 0 = 0, where Since the endemic equilibrium (x * , y * , z * ) satisfies the system (2-3), we can rewrite the above expressions as Especially, c 2 > 0, c 1 > 0, c 0 > 0, and c 2 c 1 > c 0 . By Routh-Hurwitz criterion, it follows that all eigenvalues of J * have negative real parts, which implies the local asymptotic stability of endemic equilibrium. Now, we are ready to establish the global dynamics of the system (1-3).
x * ] ≥ 0. We will prove global asymptotic stability of endemic equilibrium for each of the above three cases by constructing difference Lyapunov functions.
where K 1 is a constant such that V 1 (x * , y * , z * ) = 0. Taking the derivative of V 1 along the solution of system (1-3), we obtain Recall that the endemic equilibrium satisfies the following equations: By using x ≤ x + y, we obtain Since 3 − a − b − c ≤ 0 for any positive a, b, c such that abc = 1, it follows that The largest invariant set on which V 1 (t) = 0 is a singleton (x * , y * , z * ). We obtain from Lyapunov-LaSalle invariance principle [12, p. 30] that the endemic equilibrium is globally asymptotically stable.
Case II. y * −x * ≥ 0. We use V 1 defined in (15) to construct the Lyapunov function: where K 2 is a constant such that V 2 (x * , y * , z * ) = 0. Taking derivative along the solution of system (1-3), we obtain where we have used the equations (x + y) = b − µx − γy and b = µx * + γy * . Note that the coefficients of (x − x * )(y − y * ) in the numerator cancel out with each other. It then follows that The largest invariant set on which V 2 (t) = 0 is a singleton (x * , y * , z * ). We obtain from Lyapunov-LaSalle invariance principle [12, p. 30] that the endemic equilibrium is globally asymptotically stable.
Case III. 0 ≤ x * − y * ≤ [2µx * + 2µγ(x * + y * )/β 1 ]/(µ + γ) and (15) to construct the Lyapunov function: where and K 3 is a constant such that V 3 (x * , y * , z * ) = 0. Taking the derivative along the solution of system (1-3), we obtain where X := x − x * and Y := y − y * . To prove V 3 ≤ 0, it is sufficient to show that which is equivalent with From the choice of λ, we can rewrite the above inequality as A simple calculation shows that the above inequality is the same as the above inequality is satisfied. Again, the largest invariant set on which V 3 (t) = 0 is a singleton (x * , y * , z * ). We have the global asymptotic stability of endemic equilibrium from Lyapunov-LaSalle invariance principle [12, p. 30].
3. Discussion. In this paper, we propose a simple epidemic model coupling both direct and indirect transmission mechanisms of infectious diseases. This model has potential applications in the study of various diseases. For example, it generalizes the HIV model where x, y, z correspond to uninfected cells, infected cells and virus, respectively. In the study of Zika virus, we let x, y, z be the uninfected individuals, infected individuals, and infected mosquitoes, respectively. When we apply our model to analyze the epidemic waves of H1N1 and seasonal influenzas, x and y still denote uninfected and infected individuals, respectively, while z stands for the contaminated environment such as classrooms, buses, theaters, or other public places. We could also use our model to study the cross transmission of avian influenza among migratory birds and domestic poultry, where x, y, z denote the uninfected migratory birds, infected migratory birds, and infected domestic poultry, respectively. We would like to mention that the transmission mechanisms within and without groups may not be the same, so we assume standard incidence function for within group transmission (such as bird-to-bird transmission, or person-to-person transmission) and mass action function for transmission among different groups (such as poultry-to-bird transmission, or mosquito-to-person transmission). One may argue that there should be poultry-to-poultry transmission and bird-to-poultry transmission, so that equation (3) would have been written as where w(t) is the population of uninfected poultry. However, for simplicity, we could assume that this quantity is a large constant w(t) ≡ W with respect to the infected poultry z(t). Therefore, the above equation can be approximated by (3) with p = β 4 W and δ = σ−β 3 . A similar argument works for mosquito-borne disease models where x, y, z, w correspond to the uninfected individuals, infected individuals, infected mosquitoes, and uninfected mosquitoes, respectively. When the dynamical system on mosquito population is much faster than the disease transmission system, the population of uninfected mosquitoes reaches its ecological equilibrium w(t) ≡ W very quickly and thus its variation is negligible in the slow system. For simplicity, we further assume that the ratio of infected mosquitoes over uninfected mosquitoes is relatively small. These assumptions help us to obtain a reduced equation for the infected mosquitoes: z (t) = py(t) − δz(t). Nevertheless, a more realistic and more general model incorporating with the changes on the uninfected mosquitoes should be studied and we leave it as a future work. It is noted that the basic reproduction number of the coupled system is simply the sum of the basic reproduction numbers for the two reduced systems with only direct or indirect disease transmission, respectively. This means that, in the study of infectious disease (Zika virus, for instance) with multiple transmission mechanics, it is dangerous to ignore either direct or indirect transmission because that would underestimate the seriousness of the disease: it is possible that both R d 0 and R i 0 are less than one, but their sum exceeds the critical value one. To create a suitable policy for disease control, one should not only reduce the possibility of direct transmission, but also avoid indirect disease transmission (such as closing schools to reduce vector transmission channels of H1N1 and other pandemic influenzas [6]).
Due to the difference between direct and indirect transmission mechanisms, it becomes more difficult to analyze the model system than the one in [13] where both transmissions were chosen as mass action functions. As seen in the proof of our main theorems, the most challenging part is the construction of a suitable Lyapunov function which could nicely balance the two different nonlinear transmissions. Note that we introduce some technical conditions to prove global asymptotic stability of endemic equilibrium; see (C1)-(C3) in the proof of Theorem 2.3. From numerical simulations one could still observe global asymptotic stability of endemic equilibrium even though these technical assumptions are violated. It remains an open problem to prove global asymptotic stability of endemic equilibrium whenever it exists. Another future project would be the extension of our results to the more general situation when disease transmission delay and multi-group [16] are taken into consideration.