MODELING DISEASES WITH LATENCY AND RELAPSE

A general mathematical model for a disease with an exposed (latent) period and relapse is proposed. Such a model is appropriate for tuberculosis, including bovine tuberculosis in cattle and wildlife, and for herpes. For this model with a general probability of remaining in the exposed class, the basic reproduction number R0 is identified and its threshold property is discussed. In particular, the disease-free equilibrium is proved to be globally asymptotically stable if R0 < 1. If the probability of remaining in the exposed class is assumed to be negatively exponentially distributed, then R0 = 1 is a sharp threshold between disease extinction and endemic disease. A delay differential equation system is obtained if the probability function is assumed to be a step-function. For this system, the endemic equilibrium is locally asymptotically stable if R0 > 1, and the disease is shown to be uniformly persistent with the infective population size either approaching or oscillating about the endemic level. Numerical simulations (for parameters appropriate for bovine tuberculosis in cattle) withR0 > 1 indicate that solutions tend to this endemic state.

1. Introduction.In deterministic epidemic models, individuals are divided into a number of classes.Susceptible individuals infected with the disease but not yet infective are in the exposed (latent) class.These individuals pass into the infective class, and then on recovery (that may be natural or due to treatment) into the recovered class.For some diseases, recovered individuals may relapse with reactivation of latent infection and revert back to the infective class.This recurrence of disease is an important feature of some animal and human diseases, for example, tuberculosis, including human and bovine [1,2], and herpes [1,3].For human tuberculosis, incomplete treatment can lead to relapse, but relapse can also occur in patients who took a full course of treatment and were declared cured.Most tuberculosis in human adults (caused by Mycobacterium tuberculosis) in the USA results from reactivation of latent infection [1].In a clinical study in 1999 in Malawi, 206 P. VAN DEN DRIESSCHE, LIN WANG AND XINGFU ZOU 7.5% of patients registered with new tuberculosis were found to have had previous tuberculosis, with recurrence due to reactivation or reinfection [4].Tuberculosis patients infected with HIV are significantly more likely to relapse compared with patients uninfected with HIV [5].Specifically, it is reported that with HIV disease, individuals exposed to tuberculosis can reactivate as frequently as 10% per year, as compared with 10% per lifetime without HIV disease [6].
A model for herpes with a general relapse distribution, but ignoring the exposed class, is formulated in [7] and shown to exhibit a threshold phenomenon.No evidence of sustained oscillations is found, even for a constant relapse period, that is, a step-function relapse distribution.However, some epidemic models, formulated with step-function distributions leading to discrete delays, do exhibit sustained oscillations for some parameter values, for example, models of diseases giving temporary immunity with discrete delay in the recovered class but with no exposed class [8,9].It is thus of interest to investigate models with relapse and a general exposed distribution in order to determine whether sustained oscillations can occur.
In this paper, we formulate and analyze a model including a general exposed distribution and the possibility of relapse.This model was motivated by a formulation in [10, Section 5], which assumed a constant exposed period, for the spread of bovine tuberculosis (Mycobacterium bovis) in a cattle herd.Bovine tuberculosis is spread from animal to animal mainly by direct contact [11,12], although in some countries badgers and opossums may act as a reservoir host; see [13,14,15] and references therein.Bovine tuberculosis can be latent, may take months to develop to the infectious stage, and also can relapse.It is a chronic disease that produces lesions in lungs and causes emaciation and difficulty in breathing [12].Details of the transmission cycle of bovine tuberculosis remain unclear, and the disease is still spreading in many areas, for example, Michigan [12], Great Britain [16], New Zealand [17], and in some developing countries [18].Johne's disease in cattle, which is caused by Mycobacterium paratuberculosis, also causes animals to lose weight.Some animals appear to recover from clinical Johne's disease, but often relapse in the next stress period [19].Infection with Mycobacterium bovis in humans is rare in the USA, but is still a problem in areas where the disease in cattle is not controlled [1].
Our general model, which includes relapse in a constant population (for example, one cattle herd) and an arbitrarily distributed exposed stage, is formulated in the next section.A basic reproduction number R 0 is identified in Section 3, and the stability of the disease-free equilibrium is considered.Section 4 presents a complete analysis of the case in which the exposed period is exponentially distributed, and R 0 = 1 is proved to be a sharp threshold in the sense that the disease dies out if R 0 < 1, but approaches an endemic level if R 0 > 1.In Section 5, a constant exposed period is assumed, resulting in a delay differential equation system.The endemic equilibrium is shown to be locally stable for R 0 > 1. Persistence of the disease for R 0 > 1 is also proved in this section.In Section 6, numerical simulations, with parameters relevant for bovine tuberculosis giving R 0 > 1, indicate that solutions tend to the endemic equilibrium.We briefly summarize our results in Section 7.

2.
Formulation of the general model.Denote the size of the population by N (t).This is divided into four disjoint classes of individuals that are susceptible, exposed (latent), infective, and temporarily recovered, with class sizes denoted by S(t), E(t), I(t) and R(t), respectively.Recovery may be natural or due to treatment of infective individuals.Such treatment could entail antimicrobial drugs for tuberculosis or antiviral drugs for herpes.Depending on health, an individual in the recovered class can revert to the infective class with a constant rate α ≥ 0. For simplicity, it is assumed that the recruitment (or birth) rate and the removal (or death) rate are identical; this is denoted by b > 0. Let β > 0 be the transmission coefficient, that is, the average number of effective contacts of an infective individual per unit time, and γ ≥ 0 be the rate at which infective individuals recover.Thus the probabilities of remaining in the infective and recovered classes are assured to be exponentially distributed.It is assumed that individuals rarely die of the disease.This assumption, together with the assumption of identical recruitment and removal rates, ensures a constant population, thus N (t) = N .
Let P (t) denote the probability (without taking death into account) that an exposed individual still remains in the exposed class t time units after entering the exposed class.We assume throughout that P (t) satisfies the following reasonable properties: (A): P : [0, ∞) → [0, 1] is nonincreasing, piecewise continuous with possibly finitely many jumps and satisfies P (0+) = 1, lim t→∞ P (t) = 0 with ∞ 0 P (u)du positive and finite.Assuming that the force of infection is given by standard incidence and that initially a small number of infectives is introduced into an otherwise susceptible population and thus S(0) > 0, I(0) > 0, E(0) = R(0) = 0, the equations governing the SEIR model are Here and in the sequel, integrals are in the sense of Riemann-Stieltjes integrals.
The flow diagram for the model is depicted in Figure 1.Since the total population size N = S(t) + E(t) + I(t) + R(t) is a constant, it is convenient to work with proportions in all compartments.Rescaling gives Notice that Hence Therefore, our general model can be written as the system with S(0) + I(0) = 1.
To show basic properties of the model including well-posedness, let 3) with (1) (or ( 2)) has a unique solution with (S(t), I(t), R(t)) ∈ D and S(t) > 0 for all t ≥ 0.
Proof.The Volterra integro-differential equation system (3) with properties (A) satisfies the hypotheses stated by Miller [20, p338] that are sufficient to ensure the existence, uniqueness and continuity of solutions.It follows from the first equation in (3) that which implies that S(t) > 0 for all t > 0. To prove nonnegativity of I(t) and R(t), consider two cases: I(0) = 0 and I(0) > 0. If I(0) = 0, then S(0) = 1.By uniqueness, this implies that S(t) = 1, I(t) = 0 and R(t) = 0 for all t ≥ 0. Clearly the statement of the lemma is true.Assume that I(0) > 0. We claim both I(t) and R(t) are also nonnegative for t > 0. If not, then there must exist a finite time t 0 > 0 such that I(t) ≥ 0, R(t) ≥ 0 for t ∈ [0, t 0 ] and either (i) I(t 0 ) = 0 and I (t 0 ) ≤ 0, or (ii) R(t 0 ) = 0 and R (t 0 ) ≤ 0. But if (i) is the case, then from the second equation in (3), it follows that This implies that I(t 0 ) ≥ I(0)e −(γ+b)t 0 > 0, which shows that (i) is impossible.Next, assume that (ii) is true.As in (i), it follows that which is a contradiction.Therefore I(t) and R(t) are nonnegative for all t ≥ 0. It also follows from ( 1) that E(t) is nonnegative for t ≥ 0. The positive invariance of D now follows from the above and the fact that S(t) REMARK 1.It follows from the above proof that if (S(0), I(0), R(0)) ∈ D and I(0) > 0, then S(t) > 0, I(t) > 0, R(t) > 0 for all finite t > 0, and if  Define Here R 0 is the basic reproduction number [23], which is the product of the transmission coefficient β, with the fraction Q surviving the exposed class, and the death-adjusted mean time in I on multiple passes.This last term is given by the sum of the geometric series where 1  γ+b is the average time in the infective class on the first pass, γ γ+b is the probability of surviving this class, and α α+b is the probability of surviving the recovered class.Note that ∂R 0 ∂β , ∂R 0 ∂Q , ∂R 0 ∂α are all positive, but ∂R 0 ∂γ is negative.As expected, the parameter R 0 = 1 is a threshold value for the DFE, as shown in the following result.We use the notation: Proof.First we consider local stability of the DFE (1, 0, 0).Linearizing (3) at the DFE gives the characteristic equation as follows: where Since z = −b is a negative real root of the above equation, it suffices to show that all roots of h 1 (z) = 0 have negative real parts ([24, Chapter 5]).Suppose z = x+iy, where x ≥ 0, is a root of h 1 (z) = 0, and rewrite h 1 (z) = 0 as ( It follows from x ≥ 0 that which is equivalent to Denote the left-hand side of the above inequality by F 1 (x, y) and the right hand side by F 2 (x, y).Then F 1 (x, y) ≤ F 2 (x, y) with Notice that This contradicts (6).Therefore the real part x < 0, and this shows that the DFE is locally asymptotically stable for R 0 < 1.
Next, we show that the DFE is indeed global asymptotically stable if R 0 < 1.To this end, we only need to show it is globally attractive.Note that by Lemma 2.1, S(t), By the Fluctuation Lemma [25], there is a sequence t n with t n → ∞ as n → ∞ such that The second equation of (3) can be rewritten as It follows from the above equation that By the Lebesgue-Fatou Lemma [26, p468], this gives that From the third equation of (3), it follows that If I ∞ > 0, then the above inequality, together with (7), yields giving a contradiction.Therefore I ∞ = 0, implying R ∞ = 0 as well, and hence Applying the above and the theory of asymptotically autonomous systems (see Castillo-Chavez and Thieme [27]) to the first equation of (3), we conclude that S(t) → 1 as t → ∞.Therefore, the DFE is globally asymptotically stable if To show the DFE is unstable, it suffices to show that h 1 (z) = 0 admits a positive real root.Considering z = x > 0, it follows from (5) that Note that F 3 (x) is an increasing function of x with F 3 (0) = b(α+γ+b) α+b and F 3 (x) → ∞ as x → ∞, whereas F 4 (x) is decreasing in x with F 4 (0) = βQ.Thus R 0 > 1 implies that F 3 (0) < F 4 (0).Therefore there must be a positive x 0 such that F 3 (x 0 ) = F 4 (x 0 ) and this x 0 is a real positive root of h 1 (z).Hence the DFE is unstable if R 0 > 1.
To proceed further with the analysis, we assume particular forms for the probability of remaining in the exposed class, and investigate the dynamics.4. P (t) = e − t : an ODE system.Assume that the probability of remaining in the exposed class is negatively exponentially distributed with mean exposed time 1/ .Then P (t) = e − t with > 0 and (3) reduces by using (1) to the following system of ordinary differential equations (ODEs) In this case, Q = +b , and thus from (4) the basic reproduction number is Note that the expression for R 0 given in ( 9) can be rigorously derived by using the next generation matrix method [28] on the ODE system for E(t), I(t) and R(t) with E (t) = βS(t)I(t) − ( + b)E(t) from (1).
A complete description of the model behavior can be given in this case, with R 0 = 1 acting as a sharp threshold in the global sense.An endemic equilibrium (EE) is a steady state with disease present.THEOREM 4.1.Consider the ODE model (8).If R 0 < 1, then the DFE is globally asymptotically stable in D. If R 0 > 1, then the DFE is unstable and there is a unique EE, which is globally asymptotically stable in D \ {(1, 0, 0)}.
Proof.The global asymptotical stability of the DFE in the case that R 0 < 1 follows immediately from Theorem 3.1.This can also be proved by using the Lyapunov function defined by where k 1 = α α+b , k 2 = +b and k 3 = k 1 k 2 and by applying LaSalle's invariance principle [29,Theorem 6.4,p30].
Next we assume that R 0 > 1, then there exists a unique positive EE denoted by (S * , I * , R * ) with where R 0 is given by (9).Motivated by the Lyapunov functions and the use of the arithmetic mean-geometric mean inequality in [30,31], define Then the derivative of V 2 (t) along the solutions of ( 8) is given by
5. P (t) is a step-function: a DDE system.Assume that the probability of remaining in the exposed class P (t) is the step-function given by for finite τ .Thus all individuals remain in the exposed class for a constant period τ .For t ∈ [0, τ ], model (3) reduces to an ordinary differential equation system For t > τ , model (3) reduces to the following delay differential equation system with It is interesting to notice that once the initial condition (S(0), I(0), R(0)) ∈ D with S(0) + I(0) = 1 and E(0) = R(0) = 0 for (12) is given, the initial condition needed for the delay differential equation system ( 14) is indeed given by the solution of (12) for t ∈ [0, τ ].Note that the model given by ( 14) was formulated in [10, Section 5], where the existence of equilibria and local stability results are stated [10,Result 5.1].However, no further analysis on this model was written then, and to the best of our knowledge, no other analysis of this model has appeared in the literature.
We now consider the dynamics of (14).In this case, Q = e −bτ , and hence the basic reproduction number is THEOREM 5.1.Consider system (14).If R 0 < 1, then the DFE is globally asymptotically stable in D. If R 0 > 1, then there is a unique EE, which is locally asymptotically stable.
Next we obtain some global information about the disease in terms of persistence.Theorem 5.2 below shows that if R 0 > 1, then the disease does not go extinct, rather I(t) either approaches I * or oscillates about this value.
Proof.If the conclusion is not true, then either I ∞ = 0 or S ∞ = 1.Either case would lead to lim t→∞ I(t) = 0, lim t→∞ S(t) = 1 and lim t→∞ R(t) = 0. Thus lim t→∞ w(t) = 0, where w(t) > 0 is defined by Then, by Lemma A.25 [26, p432], we can find a sequence {t n } with t n → ∞, w(t n ) → 0 as n → ∞ and w (t n ) < 0 for all n.Since R 0 > 1, 1/R 0 < 1, there exists a positive integer K such that S(t n ) > 1/R 0 for n ≥ K.By ( 18) and (14), where η is a positive constant independent of the initial conditions.Proof.Since I(t) ≤ 1 for all t ≥ 0, the first equation of (14) gives that By the standard comparison theorem [21,Theorem B.1,p261], it follows that S(t) ≥ Ŝ(t), where Ŝ(t) is the solution of Applying the Fluctuation Lemma to the third equation of ( 14) yields Applying the Fluctuation Lemma to the second equation of ( 14), together with (19) and the fact that S * = 1/R 0 with R 0 given by ( 16), gives From Lemma 5.1 and the above inequality, it follows that S ∞ ≥ S * .By the Fluctuation Lemma as well as the first equation of ( 14), (20) indicates that X 2 is a uniform weak repeller for X 1 .Applying Theorem 1.4 of [38] to the solution semiflow Φ(t, φ) = (S t , I t , R t ), t ≥ τ, φ ∈ X 1 of ( 14), we conclude that X 2 is also a uniform strong repeller for X 1 , that is, the disease is uniformly strongly persistent, implying that there is an η > 0 such that I ∞ ≥ η with η being independent of the initial data.The Fluctuation Lemma on the second equation of ( 14) using ( 19) yields This inequality together with completing the proof of assertion (i).
Applying the Fluctuation Lemma to the first equation of ( 14) yields It then follows from (20) using (i) that Therefore, assertion (ii) holds.Assertion (iii) follows immediately from (19) and I ∞ ≥ η.
6. Numerical simulations.Theorem 5.1 gives local stability of the EE for R 0 > 1.We present some numerical simulations of (12) for t ∈ [0, τ ] and ( 14) for t > τ using Matlab.We take parameter values motivated by bovine tuberculosis in a cattle herd, with time unit of one year.It is assumed that the exposed period is 9 months (i.e., τ = 0.75).Taking b = 0.1, α = γ = 0.5 (equal periods of infection and recovery before relapse) and β = 0.6 (the same order of magnitude as reported by Barlow [14]) gives R 0 ≈ 3.04, which is in an estimated range (from 5 × 10 −1 to 2.5×10, see [10]) for bovine tuberculosis in a cattle herd.The numerical simulations shown in Figure 2   7. Discussion and Conclusions.We have proposed a general model to capture the features of some infectious disease with latency and relapse.These features are observed in diseases such as tuberculosis and herpes.
We briefly summarize our results for the disease transmission model formulated in Section 2 with the basic reproduction number identified in Section 3. If R 0 < 1, then global attractivity of the DFE is proved by using the Fluctuation Lemma.If the probability of remaining in the exposed class is assumed to be negatively exponentially distributed as in (8), then for R 0 > 1 the EE is proved (using a Lyapunov function) to be globally asymptotically stable.If this probability of remaining in the exposed class is a step-function as in (14), then for R 0 > 1, the disease persists in the population.The EE is locally asymptotically stable with I(t) either converging to the endemic level I * or oscillating about I * .No evidence of sustained oscillations about I * is found in our numerical simulations.These suggest that the EE is actually globally asymptotically stable for R 0 > 1, but we do not have an analytic proof.
To control the disease (e.g., herpes, bovine tuberculosis), a strategy should reduce the reproduction number to below one.One possibility for control is vaccination (see, for example, [18,39]), which has the effect of reducing β.Vaccination may also delay the onset of infection (i.e., increase the exposed period), which in turn decreases the reproduction number.Another possibility is treatment with drugs; see for example [2,3].However, such treatment may not be cost-effective for cattle with bovine tuberculosis.In addition, previously treated cattle may experience recrudescence of active infection; see [2] and references therein.As noted in the introduction, relapse is also reported in humans who have been treated for tuberculosis [4], especially those who also have HIV [5,6].Feng et al. [40] consider a model for tuberculosis with treatment of exposed and infective individuals, albeit with no relapse.They find that an arbitrarily distributed exposed period does not alter the qualitative dynamics of the disease, which either dies out or remains endemic depending on the value of the reproduction number.The inclusion of vaccination or treatment of exposed individuals in our model with relapse remains to be explored.It also would be interesting and worthwhile to set up a model to investigate the effect of HIV on the relapse rate in tuberculosis.

Figure 1 .
Figure 1.The flow diagram of the SEIRI model.
(we plotted the I component only) indicate that solutions tend to the EE with I * ≈ 0.34.Other simulations with different parameter values yielding R 0 > 1 show similar convergence to the EE, with no evidence of oscillations about I * .