GLOBAL DYNAMICS OF A STAGED PROGRESSION MODEL FOR INFECTIOUS DISEASES

We analyze a mathematical model for infectious diseases that progress through distinct stages within infected hosts. An example of such a disease is AIDS, which results from HIV infection. For a general n-stage stage-progression (SP) model with bilinear incidences, we prove that the global dynamics are completely determined by the basic reproduction number R0. If R0 ≤ 1, then the disease-free equilibrium P0 is globally asymptotically stable and the disease always dies out. If R0 > 1, P0 is unstable, and a unique endemic equilibrium P ∗ is globally asymptotically stable, and the disease persists at the endemic equilibrium. The basic reproduction numbers for the SP model with density dependent incidence forms are also discussed.

1. Introduction.For diseases that progress through a long infectious period, infectivity or infectiousness can vary greatly in time.The progression of a typical HIV infection can take eight to ten years before the clinical syndrome (AIDS) occurs, and the progression goes through several distinct stages, marked by drastically different CD4 + T-cell counts and viral RNA levels.HIV-infected individuals are highly infectious in the first few weeks after infection, then remain in an asymptotic stage of low infectiousness for many years and become gradually more infectious as the immune system becomes compromised and they progress to AIDS.Variability of infectiousness in time has been modelled in the literature by Markov chain models or staged-progression (SP) models (see, e.g.[1]- [10]).A deterministic SP model was proposed and analyzed in [4].To formulate an SP model, the total host population is partitioned into the following compartments: the susceptible compartment S; the infectious compartment I i , whose members are in the i-th stage of the disease progression, where i = 1, 2, • • • , n; and the terminal compartment A. Let δ i be the mean progression rate from the i-th stage to the (i + 1)-th stage, for i = 1, 2, • • • , n − 1, and δ n the mean progression rate from the n-th stage to the terminal stage of the disease.We assume that hosts in the terminal compartment are non-infectious due to inactivity.In the case of AIDS, the terminal compartment consists of people with active AIDS.AIDS patients typically either become sexually inactive or isolated from the infection process, and their infectivity is negligible.We also assume that there is no recovery from the disease, and thus the only exit from compartment A is death.Let λ i be the transmission coefficient for the infection of a susceptible from an infectious in the class I i , which takes into account average number of contact and probability of infection for each contact; then the total incidence is given by n i=1 λ i I i Sf (N ), where N = S + I 1 + • • • + I n is the total active population.Here we assume that the density dependence of the incidence is given by a function f (N ), which will be specified later.A class of special interest is f (N ) = N −α , 0 ≤ α ≤ 1, the resulting incidence term includes two of the most common incidence forms: the standard incidence form (α = 1) and the bilinear incidence (α = 0).The average death rate for the susceptible compartment is d 0 ; for the compartment I i it is d i , which may include death due to infection; and for the terminal compartment it is d A .We assume the inflow of susceptibles is a constant Λ.The population transfer among compartments is schematically depicted in the transfer diagram in Figure 1.All parameters in the model are assumed to be positive.
Based on our assumptions and the transfer diagram, the following system of differential equations can be derived for the n-stage SP model: and The incidence form is λ S, where the force of infection is density dependent.We assume that the function f (N ) satisfies the following assumptions, for N > 0, The assumptions that f (N ) > 0 and f (N ) ≤ 0 are biologically motivated.As the total population N increases, the probability of a contact with a susceptible decreases, and thus the force of infection is expected to be a decreasing function of N. The other two conditions we impose on f are needed for our analysis.It can be verified that the class f (N ) = N −α , 0 ≤ α ≤ 1, satisfies (H).This class contains the standard incidence (α = 1) and the bilinear incidence (α = 0).
Adding the equations in (1) we obtain where d = min{d 0 , d which can be verified to be positively invariant with respect to (1).We show in Section 4 that the basic reproduction number for model ( 1) is given by We prove that the global dynamics of model ( 1) are completely determined by the basic reproduction number R 0 .If R 0 ≤ 1, then the disease-free equilibrium P 0 = (Λ/d 0 , 0, • • • , 0) is globally asymptotically stable in the feasible region Γ, and the disease always dies out.If R 0 > 1, then P 0 becomes unstable, and system (1) is uniformly persistent.Moreover, a unique endemic equilibrium P * = (S * , I * 1 , • • • , I * n ) exists in the interior of Γ.For the case of bilinear incidence, namely, when f (N ) = 1, we prove that P * is globally asymptotically stable in the interior of Γ if R 0 > 1.
In the next section, we prove the uniqueness of the endemic equilibrium when R 0 > 1.In Section 3, global stability of P 0 is established for R 0 ≤ 1. Concrete expressions of R 0 are derived in Section 4. The global stability of P * for the case f (N ) = 1 is proved in Section 5.

2.
Equilibria.An equilibrium (S, where λ is given in (2).The disease-free equilibrium P 0 = (Λ/d 0 , 0, • • • , 0) exists for all positive parameter values.Next we consider the existence of endemic equilibria Then the following properties of the matrix A follow from the definition and properties of M -matrices given in the Appendix.
Proposition 2.1.The following holds for the matrix A defined above.
(2) −A −1 exists and is a non-negative matrix.
By Proposition 2.1, we know that Define the basic reproduction number of (1) as In Section 4, we will derive the explicit expression of R 0 and give its justification.
We have the following result on the number of equilibria.
Proof.The last n equations of ( 4) can be written in the form Multiplying the row vector (λ (9) Also, by (8), where, by Proposition 2.1, From the first equation of (4) we get which, together with (10), implies and thus Substitute ( 12) into ( 9), and we obtain the equation for an endemic equilibrium (S, We will show that equation ( 13) has a unique positive solution in the interval (0, Λ/d 0 ) when R 0 > 1.Let Then g(0) = 0, and Furthermore, by our assumption (H) on function f (N ), where N = pΛ + (1 − pd 0 )S.Thus the function y = g(S) is strictly monotonically increasing, and its graph has at most one intersection with the line y = 1/β.Such an intersection exists for S ∈ (0, Λ/d 0 ) if and only if This completes the proof of Theorem 2.1.

3.
Stability of the disease-free equilibrium P 0 .The following stability result can be derived from a more general result in [14] for a model with arbitrarily distributed infectious stages.The proof we presented here provides a motivation for our definition of R 0 .
Proof.The last n equations in (1) can be rewritten as where A is given as in (5).Multiplying a row vector (c 1 , c 2 , • • • , c n ) to the above equation, we obtain For the choice of c i in ( 14), define a Lyapunov function Then, using assumption (H) we obtain, along a solution of (1), Furthermore, L = 0 if and only if either (a) and S = Λ/d 0 are satisfied.In either case, the largest compact invariant subset of the set is the singleton {P 0 }.To see this, let K be the largest compact invariant subset of G.In case (a), each solution in K satisfies equation S (t) = Λ − d 0 S, and the only solution that is bounded for Therefore, all solutions in Γ converge to P 0 , by LaSalle's Invariance Principle (see [11]).The global attractivity of P 0 and the Lyapunov function L imply that P 0 is also locally stable, since otherwise P 0 will have a homoclinic orbit that is entirely contained in G, contradicting that the largest compact invariant set in G is {P 0 }.This establishes the global stability of P 0 when R 0 ≤ 1.
If R 0 > 1, L > 0 for S sufficiently close to Λ/d 0 , and thus solutions in Γ sufficiently close to P 0 move away from P 0 , except those on the invariant S-axis, along which solutions converge to P 0 .Therefore, P 0 is unstable.Furthermore, {P 0 } is the only compact invariant set on the boundary of Γ and is isolated.The local dynamics near P 0 and the boundary of Γ imply that system (1) is uniformly persistent (see [12]) with respect to Γ, when R 0 > 1.The proof of uniform persistence of ( 1) is similar to that of Proposition 3.3 in [13].
4. The basic reproduction number R 0 .Theorems 2.1 and 3.1 establish R 0 as a sharp threshold parameter.If R 0 ≤ 1, the disease dies out irrespective of the initial number of cases.If R 0 > 1, then the disease persists in the feasible region and there is a unique endemic equilibrium.Such a role of threshold parameter is expected of the basic reproduction number, the average number of infections caused by a single infective in a population at the disease-free equilibrium ( [15,16,17,18,19,20]).It is then reasonable to regard the parameter R 0 defined in (7) as the basic reproduction number.
Our derivation of R 0 is based on the stability analysis of the disease-free equilibrium P 0 using a Lyapunov function, as in [19].Other methods of deriving R 0 exist in the literature; among them are the method of second generation matrix in [17], which was later modified in [20], and the derivation based on the linear stability analysis of P 0 (see [18]).In the following, we show that the method of next generation matrix as formulated in [20] for deriving the basic reproduction number leads to the same R 0 as given in (7).
Set y = (I 1 , • • • , I n , S) T .Then model ( 1) can be written as where At the disease-free equilibrium in the new coordinates, P0 = (0, 0, where and g(N ) = N f (N ).Moreover, Here V n×n = A, where A is defined in (5).Hence, the next generation matrix (see [20]) where c i are defined in ( 14), and thus c 1 = β.The basic reproduction number is defined in [20] as the spectral radius, ρ(F V −1 ), of the matrix F V −1 .It is easy to see that which gives the same R 0 as in (7).Next, we derive the explicit expression of β defined in ( 6) and then derive R 0 according to (7).Let I n×n be the identity matrix.Solving the matrix equation −AX = I n×n , we obtain that X = −A −1 is a lower triangular matrix where c ij 's are determined by the following iterative relations: For the class f (N ) = N −α , 0 ≤ α ≤ 1, we have the corresponding basic reproduction number If α = 1, then R 0,1 gives the same expression for the basic reproduction number for the standard incidence as in [4,20], while if α = 0, R 0,0 gives the basic reproduction number for the bilinear incidence.We also note from (17) that R 0,α decreases as α increases, if Λ/d 0 > 1; the bilinear incidence gives the largest value of R 0,α while the standard incidence gives the smallest value (see Figure 2). 5. Stability of the endemic equilibrium P * .In this section, for f (N ) ≡ 1, we prove the global stability of the endemic equilibrium P * when R 0 > 1.
Theorem 5.1.Assume that f (N ) ≡ 1 and R 0 > 1.Then the endemic equilibrium P * is asymptotically stable.Furthermore, all solutions in the interior of Γ converge to P * .
Proof.Rewrite the equilibrium equations ( 4) for which leads to and , and b i > 0 is defined inductively in the following way: We note that W (x) ≥ 0, for x ∈ Int Γ, the interior of Γ, and W (x) = 0 ⇐⇒ x = x * .So function W is positive definite with respect to the endemic equilibrium x * = P * .Solving the inductive relation (21), we obtain or equivalently Computing the derivative of W along solutions of system (1), we obtain Using (1) we have In the second step of the above derivation, we substituted Λ by the right-hand side of the first equation of (18).Similarly, using (1) and ( 18), we obtain It is straightforward to verify that the last two terms of the above expression are 0, using the inductive relation (21 To see that the first term is zero, we observe by (19).Using ( 25)-(28), we can calculate and simplify (24) as From (23) we have By (22), (20), and exchange of order of summation, we have by the inequality Therefore, along any solution (S(t), I 1 (t), I 2 (t), is the endemic equilibrium P * .Therefore the directional derivative of W in the direction of the vector field of (1) is negative definite in Int Γ with respect to the endemic equilibrium P * .This implies that the basin of attraction of P * contains Int Γ.The positive definiteness of W (x) with respect to P * implies that P * is also locally stable.This completes the proof of Theorem 5.1.
Remark 1. Lyapunov functions of the type used in the proof of Theorem 5.1 have been used in the literature, for ecological models (e.g., see [21]) and more recently for epidemic models (e.g., see [22]).