Vaccination against multiple HPV types

Vaccination against multiple HPV types

Mathematical Biosciences 197 (2005) 88–117 www.elsevier.com/locate/mbs Vaccination against multiple HPV types Elamin H. Elbasha a a,* , Alison P. G...

722KB Sizes 0 Downloads 21 Views

Mathematical Biosciences 197 (2005) 88–117 www.elsevier.com/locate/mbs

Vaccination against multiple HPV types Elamin H. Elbasha a

a,*

, Alison P. Galvani

b

Merck Research Laboratories, 10 Sentry Parkway, BL 2-3, Blue Bell, PA 19422, United States b Department of Epidemiology and Public Health, Yale University School of Medicine, New Haven, CT 06520-8034, United States Available online 10 August 2005

Abstract Vaccines against the most common human papillomavirus (HPV) types are currently under development. Epidemiologic data suggest that the transmission dynamics of different HPV types are not independent. Some studies indicate that interactions among HPV types are synergistic, where infection with one type facilitates concurrent or subsequent infection with another HPV type. Other studies point to antagonistic interference among HPV types. Here we develop a mathematical model to explore how these interactions may either enhance or diminish the effectiveness of vaccination programs designed to reduce the prevalence of the HPV types associated with cervical cancer. We analyze the local stability of the infection-free and boundary equilibria and characterize the conditions leading to a coexistence equilibrium. We also illustrate the results with numerical simulations using realistic model parameters. We show that if interactions among HPV types are synergistic, mass vaccination may reduce the prevalence of types that are not even included in the vaccine. Ó 2005 Elsevier Inc. All rights reserved. Keywords: HPV types; HPV disease; Cancer; Coevolutionary interactions; Vaccine; Mathematical model; Local stability

*

Corresponding author. Tel.: +1 484 344 4045; fax: +1 484 344 3855. E-mail address: [email protected] (E.H. Elbasha).

0025-5564/$ - see front matter Ó 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2005.05.004

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

89

1. Introduction Human papillomavirus (HPV) is one of the most prevalent sexually transmitted infections (STI). Persistent infection with certain HPV types is the primary cause of cervical carcinoma and its precursor lesions [1]. Cancer of the uterine cervix is the most common malignancy among women in developing countries and a leading cause of cancer death worldwide, with an estimated 233,000 deaths in 2000 [2]. HPV infection has also been linked with other anogenital cancers, benign codylomata and genital warts among both males and females. The economic burden of HPVrelated diseases is significant. For example, healthcare costs of cervical HPV-related disease in 1998 for the US alone are estimated at $3.4 billion [3]. Over 100 HPV types have been identified, of which about 40 infect the anogenital tract. Approximately 15 anogenital HPV types are oncogenic [4]. HPV 16 is one of the most prevalent types [5,6], and is estimated to be responsible for about half of the cervical cancers worldwide [7]. Types 18, 31, and 45 are also thought to be carcinogenic [7,8], while types 33, 35, 39, 51, 52, 56, 58, 59 and 69 are less closely associated with cervical carcinoma [7]. HPV 6 and 11 are associated with more than 90% of cases of genital warts. An estimated 20 to 30% of HPV-infected women harbor multiple types that might be concurrently or sequentially acquired [6], making interactions among HPV types a potentially important public health issue. Vaccines designed to prevent specific HPV types are currently in development. In early phase human trials, administration of a monovalent HPV 16 vaccine was shown to prevent acquisition of persistent HPV 16 infection and disease in women who were HPV 16-naı¨ve at enrollment and through the vaccination period [9–11]. If these early findings are confirmed in larger studies, administration of HPV vaccine is expected to become an essential component in any public health policy aimed at reducing the prevalence of HPV and its associated morbidity and mortality. The enormous diversity of HPV types means that it is not feasible to develop a multivalent HPV vaccine that includes all HPV types implicated in the cause of cervical intraepithelial neoplasia (CIN) and cancer. Of the prophylactic vaccines currently in development, a quadrivalent HPV (against types 6, 11, 16 and 18) vaccine has the broadest coverage, while other vaccines target fewer types [12]. Since current vaccine candidates do not target all HPV types, it is paramount to determine the epidemiological implications of a vaccine that targets a subset of the HPV types in circulation. It is particularly important to consider the effect of such a vaccine on the prevalence of HPV types not included in the vaccine. Clinical and epidemiologic studies that have attempted to characterize the nature of the interactions acting among types have generated mixed results. Our model elucidates the likely impact of mass vaccination which depends on the interactions between the HPV types included in the vaccine and those against which the vaccine might offer partial protection. If infection with one type increases the probability of infection with another type, the interactions among the types are said to be synergistic. Conversely, if infection with one type reduces the likelihood of infection with another type, the interactions among the types are said to be antagonistic. Epidemiologic data have pointed to synergistic interactions among different HPV types. Thomas et al. [13] found that concurrent infection was more frequent than would be predicted if infections were independent. Similarly, Liaw et al. [14] observed that women infected initially with HPV 16 experienced an eightfold increase in the risk for subsequent infection with types 18, 39, 45, 59, and 68. Rousseau et al. [15] also found that previous HPV infection predisposed women to subsequent infection with different HPV types. HPV 16 and HPV 18 had the most pronounced

90

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

effect in terms of increasing the probability of becoming infected with further types. While persistence of infection is a risk factor for the development of cervical cancer, persistence appeared independent of infection history [15]. Immunity against HPV is thought to be primarily type-specific [16–20], with some evidence for cross-immunity between genetically related types, such as among HPV 16, 31, 33, 35, 52, and 58 [17,21,22]. For example, Ho et al. [19] found that IgG antibodies against HPV 16 reduce the hazard of subsequent infection with HPV 16, and genetically related types HPV 31, 33, 35, 52 and 58. Additionally, cross-reactivity between antibodies to HPV 16 and related types has been demonstrated [23,24]. Cross-immunity will bring types into competition for the niche of susceptible hosts. If there is competition between the two HPV types, a vaccine that reduces the prevalence of one type may promote the prevalence of other types through a process of competitive release. Conversely, if the interactions are synergistic, mass vaccination may reduce prevalence of all types, even those not included in the vaccine. Thus, in the absence of definitive evidence regarding the nature of the interaction among oncogenic HPV types, mathematical models offer useful tools for assessing the potential impact of type-specific or multivalent HPV vaccination programs. Previous mathematical models have analyzed multi-strains interactions for various diseases (e.g., HIV and dengue) with mass vaccination [25–29] or without mass vaccination [30–34]. These studies have focused on the effects of competition and cross-immunity between pathogens strains, but have not considered the role that synergistic interactions may play in the transmission dynamics of HPV and on the public health impact of HPV vaccination. The effect of cross-reactive immunological response on persistence of multiple pathogen strains has been studied in the context of dengue virus [35–37]. Few models of HPV vaccination strategies have been developed [38–41]. Even the models that have focused on HPV vaccines have not taken into account multiple types. In this paper, we compare the effects of vaccination on the prevalence of HPV types when there are synergistic versus antagonistic interactions among two HPV types. The vaccination model developed and analyzed here is based on the susceptible-infective-removed (SIR) compartmental structure [42] and features an imperfect vaccine that can protect against one or both HPV types. The model predicts three types of equilibria: elimination of infection with both types (infectionfree equilibrium), persistence of only one type (boundary equilibrium), and persistence of both types (coexistence equilibrium). We provide the conditions for the existence and local stability of the infection-free and boundary equilibria and existence of the coexistence equilibrium. In many cases, decades are required for an epidemiological system to approach equilibrium. Exploring the temporal dynamics of the system as it moves towards its equilibrium is crucial to assessing the public health impact of a vaccination program. Thus, using a realistic set of parameters from the clinical literature and fitting the model to observed prevalence levels of HPV types, we analyze the temporal dynamics and assess the impact of vaccination over the short and longer term.

2. Model description A mathematical model is developed to assess the changes in the distribution of HPV types following mass vaccination. The model analyzes the transmission of two HPV types and the effects of interactions among these HPV types. Types may interact directly in terms of synergy

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

91

or competition for colonization of the cervical epithelium, or indirectly via immunological memory or enhanced host susceptibility from previous infections. These distinctions are important when the duration of infection is relatively long, as is the case for HPV. The model stratifies the homogeneously mixing sexually active population into several classes: susceptible to both types (X), infected with type i and susceptible to the other type (Yi), infected with type 1 and type 2 (Y12), infected with type i and immune to infection with the other type because of recovery from previous infection (Ui), immune against infection with type i because of recovery from infection, but susceptible to infection with the other type (Zi), and immune against infection with either type because of recovery from infection (Z12) (Fig. 1). In the absence of vaccination, the model assumes that humans enter the sexually active population (N) into the susceptible (X) compartment at rate K and leave all compartments at rate l. Upon infection with HPV type i, the host moves into the Yi compartment. Susceptible individuals are infected with HPV type i at a rate ki. This force of infection ki depends on the probability of HPV transmission from an individual infected with type i to a susceptible, bi, the rate of partner acquisition, r, and prevalence of type i in the community. A host infected with type i can become concurrently infected with the other type and move into compartment Y12. The rate at which a host already infected with type i becomes coinfected with type j is cjkj, where kj is the rate at which a susceptible host is infected with type j, and cj captures the interaction between the HPV type i already established within the host and type j. We compare cases where infection with type 1 may reduce the probability of infection with type 2 through competitive interactions (i.e., c2 < 1), or enhance the probability of type 2 infection though synergistic interactions (i.e., c2 > 1). Infection with type i is cleared at rate ci and co-infected individuals clear type i at rate bic12, where b12 + b1 + b2 = 1. We assumed that if a host clears infection with type i, the host acquires life-long protection against that type [39]. Also, the degree of susceptibility of the host to the other type is affected (either increased or decreased) by a factor ai. Thus, when ai is zero, there is complete cross immunity whereas ai = 1 represents independence. If ai > 1, previous infection with one type facilitates infection with the other type. It is assumed that the vaccine targeted against type i offers a degree of protection (1  ki) against type i, in addition to (1  kj) protection against type j, where 0 6 ki 6 1. In our compartmental model, a fraction, f, of hosts is vaccinated and moves into compartment V. The vaccineinduced immunity is assumed, for simplicity, to be life-long. If a vaccinated host becomes infected with type i, she/he moves into compartment Wi, and then into compartment Qi, upon recovery. A host moves to compartment Pi if infected with type i, but immune against the other type. The ordinary differential equations that represent our compartmental model are dX =dt ¼ Kð1  f Þ  k1 X  k2 X  k12 X  lX ; dV =dt ¼ Kf  k 1 k1 V  k 2 k2 V  k 1 k 2 k12 V  lV ; dY i =dt ¼ ki X  cj kj Y i  ci Y i  lY i ; dY 12 =dt ¼ k12 X þ c1 k1 Y 2 þ c2 k2 Y 1  c12 Y 12  lY 12 ; dZ i =dt ¼ ci Y i  aj kj Z i  lZ i ; dU i =dt ¼ ai ki Z j þ bj c12 Y 12  ci U i  lU i ;

92

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Fig. 1. Schematical presentation of the model in which humans enter into the susceptible (X) compartment at rate K and leave all compartments at rate l. A susceptible host may be infected with either or both HPV types. Susceptible individuals acquire type i infection at rate ki. A host infected with type i can also be infected with the other type and move into compartment (Y12). This co-infection occurs at rate cj times the rate at which a susceptible individual would be infected with the same type. An infected individual clears infection with type i at rate ci. Co-infection is cleared at rate bic12. A fraction of susceptible f are vaccinated and move into compartment V. The vaccine provides complete, partial, or no protection against type i at rate 1  ki. After being infected, a vaccinated person moves into compartment Wi. If an infection with type i is cleared, a vaccinated person moves into compartment Qi and remain there until being infected with the other type and move to compartment Pj. Upon clearance of both types, the person moves to compartment Z12.

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

93

dW i =dt ¼ k i ki V  k j cj kj W i  ci W i  lW i ; dW 12 =dt ¼ k 1 c1 k1 W 2 þ k 2 c2 k2 W 1 þ k 1 k 2 k12 V  c12 W 12  lW 12 ; dQi =dt ¼ ci W i  k j aj kj Qi  lQi ; dP i =dt ¼ k i ai ki Qj þ bj c12 W 12  ci P i  lP i ; dZ 12 =dt ¼

2 X

ci ðU i þ P i Þ þ b12 c12 ðY 12 þ W 12 Þ  lZ 12 ;

ð1Þ

i¼1

ki ¼ rbi ðY i þ Y 12 þ U i þ W i þ W 12 þ P i Þ=N ; k12 ¼ rb12 ðY 12 þ W 12 Þ=N ; N ¼X þV þ

2 X ðY i þ W i þ U i þ P i þ Qi þ Z i Þ þ Y 12 þ W 12 þ Z 12 ;

i; j ¼ 1; 2; i 6¼ j;

i¼1

where b12 + b1 + b2 = 1. Definitions of the variables and parameters of the model are described in Table 1. Since at equilibrium N = K/l, we only need to analyze the limiting system where N in (1) is replaced by its equilibrium value [43,44]. Given this, the domain of biological interest is ( D¼

ðX ; V ; Y i ; W i ; U i ; P i ; Qi ; Z i ; Y 12 ; W 12 ; Z 12 Þ 2 R17 þ : X þV

) 2 X ðY i þ W i þ U i þ P i þ Qi þ Z i Þ þ Y 12 þ W 12 þ Z 12 6 K=l . þ i¼1

We will consider only the dynamics of the flow generated by (1) in this positively invariant domain D. It can be shown that unique solutions exist in D for all positive time. Thus, the model is epidemiologically and mathematically well posed [42]. We start by obtaining analytic results for some special cases and then investigate the impact of mass vaccination for more general cases. 3. Sequential acquisition of types and uniform clearance rates We first consider the simplified case where coinfections are acquired sequentially and HPV clearance rate is the same for single and multiple infections and is constant across types. Consequently, the following parameter restrictions are imposed: b12 = b1 = b2 = 0, b12 = 1, ci = c for all i = 1, 2, 12. 3.1. No vaccine The system given by Eq. (1) has three kinds of equilibria: the infection-free equilibrium, the endemic one-type (boundary) equilibrium, and the endemic (coexistence) equilibrium where both types coexist. A standard result in epidemiologic theory is that a pathogen strain or type i can invade the host population in the absence of other types, if its R0i is greater than one [45]. The

94

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Table 1 Description of variables and parameters Symbol

Description

Subscripts i, j

HPV types (1 = type 1, 2 = type 2, 12 = type 1 and 2)

Variables X(t) V(t) Yi(t) Y12(t) Ui(t) Wi(t) W12(t) Pi(t) Qi(t) Zi(t) Z12(t) ki(t) k12(t) N(t)

Susceptible to both types Vaccinated against type 1, susceptible to both types Infected with type i, susceptible to the other type Infected with both types Infected with type i, immune against the other type Vaccinated, infected with type i, susceptible to the other type Vaccinated, infected with both types Vaccinated, infected with type i, immune against the other type Vaccinated, immune against type i, susceptible to the other type Immune against type i, susceptible to the other type Immune against both types Force of infection with type i Force of infection with type 1 and 2 Size of host population

Demographic parameters K l r

New recruits into the sexually active population Death rate or duration in the sexually active population Partner acquisition rate

Biological parameters ci c12 bi b12 bi fi h(r) g ci ai

Recovery from infection with type i Recovery from infection with type 1 and 2 Relative risk of recovery from infection with type i, given coinfection Relative risk of recovery from infection with type 1 and 2 Transmission probability of type i per partner Transmission probability of type i per sexual contact Average number of contact per partnership A parameter that controls how h(r) is related to r Relative risk of infection with type i given current infection with the other type Relative risk of infection with type i given prior infection with the other type

Vaccine parameters f ki

Percentage of new recruits vaccinated, i.e., vaccination coverage Relative risk of infection with type i of a vaccinated person

basic reproduction number R0i is defined as the number of secondary infections generated by the first infection in a completely susceptible host population [46]. R0i is obtained from analyzing the stability properties of the infection-free equilibrium. Thus R0i ¼

rbi . lþc

ð2Þ

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

95

If R0i < 1 for all types, the infection-free equilibrium is locally asymptotically stable. If R0i > 1 for any type, the infection-free equilibrium is unstable. If R0i > 1 for only one type, this type will persist at equilibrium. If each R0i is greater than one and  1 ðca1 þ lc1 Þ R01 ðca2 þ lc2 Þ 1þ < <1þ ð3Þ ðR02  1Þ ðR01  1Þ; ðl þ cÞ ðl þ cÞ R02 the two types will coexist. The proofs are given in Appendix A. Similar expressions were derived in [26,36]. 3.2. Type-specific, perfect vaccine A type 1-specific vaccine with 100% efficacy is represented in this model by the case where k1 = 0 and k2 = 1. We assume R0i > 1 for i = 1, 2. 3.2.1. Competition only In the case where there is competition between the two types but no cross-immunity nor increased susceptibility to infection from prior infection, defined by ai = 1, the conditions for coexistence become   1 ðc þ lc1 Þ R01 ðc þ lc2 Þ f lð1  c2 ÞR01 ðR02  1Þ ð1  f Þ ðR01  1Þ þ . ð4Þ < <1þ 1þ ðl þ cÞ ðl þ cÞ R02 lþc The condition for elimination of the vaccine type is given by the inequality  1 R02 ðc þ lc1 Þ f >1 1þ ðR02  1Þ . ðl þ cÞ R01 This critical fraction of vaccination coverage required to eliminate the target type is less than that required without competitive interactions. Thus, the case where c1 = 1 gives a required coverage for elimination of 1  1/R01, which is greater than the critical coverage when c1 < 1. The critical vaccine coverage required to eliminate the target type increases with R01 and c1, and decreases with R02, c and l. That is, the critical coverage required to eliminate type 1 increases with R0 of type 1 and decreases with R0 of type 2, the degree of competition between the two types, and the duration of infection. 3.2.2. Competition cross-immunity or increased susceptibility We now consider the general case where there are either competitive or synergistic interactions between the two types, i.e. ci 5 1, and prior infection with one type induces cross-immunity or increases susceptibility to infection with the other type, i.e. ai 5 1. In this case, the conditions for the two types to coexist are given by the inequalities   1 ðca1 þ lc1 Þ R01 < ðR02  1Þ ð1  f Þ 1þ ðl þ cÞ R02  ðca2 þ lc2 Þ lc2 þ ca2 <1þ ðR01  1Þ þ 1  fR01 . ðl þ cÞ lþc ð5Þ

96

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

These conditions revert to those given by Eq. (3) for coexistence when there is no vaccine (f = 0). It can easily be shown that the conditions in Eq. (5) are automatically satisfied if ci > 1 and ai > 1, i = 1, 2. That is, the two types coexist whenever there are synergies. The vaccine type is eliminated if coverage satisfies the inequality  1 R02 ðca1 þ lc1 Þ ðR02  1Þ . 1þ 1>f >1 ðl þ cÞ R01 If prior infection with type 2 provides some protection against type 1 (i.e., a1 < 1), this critical fraction required to eliminate the target type is less than that in the absence of cross-immunity (i.e., a1 = 1). However, the critical vaccine coverage required to eliminate the target type is higher if prior infection with type 2 enhances the host susceptibility to infection with type 1 (i.e., a1 > 1). Fig. 2 shows the range of R0 values of type 1 and type 2 over which coexistence is possible in the absence versus presence of mass vaccination, for different levels of competition ci and cross-immunity ai. When there is complete competitive exclusion ci = 0 and perfect cross-immunity ai = 0, the type with the lower R0 is driven to extinction and hence coexistence is not possible in the absence of vaccination, unless R0 values are identical (Fig. 2(a)). Coexistence is possible if competition is less than complete ci > 0 and cross-immunity is imperfect ai > 0 (Fig. 2(c) and (e)). The region of coexistence expands as the values of ci and ai increase, corresponding to reduced competition and cross-immunity. If ci P 1 and ai P 1, the region of coexistence becomes the first quadrant with (1, 1) as its origin. Mass vaccination shifts the coexistence region to the right, indicating the diminished importance of the type targeted by the vaccine (Fig. 2). We found that a type-specific vaccine can lead to the reintroduction of a type that was previously suppressed by intense competition and crossimmunity (Fig. 2(b)). Fig. 2(d) and (f) shows the effects of mass vaccination on coexistence as levels of competition and cross-immunity fall. Fig. 3 shows the change of prevalence of type 1, type 2, and coinfection as vaccination coverage expands. Without interaction between HPV types, vaccination reduces prevalence of type 1 only (Fig. 3(a)). If there are synergistic interactions, prevalence of both types fall with more coverage (Fig. 3(b)). Competition and cross immunity mitigate the impact of vaccination and may increase prevalence of the non-target type (Fig. 3(c)). An HPV type may reemerge as a result of vaccination (Fig. 3(d)). 3.3. Bivalent, imperfect vaccine We now consider the case of a vaccine that is imperfect in its protection against one or both of the circulating types. As before, we measure the degree of protection against type i by the quantity 1  ki, where ki = 0, 0 < ki < 1, and ki = 1 represent full, partial, or no protection against type i, respectively. We first look at competitive exclusion where only type 1 prevails. Solving system (1) for k1, assuming U1 = U2 = P1 = P2 = W2 = W12 = Z2 = k2 = 0, we obtain a cubic equation in k1 k1 fk 1 k21 þ K½1  k 1 ðR01  1Þ k1  K2 ð½1  f ð1  k 1 Þ R01  1Þg ¼ 0.

ð6Þ

We have already analyzed the infection-free equilibrium k1 = 0. Thus, we assume k1 5 0. We are left with the quadratic equation within the large curly brackets. This equation has one positive solution if [1  f(1  k1)]R01 > 1, none whenever [1  f(1  k1)]R01 < 1 and 1 > k1(R01  1), and

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

97

Fig. 2. The range of R0s of type 1 and type 2 over which type 2 dominates (gray area), type 1 dominates (white area), or coexistence is possible (dark area) in the absence and presence of mass vaccination (f = 0.4), for different levels of competition (ci) and cross-reaction (ai). Vaccine protects against type 1 (k1 = 0, k2 = 1). Other parameters: l = 0.02,c = 0.5. (a) No vaccine: ai = 0, ci = 0; (b) with vaccine: ai = 0, ci = 0; (c) no vaccine: ai = 0.1, ci = 0.1; (d) with vaccine: ai = 0.1, ci = 0.1; (e) no vaccine: ai = 0.5, ci = 0.5 and (f) with vaccine: ai = 0.1, ci = 0.1

98

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Fig. 3. Prevalence of type 1 (dotted line), type 2 (solid line), and coinfection (gray line) against mass vaccination coverage (f), for various values of competition (ci) and cross-reaction (ai). Vaccine protects against type 1 (k1 = 0, k2 = 1). Other parameters: R01 = 4.0, R02 = 4.0, K = 1, l = 0.02, c = 0.2, except (d) where R02 = 2.25.

two if [1  f(1  k1)]R01 < 1 and 1 < k1(R01  1). However, the last possibility of backward bifurcation can easily be ruled out by noting that [1  f(1  k1)]R01 < 1 and 1 < k1(R01  1) leads to a contradiction. Therefore, if this boundary equilibrium exists, it is unique. The boundary equilibrium where only type i prevails can be explicitly stated as qj ¼ wj ¼ y j ¼ zj ¼ uj ¼ ui ¼ pj ¼ pi ¼ w12 ¼ y 12 ¼ 0; vi ¼

1  k i þ k i R0i  Ui ; 2ð1  k i Þk i R0i

wi ¼

l½ð1  k i Þð2fk i R0i  1Þ  k i R0i þ Ui ; 2ðl þ cÞð1  k i ÞR0i

zi ¼

xi ¼

1  k i ð1 þ R0i Þ þ Ui ; 2ð1  k i Þk i R0i

cfR0i  ð1  k i Þ½1  ð1  2f ÞR0i  k i R0i  Ui g 2ðl þ cÞð1  k i ÞR0i

qi ¼

cwi ; l

yi ¼

lzi ; c ð7Þ

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

99

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with Ui ¼ ½1 þ k i ðR0i  1Þ 2  4f ð1  k i Þk i R0i ; i; j ¼ 1; 2; i 6¼ j, where vi, xi, qi, wi, zi, pi, and yi are the respective fractions of the host population who are vaccinated (v = V/N), susceptible, vaccinated and immune against type i but susceptible to the other type, vaccinated but infected with type i and susceptible to the other type, immune against type i but susceptible to the other type, vaccinated but infected with type i and immune to the other type, and infected with type i and susceptible to the other type. Note that we added a subscript i to v and x to designate the boundary equilibrium where only type i prevails. Inspecting the above formula one can verify that a unique positive type i-only equilibrium exists if and only if R0i[1  f(1  ki)] > 1. The condition for local stability of type i-only equilibria is 1/R0j > xi + kj(vi + cjwi + aj qi) + cjyi + aj zi, where i 5 j. The equilibrium is unstable if this condition is violated (Appendix B). As shown above, the conditions for the existence of coexistence equilibrium can be completely characterized by the stability conditions of the boundary equilibria (see also, e.g., Ref. [43]). Thus, the two types will coexist if and only if R0i[1  f(1  ki)] > 1 and 1 R01 < R02 ½x2 þ k 1 ðv2 þ c1 w2 þ a1 q2 Þ þ c1 y 2 þ a1 z2 R02 < R01 ½x1 þ k 2 ðv1 þ c2 w1 þ a2 q1 Þ þ c2 y 1 þ a2 z1 . ð8Þ The conditions in Eq. (8) differ from those in the previous model, Eq. (5), by the terms involving the vaccine efficacy parameters k1 and k2. Note that because the vaccine is imperfect, elimination of the target type may not be possible even if the entire population is vaccinated (i.e., f = 1). We now assess the impact of a bivalent vaccine on prevalence of the two types under different forms of interaction between HPV types. 3.3.1. Complete independence If the types are transmitted independently, the impact of a bivalent vaccine with complete protection against type 1 and partial protection against type 2 can be assessed analytically. Setting the right-hand side of Eq. (1) to zero and solving while imposing the restrictions c1 = c2 = a1 = a2 = 1, yields k1 ¼ l½R01 ð1  f Þ  1 ; 2

k2 ¼ l½k 2 ðR02  1Þ  1 =ð2k 2 Þ þ f½k 2 ðR02  1Þ þ 1  4fk 2 ð1  k 2 ÞR02 g

1=2

.

ð9Þ

The equilibrium prevalence of type i is given by ki/bi. Thus the introduction of mass vaccination (f > 0) will reduce both the prevalence of the vaccine target type (type 1) and the prevalence of the non-target type if the vaccine offers some protection against type 2 (0 6 k2 < 1). This can be seen from the above formula (9), where k1 is inversely related to vaccine coverage f. If the vaccine does not offer protection against type 2 (k2 = 1), equilibrium prevalence of type 2 will be independent of vaccine coverage, as can be seen from Eq. (9). 3.3.2. Competition and cross-immunity If current or prior infection with one HPV type inhibits infection with the other type (i.e. interactions are antagonistic), the impact of a bivalent vaccine with complete protection against type 1

100

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

and partial protection against type 2 can also be assessed analytically. Solving for equilibrium values, while imposing the restrictions c1 = c2 = a1 = a2 = 0, yields 8 > > <

ðl½R01 ð1  f Þ  1 ; 0Þ    ðk1 ; k2 Þ ¼   > 1 fR201 1 fR01 R02 > : l R01  1 þ  ;l  þ k 2 R01  R02 k 2 R01  R02

R01  R02 ; k 2 R01 R02 R01  R02 . f > k 2 R01 R02

f 6

The two types cannot coexist in the absence of vaccination (i.e., f = 0). The HPV type with the highest R0 will eventually drive the other extinct. We assume that the surviving type is the vaccine targeted type 1. That is, R01 > R02. It can be seen that prevalence of the target HPV type falls and prevalence of the non-target type rises as a result of mass vaccination. In this case, mass vaccination has the potential to facilitate the re-emergence of an HPV type. The increase in the prevalence of one type is exactly equal to the increase in prevalence of the other type. Thus, oðk1 =b1 Þ=of ¼ lðc þ lÞR01 =ðR01  R02 Þ ¼ oðk2 =b2 Þ=of . 3.3.3. Enhanced susceptibilities and weakened immunity from infection The situation where current or prior infection with one HPV type weakens the immune system and facilitates infection with the other type occurs when the parameters c1, c2, a1, and a2 are greater than one. In the extreme case where all these parameters approach infinity, the impact of a bivalent vaccine with complete protection against type 1 and partial or full protection against type 2 can be assessed analytically as follows. First, we reduce the equilibrium system to two equations in k1 and k2 by eliminating all other variables through substitution. Second, we obtain a pre-vaccine equilibrium by solving the two equations for k1 and k2 while imposing the restrictions c1 = c2 = a1 = a2 ! 1 and f = 0. This yields  1 ; i ¼ 1; 2. ki ¼ lR0i 1  R01 þ R02 Third, to evaluate the responses of the reduced system to a change in f, we differentiate with respect to f and evaluate the result at the pre-vaccine equilibrium and f = 0. This gives the system of two equations in ok1/of and ok2/of that can be solved as   ok1 k 2 R202 ¼ lR01 1  ; of ðR01 þ R02 Þ½R01 þ ð1  k 2 þ k 2 R01 ÞR02 þ k 2 R202

ð10Þ ok2 lR02 ½R201 þ 2R01 R02 þ ð1  k 2 ÞR202

. ¼ of ðR01 þ R02 Þ½R01 þ ð1  k 2 þ k 2 R01 ÞR02 þ k 2 R202

The signs of both derivatives in (10) are negative since R01, R02 > 1 and k2 6 1. Therefore, when current or previous infection with one type facilitates infection with the other type, mass vaccination reduces prevalence of both types. When a bivalent vaccine offers only partial protection against types 1 and 2 (1 > ki > 0), the effect of mass vaccination on the prevalence of type 2 can be assessed by taking the derivative of the upper limit of the conditions for a coexistence equilibrium with respect to f and evaluating the sign

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

101

of this derivative. The effect of vaccination on the upper limit depends on the results of interactions between fractions in different compartments. As coverage expands, the susceptible fraction (x1), as defined in Eq. (7), decreases, prevalence of type 1 (y1) falls, and the fraction recovered from infection with type 1 (z1) also falls. At the same time, vaccination increases the fractions vaccinated (v1) and infected with type 1 and susceptible to type 2 (w1). The magnitude of these different effects of vaccination will depend on clinical or efficacy parameters. For example, the fall in the fraction of recovered from infection with type 1 depends, among other things, on the cross-immunity parameter a2. Formally, whether the upper limit of the expression in Eq. (8) increases or decreases as vaccination coverage f expands depends on the sign of the following expression: ðk 2  k 1 Þðl þ c  lc2  ca2 Þ  ð1  k 2 Þðlc2 þ ca2 ÞU1 ; where U1 is defined in Eq. (7) and in Appendix B. If this expression is positive, prevalence of type 2 will increase. Evaluating this expression at f = 0 (i.e., equilibrium prior to mass vaccination) gives (k2  k1)(l + c  lc2  ca2)  (1  k2)(lc2 + ca2)[1 + k1(R01  1). Note that if the vaccine is type specific (i.e., k2 = 1) and provides complete protection against type 1 (i.e., k1 = 0), prevalence of type 2 will fall only if l + c  lc2  ca2 < 0, confirming our earlier result. Fig. 4(a)–(f) shows the change in prevalence of the HPV types as vaccination coverage with a bivalent vaccine expands. If there is no interaction between types, vaccination can result in a reduction in prevalence of both HPV types (Fig. 4(a) and (b)). Under moderate levels of competition and cross immunity, prevalence of type 1 falls with more coverage whereas that of type 2 increases with low protection and decrease if protection is high. Prevalence of type 2 increases with vaccination if competition and cross-immunity are strong. The increase in type 2 prevalence can be mitigated if the vaccine provides strong protection against type 2 (Fig. 4(c) and (d)). If competition is weak and interactions are synergistic, prevalence of both types decline with increasing coverage. More pronounced reductions in prevalence can occur with increasing intensity of synergistic interactions (Fig. 4(e) and (f)). However, higher coverage may be required to achieve elimination. Because of the synergistic interactions among the two types, the required coverage to achieve elimination of type 1 is lower the more effective the vaccine in protecting against type 2.

4. Different clearance rates So far we have assumed that different HPV types have the same rate of clearance and if a host is coinfected, both HPV types are cleared simultaneously. However, empirical studies suggest that the clearance rates vary among HPV types, with high-risk types (HPV 16 in particular) being more persistent than low-risk types [47]. No significant difference between clearance rates of multiple infections or single infections has been detected. In this section, the previous model will be extended to account for differing clearance rates between the two types. Thus, the clearance rate of type 1, type 2, and co-infection are denoted by c1, c2, and c12, respectively. The model allows total clearance of coinfection at rate b12c12, type 1 clearance at rate b1c12, and type 2 clearance at rate b2c12, where b12 + b1 + b2 = 1. To simplify the analysis, we will assume that clearance of one type is independent of the presence of the other type, as is consistent with clinical data (e.g., [15]). Consequently, b1c12 = c1(1  c2), b2c12 = c2(1  c1), b12c12 = c1c2 and c12 = c1 + c2  c1c2. The

102

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117 (a) ai = 1, ci = 1

0.08

(b) ai = 1, ci = 1

0.08 Type 2 0.06 Prevalence

Prevalence

0.06

Type 1

0.04

0.02

Type 2

0.02 Coinfection

Coinfection

0.2

0.4 0.6 0.8 Coverage against type 1 ( f )

1.0

0.2

(c) ai = 0.2, ci = 0.2

0.08

0.06 Type 1

Prevalence

Type 2

0.04

0.02

0.4 0.6 0.8 Coverage against type 1 ( f )

Type 1

0.04 Type 2

0.02 Coinfection 0.2

Coinfection

0.4 0.6 0.8 Coverage against type 1 ( f )

1.0

0.2

0.4

0.6

0.8

1.0

Coverage against type 1 ( f )

(e) ai = 4, ci = 4

0.08

1.0

(d) ai = 0.2, ci = 0.2

0.08

0.06 Prevalence

Type 1

0.04

(f) ai = 4, ci = 4

0.08 Type 2

Type 1

0.06 Prevalence

Prevalence

0.06

0.04 Coinfection 0.02

Type 2 Type 1

0.04 Coinfection 0.02

0.2

0.4 0.6 0.8 Coverage against type 1 ( f )

1.0

0.2

0.4

0.6

0.8

1.0

Coverage against type 1 ( f )

Fig. 4. Prevalence of type 1 (dotted line), type 2 (solid line), and coinfection (gray line) against mass vaccination coverage (f), for various values of competition (ci), cross-reaction (ai), and degrees of protection against both types (ki). Vaccine parameters: k1 = 0.1, k2 = 0.9 (left column) and k1 = 0.1, k2 = 0.5 (right column). Other parameters: R01 = 4.0, R02 = 4.0, K = 1, l = 0.02, c = 0.2, except (c) and (d) where R02 = 2.25.

rest of the model remains the same. With these changes, the boundary equilibria and the conditions for coexistence are the same as those given in Eqs. (7) and (8), respectively, with one minor

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

103

change: c is simply replaced by ci. The upper limit of coexistence expression increases (decreases) as vaccination coverage f expands if the following expression is positive (negative): ð11Þ ðk 2  k 1 Þðl þ c1  lc2  c1 a2 Þ  ð1  k 2 Þðlc2 þ c1 a2 ÞU1 ; where U1 is defined in Appendix B. Similarly, the lower limit of coexistence expression increases (decreases) as vaccination coverage f expands if the following expression is positive (negative): ðk 1  k 2 Þðl þ c2  lc1  c2 a1 Þ  ð1  k 1 Þðlc1 þ c2 a1 ÞU2 .

ð12Þ

These expressions are obtained by differentiating the upper and lower limits with respect to f. To ease the interpretation of the above formulae, it is useful to define the following quantities: D1 = 1/(l + c1), D2 = 1/(l + c2), and L = 1/l, where D1 and D2 denotes duration of infectiousness for type 1 and type 2, respectively. L measures duration in the sexually active population. Using these definitions, we can rewrite Eqs. (11) and (12) as  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ D1 =L ð13Þ  1  ð1  k 2 Þ 1 þ k 1 ðR01  1Þ2  4f ð1  k 1 Þk 1 R01 ; ðk 2  k 1 Þ a2 þ c2 D1 =L  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ D2 =L ðk 1  k 2 Þ ð14Þ  1  ð1  k 1 Þ 1 þ k 2 ðR02  1Þ2  4f ð1  k 2 Þk 2 R02 . a1 þ c1 D2 =L Other things being equal, prevalence of the non-target type (type 2) falls after mass vaccination, the more protective the vaccine against the non-target type (lower k2), the less type 1 competitively excludes type 2 (higher c2), and the less cross-immunity to infection with type 2 is provided by previous infection with type 1 (higher a2). The persistence of the target vaccine type (higher D1) also affects prevalence of the non-target type, but the direction of change depends on the relative magnitudes of the parameters a2 and c2. The prevalence of the target vaccine type declines as vaccine coverage increases, provided that the efficacy against the target type is more than that against the non-target type (i.e., k1 < k2). This can be seen by setting f = 0 in Eq. (14) and observing that the sign of the resulting expression is always negative. 5. Temporal dynamics of the full model with vaccination So far we have assumed that coinfection occurs only through sequential infections. Appendix B considers the case where coinfections can occur simultaneously (k12 5 0). Also, most of the results obtained are for long-term outcomes that may take centuries to be realized. In this section, we explore temporal dynamics of the full model after mass vaccination under varying assumptions about antagonism and synergies among types. This is achieved by numerical analysis of the model to determine the impact of vaccination against type 1 on prevalence of the two types over time. 5.1. Parameters Estimates of parameters relevant to the model exhibit considerable variability in the literature. Here we outline the estimates employed and the sources of those estimates. The host removal rate l consists of two components: natural death rate and the rate at which individuals leave the sexually active population. Assuming a population of young adults that can live on average for 50

104

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

years and remain sexually active for 40 years, gives a removal rate of l = 0.02 + 0.025 = 0.045. To keep the total size of the sexually active population constant, we set K = l = 0.045, so that entry into and exist from the population are balanced. A few studies have estimated clearance rates ci of HPV types. Ho et al. [5] estimate a median duration of infection with HPV 16 and HPV 18 of 11 and 12 months, respectively. Two studies have estimated the median time to clearance for women infected with HPV 16 at 8.5 months [48,49]. However, the data reported in Liaw et al. [14] indicate persistence of HPV 16 infections for longer than 24 months. In general, most studies agree that infections with high-risk, oncogenic HPV types (especially the vaccine types 16 and 18) are more persistent than low-risk types (e.g., Refs. [50,51]). For example, the adjusted rate ratio of HPV 16 clearance was 0.48, relative to lowrisk types. Richardson et al. [47] provided estimates of mean time to loss of an incident infection for various HPV types. We use their estimates for HPV 16 and HPV 31 of 18 and 15 months, respectively. The ubiquitousness of HPV infection suggests high transmissibility. Frega et al. [52] reported a transmission probability of 0.53–0.55 per partner and Nebesio et al. [53] up to 0.60 per sexual contact. However, none of these two estimates was based on prospective studies of HPV-discordant sex partners. To estimate the probability of transmission per partner, we follow the approach outlined in Hyman et al. [54]: bi = 1  (1  fi)h(r), where fi is the transmission probability per contact, and the average number of contact per partnership h(r) is given by h(r) = 104rg + 1. We set f2 = 0.02 and g = 1 which correspond to 2 contacts per week (105 contacts per year) and b2 = 0.88 for r = 1 per year. For r = 5, these correspond to 21.8 contacts per year and b2 = 0.36. Since vaccine types are more prevalent, we set f1 = 0.08. As in Hyman et al. [54], the baseline value of r is set here at five partners per year. The parameter values employed produce endemic prevalence levels within the range reported in the literature [47]. 5.2. Numerical simulations The numerical simulations demonstrate how mass vaccination disturbs the endemic equilibrium and generates temporal dynamics that may be protracted before converging on the new equilibrium. With weak competition and little cross immunity, the vaccine primarily affects the targeted type 1 with little long-term impact on the non-targeted type 2 (Fig. 5(a) and (b)). Vaccination generates an immediate decline in prevalence of type 1, followed by a period of low prevalence during which susceptible individuals gradually replenish. After a sufficient susceptible population has accumulated, an epidemic occurs, followed by a smaller one, and the system approaches equilibrium within 100 years (Fig. 5(a)). The timing and size of the post-vaccination epidemics depend on vaccination coverage. Higher coverage results in greater impact on prevalence, both in the short-term and long-term. For example, type 1 is eliminated in less than 20 years at 90% coverage (Fig. 5(a)). Together with intense competition and cross immunity, high vaccination coverage produce complex dynamics (Fig. 5(a) and (b)). A previously eliminated type 2 reemerges following vaccination against the competing type 1. The higher the coverage, the sooner type 2 reemerges. For example, whereas type 2 reemerges in about 110 years when coverage is 90%, it takes close to 400 years for it to reappear when coverage is 45%. Once it reemerges, type 2 causes a large epidemic followed by smaller ones. The

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

105

Fig. 5. Prevalence of type 1 and 2 and time since mass vaccination for various values of competition (ci), cross-reaction (ai), and coverage (f). Vaccine parameters: k1 = 0.05, k2 = 0.95. Other parameters: K = 0.045, l = 0.045, c1 = 1/1.5, c2 = 1/1.25, r = 5, f1 = 0.08, f2 = 0.02, g = 1.

106

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Fig. 5 (continued)

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Fig. 5 (continued)

107

108

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

Fig. 5 (continued)

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

109

magnitude and duration of the initial epidemic depends on vaccination coverage; higher coverage generates a larger but shorter epidemic (Fig. 5(d)). The impact of these type 2 epidemics on the prevalence of type 1 is evident from Fig. 5(c) and (d). Because of competition between the two types, the large increase in type 2 prevalence is accompanied by a sharp fall in type 1 infections (Fig. 5(c) and (d)). Fig. 5(e)–(h) shows the impact of vaccination given synergies between HPV types. The existence of synergistic interactions between the two types causes prevalence of type 2 to fall in response to vaccination. For example, elimination of type 1 is possible at 90% vaccination coverage when ai = ci = 1.1 (Fig. 5(e)). Under the same level of vaccine coverage, type 1 remains endemic, albeit at low levels, if ai = ci = 4 (Fig. 5(g)). At high vaccine coverage, synergistic interactions generate complicated dynamics, especially in the behavior of type 2 (Fig. 5(f) and (h)).

6. Discussion Prophylactic vaccination against HPV holds promise for the control of HPV, cervical cancer and other HPV-related diseases [55]. Given that there are at least 35 different HPV types that infect the genital mucosa, we investigated how mathematical models can help in understanding the interactions between types and in predicting the effectiveness of vaccination against a subset of types. We compared the impact of antagonistic versus synergistic interactions between two HPV types on their prevalence in the face of mass vaccination. We found that the prevalence of the type not included in the vaccine can increase following vaccination against a competing type. If, however, there are synergistic interactions between the two, there could be an extra benefit in terms of reduction in prevalence of HPV types not targeted by the vaccine. Epidemiologic studies on the natural history of cervical HPV infections have found that current or prior infection with some HPV types is associated with an increased risk of concurrent or sequential acquisition of other HPV types. Although these results are consistent with the hypothesis of synergistic interactions among HPV types, they do not offer conclusive evidence of its validity. Furthermore, no biological mechanism through which current or previous infection with one HPV type predisposes individuals to subsequent infection with another HPV type has been suggested. The observed associations may be due to the ensuing physical lesions from HPV infection, analogous to the role played by other STIs such herpes simplex virus or Neisseria gonorrhoeae in increasing susceptibility to human immunodeficiency virus infection. Other alternative explanations for this increased risk could be presence of confounding risk factors or host susceptibility to HPV infection in general. However, it should be noted that all these studies adjusted extensively for confounding variables (e.g., age and number of sexual partners) in the statistical analysis. The efficacy of a quadrivalent vaccine is currently being assessed in Phase III trials. The candidate vaccine targets HPV 6, 11, 16, 18, the former two of which are responsible for more than 90% of genital warts, and the latter two are found in over 70% of cervical cancers. If, as the few studies suggest, there are synergistic interactions between the types included in the vaccine and other types that are also carcinogenic, such as 31 and 45, there could be an extra benefit for the vaccine

110

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

in terms of reduction in prevalence of all these HPV types and associated cancers. Targeting HPV types, such as 16 and 18, that most predispose hosts to further infection should be most effective in reducing the prevalence of HPV infections. We found that even a control program based on a highly efficacious vaccine applied with widespread coverage must be continued for many years in order to achieve elimination of HPV types. Consequently, studying temporal dynamics are crucial to understanding the impact of mass vaccination both over the short-term and the long-term. We found that the temporal dynamics of the system as it moves towards its equilibrium can follow complex trajectories, including damped limit cycles of HPV type prevalence, particularly when the interactions between the types are strong. We have established the local stability of the disease-free and boundary equilibria and proved the existence of the infection-free, boundary, and coexistence equilibria. However, we have not investigated the local stability of the coexistence equilibrium or the global stability of any of these equilibria. This task is left for future research. The model could be elaborated to take into account complications, such as gender, age and heterogeneity of mixing between different sexual activity groups. Genetic polymorphisms within the host population could also affect the degree of infectiousness and rate of recovery from infection. Inclusion of various forms of heterogeneity is likely to affect the quantitative results of the model but may not change the results qualitatively. Additionally, breakthrough infections may be less infectious and less persistent than infections among unvaccinated persons. To simplify the analysis, the model also assumes exponential recovery times and doesnt take into account any stochastic elements that may be inherent in many epidemiological processes. Further, the model traces the effects of vaccination on HPV prevalence, but does not include the relationships between HPV infection and HPV-related diseases such as cervical cancer. Also, the model assumes that the vaccine provides protection that lasts over the entire lifetime of the vaccinated individual. Current candidates appear to be highly immunogenic with titers more than 40 times higher compared with those following natural infection [56]. However, it is not currently possible to assess whether immunization confers life-long protection without conducting a longitudinal study or defining an immunological response that strongly correlates with protection against HPV infection. The model was developed specifically for HPV, but the results extend to other pathogens. The model predictions may be applicable to multi-strain pathogen populations that are subject either to competitive interactions, such as influenza, or to synergistic interactions, such as dengue. Our analysis of a pathogen population based on multiple HPV types highlights the importance of considering the population-effects of mass vaccination over both the short-term and the long-term, particularly when the types circulate in a non-independent manner.

Acknowledgments The authors wish to thank Eliav Barr, John R. Cook, and Erik J. Dasbach for helpful comments and discussions.

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

111

Appendix A In this appendix we analyze the stability of the infection-free and boundary equilibria and provide the conditions for the existence of the boundary and coexistence equilibria in the presence of mass vaccination with a perfect type-specific vaccine. As discussed in the text, we impose the following restrictions on the parameters of the model: b12 = b1 = b2 = 0, b12 = 1, ci = c for all i = 1, 2, 12. A.1. Local stability of the infection-free equilibrium, type-specific vaccine To investigate the local stability properties of the system given in Eq. (1), we compute the eigenvalues of the Jacobian matrix evaluated at the infection-free equilibrium. Note that Z12 is an absorbing state; its behavior does not affect the rest of the system. This reduces our system to 16 differential equations. At the infection-free equilibrium, X = 1  f, V = f , and all other 14 variables are zero. Evaluating the Jacobian matrix at the infection-free equilibrium it is possible to calculate all the eigenvalues as l  c (of multiplicity 8), l (of multiplicity 6), (l + c)[1  (1  f)R01], (l + c)(1  R02). Thus, all the eigenvalues are negative and the infection-free equilibrium is locally asymptotically stable if (1  f)R01 < 1 and R02 < 1. If (1  f)R01 > 1 or R02 > 1, the Jacobian has positive eigenvalues and the infection-free equilibrium is unstable. In the absence of vaccination (f = 0), the condition for stability becomes: R01 < 1 and R02 < 1. A.2. Existence of boundary equilibria, type-specific vaccine In the presence of a perfect monovalent vaccine (k1 = 0 and k2 = 1), the relevant differential equations are for the variables X, V, Y1, Z1, U1, Y2, Z2, U2, W2, and Y12 (we write the Jacobian in this order. Since the behavior of these variables is independent of Q2, its equation is not included in this analysis). Denoting the fractions of the host population by lower case letters, type-1 only endemic equilibrium can be solved for explicitly as x1 ¼

1 ; R01

v1 ¼ f ;

l y 1 ¼ z1 ; c

z1 ¼ c

R01 ð1  f Þ  1 ; ðl þ cÞR01

u1 ¼ u2 ¼ y 2 ¼ z2 ¼ w2 ¼ y 12 ¼ 0.

Clearly, there exists a unique positive equilibrium of type 1 exclusively if and only if R01(1  f) > 1. Similarly, endemic equilibrium of type 2 exclusively can be obtained as: x2 ¼

1f ; R02

v2 ¼

f ; R02

l y 2 ¼ z2 ; c

z2 ¼

cðR02  1Þð1  f Þ ; ðl þ cÞR02

w2 ¼

lðR02  1Þf ; ðl þ cÞR02

u1 ¼ u2 ¼ y 1 ¼ z1 ¼ y 12 ¼ 0. It is clear that there exists a unique positive equilibrium of type 2 exclusively if and only if R02 > 1.

112

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

A.3. Local stability of the endemic one-type equilibrium, type-specific vaccine Evaluating the Jacobian at endemic equilibrium of type 1 exclusively yields  J1 J2 J¼ ; 0 J4 where

0

lð1  f ÞR01

0 c  l 0 c  l

0

1

C B C 0 l 0 0 0 0 B C B C B l½ð1  f ÞR  1 0 0 0 cþl 0 01 C B J1 ¼B C; C B 0 0 c l 0 0 C B B 0 0 0 0 c  l la1 ½ð1  f ÞR01  1 C A @ 0 0 0 0 0 l  la1 ½ð1  f ÞR01  1

0 ðc þ lÞR02 x1 ðc þ lÞR02 x1 ðc þ lÞR02 x1 ðc þ lÞðR02 x1  c1 R01 y 1  1Þ B ðc þ lÞa2 z1 R02 ðc þ lÞða2 z1 R02  1Þ ðc þ lÞa2 z1 R02 ðc þ lÞa2 z1 R02 B J4 ¼B @ ðc þ lÞv1 R02 ðc þ lÞv1 R02 ðc þ lÞðv1 R02  1Þ ðc þ lÞv1 R02 ðc þ lÞðc1 R01 þ c2 R02 Þy 1

ðc þ lÞc2 R02 y 1

ðc þ lÞc2 R02 y 1

1 C C C. A

ðc þ lÞðc2 R02 y 1  1Þ

The eigenvalues of J are given by those of J1 and J4. The eigenvalues of J1 are l (of multiplicity 2), c  l, l(1  f)R01, l  la1[(1  f)R01  1], and the roots of the quadratic equation /2 þ lð1  f ÞR01 / þ lðl þ cÞ½ð1  f ÞR01  1 ¼ 0. Clearly all eigenvalues of J1 have negative real parts if and only if (1  f)R01 > 1. For J4, the eigenvalues are the roots of the polynomial equation: /4 þ A1 /3 þ A2 /2 þ A3 / þ A4 ¼ 0; where A1 ¼ ðc þ lÞ½4 þ c1 R01 y 1  R02 ðv1 þ x1 þ c2 y 1 þ a2 z1 Þ ; A2 ¼ ðc þ lÞ2 ½6 þ 3c1 R01 y 1  R02 ðv1 þ x1 þ c2 y 1 þ a2 z1 Þð3 þ c1 y 1 R01 Þ ; A3 ¼ ðc þ lÞ3 ½4 þ 3c1 R01 y 1  R02 ðv1 þ x1 þ c2 y 1 þ a2 z1 Þð3 þ 2c1 y 1 R01 Þ ; 4

A4 ¼ ðc þ lÞ ð1 þ c1 R01 y 1 Þ½1  R02 ðv1 þ x1 þ c2 y 1 þ a2 z1 Þ . It can be shown that all coefficients are positive, A1A2 > A3, and A1 A2 A3 > A23 þ A21 A4 if R02(v1 + x1 + c2y1 + a2z1) < 1. Therefore, according to Routh–Hurwitz criteria, all the eigenvalues of J4 have negative real parts. Hence, the endemic equilibrium of type 1 exclusively is locally asymptotically stable if (1  f)R01 > 1 and R02(v1 + x1 + c2y1 + a2z1) < 1. If R02(v1 + x1 + c2y1 + a2z1) > 1, this equilibrium is unstable. Similarly, it can be shown that endemic equilibrium of type 2 exclusively is locally asymptotically stable if R02 > 1 and R01(x2 + c1y2 + a1z2) < 1 and is unstable if R01(x2 + c1y2 + a1z2) > 1, where x2 = (1  f)/R02, y2 = lz2/c, z2 = c(R02  1)(1  f)/[(l + c)R02].

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

113

Note that these conditions can be obtained directly by evaluating the determinant of J. This strategy will be employed when analyzing the stability of the full model. In terms of parameters of the original model, the conditions where both one-type equilibria are unstable and the two types can coexist are (1  f)R01 > 1, R02 > 1, and   1 ðlc1 þ ca1 Þ R01 ðR02  1Þ ð1  f Þ 1þ < ðl þ cÞ R02  ðlc2 þ ca2 Þ lc2 þ ca2 ðR01  1Þ þ 1  <1þ fR01 . ðl þ cÞ ðl þ cÞ The conditions for coexistence in the absence of vaccination can be obtained by setting f = 0. A.4. Existence of endemic two-types equilibrium, type-specific vaccine To investigate the conditions for the existence of an endemic equilibrium where the two types can coexist in the presence of a perfect monovalent vaccine, we equate the left-hand side of Eq. (1) to zero, and solve for k1 and k2. Thus, k1 ðk2 Þ ¼

lðR02  1Þ  k2 ; 1  lR02 p1 ðk2 Þ

k2 ðk1 Þ ¼

l½ð1  f ÞR01  1  k1 ; 1  lR01 p2 ðk1 Þ

where f lðc þ l þ a2 k2 Þ þ ½lc2 þ a2 ðc þ c2 k2 Þ ½lð1  f Þ þ k2

; ðl þ k2 Þðl þ a2 k2 Þðc þ l þ c2 k2 Þ ð1  f Þ½lc1 þ a1 ðc þ c1 k1 Þ

. p2 ðk1 Þ ¼ ðl þ a1 k1 Þðc þ l þ c1 k1 Þ

p1 ðk2 Þ ¼

These functions ki(Æ) define two hyperbolas over [0, (l + c)R0j]. The term pi(kj) is strictly monotone decreasing in kj for kj P 0, i 5 j. So, ki(kj) is strictly monotone decreasing in kj for kj 2 [0,(l + c)R0j]. Inspecting the equation for k2 reveals that k2 = l(l + c)[(1  f)R01  1]/ [l + c  (1  f)(ca1 + lc1)R01] when k1 = 0, k2 = 0 when k1 = l[(1  f)R01  1], and k2 ! 1 as k1 ! +1. If l + c(1  f)(ca1 + lc1)R01 > 0, k2(k1) intersects the x-axis at k1 = l[(1  f)R01  1] and the y-axis at k2 = l(l + c)[(1  f)R01  1]/[l + c  (1  f)(ca1 + lc1)R01] > 0. If l + c  (1  f)(ca1 + lc1)R01 < 0,k2(k1) intersects the x-axis at k1 = l[(1  f)R01  1] and goes to +1 (1) as k1 approaches k1 from above (below), where k1 : 1 ¼ lR01 p2 ðk1 Þ. Similarly, k1 = l(l + c)(R02  1)/{l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02} when k2 = 0, k1 = 0 when k2 = l(R02  1), and k1 ! 1 as k2 ! +1. If l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02 > 0, k1(k2) intersects the x-axis at k1 = l[(1  f)R01  1] and the y-axis at k1 = l(l + c)(R02  1)/ {l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02} > 0 and the y-axis at k2 = l(R02  1). If l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02 < 0, k1(k2) intersects the x-axis at k2 = l(R02  1) and goes to +1 (1) as k2 approaches k2 from above (below), where k2 : 1 ¼ lR02 p1 ðk2 Þ.

114

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

If l + c  (1  f)(ca1 + lc1)R01 > 0 and l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02 > 0, the two hyperbolas will intersects in the first quadrant if and only if l(l + c)[(1  f)R01  1]/ [l + c  (1  f)(ca1 + lc1)R01] > l(R02  1) and l(l + c)(R02  1)/{l + c  [f(l + c) + (1  f) (ca2 + lc2)]R02} > l[(1  f)R01  1]. If l(l + c)[(1  f)R01  1]/[l + c  (1  f)(ca1 + lc1)R01] < 0 < l(R02  1) and l(l + c)(R02  1)/{l + c  [f(l + c) + (1  f)(ca2 + lc2)]R02} < 0 < l[(1  f) R01  1], the two hyperbolas will intersect in the first quadrant. The remaining two cases produce the same conditions. Combining these conditions for the existence of endemic equilibrium at which the two types coexist, we obtain (1  f)R01 > 1, R02 > 1, and   1 ðlc1 þ ca1 Þ R01 ðR02  1Þ ð1  f Þ 1þ < ðl þ cÞ R02  ðlc2 þ ca2 Þ lc2 þ ca2 ðR01  1Þ þ 1  <1þ fR01 . ðl þ cÞ ðl þ cÞ

Appendix B The models discussed in the text assume that coinfection with the two circulating types occurs in a sequential manner, where one type is acquired first and only after it has been established in the host can coinfection with another type occur. This assumption might be restrictive given the frequent occurrence of coinfection. In this appendix we analyze the stability of the infection-free and boundary equilibria and provide the conditions for the existence of the boundary and coexistence equilibria in the presence of an imperfect bivalent vaccine for the full model after relaxing this assumption (i.e., we allow b12 5 0). We also allow HPV clearance rates to vary by type but we will continue to assume that clearance of one type is independent of the presence of the other type. Thus, b1c12 = c1(1  c2), b2c12 = c2(1  c1), b12c12 = c1c2 and c12 = c1 + c2  c1c2. B.1. Local stability of the infection-free equilibrium, bivalent vaccine, full model Evaluating the Jacobian matrix at the infection-free equilibrium shows that the infection-free equilibrium is locally asymptotically stable if (1  f + fk1)R01 < 1, (1  f + fk2)R02 < 1, and (1  f + fk1k2)R012, where R0i = rbi/(l + ci) denotes the reproduction number of type i (12 denotes coinfection). The infection-free equilibrium is unstable if any of these conditions are violated. B.2. Existence of boundary equilibria, bivalent vaccine, full model Solving system (1) for equilibrium of type i exclusively gives qj ¼ 0; wj ¼ 0; y j ¼ 0; ui ¼ 0; uj ¼ 0; pi ¼ 0; pj ¼ 0; w12 ¼ 0; y 12 ¼ 0; 1  k i þ k i R0i  Ui 1  k i ð1 þ R0i Þ þ Ui c wi lzi vi ¼ ; xi ¼ ; qi ¼ i ; y i ¼ ; 2ð1  k i Þk i R0i 2ð1  k i Þk i R0i l ci l½ð1  k i Þð2fk i R0i  1Þ  k i R0i þ Ui

c fR0i  ð1  k i Þ½1  ð1  2f ÞR0i  Ui g wi ¼ ; zi ¼ i ; 2ðl þ ci Þð1  k i ÞR0i 2ðl þ ci Þð1  k i ÞR0i

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

115

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where Ui ¼ ½1 þ k i ðR0i  1Þ 2  4f ð1  k i Þk i R0i ; i; j ¼ 1; 2; i 6¼ j. Clearly, there exists a unique positive type i-only equilibrium if and only if R0i[1  f(1  ki)] > 1. B.3. Local stability of boundary equilibria, bivalent vaccine, full model From the determinant of the Jacobian matrix evaluated at the equilibrium of type i exclusively, the equilibrium is stable if R0j <

1  R012 ðxi þ k 1 k 2 vi Þ ; xi þ k j ðvi þ aj qi þ cj wi Þ þ cj y i þ aj zi  R012 ðxi þ k 1 k 2 vi ÞWi

where Wi ¼ aj ðk j qi þ zi Þ þ

ðl þ cj Þxi ðl þ cj Þk j vi þ ; i 6¼ j l þ cj þ ci R0i ðy i þ wi Þðl þ ci Þ l þ cj þ k i ci R0i ðy i þ wi Þðl þ ci Þ

and is unstable otherwise. The two types coexist if both boundary equilibria are unstable: R0i[1  f(1  ki)] > 1 and R0j >

1  R012 ðxi þ k 1 k 2 vi Þ ; xi þ k j ðvi þ aj qi þ cj wi Þ þ cj y i þ aj zi  R012 ðxi þ k 1 k 2 vi ÞWi

i; j ¼ 1; 2; i 6¼ j.

Note that when concurrent acquisition of infection is not allowed (i.e., R012 = 0), the above conditions become R0i[1  f(1  ki)] > 1 and 1 R01 < R02 ½x2 þ k 1 ðv2 þ a1 q2 þ c1 w2 Þ þ c1 y 2 þ a1 z2 R02 < R01 ½x1 þ k 2 ðv1 þ a2 q1 þ c2 w1 Þ þ c2 y 1 þ a2 z1 .

References [1] F.X. Bosch, A. Lorincz, N. Munoz, C.J. Meijer, K.V. Shah, The causal relation between human papillomavirus and cervical cancer, J. Clin. Pathol. 55 (2002) 244. [2] D. Parkin, F. Bray, J. Ferlay, P. Pisani, Estimating the world cancer burden: Globocan 2000, Int. J. Cancer. 94 (2001) 153. [3] R. Insinga, A. Glass, B. Rush, The healthcare costs of cervical human papillomavirus-related disease, Am. J. Obstet. Gynecol. 191 (2004) 114. [4] S. Chan, H. Delius, A. Halpern, H. Bernard, Analysis of genomic sequences of 95 papillomavirus types: uniting typing, phylogeny, and taxonomy, J. Virol. 69 (1995) 3074. [5] G. Ho, R. Bierman, L. Beardsley, C. Chang, R. Burk, Natural history of cervicovaginal papillomavirus infection in young women, N. Engl. J. Med. 338 (1998) 423. [6] R. Herrero, A. Hildesheim, C. Bratti, M. Sherman, M. Hutchinson, et al., A population-based study of human papillomavirus infection and all grades of cervical neoplasia in rural Costa Rica, J. Natl. Cancer. Inst. 92 (2000) 464. [7] F. Bosch, M. Manos, N. Munoz, M. Sherman, A. Jansen, et al., Prevalence of human papillomavirus in cervical cancer: a worldwide perspective, J. Natl. Cancer. Inst. 87 (1995) 796. [8] E.-M. DeVilliers, Heterogeneity of the human papillomavirus group, J. Virol. 63 (1989) 4898.

116

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

[9] L.A. Koutsky, K.A. Ault, C.M. Wheeler, D.R. Brown, E. Barr, et al., A controlled trial of a human papillomavirus type 16 vaccine, N. Engl. J. Med. 347 (2002) 1645. [10] D. Brown, K. Fife, C. Wheeler, L. Koutsky, L. Lupinacci, et al., Early assessment of the efficacy of a human papillomavirus type 16 L1 virus-like particle vaccine, Vaccine 22 (2004) 2936. [11] K. Fife, C. Wheeler, L. Koutsky, E. Barr, D. Brown, et al., Dose-ranging studies of the safety and immunogenicity of human papillomavirus type 11 and type 16 virus-like particle candidate vaccine in young healthy women, Vaccine 22 (2004) 2943. [12] K.U. Jansen, A.R. Shaw, Human papillomavirus vaccines and prevention of cervical cancer, Annu. Rev. Med. 55 (2004) 319. [13] K. Thomas, J. Hughes, J. Kuypers, N. Kiviat, S.-K. Lee, et al., Concurrent and sequential acquisition of different genital human papillomavirus types, J. Infect. Dis. 182 (2000) 1097. [14] K.-L. Liaw, A. Hildesheim, R. Burk, P. Gravitt, S. Wacholder, et al., A prospective study of human papillomavirus (HPV) type 16 DNA detection by polymerase chain reaction and its association with acquisition and persistence of other HPV types, J. Infect. Dis. 183 (2001) 8. [15] M.-C. Rousseau, J. Pereira, J.C. Prado, L. Villa, T. Rohan, E. Franco, Cervical coinfection with human papillomavirus (HPV) types as a predictor of acquisition and persistence of HPV infection, J. Infect. Dis. 184 (2001) 1508. [16] F. Unckell, R. Streeck, M. Sapp, Generation and neutralization of pseudovirions of human papillomavirus type 33, J. Virol. 71 (1997) 2934. [17] L.J. White, M.J. Cox, G.F. Medley, Cross immunity and vaccination against multiple microparasite strains, IMA J. Math. Appl. Med. Biol. 15 (1998) 211. [18] R. Roden, H. Greenstone, R. Kirnbauer, In vitro generation and type specific neutralization of a human papillomavirus type 16 virion pseudotype, J. Virol. 70 (1996) 5875. [19] R. Roden, N. Hubbert, R. Kirnbauer, N. Christensen, D. Lowy, J. Schiller, Assessment of the serological relatedness of genital human papillomaviruses by hemagglutination inhibition, J. Virol. 70 (1996) 3298. [20] J. Carter, L.A. Koutsky, G. Wipf, et al., The natural history of human papillomavirus type 16 capsid antibodies among a cohort of university women, J. Infect. Dis. 174 (1996) 927. [21] R. Roden, W. Yutzy, R. Fallon, S. Inglis, D. Lowy, J. Schiller, Minor capsid protein of human genital papillomaviruses contains subdominant, cross-neutralizing epitopes, Virology 270 (2000) 254. [22] G.Y.F. Ho, Y. Studentsov, C.B. Hall, R. Bierman, L. Beardsley, et al., Risk factors for subsequent cervicovaginal human papillomavirus (HPV) infection and the protective role of antibodies to HPV-16 virus-like particles, J. Infect. Dis. 186 (2002) 737. [23] R. Kirnbauer, N. Hubbert, C. Wheeler, T. Becker, D. Lowy, J. Schiller, A virus-like particle enzyme-linked immunosorbent assay detects serum antibodies in a majority of women infected with human papillomavirus type 16, J. Natl. Cancer. Inst. 86 (1994) 494. [24] T. Sasagawa, H. Yamazaki, Y. Dong, S. Satake, M. Tateno, M. Inoue, Immunoglobulin-A and -G responses against virus-like particles (VLP) of human papillomavirus type 16 in women with cervical cancer and cervical intra-epithelial lesions, Int. J. Cancer 75 (1998) 529. [25] A. McLean, Vaccination, evolution and changes in the efficacy of vaccines: a theoretical framework, Proc. R. Soc. Lond. B 261 (1995) 389. [26] M. Lipsitch, Vaccination against colonizing bacteria with multiple serotypes, Proc. Natl. Acad. Sci. 94 (1997) 6571. [27] T.C. Porco, S.M. Blower, Designing HIV vaccination policies: subtypes and cross-immunity, Interfaces 28 (1998) 167. [28] M. Lipsitch, Bacterial vaccines and serotype replacement: lessons from Haemophilus influenzae and prospects for Streptococcus pneumoniae, Emerg. Infect. Dis. 5 (1999) 336. [29] T.C. Porco, S.M. Blower, HIV vaccines: the effect of the mode of action on the coexistence of HIV subtypes, Math. Popul. Stud. 8 (2000) 205. [30] M. Nowak, R.M. May, Superinfection and the evolution of parasite virulence, Proc. R. Soc. Lond. B 255 (1994) 81. [31] R.M. May, M.A. Nowak, Coinfection and the evolution of parasite virulence, Proc. R. Soc. Lond. 261 (1995) 209. [32] S. Gupta, N.M. Ferguson, R.M. Anderson, Chaos, persistence, and evolution of strain structure in antigenically diverse infectious agents, Science 280 (1998) 912.

E.H. Elbasha, A.P. Galvani / Mathematical Biosciences 197 (2005) 88–117

117

[33] J. Gog, B. Grenfell, Dynamics and selection of many-strain pathogens, Proc. Natl. Acad. Sci. 94 (2002) 17209. [34] S. Gupta, A. Galvani, The effects of host heterogeneity on pathogen population structure, Philos. Trans. R. Soc. Lond. B 354 (1999) 711. [35] L. Esteva, C. Vargas, Analysis of a Dengue disease transmission model, Math. Biosci. 150 (1998) 131. [36] N. Ferguson, R. Anderson, S. Gupta, The effect of antibody-dependent enhancement on the transmission dynamics and persistence of multiple-strain pathogens, Proc. Natl. Acad. Sci. 96 (1999) 790. [37] L. Esteva, C. Vargas, Coexistence of different serotypes of dengue virus, J. Math. Biol. 46 (2003) 31. [38] G. Garnett, H. Waddell, Public health paradoxes and the epidemiological impact of an HPV vaccine, J. Clin. Virol. 19 (2000) 101. [39] J. Hughes, G. Garnett, L. Koutsky, The theoretical population level impact of a prophylactic human papillomavirus vaccine, Epidemiology 13 (2002) 631. [40] J. Koopman, Modeling infection transmission: the pursuit of complexities that matter, Epidemiology 13 (2002) 622. [41] R. Barnabas, G. Garnett, The Potential Impact of Vaccines against Human Papillomavirus, Imperial College of Science, Technology, and Medicine, London, UK, 2002. [42] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599. [43] C. Casillo-Chavez, W. Huang, J. Li, Competitive exclusion and coexistence of multiple strains in an SIS STD model, SIAM. J. Appl. Math. 59 (1999) 1790. [44] C. Casillo-Chavez, W. Huang, J. Li, The effects of females susceptibility on the coexistence of multiple pathogen strains of sexually transmitted diseases, J. Math. Biol. 35 (1997) 503. [45] R.M. May, S. Gupta, A.R. Mclean, Infectious disease dynamics: what characterizes a successful invader? Phil. Trans. R. Soc. Lond. B 356 (2001) 901. [46] R.M. Anderson, R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University, Oxford, 1991. [47] H. Richardson, G. Kelsall, P. Tellier, H. Voyer, M. Abrahamowicz, et al., The natural history of type-specific human papillomavirus infections in female university students, Cancer Epidemiol. Biomarkers. Prev. 12 (2003) 485. [48] A. Giuliano, R. Harris, R. Sedjo, S. Baldwin, D. Roe, et al., Incidence, prevalence, and clearance of type-specific human papillomavirus infections: the young womens health study, J. Infect. Dis. 186 (2002) 462. [49] E. Franco, L. Villa, J. Sobrinho, J. Prado, M.-C. Rousseau, et al., Epidemiology of acquisition and clearance of cervical Human Papillomavirus infection in women from a high-risk area for cervical cancer, J. Infect. Dis. 180 (1999) 1415. [50] M. Molano, A. Brule1, M. Plummer, E. Weiderpass, H. Posso, et al., Determinants of clearance of human papillomavirus infections in Colombian women with normal cytology: a population-based, 5-year follow-up study, Am. J. Epidemiol. 158 (2003) 486. [51] V. Dalstein, D. Riethmuller, J.-L. Pretet, K. Carval, J.-L. Sautiere, et al., Persistence and load of high-risk HPV are predictors for development of high-grade cervical lesions: A longitudinal French cohort study, Int. J. Cancer. 106 (2003) 396. [52] A. Frega, P. Stentella, C. Villani, D. Ruzza, G. Marcomin, et al., Correlation between cervical intraepithelial neoplasia and human papillomavirus male infections: a longitudinal study, Eur. J. Gynecol. Oncol. 20 (1999) 228. [53] C. Nebesio, G. Mirowski, T.-Y. Chuang, Human papillomavirus: clinical significance and malignant potential, Int. J. Dermatol. 40 (2001) 373. [54] J.M. Hyman, J. Li, E.A. Stanley, The initialization and sensitivity of multigroup models for the transmission of HIV, J. Theor. Biol. 208 (2001) 227. [55] W.A.A. Tjalma, M. Arbyn, J. Paavonen, T.R. Van Waes, J.J. Bogers, Prophylactic human papillomavirus vaccines: the beginning of the end of cervical cancer, Int. J. Gynecol. Cancer 14 (2004) 751. [56] D. Lowry, I. Frazer, Prophylactic human papillomavirus vaccines, J. Natl. Cancer. Inst. Mongr. 31 (2003) 111.