## Abstract

The spread of H5N1 avian influenza and the recent high numbers of confirmed human cases have raised international concern about the possibility of a new pandemic. Therefore, antiviral drugs are now being stockpiled to be used as a first line of defence. The large-scale use of antivirals will however exert a strong selection pressure on the virus, and may lead to the emergence of drug-resistant strains. A few mathematical models have been developed to assess the emergence of drug resistance during influenza pandemics. These models, however, neglected the spatial structure of large populations and the stochasticity of epidemic and demographic processes. To assess the impact of population structure and stochasticity, we modify and extend a previous model of influenza epidemics into a metapopulation model which takes into account the division of large populations into smaller units, and develop deterministic and stochastic versions of the model. We find that the dynamics in a fragmented population is less explosive, and, as a result, prophylaxis will prevent more infections and lead to fewer resistant cases in both the deterministic and stochastic model. While in the deterministic model the final level of resistance during treatment is not affected by fragmentation, in the stochastic model it is. Our results enable us to qualitatively extrapolate the prediction of deterministic, homogeneous-mixing models to more realistic scenarios.

## 1. Introduction

Recent human cases of H5N1 avian influenza have raised public concern. It is now considered an established fact that another influenza pandemic will occur sooner or later (e.g. Webby & Webster 2003). Given that our world is highly interconnected and interdependent (Hufnagel *et al.* 2004), this pandemic will certainly spread fast and widely.

Two types of antiviral drugs currently exist: the M2-inhibitors (M2I) and the neuraminidase inhibitors (NAI). M2Is target the M2 protein, an ion channel enabling the uncoating of the virion within the cell, which is a key step of the replication process. M2Is are only effective against influenza A (Stiver 2003). NAIs belong to a new class of antivirals including oseltamivir (*Tamiflu*), and zanamivir (*Relenza*). These drugs target the neuraminidase, which is required to release the budding virus from the cell and to prevent viral particles from clumping. NAIs therefore attenuate the infection.

A consequence of the use of antiviral drugs, however, is the emergence of drug resistance. Drug resistance in general is an urgent public health problem, costing human lives and money (Heymann 2006). Resistance against M2Is has risen rapidly over the past few years (Bright *et al.* 2006). The increasing incidence of M2I-resistant influenza A seems to be due to a high de novo resistance rate associated with low transmission fitness costs (Regoes & Bonhoeffer 2006). Given the chemical structure of oseltamivir, structural studies predicted some time ago that resistance to the drug could develop, allowing the NA to maintain its function (Moscona 2005). Several clinical cases have recently confirmed this prediction (de Jong *et al.* 2005). More alarmingly, transmission studies in ferrets and *in vitro* show that oseltamivir-resistant virus variants may differ in fitness and transmissibility. The *R*292*K* mutation, which also shows cross-resistance to zanamivir, hampers growth and transmission of the virus, because it reduces the sialidase activity. Viruses carrying the *E*119*V* mutation, however, have growth and transmission rates similar to wild-type viruses (Yen *et al.* 2005). In addition, it cannot be ruled out that a compensatory mutation could restore the fitness of other drug-resistant viruses.

These recent findings, combined with the fact that NAI drugs are being stockpiled world-wide underline the need for *in vivo*, *in vitro* and *in silico* studies to assess the emergence and spread of drug-resistant influenza viruses, under different public health interventions and patterns of drug use. In particular, computer simulations can provide guidelines for the most effective use of antiviral drugs (Stilianakis *et al.* 1998; Smith 2003), in terms of limiting both new infections and the emergence of drug resistance.

Recently, many complex models have been developed, that investigate the outcome of different public health interventions for controlling an influenza pandemic at the source (Ferguson *et al.* 2005; Longini *et al.* 2005), or for mitigating it once it has started to spread (Ferguson *et al.* 2006; Germann *et al.* 2006). These interventions include treatment and prophylaxis with antiviral drugs, vaccination, quarantine, restrictions on travel, school closure and other non-pharmaceutical interventions. Although these studies often mention the risk of the emergence and spread of drug-resistant viruses, none of them explicitly investigates the dynamics of resistance emergence. Some models, however, have already assessed the potential spread of drug-resistant viruses for M2Is (Stilianakis *et al.* 1998) as well as oseltamivir (Ferguson *et al.* 2003; Regoes & Bonhoeffer 2006). None of these studies took into account the spatial population structure, although Ferguson *et al.* (2003) incorporated population structure in terms of risk groups and age structure.

In this paper, we investigate the emergence and spread of oseltamivir-resistant influenza virus in a population structured into subpopulations (metapopulation), when different public-health interventions are performed. The work extends the studies of Stilianakis *et al.* (1998) and Regoes & Bonhoeffer (2006). First, we modify the initial model, and then extend it to a metapopulation. Then, we develop a stochastic version of the model. We implement spatial structure in a very abstract and generic way, rather than attempting to mimic the real world situation. This allows us to qualitatively investigate how spatial structure influences the epidemic and the efficacy of interventions.

## 2. Models and methods

We extended a previously published model describing an influenza outbreak in a school (Stilianakis *et al.* 1998) into a metapopulation model, i.e. we linked multiple epidemic foci by allowing individuals to migrate between them. We developed deterministic and stochastic versions of this metapopulation model.

### 2.1 Structure of the community

Different metapopulation structures can be studied in epidemiology, depending on what is defined as a patch (Hess *et al.* 2002). Given the organization of humans in well-defined spatial structures such as families, villages and towns (Grenfell & Harwood 1997), patches are groups of hosts in our model. Colonization of a patch corresponds to the establishment of an infection in a patch harbouring susceptible individuals. Local extinction of the disease occurs when all the infected hosts within the patch are removed (through recovery or death). It is also assumed that there is a homogeneous mixing within each subpopulation (figure 1).

We connect the disease dynamics within the patches by allowing individuals to migrate from one patch to the other. Mathematically, the migration process of individuals between patches is described as migration matrix , where *μ*_{ij} is the rate of migration of individuals from population *i* to population *j*. This matrix determines the topology and the strength of the connections between the patches.

We model different population structures. Throughout this paper, the term *population structure* refers to the number of subpopulations, and to how they are connected (i.e. which subpopulations are connected, and what is the value of the migration rate). For a given number *n* of subpopulations, two types of subpopulation connections are investigated here (figure 2*a*). In the *island*-type of connections, all the subpopulations are interconnected; in the *spider*-type, *n*−1 subpopulations are connected only to a central one, and the infection starts in a peripheral subpopulation. For a given structure, a range of migration coefficients (the nonzero coefficients in the migration matrix) and of number of subpopulations (the size of the matrix) are investigated (figure 2*b*). We assume that migration is the same among connected subpopulations, which means thatand that all subpopulations have the same size *N*/*n*, where *N* is the total number of individuals and *n* the number of subpopulations. We chose a population size of *N*=500 because larger population sizes would have rendered simulations too long to run. A population of 500 individuals might not seem realistic for an influenza pandemic. We made sure, however, that the size of the population had no impact on our conclusions. (See appendix C2 for simulations with a bigger population size.)

There is no age or social structure in the population and no seasonal variation of parameters (because we focus on the first 100 days of the epidemic only).

### 2.2 Epidemic parameters

The parameters of the model are listed in appendix A and were chosen to match with recent studies on oseltamivir-resistant viruses (Gubareva *et al.* 2001; Kiso *et al.* 2004).

*R*_{0}, the basic reproductive number, is a key concept in epidemiology (Heffernan *et al*. 2005). It is defined as the number of secondary cases produced by a single infected individual during its entire infectious period, in a fully susceptible population (Diekmann *et al*. 1990; Anderson & May 1991). The transmission rates (the *β*s) are scaled so that the basic reproductive number of the sensitive strain, , is equal to 4.5 without drugs, and , the basic reproductive number of the resistant strain, equals 0.45. The next generation method (Diekmann *et al.* 1990; Heffernan *et al.* 2005), which is described in appendix B, was used to determine the under different interventions (table 1).

Our is high when compared with the range used in recent models (always lower than 2.5 in Ferguson *et al.* (2005, 2006), Longini *et al.* (2005), Germann *et al.* (2006)). However, our value was estimated by Stilianakis *et al.* (1998) from data collected during an outbreak in a boarding school. It may be higher because a boarding school population is more homogeneously mixed than that of an entire town. If one estimated on the basis of data generated by our model (with ) using a homogeneously mixed model (e.g. Mills *et al.* 2004), the estimate of would be lower than 4.5 because the disease spreads slower in our model than in a homogeneous-mixing model. Hence, our choice of is not inconsistent with epidemiological data.

For we assumed a 10-fold resistance cost. Given the fact that oseltamivir-resistant strains found in humans to date are not transmissible (Herlocher *et al.* 2002, 2004), our may appear too high. However, because the large-scale use of oseltamivir may soon select for transmissible resistance, and compensatory mutations may reduce potential costs even further, one may argue that our is too low to describe oseltamivir-resistant influenza strains that may emerge during a future pandemic. With our choice of we try to strike a compromise between these two opposing views. Our results, however, seem to be robust with regard to the value of (see simulations performed with higher values of in appendix C1).

### 2.3 Deterministic model

Although containing many more classes, the model that we use to study the emergence of resistance is based on the classical *SIR* model (Kermack & McKendrick 1927). The set of equations describing the epidemic without any public health intervention corresponds to an *SII*_{s}*R* model. Infected people are divided into asymptomatic (*I*) and symptomatic (*I*_{s}) classes. Both of them are infectious, but their infectiousness varies, and only those who develop clinical symptoms will be treated, when treatment is provided. See appendix A for a complete description of the equations.

Two kinds of interventions will be considered.

*Treatment*. All symptomatically infected individuals receive treatment (oseltamivir).Infected people receiving treatment are less infective, develop symptoms more slowly and recover faster.

*Prophylaxis*. In addition to treatment of all symptomatic individuals, uninfected as well as asymptomatically infected individuals (as they appear to be uninfected) receive prophylaxis (oseltamivir).Prophylaxed susceptibles become less susceptible to infection, while asymptomatically infected people are less infectious and recover faster.

While treatment is a continuous intervention, prophylaxis is not. At a given time, all the susceptible and asymptomatically infected individuals receive prophylaxis, and are then moved to the corresponding treated class.

Neither treatment nor prophylaxis has any effect on people infected by a resistant virus. Note that throughout the paper, a *resistant virus* refers to a drug-resistant virus, and that a *resistant case* refers to an individual shedding drug-resistant viruses.

More classes are required when interventions are performed, in order to represent individuals receiving treatment or prophylaxis.

Infected individuals are distinguished according to whether they are symptomatically infected (s) or not, treated (tr) or not, and have developed resistance (r) or not (figure 3).

Susceptible individuals can receive prophylaxis (pr) or not.

For each class *C*, *C*∈{*S*,*S*_{pr},*I*,…,*I*_{s,}_{r,}_{tr}}, *C*_{i} refers to the number of individuals of this class present in subpopulation *i*, 1≤*i*≤*n*. Moreover, migration terms are added to the equations described in Regoes & Bonhoeffer (2006). For *C*_{i}, this migration parameter is of the following form:(2.1)The first term refers to individuals of class *C* leaving subpopulation *i* for subpopulation *j*, while the second term refers to individuals arriving (see appendix A3 for the whole equations).

The ODEs were solved numerically using R's odesolve package (R Development Core Team 2005; Setzer 2005). We also determined how many people are infected during an epidemic (i.e. ‘total number of infected cases’), among them, how many shed resistant viruses (i.e. ‘total number of resistant cases’), and among the latter, what is the proportion of resistance that appeared de novo (i.e. appeared during infection, as opposed to transmitted resistance). See appendix A3 for formulae.

### 2.4 Stochastic model

The deterministic version of the model is made stochastic by using an implementation of the Gillespie algorithm (Gillespie 1976), which numerically simulates a Markov process. The Gillespie algorithm was coded in R for the purpose of this study. For each set of parameters, at least 100 stochastic simulations are run using this algorithm.

Since computing time increases exponentially with the division of the population (i.e. the number of subpopulations, or the number of different equations), simulations were run on a cluster of computers. Parallelization was performed *by hand*, by submitting different jobs corresponding to different sets of parameters to individual nodes of the computer cluster.

## 3. Results

### 3.1 Deterministic model

In this section, we investigate the impact of population structure and of different interventions on the dynamics and the final outcome of the epidemics, focusing on the development of drug resistance. The deterministic version of the model is used here.

#### 3.1.1 Dynamics and final outcome of the epidemics after treatment

How does the structure of the population affect the dynamics of the epidemic as compared with a homogeneously mixed model? To answer this question, and for the sake of simplicity, we use the epidemic's peak time and value (i.e. the maximal instantaneous number of infected individuals in the whole population) to characterize the dynamics of the epidemic.

Figure 4 shows dynamical and final features of the epidemic for different population structures. For a larger population of 5000 individuals, figure 4 does not change. This is due to rescaling of the transmission rates such that the basic reproductive number remains the same (equation A 9). Figure 4*a*,*b*,*d*,*e* shows that the epidemic peak decreases and occurs later with decreasing migration rates. This means that the epidemic is less explosive in a fragmented population than in a homogeneous population. Treatment slows down the epidemics and reduces its peak. As expected, the level and the timing of the peak converge to the homogeneous mixing case for high migration rates.

Some results in figure 4 might appear counterintuitive: the peak time and value with 15 subpopulations (crosses) are between the values for 1 (circles) and 5 (triangles) subpopulations. As discussed later in §3, when the population is subdivided into more than one subpopulation connected according to the island model, there are two peaks, corresponding to epidemics in the focal subpopulation (first peak) and later in the peripheral subpopulations (second peak). The second peak in a model with a fragmented population is generally smaller than the only peak in a model with a homogeneous population because the total number of susceptible individuals in the peripheral subpopulations is smaller than that of the total population. This explains why the peak value in the non-fragmented population is highest. With increased numbers of subpopulations, however, the total number of susceptible individuals in the peripheral subpopulations increases. That is why the peak value is higher for 15 subpopulations than for 5 (figure 4*b*,*e*). The results with regard to the time of the peak can be explained as follows: in a single subpopulation, the speed of the epidemic decreases with its size. Hence, if the population is more fragmented the epidemics spread faster in each subpopulation. This accounts for the result concerning the peak time (figure 4*a*,*d*).

The fragmentation of the population, however, has no impact on the final outcome of the epidemic, in terms of the total number of infected and resistant cases (figure 4*c*,*f*). This is due to the fact that the total numbers of infected (*I*) and resistant (*R*) cases solely depend on the basic reproductive number, *R*_{0}, of each strain, which is independent of the migration rate. Without drugs almost everyone (98.8%) gets infected. Treatment slightly decreases this value to 91%. Nevertheless, treatment also allows the emergence of resistant cases, which represent 3.8% of the infected cases.

These results raise the following questions: can prophylaxis more efficiently decrease the total number of infected cases? Does it lead to more resistance than treatment? And lastly, when should prophylaxis be performed?

#### 3.1.2 Prophylaxis: the impact of timing

The total number of infected and resistant cases during the epidemic is lower when prophylaxis is performed (figure 5). The results are similar with a bigger population size. Unlike with treatment, the outcome of prophylaxis depends on the structure of the population. The higher the degree of fragmentation of the host population is, the fewer is the number of people getting infected during the epidemic (figure 5*a*). For instance, with our parameters, and with proph.start=15 (where proph.start denotes the time at which prophylaxis is initiated), 80% of the population finally gets infected during the epidemic in a homogeneously mixing population, while only 5% gets infected when the number of subpopulations is high or the migration is low. In both cases, around 4% of infected people shed resistant viruses. This is due to the shape of the epidemic curve: the higher the degree of fragmentation, the more peaks there are. Note that there are at most two peaks with an island shape: the first peak corresponds to the focal subpopulation, where the epidemic starts, and the second peak to all the other subpopulations, which are directly linked to the focal one. Note also that if the mean distance to the focal subpopulation is greater than 1 (e.g. with the spider shape, data not shown), there are more than two peaks. The state of the epidemic at which it is ‘hit’ by prophylaxis (peak or trough, figure 5*b*–*e*) depends on the structure of the population. After prophylaxis is initiated both and become abruptly less than 1, and the epidemic wanes, because both resistant and sensitive variants cannot sustain the epidemic anymore.

Owing to the fact that the structure of the population influences at which state the epidemic is hit by prophylaxis, the timing of prophylaxis has an impact on resistance emergence (figure 6). The earlier the prophylaxis is initiated, the fewer are the infected and resistant cases during the epidemic and the greater the proportion of infected people shedding resistant viruses. The fact that when everybody receives drugs accounts for this higher proportion of resistance when proph.start is smaller. Even though unable to trigger a real epidemic, the resistant variant outcompetes the sensitive one. We also determined the proportion of transmitted resistance versus de novo resistance. This proportion is substantially higher when proph.start is small, again because sensitive variants are outcompeted earlier.

### 3.2 Stochastic model

#### 3.2.1 The impact of population structure

Even in homogeneous populations, stochastic simulations differ from deterministic ones in some respects. First, around 40% of the simulations did not generate any secondary case (the first infected is an *I*, asymptomatic, non-treated, shedding sensitive viruses). Then, due to local extinctions, the stochastic model predicts fewer infected cases on average than the deterministic model (figures 7 and 8). Furthermore, the total number of infected people during the epidemic is lower in the stochastic than in the deterministic model (e.g. approximately 80% of the population compared with 98.8% with the deterministic version).

Unlike the deterministic model, the structure of the population has an impact on the outcome of the epidemic even without drugs or under treatment (figure 9). The higher the degree of fragmentation of the host population, the less the infection spreads. When the number of subpopulations is high, or when migration between subpopulations is low, the epidemic is often restricted to a small number of subpopulations (most of the time to only the one where the epidemic started). Comparing the effect of different interventions, we find that treatment does not substantially reduce the total number of infected cases, but leads to the emergence of resistance (figure 9). On the other hand, as in the deterministic model, prophylaxis has a strong impact on the total number of resistant cases. We also find that the proportion of resistant cases among infected people is higher with prophylaxis (between 3 and 6% when proph.start=10). This result is, as in the deterministic case, due to the change of the basic reproductive numbers when different interventions are performed.

#### 3.2.2 Prophylaxis: impact of timing

In spite of the variance in the prediction due to stochasticity, the impact of prophylaxis timing is qualitatively the same as in the deterministic model. The earlier prophylaxis starts, the smaller is the number of infected and resistant people during the epidemic (figure 10). For a larger population of *N*=5000 we obtain very similar results (see figure 13 in appendix A).

The result concerning the proportion of resistant cases is less clear in the stochastic model. Although, as in the deterministic model, the mean number of resistant cases is higher for prophylaxis starting early, the variance is also higher. In some simulations, the proportion of resistance is very low (for instance *R* : *I*=0 : 10 in one simulation), while it is very high in some others (*R* : *I*=6 : 9 in another). Moreover, among simulations that led to the emergence of resistance, the proportion of transmitted resistant strains decreases with the timing of prophylaxis: approximately 30% of transmitted resistance when prophylaxis starts at day 0, approximately 20% at day 5, but less than 5% at day 20 (data not shown). Note that there is almost no transmitted resistance when treatment is performed. In brief, if resistance appears when prophylaxis starts early, resistant strains are more likely to be transmitted than if prophylaxis is initiated later.

## 4. Discussion and conclusion

The planned, widespread use of NAI for treatment and prophylaxis during influenza pandemics is likely to lead to the emergence of drug-resistant influenza strains. In this study, we investigated the effect of population structure on the emergence of drug resistance during pandemics using mathematical models. To this end, we constructed deterministic and stochastic metapopulation models describing the spread of pandemic influenza strains.

We show that epidemic outbreaks in fragmented populations are less explosive. In a deterministic model, however, the fragmentation of the population has no influence on the final level of drug resistance in the population when nothing is done or when drugs are given to symptomatic individuals. In contrast, when antivirals are delivered prophylactically to non-symptomatic individuals in addition to treating symptomatic individuals, population structure and the timing of prophylaxis interact synergistically in resistance development and influence the proportion of resistance due to the transmission of resistant variants.

The results of the stochastic model are different in some respects. First, 40% of the epidemics go extinct very rapidly. Moreover, when the epidemic does not go extinct, the number of infected and resistant cases is still lower than in the deterministic model. This is due to an interaction of stochasticity and migration. In addition, the final level of resistance after treatment only is affected by the structure of the population, and not by stochasticity.

This study should be considered as a conceptual work, aiming to study qualitative trends, rather than to make quantitative predictions. Even though many of our model parameters are based on published studies on influenza (see Regoes & Bonhoeffer (2006) for a justification of the parameters), the topological features and parameters describing the structure of the population are unrealistic. However, the main aim of our paper was to study the effect of population fragmentation in isolation, rather than predict the level of resistance for a realistically structured population. To this end, we extended one of the simplest published influenza models (Stilianakis *et al.* 1998) into one which allowed for varying degrees of population fragmentation. This conceptual study has its justification because it allows us to extrapolate the predictions of simple models—assuming homogeneously mixed host populations and deterministic dynamics—to situations in which spatial population structure and stochasticity are important. To make this model more realistic, one could analyse more realistic network topologies and population structures. Additionally, one could also include aspects of the within-host dynamics of influenza infection and age structure, in particular because some key parameters of influenza epidemics vary with age (Ferguson *et al.* 2003; Doherty *et al.* 2006).

In spite of the lack of realism of the population structure in our model, our conclusions are similar to the ones made using very detailed models of pandemic influenza. Ferguson *et al.* (2005, 2006), Longini *et al.* (2005) and Germann *et al.* (2006) show that targeted antiviral prophylaxis is the best medical intervention, and that if *R*_{0} is high, a combination of medical and non-medical interventions (such as social distancing) have to be adopted to mitigate the pandemic. Considering social distancing as a measure that increases population fragmentation, we show that prophylaxis is the best intervention and that its effects are increased by population fragmentation (figures 5 and 9). Lipsitch *et al.* (2007) also studied the emergence of resistance, albeit without incorporating any spatial structure of the population. One of their key results is that resistance levels are higher when non-drug interventions are adopted. Hereby, they modelled non-drug interventions by simply reducing the basic reproduction number *R*_{0} of the pandemic strain. In our model, we can model non-drug interventions, such as social distancing or quarantine, by increasing the population fragmentation. In a more fragmented population, similar to Lipsitch *et al.* (2007), we observe that the prevalence of resistance is reduced, while the proportion of resistance increases (figure 11). This is due to the fact that in a fragmented population the epidemic spreads slower and prophylaxis hits the epidemic during an earlier state (figure 5*c*–*e*).

Only few studies on the containment or mitigation of an influenza pandemic investigate the emergence of resistant viruses, but most authors acknowledge the challenge that the emergence of resistant viruses represents (Ferguson *et al.* 2005, 2006; Germann *et al.* 2006). We hope that the potential emergence of drug resistance will be incorporated into more mathematical studies (especially the very detailed ones) so that we can face up to this challenge.

## Footnotes

One contribution of 20 to a Theme Issue ‘Cross-scale influences on epidemiological dynamics: from genes to ecosystems’.

- Received March 23, 2007.
- Accepted June 6, 2007.

- © 2007 The Royal Society