Diseases with chronic stage in a population with varying size

Mathematical Biosciences 182 (2003) 1–25 www.elsevier.com/locate/mbs

Diseases with chronic stage in a population with varying size Maia Martcheva

a,*

, Carlos Castillo-Chavez

b

a

b

Department of Mathematics, Polytechnic University, Brooklyn, NY 11201, USA Department of Biological Statistics and Computational Biology, Cornell University, Ithaca, NY 14853, USA Received 6 July 2002; received in revised form 13 September 2002; accepted 24 September 2002

Abstract An epidemiological model of hepatitis C with a chronic infectious stage and variable population size is introduced. A non-structured baseline ODE model which supports exponential solutions is discussed. The normalized version where the unknown functions are the proportions of the susceptible, infected, and chronic individuals in the total population is analyzed. It is shown that sustained oscillations are not possible and the endemic proportions either approach the disease-free or an endemic equilibrium. The expanded model incorporates the chronic age of the individuals. Partial analysis of this age-structured model is carried out. The global asymptotic stability of the infection-free state is established as well as local asymptotic stability of the endemic non-uniform steady state distribution under some additional conditions. A numerical method for the chronic-age-structured model is introduced. It is shown that this numerical scheme is consistent and convergent of first order. Simulations based on the numerical method suggest that in the structured case the endemic equilibrium may be unstable and sustained oscillations are possible. Closer look at the reproduction number reveals that treatment strategies directed towards speeding up the transition from acute to chronic stage in effect contribute to the eradication of the disease. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Chronic stage; Variable infectivity; Disease-age structure; Hepatitis C; Variable population size; Difference scheme; Numerical methods; Sustained oscillations

1. Introduction The impact of a chronic stage on the disease transmission and behavior in an exponentially growing or decaying population is the focus of this paper. The framework is applied to the case of *

Corresponding author. Present address: Department of Biological Statistics and Computational Biology, 434 Warren Hall, Cornell University, Ithaca, NY 14853-7801, USA. Tel.: +1-718 260 3294; fax: +1-718 260 3660. E-mail addresses: [email protected] (M. Martcheva), [email protected] (C. Castillo-Chavez). 0025-5564/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 2 5 - 5 5 6 4 ( 0 2 ) 0 0 1 8 4 - 0

2

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

hepatitis C, a disease characterized by a long chronic stage. Hepatitis C, priorly referred to as Ônon-A, non-BÕ hepatitis, is a viral infection of the liver which was first recognized as a separate disease in 1975. The causative agent was identified in 1989. The pathogen responsible for hepatitis C virus (HCV) is quite different from the causative agents for hepatitis A and B. It is an enveloped RNA virus in the flaviviridae family [1]. Researchers have classified more than 100 strains of the virus. Hepatitis C has been identified as the most common cause of post-transfusion hepatitis worldwide. It accounts for
90% of this disease in Japan, USA and Western Europe. HCV is responsible for both acute and chronic cases of hepatitis C. Data from the World Health Organization suggest that up to 3% of the world population has been infected with HCV. It is expected that more than 170 million people are chronic carriers and at risk of developing liver cirrhosis and/or liver cancer. In the United States, the annual number of acute HCV infections has declined during the past decade from 180 000 to 35 000, but an approximately 3.9 million Americans (1.8% of the population of the US) are currently infected with HCV, and an estimated 8000–10 000 deaths each year are attributed to HCV-associated chronic liver disease [2]. Hepatitis C is transmitted mainly through blood transfusions. However, in 1990 an effective screening test was developed and the transmission has begun to decrease, at least in the developed countries. HCV can also be contracted via contacts with inadequately sterilized or unsterilized equipment. Needle-sharing among drug-users is also a risk factor. Two major risk groups are the health care personnel (needle-stick injuries) as well as patients undergoing hemodialysis. Evidence of sexual and perinatal transmission exist, although these ways of transmission seem to occur less frequently and their relevancy and effectiveness are often questioned. Ear and body piercing, circumcision, tattooing and other percutaneous procedures may also be a source of HCV infections. The majority of patients with acute hepatitis C develop a chronic infection which is characterized by detection of HCV RNA for a period of at least six months after a newly acquired infection. Data suggest that at least 85% of patients with acute HCV infection eventually progress to the chronic stage. Chronic hepatitis C can lead to cirrhosis and end-stage liver disease (HCC). Cirrhosis can develop either within 1–2 years or after 2–3 decades following exposure. Some long term studies suggest that cirrhosis develops in 20–30% of patients. A fraction (1–5%) of cirrhotic individuals develop cancer of the liver in about a decade after developing cirrhosis. The disease has extremely slow average rate of progression, ranging from 2 to 4 or even more decades. Using data collected in Japan, investigators estimated that, following acute infection, chronic hepatitis could be identified 13:7  10:9 years later, chronic active hepatitis 18:4  11:2 years later, cirrhosis 20:6  10:1 years later, and HCC, 28:3  11:5 years later [3]. The most common symptoms of acute hepatitis C are fatigue and jaundice. However, the vast majority of cases (up to 90%), including those with chronic disease, are asymptomatic. This makes the diagnosis of hepatitis C very difficult and is the reason why the HCV epidemic is often called Ôthe silent epidemicÕ [4]. The enzyme immunosorbent assay is used for the detection of specific antibodies which confirm the presence of HCV. The presence of HCV RNA in serum indicates the presence of active infection and a potential for transmission of the infection and/or the development of chronic liver disease. Hepatitis C can be treated with interferon. The treatment is effective in about 20% of the cases. Ribavirin in combination with interferon is effective as an antiviral agent against HCV in approximately 40% of previously untreated patients. The combination therapy has been approved by

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

3

FDA in 1998. However, it also appears that the cost of such combined treatment is high – about $20 000 for 48-week treatment. No vaccine is available for hepatitis C. The high mutability of hepatitis C genome [5] complicates its development. There is no evidence that the successful treatment of HCV gives any kind of partial or temporary immunity. Hence the models developed fall within the class of models that treated or recovered individuals move back to the susceptible class. Several studies exist treating epidemics in populations with active demographic patterns, particularly in the ODE case, [6–10]. The epidemic contact process or the force of infection is typically assumed to be modeled by proportionate mixing. Derrick and van den Driessche [11] discuss a general SIRS disease transmission model in which the population varies and the rate at which susceptibles become infective is given by a rather general non-linear function. Another model with general contact function is treated in [12] where the effective contact rate is an arbitrary function of the total population size. The rate of transmission used here allows for the possibility of generation of secondary infections via contacts with two types of infectious individuals: those with acute and those with chronic infections. Very little analysis of diseases with chronic stage has been performed. In fact, the only two works known to the authors are [13,14]. A model of Cytomegalovirus infection with disease-age structured chronic stage is considered in [13]. A model structured by age-since-infection has also been discussed in relation to HIV in [15,16]. Reade et al. consider an ODE model for infections with acute and chronic phases with feline calicivirus [14]. Their work is mostly numerical and concentrates on the impact of vaccination on the acute and chronic phases. A mathematical model of hepatitis C has been considered in [17]. The authors use a backcalculation approach to reconstruct the history of HCV infection. Their parameters are estimated with data from France. The model is structured by age and sex since there is evidence that progression to cirrhosis depends strongly on these two characteristics. The model traces the hepatitis C epidemic in France back to the 1940s and makes projections until the year 2025. This paper is structured as follows. Section 2 presents and analyzes the baseline ODE model of hepatitis C infection. Section 3 extends the model of Section 2 to the case where the chronic class is structured by age-since-infection. The basic reproduction number is computed and used to show that the disease-free equilibrium is locally stable and to establish the existence of a unique endemic state. Section 4 discusses the global stability of the disease-free equilibrium and the local stability of the unique endemic equilibrium. In Section 5 we present a numerical method and discuss the results of our simulation. Section 6 summarizes our findings and conclusions. The Appendix A contains the proof of Theorem 5.2.

2. A simple ODE model The asymptomatic nature of hepatitis C and its slow progression make it difficult to characterize the natural history of the disease. The following assumptions are used in the construction of the models. (1) Only the acute and chronic stages are differentiated. Patients with either acute or chronic infections are capable of transmitting the disease. (2) All infected individuals develop an acute hepatitis C first.

4

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

(3) Some individuals with an acute infection progress towards the chronic state and later on develop cirrhosis. (4) Since the disease-induced death rate is relatively low, it is ignored. (5) The acute stage of infection is short and often asymptomatic and there is no possibility for treatment during this state. Since the disease progression is slow the model should also incorporate demographic changes in the population. It is assumed that the total population is growing exponentially ([18], Fig. 1.8). We divide the population in three subclasses: S––susceptible; I––infected with acute hepatitis C; V ––infected with chronic hepatitis C (with or without cirrhosis); N ¼ S þ I þ V (total number of people). The model introduces the following rates of transition from one class to another (parameters): b––birth/recruitment rate; l––death rate; a––treatment/recovery rate for the chronic state; k––rate of progression to chronic stage; c––effective contact rate of individuals with acute hepatitis C; d–– effective contact rate of individuals with chronic hepatitis C. All parameters are assumed positive. The model is homogeneous (of degree one). 8 < S 0 ¼ bN  ðcI þ dV Þ NS  lS þ aV : ð2:1Þ I 0 ¼ ðcI þ dV Þ S  ðl þ kÞI : V 0 ¼ kI  ðl þNaÞV The demographic equation for the dynamics of the total population size is given by: N 0 ¼ bN  lN: We set r ¼ b  l. We obtain N 0 ¼ rN and therefore, N ¼ N0 ert . Hence, r gives the growth rate of the population. If r > 0, that is b > l, the population grows exponentially, if r < 0, that is b < l, the population decreases exponentially. The case r ¼ 0 or b ¼ l implies that the population is stationary. These thresholds are often interpreted in terms of the demographic reproduction number. In particular, in [19] a net reproduction number of the population is thus defined as b R¼ : l Clearly, if R > 1 the population is exponentially increasing, if R < 1 the population is exponentially decreasing, and if R ¼ 1, the population is stationary. Since the model (2.1) is homogeneous of degree one, we consider also the equations for the normalized quantities. Setting s ¼ S=N , i ¼ I=N , v ¼ V =N leads to the following equivalent non-homogeneous system 8 < s0 ¼ bð1  sÞ  ðci þ dvÞs þ av i0 ¼ ðci þ dvÞs  ðb þ kÞi ð2:2Þ : v0 ¼ ki  ðb þ aÞv where s þ i þ v ¼ 1:

ð2:3Þ

In a sequence of papers Hadeler [20–22] develops a theory of ODEs with homogeneous nonlinearities. These models, in general, do not have time independent solutions. Instead, they support solutions of the form S ¼ ek t S , I ¼ ek t I , V ¼ ek t V . The normalized model (2.2), on the other hand, is non-linear, non-homogeneous and it has steady states in the classical sense. Hadeler calls the exponential solutions of (2.1) stable (unstable), if the corresponding steady states of (2.2) are

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

5

stable (unstable) in the classical sense. Under these conditions investigating model (2.1) is to a large extent reduced to investigating model (2.2). The disease-free state of (2.2) is (1,0,0) and the epidemiological reproduction number for the proportions is c kd þ : ð2:4Þ R0 ¼ b þ k ðb þ kÞðb þ aÞ The first term c=ðb þ kÞ can be interpreted as the contribution to the reproduction number due to secondary infections generated by an infective with acute hepatitis C. Naturally, it increases as the effective contact rate, c, of individuals with acute hepatitis C increases. The second term kd=ðb þ kÞðb þ aÞ is the contribution to the reproduction number due to secondary infections generated by an infective with chronic hepatitis C. It increases with the increase of effective contact rate of chronic individuals, d, as well as with the rate of progression to chronic stage, k. Overall, the reproduction number R0 has a more complicated response to variations of the rate of progression to chronic stage, k. In particular, it decreases when bðd  cÞ  ac < 0 and increases when the opposite inequality is valid. However, the probability of transmitting the disease from an individual with acute infection is larger than that from an individual with chronic infection, that is, d < c. Therefore, we expect that for realistic values of the parameters R0 will decrease as the rate progression to chronic stage increases. Lemma 2.1. If R0 < 1 then the disease-free equilibrium (1,0,0) is globally stable. If R0 > 1 the disease-free equilibrium is unstable and there is a unique endemic equilibrium ðs ; i ; v Þ which is locally stable. Proof. The Jacobian of (2.2) is 0 b c Jð1; 0; 0Þ ¼ @ 0 c  ðb þ kÞ 0 k

1 d þ a A: d ðb þ aÞ

This is a block-triangular matrix. It can be seen that one of the eigenvalues is b. The other two have negative real part if and only if c  ðb þ kÞ  ðb þ aÞ < 0 and ½c  ðb þ kÞ½ðb þ aÞ  kd > 0. The second inequality is equivalent to R0 < 1 (in fact, this can be used for the computation of R0 ). The condition R0 < 1 guarantees both inequalities. If R0 > 1 the second inequality is not satisfied and the disease-free equilibrium is unstable. In this case an endemic equilibrium exists and is given by ðs ; i ; v Þ, where 1 ; s ¼ R0  ðb þ aÞ 1 i ¼ ; 1 ð2:5Þ bþaþk R0 k 1 : v ¼ 1 bþaþk R0 Since s < 1, i < 1, v < 1 the quantities in (2.5) are well defined. To see the global stability of the disease-free equilibrium, we integrate the second and third equation in (2.2) to obtain the following pair of inequalities:

6

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

iðtÞ 6 e

ðbþkÞt

i0 þ

Z

t

eðbþkÞs ½ciðt  sÞ þ dvðt  sÞ ds; 0

vðtÞ 6 e

ðbþaÞt

v0 þ k

Z

t

eðbþaÞs iðt  sÞ ds:

0

We have lim supt i < 1 and lim supt v < 1. Taking the lim sup in the above inequalities leads to c d lim sup i þ lim sup v; bþk bþk t t t k lim sup v 6 lim sup i: bþa t t lim sup i 6

From here we obtain lim sup i 6 R0 lim sup i t

t

and hence, lim supt i ¼ 0 and lim supt v ¼ 0 if R0 < 1. Let R0 > 1. To see the local stability of the endemic state we consider the equations for i and v only. We use s ¼ 1  i  v to eliminate s. We obtain  0 i ðtÞ ¼ ½c  ðb þ kÞi þ dv  ðc þ dÞiv  ci2  dv2 ; ð2:6Þ v0 ðtÞ ¼ ki  ðb þ aÞv: The Jacobian at the endemic state is  c  ðb þ kÞ  ðc þ dÞv  2ci Jði ; v Þ ¼ k

d  ðc þ dÞi  2dv : ðb þ aÞ

Since the endemic state satisfies the system obtained from (2.6) by putting the derivatives equal to zero, the (1,1) entry in the Jacobian can be written as c  ðb þ kÞ  ðc þ dÞv  2ci ¼ 

dk ð1  v Þ  ci : ðb þ aÞ

It can be easily seen that the trace of the Jacobian is negative and the determinant is positive. Therefore, ðs ; i ; v Þ is locally asymptotically stable.  System (2.2) is two dimensional (see (2.6)) and direct application of the DulacÕs criterion is possible. We define the relevant region D ¼ fi P 0; v P 0; i þ v 6 1g: Proposition 2.2. The system (2.2) has no periodic solutions, homoclinic loops, or oriented phase polygons inside the region D. Proof. We begin by denoting F ði; vÞ ¼ ½c  ðb þ kÞi þ dv  ðc þ dÞiv  ci2  dv2 ; Gði; vÞ ¼ ki  ðb þ aÞv:

ð2:7Þ

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

7

As a Dulac multiplier we use 1=i. We have oF ði; vÞ=i oGði; vÞ=i dv ðb þ aÞ þ ¼  2 ð1  vÞ  c  <0 oi ov i i and therefore there are no closed orbits in the region D.



Since from Lemma 2.1 it follows that the endemic equilibrium is locally stable, then the result in Proposition 2.2 rules out sustained oscillations and guarantees that it is actually globally stable. The endemic proportions might approach the endemic equilibrium or (1,0,0) and in either case the total population size may be increasing, decreasing or constant, depending on the growth rate r. Similarly the behavior of the proportions does not give us much insight on the behavior of the total number IðtÞ of infecteds. In particular, if the number of infecteds IðtÞ increases exponentially but at a lower rate than the total population size, then the proportion of the two tends to zero, that is, the growth of the population outperforms the growth of the number of infecteds and, consequently, they become more and more diluted. We can also imagine the reversed situation––both the number of infecteds and the total population decline to zero but the proportion remains almost constant (and non-zero) at all times. In this case the disease remains in the population as long as the population exists. To describe the behavior of IðtÞ one needs an additional threshold parameter. Following the terminology introduced by Thieme [12] we define the basic replacement ratio: 8 c dx > > þ R0 < 1 < lþk lþk R1 ¼ cs kds > > þ R0 > 1 : l þ k ðl þ kÞðb þ aÞ where x is the unique positive solution of the equation dx2 þ ðc þ a  kÞx  k ¼ 0: The validity of the basic replacement ratio is established in Lemma 2.3. Lemma 2.3 also gives the exponential rate of increase (decrease) for SðtÞ, IðtÞ, V ðtÞ. Lemma 2.3 (1) The number of susceptible individuals SðtÞ increases (decreases) with the total population. In particular, the exponential rate of increase (decrease) is lim

t!1

S0 ¼ b  l: S

ð2:8Þ

(2) If R0 > 1 the number of individuals with acute hepatitis C, IðtÞ, decreases if R1 < 1 and increases if R1 > 1. Moreover, the exponential rate of increase (decrease) is lim

t!1

I0 ¼ ðl þ kÞðR1  1Þ ¼ b  l: I

ð2:9Þ

8

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

The exponential rate of increase (decrease) of the chronic class is equal to the exponential rate of increase (decrease) of the total population, that is, lim

t!1

V0 ¼ b  l: V

ð2:10Þ

(3) If R0 < 1 the number of individuals with acute hepatitis C, IðtÞ, and the number of chronic individuals, V ðtÞ, decrease if R1 < 1 and increase if R1 > 1. Moreover, the exponential rates of increase (decrease) are I0 lim ¼ ðl þ kÞðR1  1Þ; t!1 I ð2:11Þ V0 lim ¼ ðl þ kÞðR1  1Þ: t!1 V Proof (1) The cases R0 > 1 and R0 < 1 must be considered separately to establish the limit in (2.8). The case R0 < 1 is straight forward. To see (2.8) in the case R0 > 1 we divide by S in the first equation in (2.1) and take the limit:  S0 b i v v 1 lim ¼  c þ d s þ a  l ¼ ðb  ðci þ dv Þs þ av  ls Þ; t!1 S s s s s s where s , i and v are given by (2.5). From the first equation in (2.2) it follows that the equilibrium satisfies b  ðci þ dv Þs þ av ¼ bs : Thus limit (2.8) is established. (2) If R0 > 1 then the proportions of infected and chronic individuals tend to i and v respectively. Hence, we see that lim

t!1

V v v k ¼ lim ¼ ¼ : t!1 i I bþa i

ð2:12Þ

From the second equation in (2.1) and (2.12) we obtain (2.9). From the third equation in (2.1) and (2.12) we also obtain (2.10). (3) The key element in establishing the limits (2.11) is the limit of the quotient V =I. Since in this case the proportions i ! 0 and v ! 0, LÕHospitalÕs rule applies. We have V v v0 ki  ðb þ aÞv : lim ¼ lim ¼ lim 0 ¼ lim t!1 I t!1 i t!1 i t!1 ðci þ dvÞs  ðb þ kÞi If we denote v lim ¼ x t!1 i then x satisfies the equation x ¼

k  ðb þ aÞx ðc þ dx Þ  ðb þ kÞ

ð2:13Þ

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

9

which implies that x is finite and non-zero and is the positive solution of the equation dx2 þ ðc þ a  kÞx  k ¼ 0: Consequently, lim

t!1

V v ¼ lim ¼ x : t!1 I i

From the second equation in (2.1) we have lim

t!1

I0 ¼ ðc þ dx Þ  ðl þ kÞ ¼ ðl þ kÞðR1  1Þ: I

From the third equation in (2.1) and (2.13) we have V0 k ¼  ðl þ aÞ ¼ ðc þ dx Þ  ðl þ kÞ ¼ ðl þ kÞðR1  1Þ: t!1 V x This completes the proof.  lim

As we expect, in the case when the proportions approach the Ôendemic equilibriumÕ (R0 > 1) the exponential rate of increase (decrease) of the susceptible individuals, infected individuals with acute hepatitis and chronic individuals is the same as the exponential rate of increase (decrease) of the total population, that is, it is b  l. In the case when the disease is dying out in the sense that the endemic proportions are approaching the Ôdisease-free equilibriumÕ (R0 < 1), the exponential rate of increase of the infected and chronic class is smaller than the exponential rate of increase of the total population, that is, ðl þ kÞðR1  1Þ < b  l if R0 < 1:

ð2:14Þ



To see this, we notice that x > 0 and hence the numerator and the denominator of the fraction on the right-hand side in (2.13) have the same signs. If both are positive, then x <

k bþa

and (2.14) is valid. If both are negative, then bþkc d and (2.14) is again valid. We note that since R0 < 1, the expression b þ k  c is positive. x <

3. A model with age of infection The chronic stage in hepatitis C consists of several substages. In this section we consider an extension of the ODE model from the previous section in which the chronic stage is stratified by the age-since-infection, that is, the time spent in the chronic stage.

10

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

Let h be the age-since-infection. With the notation from the previous section we consider R1 R1 8 0 S ¼ bN  cI NS  NS 0 dðhÞvðh; tÞ dh  lS þ 0 aðhÞvðh; tÞ dh; > > > R1 > 0 S S > > < I ¼ cI N þ N 0 dðhÞvðh; tÞ dh  ðl þ kÞI; ð3:1Þ vt þ vh ¼ ðl þ aðhÞÞvðh; tÞ; > > > > vð0; tÞ ¼ kI; > > : Sð0Þ ¼ S0 ; Ið0Þ ¼ I0 ; vðh; 0Þ ¼ v0 ðhÞ; where aðhÞ and dðhÞ are non-negative functions of h. As before, aðhÞ is the age structured treatment/recovery rate for the chronic state, and dðhÞ is the age structured effective contact rate of individuals with chronic hepatitis C. We will use the following notation: Z 1 vðh; tÞ dh; V ðtÞ ¼ 0 Rh  aðsÞ ds CðhÞ ¼ e 0 ; where V ðtÞ is the number of chronic individuals and CðhÞ is the probability of being still chronic h time units after becoming chronic. In addition, for an arbitrary function gðh; tÞ we introduce the quantities: Z 1 dðhÞgðh; tÞ dh; Dg ðtÞ ¼ 0 Z 1 ð3:2Þ aðhÞgðh; tÞ dh: Ag ðtÞ ¼ 0

We implicitly assume the integrals exist. Integrating the third equation in (3.1) and assuming that there are no individuals with infinite age-since-infection, that is limh7!1 vðh; tÞ ¼ 0 for all t, we get V 0 ðtÞ ¼ kI  lV  Av ðtÞ:

ð3:3Þ

For the total population size N ¼ S þ I þ V we obtain a Malthus equation of exponential growth: N 0 ¼ ðb  lÞN . We also note that if d and a are constant we obtain the model (2.1). Introducing the proportions s ¼ S=N , i ¼ I=N, u ¼ v=N and using the notation in (3.2) leads to the normalized system: 8 0 s ¼ bð1  sÞ  cis  sDu ðtÞ þ Au ðtÞ; > > > 0 > < i ¼ cis þ sDu ðtÞ  ðb þ kÞi; ut þ uh ¼ ðb þ aðhÞÞuðh; tÞ; ð3:4Þ > > uð0; tÞ ¼ ki; > > : sð0Þ ¼ s0 ; ið0Þ ¼ i0 ; uðh; 0Þ ¼ u0 ðhÞ and the relation Z sðtÞ þ iðtÞ þ

1

uðh; tÞ dh ¼ 1:

ð3:5Þ

0

Two steady states are computed: the disease-free equilibrium with i0 ¼ 0, u0 ¼ 0 and s0 ¼ 1 and an unique endemic state (a non-uniform steady state distribution)

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

bþk ; c þ 0 dðhÞkebh CðhÞ dh 1  s R1 i ¼ ; 1 þ 0 kebh CðhÞ dh s ¼

11

R1

ð3:6Þ

u ðhÞ ¼ ki ebh CðhÞ: We note that although these expressions look very different from the ones in (2.5), it can be easily shown that in the case of constant d and a they coincide. To understand the local dynamics around the disease-free equilibrium, we linearize around it. Set s ¼ s0 þ x, i ¼ y, u ¼ z. Substituting in (3.4) and taking into account the relations satisfied by the disease-free equilibrium we obtain 8 0 x ¼ bx  cy  Dz ðtÞ þ Az ðtÞ; > > < 0 y ¼ cy þ Dz ðtÞ  ðb þ kÞy; ð3:7Þ zt þ zh ¼ ðb þ aðhÞÞzðh; tÞ; > > : zð0; tÞ ¼ ky: Looking for exponential solutions in (3.7), that is, solutions of the form x ¼ ept X ;

y ¼ ept Y ;

z ¼ ept ZðhÞ:

where p is a constant, we obtain 8 pX ¼ bX  cY  DZ þ AZ ; > > < pY ¼ cY þ DZ  ðb þ kÞY ; pZ þ Zh ¼ ðb þ aðhÞÞZðhÞ; > > : Zð0Þ ¼ kY :

ð3:8Þ

ð3:9Þ

We solve the last differential equation and obtain: ZðhÞ ¼ kY eðpþbÞh CðhÞ. Substituting in the second equation and canceling Y we get the characteristic equation: Z 1 dðhÞkeðpþbÞh CðhÞ dh: ð3:10Þ pþbþk ¼cþ 0

We can rewrite the characteristic equation in the form R1 c þ 0 dðhÞkeðpþbÞh CðhÞ dh ¼1 pþbþk and define the basic reproduction number R1 c þ 0 dðhÞkebh CðhÞ dh R0 ¼ : bþk

ð3:11Þ

ð3:12Þ

For p < ðb þ kÞ the expression in the left hand side of (3.11) is negative and the equation has no solution. For p > ðb þ kÞ the expression on the left hand side of (3.11) is a decreasing function of p which approaches 1 as p ! ðb þ kÞ and zero as p ! 1. Therefore, the characteristic equation always has a unique real solution p . If R0 < 1, then for p ¼ 0 the left hand side in (3.11) is smaller than one. Hence, the real solution satisfies p < 0. In addition every complex solution p has as a real part, Rp, that satisfies

12

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

cþ

R1 0

dðhÞkeðRpþbÞh CðhÞ dh P 1: Rp þ b þ k

Since the function on the left hand side of (3.11) is a decreasing function in Rp and p satisfies (3.11), we have that Rp 6 p . Hence, any complex solution of (3.11) has a real part smaller than the unique real solution of (3.11). Therefore, if R0 < 1, then the disease-free equilibrium is locally stable. We note that in this case s > 1 (see (3.6)) and an endemic state does not exist. If R0 > 1, then s as given by (3.6) is smaller than one (since it is the reciprocal of R0 ) and an endemic state exists and is given by (3.6). Moreover, the real solution of the characteristic equation satisfies p > 0 and therefore the disease-free equilibrium is locally unstable. We have proved the following theorem. Theorem 3.1. If R0 < 1 then the disease-free equilibrium (1,0,0) is locally stable. If R0 > 1 the disease-free equilibrium is unstable and there is a unique endemic equilibrium ðs ; i ; v Þ. In the next section we discuss the global stability of the disease-free equilibrium and, as in the unstructured case, it is shown that if R0 < 1 the disease-free equilibrium is globally stable and there is a unique endemic equilibrium. We also derive and investigate the characteristic equation obtained from the linearization around the endemic equilibrium. We conclude that the endemic equilibrium is locally stable under appropriate conditions. 4. Global stability of the disease-free equilibrium and local stability of the endemic equilibrium We begin by showing global stability of the disease-free proportions. We use a standard argument discussed in [23]. Lemma 4.1. Assume dðhÞ is a bounded function and d ¼ sup dðhÞ:

ð4:1Þ

h2½0;1Þ

If R0 < 1 then the disease-free equilibrium (1,0,0) is globally asymptotically stable. Proof. Integrating the third equation of (3.4) along the characteristic lines we obtain:  u0 ðh  tÞebt Cðh  t; hÞ h P t; uðh; tÞ ¼ ð4:2Þ h < t; kiðt  hÞebh CðhÞ Rh  aðsÞ ds is the probability of an individual who has been chronic for h  t where Cðh  t; hÞ ¼ e ht units after infection to remain chronic until h units after infection. We can estimate the behavior of Du ðtÞ at infinity: Z 1 Du ðtÞ ¼ dðhÞuðh; tÞ dh 0 Z t Z 1 bh ¼ dðhÞkiðt  hÞe CðhÞ dh þ dðhÞu0 ðh  tÞebt Cðh  t; hÞ dh 0 t Z t 6 dðhÞkiðt  hÞebh CðhÞ dh þ dku0 kL1 ebt : ð4:3Þ 0

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

Therefore, considering the lim sup as t ! 1, Z 1 dðhÞkebh CðhÞ dh lim sup iðtÞ: lim sup Du ðtÞ 6 t

13

ð4:4Þ

t

0

To estimate lim supt iðtÞ we integrate the equation for the proportion of the individuals with acute hepatitis C (see (3.4)) and use the fact that s 6 1 to obtain Z t ðbþkÞt i0 þ eðbþkÞs ½ciðt  sÞ þ Du ðt  sÞ ds; ð4:5Þ iðtÞ 6 e 0

where i0 is the initial value for i. Hence, lim sup iðtÞ 6 t

1 lim sup½ciðtÞ þ Du ðtÞ: bþk t

Combining (4.4) and (4.6) we have   Z 1 1 bh cþ lim sup iðtÞ 6 dðhÞke CðhÞ dh lim sup iðtÞ: bþk t t 0

ð4:6Þ

ð4:7Þ

The coefficient on the right hand side of this inequality is exactly R0 (see (3.12)) which we have assumed strictly smaller than one. In this case, the only way inequality (4.7) can hold is if lim sup iðtÞ ¼ 0: t

From (4.2) we also have lim sup uðh; tÞ ¼ 0 for all h fixed:



t

Lemma 4.1 guarantees that if R0 < 1 the proportions of infected and chronic individuals tend to zero. As we mentioned before this does not imply that the number of infecteds is going to zero or even that the number of infecteds is constant or decreasing. Hence, this result does not tell us that the disease is disappearing from the population but only that its basis is becoming more diluted. In the rest of this section we compose and investigate the characteristic equation. As a result we establish local stability of the endemic equilibrium under some appropriate conditions. From the simulations carried out with the numerical method in Section 5, it appears that the endemic equilibrium may be unstable and sustained oscillations are possible. Thus, the presence of host age-structure leads to much more complex dynamical behavior than the one exhibited by the ODE model where age-structure has been ignored. We begin by linearizing (3.4) around the endemic state. We denote by x, y and z the perturbations as follows: s ¼ s þ x;

i ¼ i þ y;

u ¼ u þ z;

where s ; i ; u are given in (3.6). The perturbations can take both positive and negative values. We obtain the following system of equations:

14

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

8 0 x ¼ bx  ½ci þ Du x  ½cy þ Dz ðtÞs þ Az ðtÞ; > > > 0 > > < y ¼ ½ci þ Du x þ ½cy þ Dz ðtÞs  ðb þ kÞy; zh þ zt ¼ ðb þ aðhÞÞzðh; tÞ; > > > zð0; tÞ ¼ ky; > > : x0 ¼ sð0Þ  s ; y0 ¼ ið0Þ  i ; z0 ðhÞ ¼ u0 ðhÞ  u ðhÞ: From equation (3.5) it follows that the following algebraic equation is also satisfied: Z 1 zðh; tÞ dh ¼ 0: xðtÞ þ yðtÞ þ

ð4:8Þ

ð4:9Þ

0

We will use the following notation ha ðhÞ ¼ aðhÞkebh CðhÞ; hd ðhÞ ¼ dðhÞkebh CðhÞ;

ð4:10Þ

hðhÞ ¼ kebh CðhÞ: For any function gðtÞ we use g^ðkÞ to denote the Laplace transform of g defined as follows: Z 1 g^ðkÞ ¼ ekt gðtÞ dt: ð4:11Þ 0

The stability of the endemic equilibrium is determined by the eigenvalues of the linear operator in (4.8). To obtain the system for the eigenvalues we formally look for separable solutions of the system (4.8) in the form xðtÞ ¼ ekt x;

yðtÞ ¼ ekt y ;

zðh; tÞ ¼ ektzðhÞ:

To avoid cumbersome notation we will omit the bars and denote both the time dependent perturbations and the time independent perturbations by x; y; zðhÞ. Hence, we arrive at the following linear eigenvalue problem (where k denotes an eigenvalue): 8 kx ¼ bx  ½ci þ Du x  ½cy þ Dz ðkÞs þ Az ðkÞ; > > > < ky ¼ ½ci þ D x þ ½cy þ D ðkÞs  ðb þ kÞy; u z ð4:12Þ > zh ðhÞ ¼ ðk þ b þ aðhÞÞzðhÞ; > > : zð0Þ ¼ ky with the algebraic condition Z 1 zðhÞ dh ¼ 0: xþyþ

ð4:13Þ

0

We solve the differential equation for z: zðhÞ ¼ kyeðkþbÞh CðhÞ

ð4:14Þ

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

and eliminate z from the remaining equations. We have Z 1 dðhÞeðkþbÞh CðhÞ dh ¼ h^d ðkÞy; Dz ðkÞ ¼ ky 0 Z 1 Az ðkÞ ¼ ky aðhÞeðkþbÞh CðhÞ dh ¼ h^a ðkÞy; 0 Z 1 Ez ðkÞ ¼ ky eðkþbÞh CðhÞ dh ¼ h^ðkÞy:

15

ð4:15Þ

0

With the notation in (4.15) the equations (4.12) can be rewritten as  ðk þ b þ ½ci þ Du Þx ¼ ½cs þ s h^d ðkÞy þ h^a ðkÞy; ðk þ b þ kÞy ¼ ½ci þ Du x þ ½cs þ s h^d ðkÞy:

ð4:16Þ

Condition (4.13) becomes x þ y þ h^ðkÞy ¼ 0:

ð4:17Þ

Eqs. (4.16) and (4.17) give three conditions for the unknowns x and y. However, one of the equations is a consequence of the other two. In particular, adding the equations in (4.16) and using the fact that h^a ðkÞ ¼ k  ðk þ bÞh^ðkÞ leads to Eq. (4.17). Eliminating x and y from the second equation in (4.16) and (4.17) we obtain an equation for k, referred to as the characteristic equation. k þ b þ k þ ½ci þ Du ð1 þ h^ðkÞÞ ¼ cs þ s h^d ðkÞ:

ð4:18Þ

Let k have real part n and imaginary part g. The equations for the real and the imaginary part become n þ b þ k þ ½ci þ Du ð1 þ Rh^ðkÞÞ ¼ cs þ s Rh^d ðkÞ; g  ½ci þ Du Ih^ðkÞ ¼ s Ih^d ðkÞ;

ð4:19Þ

where Rh^d ðkÞ ¼

Z

1

enh cosðghÞdðhÞkebh CðhÞ dh; 0

Ih^d ðkÞ ¼

Z

1

enh sinðghÞdðhÞkebh CðhÞ dh; 0

Rh^ðkÞ ¼

Z

1

enh cosðghÞkebh CðhÞ dh;

0

Ih^ðkÞ ¼

Z

1

enh sinðghÞkebh CðhÞ dh:

0

In the expressions above Rh^d and Ih^d are the real part and the imaginary part of h^d , while Rh^ and Ih^ are the real part and the imaginary of h^.

16

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

Proposition 4.2. If R0 > 1 and if Z 1 kebh CðhÞ dh < 1

ð4:20Þ

0

then the endemic equilibrium (3.6) is locally asymptotically stable. Proof. We note that for n P 0, the real part Rh^ðkÞ < h^ð0Þ and hence, 1 þ Rh^ðkÞ > 1  h^ð0Þ > 0: The last inequality follows from condition (4.20). Consequently, for n P 0 the left hand side of the first equation in (4.19) is larger than b þ k. On the other hand for n P 0 the right hand side of the same equation satisfies cs þ s Rh^d ðkÞ 6 cs þ s jRh^d ðkÞj 6 cs þ s h^d ð0Þ ¼ b þ k: Thus, the system (4.19) does not have a solution with n P 0, or equivalently, there are no eigenvalues with non-negative real part. This completes the proof.  In interpreting inequality (4.20) we notice that Z 1 Z 1 1 bh ke CðhÞ dh ¼ u ðhÞ dh: i 0 0 Thus, inequality (4.20) is a requirement that the ratio of the proportion of chronic individuals to the proportion of acutely infected individuals at equilibrium is smaller than one.

5. A numerical method Homogeneous models have routinely become present among epidemic models. Most of the interest is directed towards ODEs; however, in the recent years some work on PDEs has also appeared [24–27]. Yet no numerical methods specifically designed for homogeneous PDEs have been developed. We propose here a first order numerical method for the infection-age structured model of hepatitis C. To avoid problems with the exponential growth and decay of the solutions, we construct the numerical scheme for the normalized system of equations. The difficulty here is that the system in ðN þ 3Þ equations in ðN þ 2Þ (for each time step) unknowns has to be consistent. Consistency means that one of the equations has to follow from the others at each time level. Although most of the numerical methods for epidemic models are implicit [28–30], here, the purely implicit schemes do not satisfy the consistency condition. A mixed implicit–explicit method is therefore developed for this system. We discretize the system (3.4) using a finite difference method along the characteristic lines for the PDE. Since the characteristic lines have slope one, it is convenient to take the step in time and age to be the same. We denote the step by h. Let T ¼ Mh be the final time of the simulation and A ¼ Nh be the final age. The maximum possible human age is now assumed to be 125 years. We take the age-since-infection A smaller than 125. In time, one can simulate (at least in theory) for

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

17

very long periods of time. Hence, T can be larger than A. We denote by S n , I n , Ujn the approximation of sðnhÞ, iðnhÞ, and uðjh; nhÞ respectively. We have also used the notation kU n k ¼

N X

jUkn jh;

1

DnU ¼

N X

dk jUkn jh;

ð5:1Þ

1

AnU ¼

N X

ak jUkn jh:

1

where dk ¼ dðkhÞ and ak ¼ aðkhÞ for k ¼ 0; . . . ; N . The method is given by 8 n n1 S S n n > ¼ bð1  S n Þ  cI n1 S n  Dn1 nP1 U S þ AU > h > > > I n I n1 n n n1 > > ¼ cI n1 S n þ Dn1 nP1 > U S  bI  kI < h U n U n1

j j1 ¼ ðb þ aj ÞUjn > h > > > > U0n ¼ kI n > > > : n S þ I n þ kU n k ¼ 1

j; n P 1

ð5:2Þ

nP1 nP0

We note that the first and the third equations are implicit but the second is not since the righthand side of the second equation depends not only on I n but also on I n1 . If we know S n1 , I n1 , Ujn1 for j ¼ 0; . . . ; N we can first compute Ujn from the third equation, next S n from the first equation, then I n from the second equation and finally U0n from the boundary condition which goes with the third equation. We show that the system of ðN þ 3Þ equations in ðN þ 2Þ unknowns (for each time step) has a solution. We do that by showing that the fourth equation is a consequence of the rest. Multiplying by h and summing the third equation, we get kU n k  kU n1 k ¼ kI n1  bkU n k  AnU : h

ð5:3Þ

Adding this equation to the sum of the first and second equation in (5.2) we obtain S n þ I n þ kU n k ¼

bh þ ðS n1 þ I n1 þ kU n1 kÞ : bh þ 1

Using this relation, one can see inductively that if S 0 þ I 0 þ kU 0 k ¼ 1 then S n þ I n þ kU n k ¼ 1 for every n. Analogously, we can derive any of the equations from the other three. The price we pay for using mixed implicit–explicit method is that the scheme may not generate physically meaningful, that is non-negative, solutions for all values of the step h. In particular, we have the following result. Proposition 5.1. Assume the step h satisfies h < 1=k. Then the solution of the system (5.2) is nonnegative and satisfies 0 6 S n 6 1, 0 6 I n 6 1, 0 6 Ujn 6 1 for all j ¼ 0; . . . ; N and n P 1.

18

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

Proof. Solving the system (5.2) we obtain 8 S n1 þ bh þ hAnU > n > S ¼ > > > 1 þ bh þ chI n1 þ hDn1 U > > > n1 n1 n n > ð1  khÞI þ chI S þ hDn1 > n U S > > > > Ujn ¼ > > 1 þ hðb þ aj Þ > > > > n > U ¼ kI n > > : n0 S þ I n þ kU n k ¼ 1

nP1 nP1 ð5:4Þ j; n P 1 nP1 nP0

We know that S 0 P 0, I 0 P 0 and Uj0 P 0 for j ¼ 0; . . . ; N . Assume S n1 P 0, I n1 P 0 and Ujn1 P 0 for j ¼ 0; . . . ; N . Clearly, in this case we have Ujn P 0 for j ¼ 1; . . . ; N . This implies AnU P 0 and hence we have S n P 0. The restriction on the step guarantees that I n P 0. Finally, we have U0n P 0.  We now demonstrate that the numerical method (5.2) converges at a first order rate. We first introduce the following notation for the evaluations of the solution of (3.4): sn ¼ sðtn Þ;

in ¼ iðtn Þ;

unj ¼ uðhj ; tn Þ:

ð5:5Þ

To compute the errors at each step we evaluate (3.4) at time level tn and approximate the derivatives with backward Euler finite differences: 8 sn sn1 ¼ bð1  sn Þ  cin sn  Dnu sn þ Anu þ OðhÞ n P 1 > > > in ihn1 > > nP1 < h ¼ cin sn þ Dnu sn  ðb þ kÞin þ OðhÞ unj un1 ð5:6Þ j1 n ¼ ðb þ aj Þuj þ OðhÞ j; n P 1 > h > > n n > nP1 > : un0 ¼ ki n nP0 s þ i þ kun k ¼ 1 þ OðhÞ We introduce notation for the errors in the approximation: nn ¼ sn  S n ;

gn ¼ i n  I n ;

rnj ¼ unj  Ujn :

ð5:7Þ

With this notation we have the following result about the convergence of the numerical scheme (5.2): Theorem 5.2. Assume the step h satisfies h < 1=k. Assume that the solution of (3.4) is twice continuously differentiable and that a ¼ sup aðhÞ; d ¼ sup dðhÞ: h2½0;1Þ

h2½0;1Þ

Then, for n ¼ 0; . . . ; M, we have jnn j þ jgn j þ krn k 6 Qh; where Q is a constant independent of h.

ð5:8Þ

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

19

The proof of Theorem 5.2 is given in Appendix A. We have used the numerical method (5.2) to perform simulations with s, i, kuk. In a sequence of runs we varied the parameters of the model, and specifically the functions aðhÞ and dðhÞ, as well as the initial conditions. The results of our simulations in terms of the long term behavior of the functions s, i, kuk varied significantly. They reveal that the endemic equilibrium of (3.4) may be locally stable with either asymptotically approaching solutions or solutions approaching through damped oscillations, or it may be unstable and solutions approaching sustained oscillations. We include here the results of two examples which are representative of the more interesting results that we obtained from all simulations. The examples we present here have parameters that do not satisfy the condition in Proposition 4.2. Nonetheless, in our first example the infected classes i, kuk tend towards equilibrium values, but they do so via damped oscillation (see Figs. 1 and 2). The parameters of the model are b ¼ 1, k ¼ 25, c ¼ 87:4075053,  ð7:2125Þhð0:5  hÞ 0 6 h 6 0:5; dðhÞ ¼ 0 h > 0:5;  0 h 6 0:5; aðhÞ ¼ 10 000 h > 0:5: In this example we begin from the following initial conditions: s0 ¼ 0:29, i0 ¼ 0:07 and u0 ðhÞ ¼ 1:28e2h . In Fig. 1 the proportion of acutely infected individuals tends to the steady state i
0:0657 via damped oscillations. In Fig. 2 the integral of u, denoted by U ðtÞ Z 1 uðh; tÞ dh; U ðtÞ ¼ kuk ¼ 0

tends to the steady state 0:64652 also via damped oscillations. The proportion of susceptible individuals tends to a steady state of s
0:28777 which gives a reproduction number R0
3:475. The second example has parameters similar to the first example, however, the endemic equilibrium in this case is unstable and the proportions stabilize in an oscillatory solution. Consequently, the endemic equilibrium in the structured case cannot be globally stable. Hence, the incorporation of host age-structure into the model leads to a significantly different and more complex dynamics. The parameters of the model in this case are b ¼ 0:5, k ¼ 25, c ¼ 65:817518, aðhÞ is the same as in the first example and dðhÞ is given below:  ð2:543688Þhð0:5  hÞ 0 6 h 6 0:5; dðhÞ ¼ 0 h > 0:5: The initial conditions are the same as in the previous example. In Fig. 3 the proportions of infected individuals clearly exhibit sustained oscillations. Note that they become very close to zero as they oscillate. In Fig. 4 a plot of U ðtÞ is given which also clearly stabilizes in sustained oscillations. Comparing Figs. 3 and 4 we see that inequality (4.20) is violated. The mechanism of constructing an oscillatory solution is similar to the one described in [28]. The simulations are performed with step 0.001. The final infection-age is taken to be 10. The motivation for this choice is that u ðhÞ and u0 ðhÞ are of order of magnitude of e2h and, therefore, the integral from 0 to 10 of u ðhÞ is almost the same as the integral from zero to infinity.

20

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

Fig. 1. Plot of the proportions of infecteds as a function of time with parameters as given in the text (Example 1). The proportion of infecteds tends to an equilibrium via damped oscillations.

Fig. 2. Plot of the total proportion (the integral of the age-density) of chronic individuals as a function of time with parameters as given in the text (Example 1). The total proportion of chronic individuals tends to an equilibrium via damped oscillations.

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

21

Fig. 3. Plot of the proportion of infecteds as a function of time with parameters as given in the text (Example 2). The proportion of infecteds exhibits sustained oscillations.

Fig. 4. Plot of the total proportion (the integral of the age-density) of chronic individuals as a function of time with parameters as given in the text (Example 2). The total proportion of chronic individuals exhibits sustained oscillations.

22

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

6. Discussion Hepatitis C is a disease of the liver characterized by its long chronic stage. Since it is vastly asymptomatic and the acute stage is relatively short, patients are most often identified as having hepatitis C during the infectious chronic stage. The focus of this paper has been on the study of the role of the chronic infectious stage on the long-term dynamics of HCV. First we look at the simplest setting, that is, it is assumed that the chronic infectious stage is exponentially distributed while the infectivity of each chronically infectious individual is the same. Second, the assumptions on the distribution of the chronically infectious period and the infectivity are weakened. It is assumed that both depend on the age of the chronic infection. It is shown that the local and global results concerning the disease-free equilibrium remain the same. It is established that the unique endemic is locally stable under some additional conditions. Numerical simulations demonstrate that the weakening of the assumptions on the infectivity and the chronic-infectious-period distribution, however, can lead to sustained oscillations and, consequently, do generate qualitatively distinct dynamics than those exhibited by the ODE model. This result is in agreement with results from other models with variable infectivity and age-of-infection which are known to become unstable and have sustained oscillations [15,16,28]. In summary, we have completely investigated the long term dynamical behavior of the ODE model. In particular, we have shown that there exists a threshold quantity, which we have denoted by R0 , such that if it is below one the proportions of infected and chronic individuals decrease in time to zero. In the case of hepatitis C the reproduction number depends on several parameters all of which can influence its magnitude although in a different way. With all other parameters fixed R0 increases linearly with the effective contact rate of susceptible individuals with individuals who have acute hepatitis. The larger the slope of the line, the bigger the effect of the change of c on R0 . The slope of the line depends in an inversely proportional manner on the birth/immigration rate and the rate of progression to the chronic stage. Similarly, the reproduction number depends linearly on the effective contact rate of susceptible individuals with chronic individuals. In this case the slope of the line increases as the rate of progression to chronic stage increases and as the birth/immigration rate decreases. Hence, the reproduction number can be reduced by decreasing the transmission rates and this decrease will be more substantial if the birth/immigration rate is as small as possible. Consequently, in countries with small birth/immigration rates the control measures should be directed towards reducing the effective transmission rates from infected and chronic individuals. The transmission rates were substantially reduced in 1992 when effective screening of the blood for HCV was implemented [31, p. 31]. Secondary safety procedures can also be implemented to further decrease the rates of transmission. It is easy to see that R0 decreases if the rate at which new healthy individuals enter the system increases. This effect is a consequence of the fact that the proportions of infecteds and carriers decrease which hampers the distribution of the disease. The dependence of the reproduction number on the rate of progression to chronic stage is more complicated. We expect that for all realistic values of the parameters R0 will decrease as the rate of progression to chronic stage increases. This mechanism will have a suppressed impact in the case of hepatitis C since the acute stage is in most cases asymptomatic and, on a personal level, transition to the chronic stage will not be viewed as improvement. However, for diseases where the

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

23

acute stage is characterized by severe symptoms, it is likely that a drug which speeds the transition to a milder chronic stage will be better accepted. Our observation implies that such speeding of the transition to chronic stage (possibly with a drug) in effect contributes to the eradication of the disease. We would like to note that in the unlikely case when the transmission rate of the chronic stage is much larger than the transmission rate of the acute stage, it is possible that R0 increases with the rate of progression to the chronic stage. Treatment which leads to a healthy state is, of course, an obvious and effective way of reducing the reproduction number and eradicating the disease. There are many interesting and unanswered questions associated with HCV treatment and therapies which may impact the epidemiology of the disease. Our model is a macro (population) model and therefore it cannot give answers to most of the above questions. Models of hepatitis C at the immune system level should be used. In [32] the authors address a possible mechanism of action of interferon-a on HCV and the efficacy of the therapy with it. Acknowledgements This research in its larger part was performed while M.M. was a postdoctoral associate at IMA, University of Minnesota. M.M. is partially supported by NSF grant DMS-0137687. The authors thank the reviewers for their valuable comments. The authors also thank Horst Thieme for a helpful discussion. Appendix A Proof of Theorem 5.2. We subtract Eq. (5.2) from Eq. (5.6) to obtain the equations for the errors: 8 nn nn1 n n1 n n ¼ bnn  cS n gn1  cin nn  Dn1 > U n  Dr s þ Ar þ OðhÞ; > > gn gh n1 < n n1 n n n1 ¼ cS n gn1 þ cin nn þ Dn1 þ OðhÞ; U n þ Dr s  bg  kg h ðA:1Þ n n1 rj rj1 > n > > : n h ¼n ðb þ aj Þrj þ OðhÞ; r0 ¼ kg : From (3.4) it follows that sn 6 1, in 6 1, kun k 6 1. We also know that S n 6 1, I n 6 1, kU n k 6 1. Multiplying by h the third equation and adding, we obtain krn k 6 khjgn1 j þ krn1 k þ Oðh2 Þ: Solving for nn , gn , and taking absolute value, we obtain: 8 n < jn j 6 jnn1 j þ C1 hjgn1 j þ C2 hkrn1 k þ Oðh2 Þ; jgn j 6 C3 hjnn1 j þ ð1 þ C4 hÞjgn1 j þ C5 hkrn1 k þ Oðh2 Þ; : n kr k 6 khjgn1 j þ krn1 k þ Oðh2 Þ;

ðA:2Þ

ðA:3Þ

where C1 ; . . . ; C5 are positive constants, independent of h. Adding the three inequalities above, we get: jnn j þ jgn j þ krn k 6 ð1 þ PhÞjnn1 j þ ð1 þ RhÞjgn1 j þ ð1 þ VhÞkrn1 k þ Oðh2 Þ;

ðA:4Þ

24

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

where P ; R; V are constants independent of h. Set q ¼ maxðP ; R; V Þ. Then (A.4) becomes jnn j þ jgn j þ krn k 6 ð1 þ qhÞðjnn1 j þ jgn1 j þ krn1 kÞ þ Oðh2 Þ:

ðA:5Þ

Iterating this inequality in itself, we obtain n

n1

jnn j þ jgn j þ krn k 6 ð1 þ qhÞ ðjn0 j þ jg0 j þ kr0 kÞ þ ðð1 þ qhÞ

n2

þ ð1 þ qhÞ

þ    þ 1ÞOðh2 Þ:

Taking into account that ð1 þ qhÞn 6 eqT we get jnn j þ jgn j þ krn k 6 Qh; where Q is a constant independent of h but dependent on T . 

References [1] R. Purcell, Hepatitis C virus: an introduction, NIH Consensus Development Conference on Management of Hepatitis C, March 24–26, Bethesda, MD, 1997. [2] M. Alter, Hepatitis C: the clinical spectrum of the disease, NIH Consensus Development Conference on Management of Hepatitis C, March 24–26, Bethesda, MD, 1997. [3] L. Seeff, Natural history of hepatitis C, NIH Consensus Development Conference on Management of Hepatitis C, March 24–26, Bethesda, MD, 1997. [4] U. Bandy, Hepatitis C virus (HCV): a silent epidemic, Med. Health R.I. Jun 82(6) (1999) 223. [5] R. Colina, C. Azambuja, R. Uriarte, C. Mogdasy, J. Cristina, Evidence of increasing diversification of hepatitis C viruses, J. General Virol. 80 (1999) 1377. [6] S. Busenberg, P. van den driessche, Analysis of a disease transmission model in a population with varying size, J. Math. Biol. 28 (1990) 257. [7] S. Busenberg, K. Hadeler, Demography and epidemics, Math. Biosci. 101 (1990) 63. [8] K. Hadeler, C. Castillo-Chavez, A core group model for disease transmission, Math. Biosci. 128 (1995) 41. [9] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol. 38 (1999) 22. [10] W.-M. Liu, P. van den Driessche, Epidemiological models with varying population size and dose-dependent latent period, Math. Biosci. 128 (1995) 57. [11] W. Derrick, P. van den Driessche, A disease transmission model in a nonconstant population, J. Math. Biol. 31 (1993) 495. [12] H. Thieme, Epidemic and demographic interaction in the spread of potentially fatal diseases in a growing population, Math. Biosci. 111 (1992) 99. [13] R. Saenz, J. Velasco-Hernandez, C. Castillo-Chavez, Analysis of an age-structured epidemic model with a chronic state (manuscript). [14] B. Reade, R. Bowers, M. Begon, R. Gaskell, A Model of disease and vaccination for infections with acute and chronic phases, J. Theor. Biol. 190 (1998) 355. [15] H. Thieme, C. Castillo-Chavez, How may infection-age-dependent infectivity affect the dynamics of HIV/AIDS? SIAM J. Appl. Math. 53 (5) (1993) 1447. [16] H. Thieme, C. Castillo-Chavez, On the role of variable infectivity in the dynamics of HIV epidemic, Mathematical and statistical approaches to AIDS epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer, Berlin, 1993, p. 157. [17] S. Deuffic, L. Buffat, T. Poynard, A.J. Valleron, Modelling hepatitis C virus epidemic in France, Hepatology 29 (5) (1999) 1596. [18] T.J. Case, An Illustrated Guide to Theoretical Ecology, Oxford University, New York, 2000. [19] S. Busenberg, K. Cooke, H. Theime, Demographic change and persistence of HIV/AIDS in a heterogeneous population, SIAM J. Appl. Math. 51 (4) (1990) 1030.

M. Martcheva, C. Castillo-Chavez / Mathematical Biosciences 182 (2003) 1–25

25

[20] K. Hadeler, R. Waldst€atter, A. W€ orz-Busekros, Models for pair formation in bisexual populations, J. Math. Biol. 26 (1988) 635. [21] K. Hadeler, Pair formation with maturation period, J. Math. Biol. 32 (1993) 1. [22] K. Hadeler, Periodic solutions of homogeneous equations, J. Differ. Eq. 95 (1992) 183. [23] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, Comitato Nazionale per le Scienze Matematiche, C.N.R., 7 addr Gardini Editori e Stampatori in Pisa, 1995. [24] K. Hadeler, J. M€ uller, in: V. Isham, G. Medley (Eds.), Epidemic Models: Their Structure and Relation to Data, Cambridge University, Cambridge, UK, 1996. [25] K. Hadeler, J. M€ uller, in: V. Isham, G. Medley (Eds.) Epidemic Models: Their Structure and Relation to Data, Cambridge University, Cambridge, UK, 1996. [26] G. Webb, Asynchronous exponential growth in differential equations with homogeneous nonlinearities. Differential equations in Banach spaces (Bologna, 1991), Lecture Notes in Pure and Applied Mathematics, vol. 148, Dekker, New York, 1993, p. 225. [27] J. M€ uller, Optimal vaccination patterns in age-structured populations, SIAM J. Appl. Math. 59 (1) (1998) 222. [28] F. Milner, A. Pugliese, Periodic solutions: a robust numerical method for an S-I-R model of epidemics, J. Math. Biol. 39 (6) (1999) 471. [29] M. Iannelli, M. Kim, E. Park, Splitting methods for the numerical approximation of some models of age-structured population dynamics and epidemiology, Appl. Math. Comp. 87 (1997) 69. [30] M. Iannelli, F. Milner, A. Pugliese, Analytical and numerical results for the age-structured S-I-S epidemic model with mixed inter–intracohort transmission, SIAM J. Math. Anal. 23 (3) (1992) 662. [31] C. Turkington, Hepatitis C: The Silent Killer, Contemporary Books, 1998. [32] A. Neumann, N. Lam, H. Dahari, D. Gretch, T. Wiley, T. Layden, A. Perelson, Hepatitis C viral dynamics in vivo and the antiviral efficacy of Interferon-a therapy, Science 282 (1998) 103.