A geographical spread of vaccine-resistance in avian influenza epidemics

ARTICLE IN PRESS Journal of Theoretical Biology 259 (2009) 219–228

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

A geographical spread of vaccine-resistance in avian influenza epidemics Shingo Iwami a,, Yasuhiro Takeuchi a, Xianning Liu b, Shinji Nakaoka c a b c

Graduate School of Science and Technology, Shizuoka University, Japan School of Mathematics and Statistics, Southwest University, China Aihara Complexity Modelling Project, ERATO, JST, The University of Tokyo, Japan

a r t i c l e in fo

abstract

Article history: Received 13 August 2008 Received in revised form 20 February 2009 Accepted 31 March 2009 Available online 8 April 2009

Vaccination can be a useful tool for control of avian influenza outbreaks in poultry, but its use is reconsidered in most of the countries worldwide because of its negative effects on the disease control. One of the most important negative effects is the potential for emergence of vaccine-resistant viruses. Actually, in the vaccination program in China and Mexico, several vaccine-resistant strains were confirmed. Vaccine-resistant strains usually cause a loss of the protection effectiveness of vaccination. Therefore, a vaccination program that engenders the emergence of the resistant strain might promote the spread of the resistant strain and undermine the control of the infectious disease, even if the vaccination protects against the transmission of a vaccine-sensitive strain. We designed and analyzed a deterministic patch-structured model in heterogeneous areas (with or without vaccination) illustrating transmission of vaccine-sensitive and vaccine-resistant strains during a vaccination program. We found that the vaccination program can eradicate the vaccine-sensitive strain but lead to a prevalence of vaccine-resistant strain. Further, interestingly, the replacement of viral strain could occur in another area without vaccination through a migration of non-infectious individuals due to an illegal trade of poultry. It is also a novel result that only a complete eradication of both strains in vaccination area can achieve the complete eradication in another areas. Thus we can obtain deeper understanding of an effect of vaccination for better development of vaccination strategies to control avian influenza spread. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Epidemic model Patch-structured model Avian influenza Vaccination program Geographical spread Illegal poultry trade

1. Introduction It is likely that H5N1 virus infections among domestic poultry have become endemic in certain areas such as China, Thailand and Vietnam (Chen et al., 2004; Guan et al., 2004; Li et al., 2004; Webster et al., 2006). Vaccination of domestic poultry against the H5N1 subtype of avian influenza has been used in several countries such as Pakistan, Hong Kong, Indonesia, China, and Vietnam (Capua, 2007; Marangon et al., 2008; Peyre et al., 2008; Tiensin et al., 2005). Using vaccination to reduce the transmission rate might provide an alternative to mass culling, by reducing both the susceptibility of healthy birds and the infectiousness of infected birds (Capua and Marangon, 2004, 2006; Capua, 2007). However, it is known that several negative effects exist for using vaccination: a vaccination sometimes engenders vaccine-resistances, silent carriers, silent spreads, and antigen drifted or shifted variants (Gambotto et al., 2008; Lee et al., 2004; Pasquato and Seidah, 2008; Peyre et al., 2008; Savill et al., 2006; Smith et al., 2006). Therefore, if not used appropriately, vaccination might fail in avian influenza control among poultry.

 Corresponding author. Tel.: +81 53 478 1200.

E-mail address: [email protected] (S. Iwami). 0022-5193/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2009.03.040

The most important issue related to those negative effects is the potential for the emergence of vaccine-resistant influenza viruses. The vaccine-resistant strain, in general, causes a loss of the protection effectiveness of vaccination (Lee et al., 2004; Pasquato and Seidah, 2008; Peyre et al., 2008; Smith et al., 2006). Consequently, a vaccination program that engenders the emergence of the resistant strain might promote the spread of the resistant strain, even if the vaccination protects against the transmission of a vaccine-sensitive strain. Actually, it was reported that vaccination programs at Mexico and China led to vaccine-resistant strains (Lee et al., 2004; Pasquato and Seidah, 2008; Peyre et al., 2008; Smith et al., 2006). In Mexico, the H5N2 vaccines had been used since 1995, and phylogenetic analysis suggests the presence of (previously uncharacterized) multiple sublineages of Mexican lineage isolates. And also, in China, the H5N1 vaccines have been used since September 2005, and genetic analysis revealed that an H5N1 influenza variant (Fujian-like, FJlike), which is a previously uncharacterized H5N1 virus sublineage, had emerged (Pasquato and Seidah, 2008; Smith et al., 2006) (some researchers consider that the emergence and replacement of FJ-like virus are questionable, Guan et al., 2004; Leung, 2007). Further, interestingly, this virus had also transmitted to nonvaccinated area such as Hong Kong, Laos, Malaysia, and Thailand, resulting in a new transmission and outbreak wave in Southeast

ARTICLE IN PRESS 220

S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

Asia (Smith et al., 2006). The emergence of vaccine-resistant strains presents the risk of generating a new pandemic virus that is dangerous for humans through an avian-human link because of the spread of vaccine-resistant strain (Capua and Marangon, 2004; Lee et al., 2004; Peyre et al., 2008). We must understand, in order to develop effective vaccination strategies, why the resistant strain emerges and spreads over the large geographical region with or without vaccination program. In this study, we developed a deterministic patch-structured model in heterogeneous areas (with or without vaccination) illustrating transmission of vaccine-sensitive and vaccine-resistant strains during the vaccination program. In the model, we considered risk factors such as illegal domestic poultry and exotic bird trade to explain a large geographical spread of the resistant strain (Gilbert et al., 2006; Kilpatrick et al., 2006; Yee et al., 2009). Some researchers suspect illegal trade of poultry or poultry products as a source for H5N1 outbreaks (Kilpatrick et al., 2006). Actually, in some counties with H5N1 cases, where the demand for poultry is high, despite known risks of H5N1 transmission, poultry is transported illegally (for example, authorities in Vietnam estimated up to 70% of poultry that are illegally transported from China, go undetected, Yee et al., 2009). Our findings are that the vaccination program can eradicate the vaccine-sensitive strain but lead to a prevalence of vaccineresistant strain in the both heterogeneous areas (i.e., the large geographical spread of the resistance can occur). Further, from our exhaustive study, we proposed some possibility of a complete eradication of both strains from all areas.

2. Materials and methods Herein, we describe a deterministic patch-structured model of avian influenza in heterogeneous areas and its control using a vaccination program in the context of a pre-existence of a vaccineresistant strain (Fig. 1). Although, actually, an occurrence of the vaccine-resistant strain might be caused by immunological pressures of the vaccination, the vaccine-resistant strain is assumed to be present at low levels in both areas before the

X1

(1-p)c

c

ω1

eX1

b + my

b

V1

σφ1

pc b + mz Area 1

ω2

Y2 b + my

b

e(X1 + V1)

φ1

b

X2

Y1

φ2

V2

Z1 eV1

b

σφ2

Z2 b + mz

Area 2

Fig. 1. Model schematic showing a vaccination program and illegal trades. We consider, in the context of a pre-existence of the vaccine-resistant strain, whether the resistance is selected by the program in each area. Note that the vaccineresistant strain is assumed to be present at low levels in both areas before the program. The vaccination is executed in Area 1 and not in Area 2, but these areas are combined by illegal trades of poultry from Area 1 to Area 2. We mention that only susceptible and vaccinated birds can be exported because those strains can cause severe illness and high mortality for birds. Therefore the migration of susceptible and vaccinated birds affects a balance of prevalence between those strains in Area 2. In each area, susceptible birds ðX 1 ; X 2 Þ become infected with vaccine-sensitive ðY 1 ; Y 2 Þ and vaccine-resistant ðZ 1 ; Z 2 Þ strains at rates in direct relation to the number of respective infectious birds. We consider that vaccinated birds ðV 1 ; V 2 Þ can be protected completely from the vaccine-sensitive strain, but are partially protected from vaccine-resistant strains. See the Mathematical model section for corresponding equations.

program. We consider that the vaccination is executed in Area 1 and not in Area 2, but these areas are combined by illegal trades of poultry from Area 1 to Area 2. We regard that only susceptible and vaccinated birds can be exported at the rate e because those strains can cause severe illness and high mortality for birds (we can expect that the migration can affect a balance of prevalence between those strains in Area 2). We investigate whether the resistance is eventually selected by the program in each area. All birds in the effective population are divided into several compartments, respectively, including susceptible birds ðX i Þ, vaccinated birds ðV i Þ, birds infected with vaccine-sensitive strain ðY i Þ, and birds infected with vaccine-resistant strain ðZ i Þ in Area i ði ¼ 1; 2Þ. We assume that susceptible birds are born or restocked at a rate of c per day and that all birds are naturally dead or removed from the effective population at a rate of b per day in each area. In the absence of vaccination, transmission occurs at a rate that is directly related to the number of infectious birds, with respective transmission rate constants oi and fi from infected birds with the vaccine-sensitive strain and with the vaccineresistant strain in Area i. The infectious periods of vaccinesensitive and vaccine-resistant strain are assumed to be exponentially distributed, respectively, with mean durations of 1=ðb þ my Þ and 1=ðb þ mz Þ days. We implicitly assume that the infected bird with one strain cannot be infected with other strain. Actually, the infected birds rapidly die before the infection of other strains because the mean infectious period of infected birds is very short (Stegeman et al., 2004; Tiensin et al., 2005, 2007). At the beginning of the vaccination program, X 1 moves directly to V 1 by the vaccination. However, after some period after the initial vaccination, the direct movement might vanish because almost all birds are vaccinated. Therefore, we can assume that vaccination is only administered to the newly hatched birds. The newly hatched birds are vaccinated at the rate 0ppp1 (more appropriately, p is proportional). We remark that p represents a vaccination coverage among poultry in Area 1. To simplify the theoretical treatment, as described in Savill et al. (2006), we assume that the vaccinated birds can be protected completely from the vaccine-sensitive strain. Actually, in laboratory experience, many avian influenza vaccines confer a very high level of protection against clinical signs and mortality (90–100% protected birds) (Peyre et al., 2008). However, many factors determine whether a vaccinated bird becomes infected, including age, species, challenge dose, health, antibody titre, infections of immunosuppressive diseases, and crossreactivity of other avian influenza serotypes (Savill et al., 2006; Seo and Webster, 2001; Swayne et al., 1999; van den Berga et al., 2007). On the other hand, we assume that the vaccinated birds are partially protected from the vaccine-resistant strain at the rate (proportion) 0p1  sp1 because of cross-reactivity of immune systems (Hayden, 2001; Lee et al., 2004; Pasquato and Seidah, 2008; Smith et al., 2006; van den Berga et al., 2007) (e.g., s ¼ 0 represents complete cross-immunity against vaccine-resistant strains). Note that s represents a loss of protection effectiveness of the vaccination caused by a vaccine-resistant strain.

2.1. Mathematical model We extended the standard susceptible–infective model (Anderson and May, 1991) including the effect of the vaccination program and the illegal trade of poultry in heterogeneous areas. Our patch-structured mathematical model is given by the

ARTICLE IN PRESS S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

The epidemiological and biological feature of antiviral drugresistance is well reported in Hayden (2001). The transmissibility and virulence of drug-resistant strains are usually lower than those of the wild strain because of its mutation cost (Lipsitch et al., 2007; Regoes and Bonhoeffer, 2006; Hayden, 2001; Handel et al., 2006). Actually, antiviral drugs are also used for prophylaxis drug intervention as vaccination (Lipsitch et al., 2007; Regoes and Bonhoeffer, 2006; Stilianakis et al., 1998). Herein, we use some reduced values of transmissibility (f=o ¼ 0:58 where f ¼ f1 ¼ f2 ) and the infectious period of the vaccine-sensitive strain ððb þ my Þ=ðb þ mz Þ ¼ 1:32Þ for parameters of vaccine-resistant strain in both areas. Further, we, respectively, sample the vaccination coverage p from the range of [0,1] (p is proportion), the loss of protection effectiveness of the vaccination s from the range of [0,1] (s is proportion), and the export rate e from the range of [0,0.03]. We assume that minimum mean duration that birds in Area 1 are not exported to Area 2 is about a month (i.e., the rate of the birds in 1 Area 1 exported to Area 2 is e ¼ 0:03 day ). Actually we do not have any justification for e, because there are not any data or paper about it, but the range is not biologically unreasonable and can be chosen arbitrarily (if it is relatively large) with little effect on the meaning of the results.

following equations: X 01 V 01 Y 01 Z 01 X 02 V 02 Y 02 Z 02

221

¼ ð1  pÞc  ðb þ eÞX 1  ðo1 Y 1 þ f1 Z 1 ÞX 1 , ¼ pc  ðb þ eÞV 1  sf1 Z 1 V 1 , ¼ o1 Y 1 X 1  ðb þ my ÞY 1 , ¼ f1 Z 1 ðX 1 þ sV 1 Þ  ðb þ mz ÞZ 1 , ¼ c þ eX 1  bX 2  ðo2 Y 2 þ f2 Z 2 ÞX 2 , ¼ eV 1  bV 2  sf2 Z 2 V 2 , ¼ o2 Y 2 X 2  ðb þ my ÞY 2 , ¼ f2 Z 2 ðX 2 þ sV 2 Þ  ðb þ mz ÞZ 2 .

We do focus on the illegal trade of poultry but do not focus on a migration of wild birds in the model. Note that both infected birds with vaccine-sensitives strains ðY 1 Þ and those with vaccineresistant strains ðZ 1 Þ must directly move from Area 1 to Area 2, if we consider the migratory birds movement (only susceptible and vaccinated birds can move here). Certainly an effect of migratory birds movement is also considered as one of important risk factors for the spread of avian influenza strains (Kilpatrick et al., 2006; Normile, 2005; Russell et al., 2008; Yee et al., 2009). However, there are countries that have reported H5N1 infection in poultry in which infections are not associated with migratory bird movements and that did not report poultry trade with other reported infected countries (Kilpatrick et al., 2006; Yee et al., 2009). Actually, in South Asia such as Vietnam, Thailand and Malaysia, these illegal or improper trades are common and persistent and these birds were not vaccinated legally against H5N1 because of their illegal status (Gilbert et al., 2006; Yee et al., 2009). This implies that these trades might be a source for H5N1 outbreaks and make an avian influenza control by a vaccination difficult. Here we are interested in the effect of illegal poultry trade on the avian influenza control.

2.3. Reproduction numbers A measure of transmissibility and of the stringency of control policies required to stop an epidemic is the basic reproduction number, which is the number of secondary cases produced by each primary case during his/her period of infectiousness (Anderson and May, 1991). We obtain basic reproduction quantities of vaccine-sensitive strain Rsi ðpÞ and vaccine-resistant strain Rri ðpÞ in Area i ði ¼ 1; 2Þ. During the vaccination program, the basic reproduction numbers depend on the vaccination coverage ð0ppp1Þ. Here we remark that Rs1 ð1Þ is always 0. We derived these basic reproduction numbers in Appendix. With the estimated parameters in Table 1 the basic reproduction number of vaccine-sensitive and vaccine-resistant strain before the program is Rs1 ð0Þje¼0 ¼ Rs2 ð0Þje¼0 ¼ 6:53 and Rr1 ð0Þje¼0 ¼ Rr2 ð0Þje¼0 ¼ 4:96, respectively, if we do not consider any illegal trades (i.e., e ¼ 0). Note that Rs1 ð0Þje¼0 and Rs2 ð0Þje¼0 , respectively, correspond to estimated values in Stegeman et al. (2004). Furthermore, to clarify the concept of competition among strains simply, we introduce an invasive reproduction number for r the vaccine-resistant strain R¯ i ðpÞ in each area, which signifies an expected number of new infectious cases with the vaccineresistant strain after a spread of the vaccine-sensitive strain among birds (the detailed derivations are referred to Appendix). Similarly, we define the invasive reproduction number for the s vaccine-sensitive strain R¯ i ðpÞ. During the vaccination program, the invasive numbers also depend on the vaccination coverage. Here s r we remark that R¯ 1 ð1Þ and R¯ 1 ð1Þ cannot be defined (see Appendix). The invasive reproduction number is considered as a competitive

2.2. Estimation of epidemiological parameters Baseline values of model parameters used for simulations are presented in Table 1. These parameters are based on avian influenza epidemics among poultry in The Netherlands in 2003 (Elbers et al., 2004, 2005; Stegeman et al., 2004). The initial population size was c=b ¼ 984 birds at the 2003 epidemic (Stegeman et al., 2004). Usually, the mean lifespan of poultry is about 2 years. However, we assume that the mean duration of a bird being in effective population is about 1=b ¼ 100 days in each area because of migration and marketing. Therefore, the birth or restocking rate of birds is c ¼ 9:84 birds per day in each area. Estimated infectious period and transmission parameters are 1 1 13.8 days and 4:78  104 day individual , respectively (Stegeman et al., 2004). These pathogenic characteristics such as infectious and transmission parameters are used in our model as parameters of the vaccine-sensitive strain in both areas. That is, we assume that 1=ðb þ my Þ ¼ 13:8 days and o ¼ o1 ¼ o2 ¼ 1 1 4:78 104 day individual .

Table 1 Description of physical characteristics, transmission, and infectious parameters of the model with their baseline values used for simulations. Symbol

Description

Value (range)

Reference

c=b 1=ðb þ my Þ

Initial bird population size Mean infectious period of V-S strain Transmissibility of V-S strain

984 birds 13.8 days

Stegeman et al. (2004) Capua and Alexander (2004), Stegeman et al. (2004) Stegeman et al. (2004)

o

4:78  104 day ðb þ my Þ=ðb þ mz Þ Relative mean infectious period of V-R strain 1.32 Relative transmissibility of V-R strain 0.58 f=o

1

individual

1

Hayden (2001), Regoes and Bonhoeffer (2006), Stilianakis et al. (1998) Hayden (2001), Regoes and Bonhoeffer (2006), Stilianakis et al. (1998)

These parameters are based on avian influenza epidemics in The Netherlands in 2003 (Elbers et al., 2004, 2005; Stegeman et al., 2004). Actually, V-S and V-R represent ‘‘vaccine-sensitive’’ and ‘‘vaccine-resistant’’, respectively.

ARTICLE IN PRESS 222

S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

Table 2 Exhaustive study of the selection and eradication of prevalent strains during the vaccination program: we investigated what strains are eventually selected by the program in each area comparing the basic and invasion reproduction numbers. (I)

Area 1

14Rs1 ; 14Rr1     1

Area 2

(i) 14Rs2 ; 14Rr2

V-S and V-R are eradicated

r

(ii) f1oRs2 ; 14Rr2 g or f1oRs2 ; 1oRr2 ; R¯ 2 o1g s (iii) f14Rs2 ; 1oRr2 g or f1oRs2 ; 1oRr2 ; R¯ 2 o1g s

r

r

Area 1

f1oRs1 ; 14Rr1 g or f1oRs1 ; 1oRr1 ; R¯ 1 o1g

V-S is selected

Area 2

r (i) 1oRs2 , 1oRr2 , R¯ 2 o1 s (ii) 1oRs2 , 1oRr2 , R¯ 2 o1

V-S is selected

s

V-R is selected

r

V-S and V-R are selected

(iii) 1oR¯ 2 , 1oR¯ 2 (III)

(IV)

V-R is selected V-S and V-R are selected

(iv) 1oR¯ 2 , 1oR¯ 2 (II)

V-S and V-R are eradicated V-S is selected

s

Area 1

f14Rs1 ; 1oRr1     2g or f1oRs1 ; 1oRr1 ; R¯ 1 o1g

V-R is selected

Area 2

r (i) 1oRs2 , 1oRr2 , R¯ 2 o1

V-S is selected

s (ii) f14Rs2 ; 1oRr2 g or f1oRs2 ; 1oRr2 ; R¯ 2 o1g s r (iii) 1oR¯ 2 , 1oR¯ 2

V-R is selected

s

r

V-S and V-R are selected

Area 1

1oR¯ 1 , 1oR¯ 1

Area 2

(i) 1oRs2 ; 1oRr2 ; R¯ 2 o1 s (ii) 1oRs2 , 1oRr2 , R¯ 2 o1

r

s

V-S and V-R are selected

r

V-S is selected V-R is selected V-S and V-R are selected

(iii) 1oR¯ 2 , 1oR¯ 2

Our results imply that the emergence and spread of the resistance over the large geographical region is a possible phenomenon. Further, interestingly, we found a possibility that the program can eradicate both strains in both areas. Note that the cases Rr1 ð1Þo1 and Rr1 ð1Þ41 are, respectively, included in 1 and 2.

condition (relative fitness), which represents some advantage measure of the vaccine-resistant (vaccine-sensitive) strain against the vaccine-sensitive (vaccine-resistant) strain. If one of the invasive numbers is less than one, then the corresponding strain dies out and only the other strain survives because the latter strain overcomes the former strain. On the other hand, if both invasive numbers are greater than one, then both strains can coexist. Thus these concepts make it easy for us to understand the competitive exclusion or coexistence of these strains. The estimated invasive reproduction number of the vacciner r resistant strain (vaccine-sensitive) is R¯ 1 ð0Þje¼0 ¼ R¯ 2 ð0Þje¼0 ¼ 0:76 s s ðR¯ 1 ð0Þje¼0 ¼ R¯ 2 ð0Þje¼0 ¼ 1:32Þ if we do not consider any illegal trades.

2.4. Epidemiological scenarios We consider a scenario in which the vaccine-resistant strain can prevail during the vaccination program designed to be effective against the spread of a vaccine-sensitive strain. In order to investigate whether the resistance can be eventually (i.e., at a final phase of the epidemics) selected over a large geographical region, we assume that Rri ð0Þ41: otherwise the vaccine-resistant strain cannot emerge at all because Rri ðpÞ is a monotonically decreasing function of the vaccination coverage p (see Appendix). Acquisition of resistance ability usually engenders a strain which, in the absence of a pharmaceutical intervention, is less fit than the sensitive strain (Handel et al., 2006; Lipsitch et al., 2007; Moghadas et al., 2008; Stilianakis et al., 1998). Therefore, Rri ð0ÞoRsi ð0Þ in each area. We generally assume the following conditions for reproduction numbers before the vaccination program (our baseline parameter values are satisfied with these assumptions for e 2 ½0; 0:03): Rsi ð0Þ41;

Rri ð0Þ41;

r R¯ i ð0Þo1.

The assumption precludes the possibility that the pre-existing vaccine-resistant strain beats the vaccine-sensitive strain before r the vaccination program because R¯ i ð0Þo1. Further, since the resistance presents at low levels in both areas and the sensitive strain has already spread widely through the populations before the program, we assume that Z i ð0Þ40 and Y i ð0Þ is near some steady state in each area.

3. Results We explain how the vaccination program and the illegal trade affect a selection of prevalent strains over the large geographical region in the context of the pre-existence of the vaccine-resistant strain. Because we revealed how the final phase of the epidemic is composed (i.e., we proved existence and stability conditions of equilibria of the model in Appendix) by strictly mathematical analysis, we could consider the selection in each area by comparing the basic and invasive reproduction numbers.

3.1. Selection and eradication of prevalent strains We exhaustively investigated what strains are eventually selected by the program in each area as follows (see Table 2). When the program is executed, the patterns of the selection and eradication of prevalent strains in Area 1 are divisible into four cases (see Appendix). (I) both the vaccine-sensitive and vaccineresistant strains are eradicated, (II) the vaccine-sensitive strain is selected, (III) the vaccine-resistant strain is selected, (IV) both the vaccine-sensitive and vaccine-resistant strains are selected ((III) and (IV), respectively, represent a complete and partial selection of the resistance in Area 1). For each case, we evaluated the selection and eradication in Area 2 where is not vaccinated but

ARTICLE IN PRESS S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

affected on the vaccination program through the illegal trade (see Appendix). Results in Table 2 show that the replacement and spread of the resistance over the large geographical region is a possible phenomenon. In the case of (III-ii), (III-iii), (IV-ii), and (IV-iii), the resistance eventually spreads in both areas. We assume, for example, that the vaccination coverage is p ¼ 80%, the export rate is e ¼ 1%, and the loss of the protection effectiveness is s ¼ 80%, which correspond to the case (III-iii). We calculate epidemic curves with the vaccination program for 365 days in Fig. 2. The blue and red curves, respectively, represent the number of infected individuals with vaccine-sensitive and vaccine-resistant strains. The top and bottom figures, respectively, depict time courses of infection in Area 1 and in Area 2. The resistant strain excludes the sensitive strain in Area 1 and invades into Area 2. The program changes the prevalent strain over the large geographical region. And also the program seems to promote a coexistence of multiple strains. Further, we can observe that it takes about several months for the resistant strain to spread in each area. Actually, the replacement time of the resistant strain was reported as several months in the China and Mexico epidemics (Lee et al., 2004; Peyre et al., 2008; Smith et al., 2006). Thus, these results

Vaccinated Area 1 300

Resistantstrain Red

Sensitivestrain Blue

250 200 150 100 50 50

100

150 200 time days Non

250

300

350

300

350

vaccinated Area 2

300

Resistantstrain Red

Sensitivestrain Blue

250 200 150 100 50 50

100

150

200

250

time days Fig. 2. Time-course of the spread of the disease with the vaccination program. We assume that the vaccination coverage is p ¼ 0:8, the export rate is e ¼ 0:01, and the loss of the protection effectiveness is s ¼ 0:8. We calculate epidemic curves with the vaccination program for 365 days. The blue and red curves, respectively, represent the number of infected individuals with vaccine-sensitive and vaccineresistant strains. The top and bottom figures, respectively, depict time courses of infection in Area 1 and in Area 2. The program completely changes the prevalent strain in Area 1 (the resistant strain excludes the sensitive strain) and partially changes one in Area 2 (the both strains coexist). And also we can observe that it takes about several months for the resistant strain to spread in each area. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

223

show the possibility that the emergence and replacement of the resistant strain can be facilitated by the vaccination program, as in some vaccination programs (Lee et al., 2004; Peyre et al., 2008; Smith et al., 2006). Note that some sensitivity analyses concerned about the change of prevalent strains for p, e, and s are referred to Fig. 3. Furthermore, we can find a possibility that the program can eradicate both strains in both areas. That is, only a complete eradication of both strains in vaccination area can achieve the complete eradication in another area, which correspond to the case (I-i). The eradication cannot occur in the other situations. This is a very important information for disease control to prevent and eradicate some disease spread. 3.2. Impact of vaccination program and illegal trade We conducted simulations using our baseline parameters to elucidate how the vaccination program and the illegal trade affect the selection of the resistant strain at the final phase of the epidemic with the loss of the protection effectiveness s ¼ 40%, 60%, and 80% in the left, middle, and right figures in Fig. 3, respectively. The top and bottom figures, respectively, represent the outcomes of the vaccination program in the area with the vaccination program (i.e., Area 1) and in the area without the program (i.e., Area 2). The blue, green, red, and pink regions, respectively, correspond to the situation in which only the vaccine-sensitive strain is selected, both the vaccine-sensitive and the vaccine-resistant strains are selected, only the vaccineresistant strain is selected, and both strains are eradicated. Note that the resistance is partially and completely selected in the green and red regions, respectively. Results showed that the final phase is significantly affected on by the vaccination program and the illegal trade. In general, the high vaccination coverage leads to a spread of the vaccineresistant strain at the final phase. However, as the export rate e increases, the resistance becomes difficult to be selected in Area 1. Further, when the loss of the protection rate is relatively small, the high coverage can eradicate both strains in Area 1. On the other hand, in Area 2, as the export rate increases, the resistance tends to be easily selected. From these asymmetrical effects of the program and the trade, we could observe non-synchronized changes of the prevalent strain over the large geographical region. For example, if the export rate is relatively high (e.g., e ¼ 0:02), the resistant strain is partially selected in Area 2 before the selection in Area 1, but the sensitive strain is eradicated in Area 1 before the eradication in Area 2, as the coverage increases. Thus the illegal trade can affect a balance law of the prevalence strain in nonvaccinated area and make the avian influenza control difficult and complex over the large geographical region (Capua and Marangon, 2004; Lee et al., 2004; Pasquato and Seidah, 2008; Peyre et al., 2008; Smith et al., 2006). Illegal trades in poultry are a serious social behavior in order to evaluate the effect of vaccination programs more precisely. Note that the above results are significantly dependent on the loss of the protection effectiveness s. It is expected that, as the loss decreases, the possibility of the eradication of both strains increases (i.e., the pink regions expand). Although the exact value of the loss is not clarified, we investigated an impact of the loss on the outcome of the program, the final size of the epidemic, and the optimal vaccination coverage during the program in Iwami et al. (2009).

4. Discussion A reason of the emergence of new H5N1 variants might be a natural evolution of avian viruses (for example, adaptation

ARTICLE IN PRESS 224

S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

e

e

Area 1:  =0.4

e

Area 1:  =0.6

0.03

0.03

0.03

0.02

0.02

0.02

0.01

0.01

0.01

0 0

0.2

0.4

e

0.6

0.8

1

p

0 0

0.2

0.4

e

Area 2:  =0.4

0.6

0.8

1

p

0 0

0.03

0.03

0.02

0.02

0.02

0.01

0.01

0.01

0

0.2

0.4

0.6

0.8

1

p

0 0

0.2

0.4

0.6

0.8

0.2

0.4

e

Area 2:  =0.6

0.03

0

Area 1:  =0.8

1

p

0.6

0.8

1

Area 2:  =0.8

0 0

p

0.2

0.4

0.6

0.8

1

p

Fig. 3. The outcomes of the vaccination program over the large geographical region: we assumed that the loss of the protection effectiveness in the left, middle, right figures is s ¼ 40%, 60%, and 80%, respectively. The blue, green, red, and pink region, respectively, corresponds to the situation in which only the vaccine-sensitive strain is selected, both the vaccine-sensitive and the vaccine-resistant strains are selected, only the vaccine-resistant strain is selected, and both strains are eradicated. The top and bottom figures, respectively, represent which strain is selected in Area 1 and Area 2. Although the selection significantly depends on the vaccination coverage and the export rate, the high vaccination coverage generally leads to a spread of the vaccine-resistant strain at the final phase. Further, we could observe non-synchronized changes of the prevalent strain in both areas. Thus the illegal trade can affect a balance law of the prevalence strain in non-vaccinated area and make the avian influenza control difficult and complex over the large geographical region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

following a cross-species jump from wild birds to domesticated birds) (Skeik and Jabr, 2007). Actually, we now have a diversity of avian H5N1 lineages co-circulating, even in countries which did not use vaccines (Lam et al., 2008; Vijaykrishna et al., 2008; Yee et al., 2009). The natural evolution of the virus might simply lead to the change of prevalent strains over a large geographical region. However, our theoretical study shows the possibility that the emergence might be caused by a vaccination program. In the context of a pre-existence of the vaccine-resistant strain, the program can change a balance of prevalence between vaccinesensitive and vaccine-resistant strains in both vaccinated and non-vaccinated areas through an illegal trade of poultry. The case (III-ii) in Table 2 represents that the resistant strain excludes the sensitive strain and spreads in both areas. The cases (III-iii), (IV-ii), and (IV-iii) indicate that the resistance invades into both areas. And also Fig. 3 shows which strain is selected by the program in each area using our baseline parameters in Table 1. The program seems to be able to promote the spread of the resistant strain. Further, interestingly, if the illegal export rate is low, both strains can be selected in both areas (green region), but if the export rate is high, the resistance cannot be selected in the vaccinated area (blue region) but be partially selected in the non-vaccinated area (green region) in Fig. 3. This non-synchronized changes of the prevalent strain can be explained as follows. As the export rate increases, the susceptible and vaccinated birds move from Area 1 to Area 2, which can increase some herd immunity of bird population against the resistant strain in Area 1 but decrease that in Area 2. This is because, although the sensitive strain can be maintained by infections of only the susceptible birds, the resistance must be maintained by infections of both the susceptible and vaccinated birds. Actually, we assumed that a fitness of the resistance strain is less than one of the sensitive r strain at the beginning of the program ðR¯ i ð0Þo1Þ. The poor resource, because of high export, leads to some advantage of the sensitive strain in vaccinated area and the rich resource leads to

some advantage of the resistant strain in non-vaccinated area in the context of the lower fitness of the resistance. Thus the program can affect the balance of prevalent strains in both vaccinated and non-vaccinated areas asymmetrically. These results imply that the emergence of the vaccine resistant strain such as Fujian-like H5N1 and Mexican lineage H5N2 can be caused by the vaccination programs(Lee et al., 2004; Pasquato and Seidah, 2008; Peyre et al., 2008; Smith et al., 2006). Many experts for the avian influenza worry about the change of prevalent avian influenza strain in Asia (Chen et al., 2004; Lam et al., 2008; Smith et al., 2006; Vijaykrishna et al., 2008) because the new strain may easily mutate and obtain a sustained humanto-human transmission ability than current strains. For example, since November 2005, 22 H5N1 human infection cases from 14 provinces of China have been reported. Genetic findings also revealed that FJ-like viruses also were responsible for all recently reported human infection cases in China (Smith et al., 2006). This implies that the prevalent strain in poultry directly affects the human infection because almost human infections result from contacts with the infected poultry, or with the surfaces contaminated with secretion/excretions of the infected birds (Chen et al., 2004; Gambotto et al., 2008). That is, if the mutation acquiring the transmission ability among humans occurs, global pandemic such as the 1918 ‘‘Spanish influenza’’ might occur and have devastating consequences in human. Therefore we must manage carefully the vaccination program which might select a vaccine-resistant strain to control avian influenza in poultry. Fortunately, we could find the possibility that the program can eradicate both strains in both areas. Our deeper understanding of an effect of vaccination is one of the important priority in development of vaccination strategies to control avian influenza spread. Although we can obtain the novel conclusions noted above, our model may have some limitations that it includes no capacity for the virus to evolve more antigenic novelty, either through mutation accumulation or, perhaps more likely, through

ARTICLE IN PRESS S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

reassortment. To ignore the evolutionary dynamics of the virus means that our model may only ever have limited applicability. Further, we only consider a uni-directional migration of birds from the vaccinated area to the unvaccinated area. A bi-directional migration between the areas is probably more realistic for birds. These additional effects such as viral evolution and bi-directional migration may impact on our results but we leave the inclusion of the effects as future works.

Acknowledgments The authors would like to thank anonymous referees for very helpful suggestions and comments which improved the quality of this paper and study. Shingo is supported by Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists, Xianning was supported by the National Natural Science Foundation of China (10571143), the Science Foundation of Southwest University (SWNUB2004001) and the Japanese Government (Monbukagakusho:MEXT) Scholarship, and Shinji was partly supported by the Sasakawa Scientific Research Grant from The Japan Science Society.

Appendix A In order to understand the dynamical behavior of our mathematical model, we perform mathematical analysis. Note that the dynamics of X 1 , V 1 , Y 1 , and Z 1 are independent on that of X 2 , V 2 , Y 2 , and Z 2 . Therefore, we can consider Area 1 with separating Area 2 (but Area 2 is always dependent on Area 1).

The following are basic and invasive reproduction numbers which are functions of p: Rs1 ðpÞ ¼ Rr1 ðpÞ ¼

o1 b þ my

f1 b þ mz

V 01 ¼ pc  ðb þ eÞV 1  sf1 Z 1 V 1 , Y 01 ¼ o1 Y 1 X 1  ðb þ my ÞY 1 ,

X 01 þ

s R¯ 1 ðpÞ ¼

sf1 b þ mz

o1 b þ my

V 01 ;

ð1  pÞc ¼ where ¼ ; bþe b þ my ; Es1 ¼ ðX s1 ; V s1 ; Y s1 ; 0Þ where X s1 ¼

V 01

(i) (ii) (iii) (iv)

Proof. (i) It is clear that E01 always exists. (ii) We can easily show that Es1 exists iff Rs1 ðpÞ41. (iii) From Eq. (1), the existence condition of Er1 is given by

f1 ð1  pÞc þ sf1 pc bþe

ð1  pÞc ; b þ e þ f1 Z r1

f1 ð1  pÞc sf1 pc þ ¼ b þ mz , b þ e þ f1 Z 1 b þ e þ sf1 Z 1

Yþ 1

¼

1

o1



where X þ 1 ¼

b þ my

o1

;

() 1oRr1 ðpÞ.

(iv) Let F be the following function of X 1 :   1 2 X FðX 1 Þ ¼ ðb þ eÞf1 1  s 1     1 þ f1 c X 1 þ ð1  pÞcðb þ mz Þ.  ðb þ eÞðb þ mz Þ 1  Then we obtain the following existence condition of Eþ 1: b þ mz

Vþ 1 ¼

 ð1  pÞc  ðb þ eÞX þ 1  f1 Z þ 1 ; þ X1



V r1 ¼

 1 b þ mz

s

f1

Zþ 1



pc . b þ e þ sf1 Z r1

b þ my

o1

pcs b þ mz oX þ ; 1o bþe f1

0oFðX þ 1 Þ.

pcs b þ mz oX þ ; 0oFðX þ 1o 1Þ bþe f1   b þ mz pcs  oX þ  () max 0; 1 oX 1 , f1 bþe

b þ mz

 ,

pc  ðb þ eÞV þ 1 ¼ . sf1 V þ1



where X 1 is the larger root of FðX 1 Þ ¼ 0. From straightforward but tedious calculations, we can evaluate X 1 ¼ X r1 . This implies that   b þ mz pcs r  ¯s max 0; oX þ  1oR¯ 1 ðpÞ: & 1 oX 1 () 1oR1 ðpÞ; f1 bþe

In addition, Z r1 is the unique root of the following equation:

þ þ þ þ Eþ 1 ¼ ðX 1 ; V 1 ; Y 1 ; Z 1 Þ

V s1 .

E01 always exists in R4þ . Es1 exists in R4þ iff 1oRs1 ðpÞ. Er1 exists in R4þ iff 1oRr1 ðpÞ. 4 ¯s ¯r Eþ 1 exists in Rþ iff 1oR1 ðpÞ and 1oR1 ðpÞ.

f1

where X r1 ¼

sf1 b þ mz

Since 0oFð0Þ, Fðb þ mz =f1 Þo0 and F 00 ðX 1 Þo0, we can obtain the following relation:

ð1  pÞc  ðb þ eÞX s1 ¼ , o1 X s1

Er1 ¼ ðX r1 ; V r1 ; 0; Z r1 Þ

X s1 þ

Theorem A.1.

f1

pc ¼ , bþe pc V s1 ¼ , bþe

X 01

o1

Y s1

f1 b þ mz

s

This model has the following four possible equilibria which are functions of p: ðX 01 ; V 01 ; 0; 0Þ

r R¯ 1 ðpÞ ¼

The existence and stability conditions of E01 , Es1 , Er1 , and Eþ 1 are given in the following theorems:

Z 01 ¼ f1 Z 1 ðX 1 þ sV 1 Þ  ðb þ mz ÞZ 1 .

E01

X r1 ,

s r Lemma A.1. R¯ 1 ðpÞo1oRs1 ðpÞ and R¯ 1 ðpÞo1oRr1 ðpÞ cannot hold simultaneously.

b þ mz o

X 01 ¼ ð1  pÞc  ðb þ eÞX 1  ðo1 Y 1 þ f1 Z 1 ÞX 1 ,

X 01 ;

s r Note that Rs1 ð1Þ is always 0, R¯ 1 ð1Þ and R¯ 1 ð1Þ cannot be defined. Therein, the subscript ‘‘1’’ means Area 1 and the superscripts ‘‘0’’, ‘‘s’’, ‘‘r’’, and ‘‘þ’’, respectively, signify the disease-free equilibrium, vaccine-sensitive strain existing equilibrium, vaccine-resistant strain existing equilibrium, and the both-strains existing equilibrium. Here we have to note the relation between these reproduction numbers in the following lemma. The lemma can be proved directly by tedious and complex analysis but it will be clear in the later Theorem A.2.

A.1. Mathematical analysis of Area 1 The dynamics of X 1 , V 1 , Y 1 , and Z 1 are determined by the following equations:

225

Theorem A.2. (1)

Note that, if p ¼ 0, limt!1 V 1 ðtÞ ¼ 0 and then Eþ 1 disappears (except for a special case: Rs1 ð0Þ ¼ Rr1 ð0Þ) because of competitive exclusion principle (Bremermann and Thieme, 1989), and limt!1 X 1 ðtÞ ¼ 0, then Es1 , Eþ 1 clearly disappear if p ¼ 1.

(i) If Rs1 ðpÞp1 and Rr1 ðpÞp1, then E01 is globally asymptotically stable (GAS) which means that the orbit converges to the equilibrium as t ! 1 for arbitrary initial point. r (ii) If Rs1 ðpÞ41 and R¯ 1 ðpÞp1, then Es1 is GAS. s r (iii) If R1 ðpÞ41 and R¯ 1 ðpÞp1, then Er1 is GAS. s r ¯ ¯ (iv) If R1 ðpÞ41 and R1 ðpÞ41, then Eþ 1 is GAS.

ARTICLE IN PRESS 226

S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228 s We remark that o1 X r1  ðb þ my Þp0 iff R¯ 1 ðpÞp1. In the similar s _ manner, we can show that V r p0 because R¯ 1 ðpÞp1. This completes the proof. (iv) Let us consider the Lyapunov function

Proof. (i) Let us consider the Lyapunov function V 0 ¼ X 1  X 01 log X 1 þ V 1  V 01 log V 1 þ Y 1 þ Z 1 . We have   ð1  pÞc V_ 0 ¼ ðX 1  X 01 Þ  ðb þ eÞ  o1 Y 1  f1 Z 1 X1   pc  ðb þ eÞ  sf1 Z 1 þ ðV 1  V 01 Þ V1

þ V þ ¼ X1  Xþ 1 log X 1 þ V 1  V 1 log V 1 þ Y 1 þ  Yþ 1 log Y 1 þ Z 1  Z 1 log Z 1 .

Then   ð1  pÞc  ðb þ eÞ  o1 Y 1  f1 Z 1 V_ þ ¼ ðX 1  X þ 1Þ X1   pc þ þ ðV 1  V 1 Þ  ðb þ eÞ  sf1 Z 1 V1

þ Y 1 fo1 X 1  ðb þ my Þg þ Z 1 ff1 ðX 1 þ sV 1 Þ  ðb þ mz Þg ! ! X 01 X 1 V 01 V 1 þ pc 2    ¼ ð1  pÞc 2  X 1 X 01 V 1 V 01     b þ my b þ mz þ f1 Z 1 X 01 þ sV 01  . þ o1 Y 1 X 01 

o1

þ þ ðY 1  Y þ 1 Þfo1 X 1  ðb þ my Þg þ ðZ 1  Z 1 Þ

f1

We remark that X 01  ðb þ my Þ=o1 p0 iff Rs1 ðpÞp1 and X 01 þ sV 01  ðb þ mz Þ=f1 p0 iff Rr1 ðpÞp1. Further it is clear that 2

X 01 X 1  p0; X 1 X 01

2

V 01 V 1  p0, V 1 V 01

because the arithmetic mean is larger than, or equals to the geometric mean. Therefore V_ 0 p0 because Rs1 ðpÞp1 and Rr1 ðpÞp1, and we can conclude that by the Lyapunov–LaSalle’s invariance principle, all the trajectories of the model converges to E01 . (ii) Let us consider the Lyapunov function

ff1 ðX 1 þ sV 1 Þ  ðb þ mz Þg. þ þ þ þ Since b þ e ¼ ð1  pÞc=X þ 1  o1 Y 1  f1 Z 1 ¼ pc=V 1  sf1 Z 1 , b þ þ þ þ my ¼ o1 X 1 and b þ mz ¼ f1 ðX 1 þ sV 1 Þ, we can evaluate     Xþ X1 Vþ V1 V_ þ ¼ ð1  pÞc 2  1  þ þ pc 2  1  þ p0. X1 X1 V1 V1

This completes the proof.

&

Note that the cases Rr1 ð1Þo1 and Rr1 ð1Þ41 are, respectively, included in cases (i) and (iii).

V s ¼ X 1  X s1 log X 1 þ V 1  V s1 log V 1 þ Y 1  Y s1 log Y 1 þ Z 1 .

A.2. Mathematical analysis of Area 2

Then

Using a method of dynamical system theory in Castillo-Chavez and Thieme (1995), the dynamics of X 2 , V 2 , Y 2 , and Z 2 are determined by the following equations:



ð1  pÞc V_ s ¼ ðX 1   ðb þ eÞ  o1 Y 1  f1 Z 1 X1   pc  ðb þ eÞ  sf1 Z 1 þ ðV 1  V s1 Þ V1 X s1 Þ



X 02 ¼ c þ eX 1  bX 2  ðo2 Y 2 þ f2 Z 2 ÞX 2 ,

þ ðY 1  Y s1 Þfo1 X 1  ðb þ my Þg þ Z 1 ff1 ðX 1 þ sV 1 Þ  ðb þ mz Þg. Since b þ e ¼ ð1  pÞc=X s1  o1 Y s1 ¼ pc=V s1 and b þ my ¼ o1 X s1 , we can evaluate     Xs X1 Vs V1 V_ s ¼ ð1  pÞc 2  1  s þ pc 2  1  s X1 X1 V1 V1   b þ m z . þ f1 Z 1 X s1 þ sV s1 

f1

r We remark that þ sV s1 Þ  ðb þ mz Þp0 iff R¯ 1 ðpÞp1. In the r similar manner, we can show that V_ s p0 because R¯ 1 ðpÞp1. This completes the proof. (iii) Let us consider the Lyapunov function

f1 ðX s1

V 02 ¼ eV 1  bV 2  sf2 Z 2 V 2 , Y 02 ¼ o2 Y 2 X 2  ðb þ my ÞY 2 , Z 02 ¼ f2 Z 2 ðX 2 þ sV 2 Þ  ðb þ mz ÞZ 2 . Here X 1 and V 1 represent a corresponding pair of X 01 and V 01 (if r Rs1 ðpÞp1 and Rr1 ðpÞp1), X s1 and V s1 (if Rs1 ðpÞ41 and R¯ 1 ðpÞp1), X r1 s s r r þ þ and V 1 (if R1 ðpÞ41 and R¯ 1 ðpÞp1), or X 1 and V 1 (if R¯ 1 ðpÞ41 and r R¯ 1 ðpÞ41), respectively. This model has the following four possible equilibria which are functions of p: E02 ¼ ðX 02 ; V 02 ; 0; 0Þ Es2 ¼ ðX s2 ; V s2 ; Y s2 ; 0Þ s

c þ eX 1  bX 2 , o2 X s2

V r ¼ X 1  X r1 log X 1 þ V 1  V r1 log V 1 þ Y 1 þ Z 1  Z r1 log Z 1 .

Y s2 ¼

We have

Er2 ¼ ðX r2 ; V r2 ; 0; Z r2 Þ

  ð1  pÞc V_ r ¼ ðX 1  X r1 Þ  ðb þ eÞ  o1 Y 1  f1 Z 1 X1   pc r þ ðV 1  V 1 Þ  ðb þ eÞ  sf1 Z 1 V1

f1 Z r1

pc=V r1

f1 Z r1

Since b þ e ¼ ð1   ¼ s b þ mz ¼ f1 ðX r1 þ sV r1 Þ, we can evaluate     Xr X1 Vr V1 V_ r ¼ ð1  pÞc 2  1  r þ pc 2  1  r X1 X1 V1 V1   b þ my r . þ o1 Y 1 X 1 

o1

where X r2 ¼

c þ eX 1 ; b þ f2 Z r2

V r2 ¼

eV 1 . b þ sf2 Z r2

In addition, Z r2 is the unique root of the following equation:

þ Y 1 fo1 X 1  ðb þ my Þg þ ðZ 1  Z r1 Þff1 ðX 1 þ sV 1 Þ  ðb þ mz Þg. pÞc=X r1

c þ eX 1 eV 1 ; V 02 ¼ , b b  b þ my eV 1 , where X s2 ¼ ; V s2 ¼ o2 b

where X 02 ¼

and

f2 ðc þ eX 1 Þ sf2 eV 1 þ ¼ b þ mz , b þ sf2 Z 2 b þ f2 Z 2 þ þ þ þ Eþ 2 ¼ ðX 2 ; V 2 ; Y 2 ; Z 2 Þ

Yþ 2 ¼

where X þ 2 ¼

b þ my

o2

  þ c þ eX 1  bX 2  f2 Z þ 2 ; þ o2 X2 1

;

Vþ 2 ¼

 1 b þ mz

s

f2



b þ my

o2

 ,

þ

Zþ 2 ¼

eV 1  bV 2 . sf2 V þ2

Note that, if p ¼ 0, limt!1 V 2 ðtÞ ¼ 0 and then Eþ 2 disappears (except for a special case: Rs2 ð0Þ ¼ Rr2 ð0Þ).

ARTICLE IN PRESS S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

The following are basic and invasive reproduction numbers which are functions of p: Rs2 ðpÞ ¼ Rr2 ðpÞ ¼

o2 b þ my

f2 b þ mz

X 02 ; X 02 þ

s R¯ 2 ðpÞ ¼

sf2 b þ mz

o2 b þ my

V 02 ;

X r2 ,

r R¯ 2 ðpÞ ¼

f2 b þ mz

X s2 þ

sf2 b þ mz

V s2 ,

where the subscript ‘‘2’’ means Area 2. Here we also remark the relation between these reproduction numbers in the following lemma. s r Lemma A.2. R¯ 2 ðpÞo1oRs2 ðpÞ and R¯ 2 ðpÞo1oRr2 ðpÞ cannot hold simultaneously.

The existence and stability conditions of E02 , Es2 , Er2 , and Eþ 2 are given in the following theorems: Theorem A.3. (i) (ii) (iii) (iv)

E02 always exists in R4þ . Es2 exists in R4þ iff 1oRs2 ðpÞ. Er2 exists in R4þ iff 1oRr2 ðpÞ. 4 ¯s ¯r Eþ 2 exists in Rþ iff 1oR2 ðpÞ and 1oR2 ðpÞ.

(i) (ii) (iii) (iv)

If If If If

Rs2 ðpÞp1 Rs2 ðpÞ41 Rr2 ðpÞ41 s R¯ 2 ðpÞ41

(II) The vaccine-sensitive strain is selected in Area 1: When r Rs1 ðpÞ41 and Rr1 ðpÞp1, or Rs1 ðpÞ41, Rr1 ðpÞ41, R¯ 1 ðpÞo1, then the vaccine-sensitive strain is selected in Area 1 by Theorem A.2 (Es1 is GAS). Let 0opp1. It is clear that Rs2 ð0Þ ¼ Rs2 ðpÞ, Rr2 ð0ÞoRr2 ðpÞ and r r s s R¯ 2 ð0ÞoR¯ 2 ðpÞ. Further, since X r2 oX r2 jp¼0 , we can show R¯ 2 ð0Þ4R¯ 2 ðpÞ. We can conclude that the vaccination program leads to the following possibilities of the change in Area 2 by Theorem A.4: r (i) 1oRs2 ðpÞ, 1oRr2 ðpÞ, R¯ 2 ðpÞo1; Es2 is GAS. s s r (ii) 1oR2 ðpÞ, 1oR2 ðpÞ, R¯ 2 ðpÞo1; Er2 is GAS. s r (iii) 1oR¯ 2 ðpÞ, 1oR¯ 2 ðpÞ; Eþ 2 is GAS.

(III) The vaccine-resistant strain is selected in Area 1: When s Rr1 ðpÞ41 and Rs1 ðpÞp1, or Rs1 ðpÞ41, Rr1 ðpÞ41, R¯ 1 ðpÞo1, then the vaccine-resistant strain is selected in Area 1 by Theorem A.2 (Er1 is GAS). s s r r Let 0opp1. It is clear that R¯ 2 ð0Þ4R¯ 2 ðpÞ and R¯ 2 ð0ÞoR¯ 2 ðpÞ. We s r s ¯ remark that R1 ðpÞo1 is equivalent to X 1 oX 1 . This implies that r Rs2 ð0Þ4Rs2 ðpÞ. Remember that Rs1 ð0Þ41, Rr1 ð0Þ41, R¯ 1 ð0Þo1, s s r R1 ðpÞ41, R1 ðpÞ41 and R¯ 1 ðpÞo1. In this case, we can show c þ eX r1 jp¼0 f2 c þ eX r1 sf2 eV r1 o þ b þ mz b b þ mz b b þ mz b s c eX 1 jp¼0 c eX r1 seV r1 o þ þ () þ b b b b b s r r ()X 1 jp¼0 oX 1 þ sV 1

Rr2 ð0ÞoRr2 ðpÞ()

Theorem A.4. and and and and

Rr2 ðpÞp1, r R¯ 2 ðpÞp1, s R¯ 2 ðpÞp1, ¯Rr2 ðpÞ41,

then then then then

E02 is GAS. Es2 is GAS. Er2 is GAS. Eþ 2 is GAS.

The proofs of the above theorems are essentially similar to the case of Area 1.

A.3. Evaluation of vaccination program When the vaccination program is executed in Area 1, the outcome of the selection and eradication of prevalent strains during the vaccination program in Area 1 are divisible into four cases. For each case, we investigate possibilities of the outcomes in Area 2 during the program by comparing the basic and invasive reproduction numbers. (I) Both the vaccine-sensitive and vaccine-resistant strains are eradicated in Area 1: When Rs1 ðpÞp1 and Rr1 ðpÞp1, then both the vaccine-sensitive and vaccine-resistant strains are eradicated in Area 1 by Theorem A.2 (E01 is GAS). r r Let 0opp1. It is clear that R¯ 2 ð0ÞoR¯ 2 ðpÞ. Since Rs1 ðpÞo1, we can s s s s obtain R2 ð0Þ4R2 ðpÞ. We can also evaluate R¯ 2 ð0Þ4R¯ 2 ðpÞ because r r r r X 2 oX 2 jp¼0 . Further we can show that R2 ð0Þ4R2 ðpÞ is equivalent to Rs1 ðpÞ þ so1 pc=ðb þ my Þðb þ eÞo1. The last inequality can be satisfied if s is small because Rs1 ðpÞo1. Remember that 1oRs2 ð0Þ, r 1oRr2 ð0Þ and R¯ 2 ð0Þo1. Therefore we can conclude that the vaccination program leads to the following possibilities of change of prevalent strains in Area 2 by Theorem A.4: (i) Rs2 ðpÞo1, Rr2 ðpÞo1; E02 is GAS. r (ii) (1oRs2 ðpÞ,Rr2 ðpÞo1) or (1oRs2 ðpÞ, 1oRr2 ðpÞ, R¯ 2 ðpÞo1); Es2 is GAS. s (iii) (Rs2 ðpÞo1, 1oRr2 ðpÞ) or (1oRs2 ðpÞ, 1oRr2 ðpÞ, R¯ 2 ðpÞo1); Er2 is GAS. s r (iv) 1oR¯ 2 ðpÞ, 1oR¯ 2 ðpÞ; Eþ 2 is GAS. Note that (i) both the vaccine-sensitive and vaccine-resistant strains are eradicated, (ii) the vaccine-sensitive strain is selected, (iii) the vaccine-resistant strain is selected, and (iv) both the vaccine-sensitive and vaccine-resistant strains are selected in Area 2, respectively.

227

f2

()X s1 jp¼0 oX s1 þ sV r1 ()0osV r1 . We have that Rr2 ð0ÞoRr2 ðpÞ. The vaccination program reduces the infection force of vaccine-sensitive strain and intensifies one of vaccine-resistant strain in Area 2. We can conclude that the vaccination program leads to the following possibilities of the change in Area 2 by Theorem A.4: r (i) 1oRs2 , 1oRr2 , R¯ 2 o1; Es2 is GAS. s s r (ii) (R2 o1, 1oR2 ) or (1oRs2 , 1oRr2 , R¯ 2 o1); Er2 is GAS. s r is GAS. (iii) 1oR¯ 2 , 1oR¯ 2 ; Eþ 2

(IV) Both the vaccine-sensitive and vaccine-resistant strains are s r selected in Area 1: When R¯ 1 ðpÞ41 and R¯ 1 ðpÞ41, then both the vaccine-sensitive and vaccine-resistant strains are selected in Area 1 by Theorem A.2 (Eþ 1 is GAS). s s It is clear that Rs2 ð0Þ ¼ Rs2 ðpÞ, Rr2 ð0ÞoRr2 ðpÞ, R¯ 2 ð0Þ4R¯ 2 ðpÞ and r r R¯ 2 ð0ÞoR¯ 2 ðpÞ. We can conclude that the vaccination program leads to the following possibilities of the change in Area 2 by Theorem A.4: r

(i) 1oRs2 ðpÞ, 1oRr2 ðpÞ, R¯ 2 ðpÞo1; Es2 is GAS. s (ii) 1oRs2 ðpÞ, 1oRr2 ðpÞ, R¯ 2 ðpÞo1; Er2 is GAS. s r (iii) 1oR¯ 2 ðpÞ, 1oR¯ 2 ðpÞ; Eþ 2 is GAS.

References Anderson, R.M., May, R.M., 1991. Infectious Disease of Humans: Dynamics and Control. Oxford University Press, Oxford. Bremermann, H.J., Thieme, H.R., 1989. A competitive exclusion principle for pathogen virulence. J. Math. Biol. 27, 179–190. Capua, I., Alexander, D.J., 2004. Human health implications of avian influenza viruses and paramyxoviruses. Eur. J. Clin. Microbiol. Infect. Dis. 23, 1–6. Capua, I., Marangon, S., 2004. Vaccination for avian influenza in Asia. Vaccine 22, 4137–4138. Capua, I., Marangon, S., 2006. Control of avian influenza in poultry. Emerg. Infect. Dis. 12, 1319–1324.

ARTICLE IN PRESS 228

S. Iwami et al. / Journal of Theoretical Biology 259 (2009) 219–228

Capua, I., 2007. Vaccination for notifiable avian influenza in poultry. Rev. Sci. Tech. Off. Int. Epiz. 26, 217–227. Castillo-Chavez, C., Thieme, H.R., 1995. Asymptotically autonomous epidemic models. In: Arino, O., Axelrod, D., Kimmel, M., Langlais, M. (Eds.), Mathematical Population Dynamics: Analysis of Heterogeneity. Theory of Epidemics, vol. 1. Wuerz, pp. 33–50. Chen, H., Deng, G., Li, Z., Tian, G., Li, Y., Jiao, P., Zhang, L., Liu, Z., Webster, R.G., Yu, K., 2004. The evolution of H5N1 influenza viruses in ducks in southern China. Proc. Natl. Acad. Sci. USA 101, 10452–10457. Elbers, A.R., Fabri, T.H., de Vries, T.S., de Wit, J.J., Pijpers, A., Koch, G., 2004. The highly pathogenic avian influenza A (H7N7) virus epidemic in The Netherlands in 2003—lessons learned from the first five outbreaks. Avian Dis. 48, 691–705. Elbers, A.R., Koch, G., Bouma, A., 2005. Performance of clinical signs in poultry for the detection of outbreaks during the avian influenza A (H7N7) epidemic in The Netherlands in 2003. Avian Pathol. 34, 181–187. Gambotto, A., Barratt-Boyes, S.M., de Jong, M.D., Neumann, G., Kawaoka, Y., 2008. Human infection with highly pathogenic H5N1 influenza virus. Lancet 371, 1464–1475. Gilbert, M., Xiao, X., Domenech, J., Lubroth, J., Martin, V., Slingenbergh, J., 2006. Anatidae migration in the Western Palearctic and spread of highly pathogenic avian influenza H5N1 virus. Emerg. Infect. Dis. 12, 1650–1656. Guan, Y., Poon, L.L.M., Cheung, C.Y., Ellis, T.M., Lim, W., Lipatov, A.S., Chan, K.H., Sturm-Ramirez, K.M., Cheung, C.L., Leung, Y.H.C., Yuen, K.Y., Webster, R.G., Peiris, J.S.M., 2004. H5N1 influenza: a protean pandemic threat. Proc. Natl. Acad. Sci. USA 101, 8156–8161. Handel, A., Regoes, R.R., Antia, R., 2006. The role of compensatory mutations in the emergence of drug resistance. PLoS Comput. Biol. 2 (10), e137. Hayden, F.G., 2001. Perspectives on antiviral use during pandemic influenza. Philos. Trans. R. Soc. London B 356, 1877–1884. Iwami, S., Suzuki, T., Takeuchi, Y., 2009. Paradox of vaccination: Is vaccination really effective against avian flu epidemics? PLoS ONE 4 (3), e4915. Kilpatrick, A.M., Chmura, A.A., Gibbons, D.W., Fleischer, R.C., Marra, P.P., Daszak, P., 2006. Predicting the global spread of H5N1 avian influenza. Proc. Natl. Acad. Sci. USA 103, 19368–19373. Lam, T.T., Hon, C.C., Pybus, O.G., Kosakovsky, S.L., Wong, R.T., Yip, C.W., Zeng, F., Leung, F.C., 2008. Evolutionary and transmission dynamics of reassortant H5N1 influenza virus in Indonesia. PLoS Pathog. 4 (8), e1000130. Lee, C.W., Senne, D.A., Suarez, D.L., 2004. Effect of vaccine use in the evolution of Mexican lineage H5N2 avian influenza virus. J. Virol. 78, 8372–8381. Leung, F.C., 2007. Comments on the Fujian-Like strain of avian influenza H5N1. Poult. Sci. 86, 435–436. Li, K.S., Guan, Y., Wang, J., Smith, G.J.D., Xu, K.M., Duan, L., Rahardjo, A.P., Puthavathana, P., Buranathai, C., Nguyen, T.D., Estoepangestie, A.T.S., Chaisingh, A., Auewarakul, P., Long, H.T., Hanh, N.T.H., Webby, R.J., Poon, L.L.M., Chen, H., Shortridge, K.F., Yuen, K.Y., Webster, R.G., Peiris, J.S.M., 2004. Genesis of a highly pathogenic and potentially pandemic H5N1 influenza virus in eastern Asia. Nature 430, 209–213. Lipsitch, M., Cohen, T., Murray, M., Levin, B.R., 2007. Antiviral resistance and the control of pandemic influenza. PLoS Med. 4 (1), e15. Marangon, S., Cecchinato, M., Capua, I., 2008. Use of vaccination in avian influenza control and eradication. Zoo. Public Health 55, 65–72. Moghadas, S.M., Bowman, C.S., Rost, G., Wu, J., 2008. Population-wide emergence of antiviral resistance during pandemic influenza. PLoS ONE 3 (3), e1839. Normile, D., 2005. Are wild birds to blame? Science 310, 426–428.

Pasquato, A., Seidah, N.G., 2008. The H5N1 influenza variant Fujian-like hemagglutinin selected following vaccination exhibits a compromised furin cleavage: neurological consequences of highly pathogenic Fujian H5N1 strains. J. Mol. Neurosci. 35, 339–343. Peyre, M., Fusheng, G., Desvaux, S., Roger, F., 2008. Avian influenza vaccines: a practical review in relation to their application in the field with a focus on the Asian experience. Epidemiol. Infect. 14, 1–21. Regoes, R.R., Bonhoeffer, S., 2006. Emergence of drug-resistant influenza virus: population dynamical considerations. Nature 312, 389–391. Russell, C.A., Jones, T.C., Barr, I.G., Cox, N.J., Garten, R.J., Gregory, V., Gust, I.D., Hampson, A.W., Hay, A.J., Hurt, A.C., de Jong, J.C., Kelso, A., Klimov, A.I., Kageyama, T., Komadina, N., Lapedes, A.S., Lin, Y.P., Mosterin, A., Obuchi, M., Odagiri, T., Osterhaus, A.D.M.E., Rimmelzwaan, G.F., Shaw, M.W., Skepner, E., Stohr, K., Tashiro, M., Fouchier, R.A.M., Smith, D.J., 2008. The global circulation of seasonal influenza A (H3N2) viruses. Science 320, 340–346. Savill, N.J., St Rose, S.G., Keeling, M.J., Woolhouse, M.E., 2006. Silent spread of H5N1 in vaccinated poultry. Nature 442, 757. Seo, S.H., Webster, R.G., 2001. Cross-reactive, cell-mediated immunity and protection of chickens from lethal H5N1 influenza virus infection in Hong Kong poultry markets. J. Virol. 75, 2516–2525. Skeik, N., Jabr, F.I., 2007. Influenza viruses and the evolution of avian influenza virus H5N1. Int. J. Infect. Dis. 12, 233–238. Smith, G.J.D., Fan, X.H., Wang, J., Li, K.S., Qin, K., Zhang, J.X., Vijaykrishna, D., Cheung, C.L., Huang, K., Rayner, J.M., Peiris, J.S.M., Chen, H., Webster, R.G., Guan, Y., 2006. Emergence and predominance of an H5N1 influenza variant in China. Proc. Natl. Acad. Sci. USA 103, 16936–16941. Stegeman, A., Bouma, A., Elbers, A.R.W., de Jong, M.C.M., Nodelijk, G., de Klerk, F., Koch, G., van Boven, M.J., 2004. Avian influenza A virus (H7N7) epidemic in The Netherlands in 2003: course of the epidemic and effectiveness of control measures. J. Infect. Dis. 190, 2088–2095. Stilianakis, N.I., Perelson, A.S., Hayden, F.G., 1998. Emergence of drug resistance during an influenza epidemic: insights from a mathematical model. J. Infect. Dis. 177, 863–873. Swayne, D.E., Beck, J.R., Garcia, M., Stone, H.D., 1999. Influence of virus strain and antigen mass on efficacy of H5 avian influenza inactivated vaccines. Avian Pathol. 28, 245–255. Tiensin, T., Chaitaweesub, P., Songserm, T., Chaisingh, A., Hoonsuwan, W., Buranathai, C., Parakamawongsa, T., Premashthira, S., Amonsin, A., Gilbert, M., Nielen, M., Stegeman, A., 2005. Highly pathogenic avian influenza H5N1, Thailand, 2004. Emerg. Infect. Dis. 11, 1664–1672. Tiensin, T., Nielen, M., Vernooij, H., Songserm, T., Kalpravidh, W., Chotiprasatintara, S., Chaisingh, A., Wongkasemjit, S., Chanachai, K., Thanapongtham, W., Srisuvan, T., Stegeman, A., 2007. Transmission of the highly pathogenic avian influenza virus H5N1 within flocks during the 2004 epidemic in Thailand. J. Infect. Dis. 196, 1679–1684. van den Berga, T., Lambrechta, B., Marche, S., Steenselsa, M., Borma, S.V., Bublot, M., 2007. Influenza vaccines and vaccination strategies in birds. Comp. Immunol. Microbiol. Infect. Dis. 31, 121–165. Vijaykrishna, D., Bahl, J., Riley, S., Duan, L., Zhang, J.X., Chen, H., Peiris, J.S.M., Smith, G.J.D., Guan, Y., 2008. Evolutionary dynamics and emergence of panzootic H5N1 influenza viruses. PLoS Pathog. 4 (9), e1000161. Webster, R.G., Peiris, M., Chen, H., Guan, Y., 2006. H5N1 outbreaks and enzootic influenza. Emerg. Infect. Dis. 12, 3–8. Yee, K.S., Carpenter, T.E., Cardona, C.J., 2009. Epidemiology of H5N1 avian influenza. Comp. Immunol. Microbiol. Infect. Dis. 32, 325–340.