Impact of Emerging Antiviral Drug Resistance on Influenza Containment and Spread: Influence of Subclinical Infection and Strategic Use of a Stockpile Containing One or Two Drugs

Background Wide-scale use of antiviral agents in the event of an influenza pandemic is likely to promote the emergence of drug resistance, with potentially deleterious effects for outbreak control. We explored factors promoting resistance within a dynamic infection model, and considered ways in which one or two drugs might be distributed to delay the spread of resistant strains or mitigate their impact. Methods and Findings We have previously developed a novel deterministic model of influenza transmission that simulates treatment and targeted contact prophylaxis, using a limited stockpile of antiviral agents. This model was extended to incorporate subclinical infections, and the emergence of resistant virus strains under the selective pressure imposed by various uses of one or two antiviral agents. For a fixed clinical attack rate, R 0 rises with the proportion of subclinical infections thus reducing the number of infections amenable to treatment or prophylaxis. In consequence, outbreak control is more difficult, but emergence of drug resistance is relatively uncommon. Where an epidemic may be constrained by use of a single antiviral agent, strategies that combine treatment and prophylaxis are most effective at controlling transmission, at the cost of facilitating the spread of resistant viruses. If two drugs are available, using one drug for treatment and the other for prophylaxis is more effective at preventing propagation of mutant strains than either random allocation or drug cycling strategies. Our model is relatively straightforward, and of necessity makes a number of simplifying assumptions. Our results are, however, consistent with the wider body of work in this area and are able to place related research in context while extending the analysis of resistance emergence and optimal drug use within the constraints of a finite drug stockpile. Conclusions Combined treatment and prophylaxis represents optimal use of antiviral agents to control transmission, at the cost of drug resistance. Where two drugs are available, allocating different drugs to cases and contacts is likely to be most effective at constraining resistance emergence in a pandemic scenario.


S1.1 The single-drug model
We have states I y p with y ∈ {w, r} representing symptomatic infections for those taking prophylaxis (p) and infected with the wild-type (y = w) or resistant-type (y = r) strain. States A y p represent equivalent asymptomatic infections.
We also have the state I w np,t representing symptomatic infections with the wild-type strain for those not taking prophylaxis (np), but being provided with antiviral drugs as treatment (t). For symptomatic infections with the wild-type strain that are not provided with antiviral drugs at all, we have the state I w np,nt . Symptomatic infections with the resistant-type strain may or may not be provided with treatment but it has no effect so we require just one state I r np . Asymptomatic infections are not able to be identified and provided with antiviral drugs and so we only require two states A y np with y ∈ {w, r}. The force of infection arises from each of these nine infectious classes: λ w p = βe i I w p + χA w p (S1) λ w np,nt = β I w np,nt + χA w np (S2) λ w np,t = βe t I w np,t (S3) λ r = βφ I r p + I r np + χ A r p + A r np , where φ is the relative transmissibility of the resistant strain relative to uncontrolled wild-type and χ is the relative transmissibility of asymptomatic infections. Other parameters are defined as in McCaw and McVernon (2007). Infections arising from these nine states cause susceptible individuals (S) to become exposed (E y x with x ∈ {p, np} and y ∈ {w, r}). Exposed individuals will become infectious after a latent period of mean duration 1/ω days. New infections (I and A states) have contacts (C y x states) which may or may not be exposed. The mean infectious period is 1/γ days.
A susceptible, upon exposure, can proceed down a number of paths, dependent upon the characteristics of the infectee (wild-type, resistant-type) and whether or not they are provided with post-exposure prophylaxis. A proportion of exposures, α, are symptomatic. Our model allows for this proportion to depend upon the nature of the strain (w or r) and the provision of prophylaxis (p or np).
Resistance development ("seeding") is captured by two parameters. A proportion, ρ t , of those exposed to wild-type virus and treated (α w np ψE w np ) convert to become symptomatic resistant-strain infectors (I r np ). A proportion, ρ p , of those exposed to wild-type virus and taking prophylaxis (E w p ) convert to become either symptomatic or asymptomatic resistant-strain infectors (I r p or A r p ). Figures S1 and S2 show the possible paths a susceptible can take on becoming infectious. These diagrams are used to construct the model equations that follow. (1 − ψ) We introduce four contact states C y x (x ∈ {p, np}, y ∈ {w, r}) which classify a contact by their own type if they were to become infectious. We introduce: Contacts of wild-type infectives who are provided with AV drugs (C w p ) have a reduced susceptibility accounted for by the factor e s in the equation for Θ w p . All other contacts (those not on prophylaxis and those on prophylaxis but in contact with resistant-strain infectives) remain fully susceptible to infection.

S1.1.1 ODEs
The model equations are: and, for x ∈ {p, np}, y ∈ {w, r} and, for symptomatic infectious states and, for asymptomatic infectious states and, for x ∈ {p, np}, y ∈ {w, r} dR y and finally, for the depletion of the finite stockpile of AV drugs, Define the column vector ( is transpose) We may ignore the E states by factoring them into the right hand side. Defining this construction (˙ Y ≡˙ X + appropriate parts of˙ E) allows us to write: and R 0 is the maximum eigenvalue of the 9 × 9 matrix. We construct this matrix as a linear combination of tensor products of force of infection row vectors, λ w and λ r , with suitable column vectors constructed from the differential equations, F w and F r : where and where K w and K r are the steady state proportions of C x p / C x p + C x np (x ∈ {w, r}), which, in the case α w p = α w np , are simply α x np (x ∈ {w, r}). The matrix F w λ w + F r λ r can be put into block-triangular form, effectively breaking the eignevalue calculation into two separate calculations, one from each of F w λ w and F r λ r . As each of these matrices is, by construction, rank-1 the eigenvalues of the system are given by the scalar product of F x and λ x (x ∈ {w, r}). From x = w, we obtain the reproduction number for the controlled wild strain; from x = r, the reproduction number for the reduced-fitness resistant strain: This expression deserves some explanation. The first eigenvalue is for the wild-type strain. It's form is the appropriate linear combination of prophylaxis (K w ) and non-prophylaxis (1 − K w ) sectors. The prophylaxis sector is itself a linear combination of symptomatic and asymptomatic parts, capturing the reduced susceptibility (e s ) and reduced infectivity (e i ) due to prophylaxis, as well as the leakage to the resistant sector. The non-prophylaxis sector is, again, a linear combination of symptomatic and asymptomatic parts. Furthermore, the symptomatic part is itself a linear combination of treated (ψ) and untreated (1 − ψ) components. The treatment part captures leakage to the resistant sector.
The second reproduction number is for the resistant strain. It is the appropriate linear combination of reproduction numbers for the symptomatic resistant-strain infections (βφ/γ) and asymptomatic resistant-strain infections (βφχ/γ). The simple form arises due to the fact that antiviral drugs have no impact on the resistant strain.

S1.2 Multi-drug models
The multi-drug models build upon the single-drug model just detailed. More states are included to account for the four possible strains in cicurlation (wild-type, two single-drug resistant types and multi-drug resistant type). Parameters relating to the resistant strain (φ, ρ t and ρ p as well as α r ) pick up labels to tie them to the drug in use and the strain in circulation.

S1.2.1 Random allocation and cycling models
We have four strains in circulation: with wild-type virus (w label), drug 1 resistant virus (r 1 label), drug 2 resistant virus (r 2 label) and multiple resistant virus (r 12 label). We introduce new parameters: • φ r 1 ,r 2 ,r 12 -the transmissibility of the (drug 1, drug 2, multiple)-resistant strain relative to the transmissibility of the uncontrolled wild-type strain.
• ρ t1,t2 -the proportion of offspring (that is, secondary infectives) of treated wild-type infectives who carry the drug 1 (drug 2) resistant strain.
We assume that all offspring of r 1 -resistant infectives are either r 1 -resistant or r 12 -resistant. Similarly, all offspring of r 2 -resistant infectives are either r 2 -resistant or r 12 -resistant. All offspring of r 12 -resistant infectives are r 12 -resistant.
The provision of AV drug 1 (drug 2) as prophylaxis to contacts of drug 2 (drug 1) resistant infectives is assumed to have the same protective effect as provision of prophylaxis to contacts of wild-type infectives. Prophylaxis with drug 1 (drug 2) of contacts of drug 1 (drug 2) resistant infectives has no protective effect. Similarly, the provision of AV drug 1 (drug 2) as treatment to drug 2 (drug 1) resistant infectives is assumed to have the same effect as provision of treatment to wild-type infectives. Treatment with drug 1 (drug 2) has no effect on infectives carrying the drug 1 (drug 2) resistant strain. Provision of AV drugs as prophylaxis to contacts or as treatment to infectives carrying the multiple-resistant strain is assumed to have no effect on either infectiousness or susceptibility.
In the random allocation model, which drug is used for treatment or prophylaxis is chosen at random, assumed to be a function of the stockpile size of each drug. For treatment or prophylaxis, drug 1 will be provided a proportion f 1 (O 1 , O 2 ) of the time and drug 2 a proportion f 2 (O 1 , O 2 ) of the time. In the simplest case, f 1 and f 2 are given by the proportion of the stockpile that is of type 1 or 2: In the cycling model, when using drug 1 we set f 1 = 1, f 2 = 0 and vice-versa when using drug 2.
The differential equations that describe the dynamics are: and for the wild-type contact classes and for the drug 1 resistant contact classes and for the drug 2 resistant contact classes and for the multi-drug resistant contact classes and, for x ∈ {p1, p2, np}, y ∈ {w, r 1 , r 2 , r 12 }, and, for symptomatic infectious states and, for asymptomatic infectious states and, for x ∈ {p1, p2, np}, y ∈ {w, r 1 , r 2 , r 12 }, and finally, for the depletion of the finite stockpiles of AV drugs, R 0 is the maximum eigenvalue of a 30 × 30 matrix. As in the single-drug model, each of the four non-zero eigenvalues corresponds to the reproduction number for one of the four strains in circulation:

S1.2.2 Treatment and Prophylaxis model
We have four strains in circulation:wild-type virus (w label), treatment-resistant virus (r t label), prophylaxis-resistant virus (r p label) and multi-resistant virus (r tp label). We introduce new parameters: • φ rt,rp,rtp -the transmissibility of the (treatment, prophylaxis, multiple)-resistant strain relative to the transmissibility of the uncontrolled wild-type strain.
• ρ t -the proportion of offspring (that is, secondary infectives) of treated wild-type infectives who carry the resistant strain.
• ρ p -the proportion of offspring of breakthrough (prophylaxis) wild-type infectives who carry the resistant strain.
We assume that all offspring of r t -resistant infectives are either r t -resistant or r tp -resistant. Similarly, all offspring of r p -resistant infectives are either r p -resistant or r tp -resistant. All offspring of r tp -resistant infectives are r tp -resistant.
The provision of AV drugs as prophylaxis to contacts of treatment-resistant infectives is assumed to have the same protective effect as provision of prophylaxis to contacts of wild-type infectives. Treatment has no effect on infectives carrying the treatment-resistant strain. Similarly, the provision of AV drugs as treatment to prophylaxis-resistant infectives is assumed to have the same effect as provision of treatment to wild-type infectives. Prophylaxis of contacts of prophylaxisresistant infectives has no protective effect.
Provision of AV drugs as prophylaxis to contacts or as treatment to infectives carrying the multiple-resistant strain is assumed to have no effect on either infectiousness or susceptibility.
The force of infections arises from 18 infectors: (S122) We have eight contact states: C y x (x ∈ {p, np}, y ∈ {w, r t , r p , r tp }): .
The differential equations that describe the dynamics are: and for the wild-type contact classes and for the prophylaxis-strain resistant contact classes and for the treatment-strain resistant contact classes and for the multi-drug resistant contact classes and, for x ∈ {p, np}, y ∈ {w, r p , r t , r tp }, and, for symptomatic infectious states and, for asymptomatic infectious states and, for x ∈ {p, np}, y ∈ {w, r p , r t , r tp } dR y Ix dt = γ I I y x (S159) and finally, for the depletion of the finite stockpiles of AV drugs, When the stockpile of AV drugs for prophylaxis, O p , or treatment, O t , runs out that form of intervention is no longer possible. Within the simulations, we set or ψ to zero as appropriate from that point onwards.
R 0 is, again, the maximum of four eigenvalues (one for each of the strains):

S2 Results
We consider a number of different scenarios to demonstrate the general nature of the results presented in the main text. We also provide a sensitivity analysis for the single-drug, random allocation and treatment and prophylaxis models.

S2.1 Alternative scenarios
In this section we examine two plausible alternative scenarios where it would be expected that a successful intervention is possible: 1. The resistant strain has a higher relative fitness than φ = 0.8, but a lower seeding rate. We examine a 90% fit resistant strain (φ = 0.9) with a seeding rate of 10 or 100 times lower (ρ t = 10 −2 or ρ t = 10 −3 ). We keep ρ p = ρ t /10.

The symptomatic proportion is somewhat lower. With a fixed clinical attack rate of 40%
this implies a higher R 0 . We examine a scenario with α = 0.6 and an increased (but still realistic) intervention of providing prophylaxis to 40% of contacts (up from 30% in the main text) and treatment to 50% of symptomatic cases (up from 40% in the main text).

S2.1.1 Extremely high fitness and lower seeding
We keep the symptomatic proportion unchanged at α w np = 0.7, and the intervention unchanged at = 0.3, ψ = 0.4. We restrict ourselves to a 90-10 stockpile. Figures S3-S5 are the equivalent of Figures 5-7 in the main text. While each strategy is less capable of controlling the epidemic compared to the main text, the treatment and prophylaxis strategy remains the most effective choice. In fact, especially in the case of lower seeding, it is more clear cut that treatment and prophylaxis is a better strategy compared to either single-drug usage or a random allocation strategy. The reason is clear: with higher fitness (φ = 0.9) it becomes even more important to delay emergence of the multi-drug resistant strain.  Figure S3: a. ρ t = 0.01 (one order of magnitude less than usual) b. ρ t = 0.001 (two orders of magnitude less than usual).  Figure S4: a. ρ t = 0.01 (one order of magnitude less than usual) b. ρ t = 0.001 (two orders of magnitude less than usual).  Figure S5: a. ρ t = 0.01 (one order of magnitude less than usual) b. ρ t = 0.001 (two orders of magnitude less than usual).

S2.1.2 A reduced symptomatic proportion
Consider a scenario where α = 0.6 (R 0 = 1.x). Without an increase in the intervention the impact of treatment and/or prophylaxis is minimal (see Figure 2b in the main text). Here, we examine the ability of an increased intervention to delay the onset of the epidemic. Because the baseline epidemic duration is shorter, an increased intervention can not produce delays to median infection as great as for α = 0.7. However, a similar multiplicative increase in the time to median infection can be achieved and the different strategies (single-drug, random allocation, treatment and prophylaxis) have the same qualitative impacts. Figure S6 shows the epidemic curves, cumulative curves and cumulative resistant proportion for a 90-10 stockpile scenario. We explore this scenario (α = 0.6, = 0.4, ψ = 0.5) a little more below when performing a sensitivity analysis for multi-drug strategies.

S2.2 Sensitivity Analysis
Throughout the main text we assumed that: 1. antiviral drugs had a fixed efficacy. For treatment, we set e t = 0.7 (relative infectiousness).
For prophylaxis, we set e s = 0.5 (relative susceptibility) and e i = 0.4 (relative infectiousness if nevertheless infected).
Here, we examine sensitivity to these parameter choices. We also compare conclusions regarding optimal strategy across these same parameter sensitivity analyses. Throughout, we restrict ourselves to a high-fitness (φ = 0.8) resistant strain scenario. In the multi-drug strategy section we further restrict the analysis to high seeding only (ρ t = 10 −1 , ρ p = 10 −2 ). We perform this analysis with the symptomatic proportion fixed as in the main text (α w np = 0.7). The intervention also remains as in the main text ( = 0.3, ψ = 0.4).

S2.2.1 Single drug model
With a single drug for treatment and prophylaxis, it is implausible that e t and e i vary independently. Thus, we tie e t to e i : as e i varies from 0.4 to 1, e t varies from 0.7 to 1. Figure S7 shows the singledrug model's sensitivity to variation in antiviral drug efficacy.  Figure S8 shows the single-drug model's sensitivity to variation in the seeding rate of the resistant strain. We always keep ρ p = ρ t /10.  Figure S8: Impact of varying ρ t and ρ p (logscale) for φ above the threshold (solid line) and φ below the threshold (dashed line). a. Above the threshold, increasing the seeding rate increases the proportion resistant. Below the threshold, the proportion resistant is essentially independent of ρ t and ρ p . b. Time to median infection decreases as seeding rate increases above the threshold and increases above the threshold (the "immunising" effect). Figure S9 shows the single-drug model's sensitivity to variaiton in the assumed relative transmissibility of asymptomatic infections.  Figure S9: Impact of varying χ. Colour is seeding rate (ρ t = 0.1, ρ p = 0.01 (blue), ρ t = 0.01, ρ p = 0.001 (red), ρ t = 0.001, ρ p = 0.0001 (green)). Solid lines are high fitness (φ = 0.8), dashed lines are low fitness (φ = 0.3). a. For low fitness strains there is no ability to generate significant proportions of resistance. For high fitness strains, and for high seeding (blue) we have a resistantstrain dominated epidemic. Therefore, relative transmissibility of asymptomatics is not important. For lower seeding (red and green), the epidemic is a mixture of wild-type and resistant-type strains and thus, as χ increases and the wild-type is more able to continue transmitting, the proportion resistant is reduced. b. For transmissible resistant strains (solid lines) we generally have significant proportions resistant so the time to median infection is fairly independent of χ. For untransmissible resistant strains (dotted lines) the time to median infection is reduced as asymptomatics are more capable of transmitting.

S2.2.2 Random allocation model
In a random allocation strategy, both drugs 1 and 2 are used for treatment and prophylaxis. Therefore, for each drug we must tie e i , e t and e s as in the single drug strategy. However, we can vary the efficacy of the two drugs. For a fixed fitness and seeding (φ = 0.8, ρ t = 10 −1 , ρ p = 10 −2 ), we vary e i1 and e i2 independently. e t1 and e s1 are locked to e i1 and e t2 and e s2 are locked to e i2 as in the single drug strategy sensitivity analysis. Figure S10 shows the random allocation model's sensitivity to variation in antiviral drug efficacy for a 90-10 stockpile split. Figure S11 shows the result for a 50-50 stockpile split.
For the 90-10 split, we see, unsurprisingly, little sensitivity to the properties of drug 2. For the 50-50 split the system is symmetric in drug 1 and drug 2. Figures S12 and S13 show the random allocation model's sensitivity to variation in the assumed relative transmissibility of asymptomatic infections with a 90-10 and 50-50 stockpile split respectively. We show two scenarios: as in the main text (α w np = 0.7, = 0.3, ψ = 0.4, blue) and the alternative scenario presented above in Section S2.1.2 (α w np = 0.6, = 0.4, ψ = 0.5, red). For the 90-10 stockpile split, most resistance (that is, more than half) is single-drug resistance. The results are largely insensitive to χ While the proportion resistant at stockpile expiry is similar for the α = 0.6 scenario (compared to the α = 0.7 scenario), the time to median infection is much shorter. This reflects the naturally faster dynamics due to the higher R 0 .
For a 50-50 stockpile, fewer of the infections are resistant in total compared to in the 90-10 case, but of those that are, multi-strain resistance is dominant.  Figure S13: Details as in Figure S12 but for a 50-50 stockpile.

S2.2.3 Treatment and Prophylaxis model
In a treatment and prophylaxis strategy we have one drug used solely for treatment and the other solely for prophylaxis. It follows that we examine sensitivity to antiviral drug efficacy in the (e i , e t )-plane. We tie e s to e i as in the single-drug sensitivity analysis. Figure S14 shows the proportion resistant at stockpile expiry and the time of median infection as a function of e i and e t .  Figure S15 shows the variation in output as a function of relative transmissibility of asymptomatic infections, for both the main text scenario (α = 0.7, blue) and the alternative scenario considered earlier in this document (α = 0.6, red).  Figure S15: Details as in Figure S12 but for the treatment and prophylaxis strategy.

S2.2.4 Comparing models
We have explored the sensitivity of each model to key parameter assignments. Here, as in the main text, we make a comparison between models.
A comparison of Figures S12 and S15 demonstrates that across the range of χ values considered, the treatment and prophylaxis strategy provides better outcomes, in terms of both proportion of infections resistant and time of median infection, than does the random allocation strategy. The single-drug strategy result ( Figure S9, solid blue line) is inferior to both multi-drug strategies.
To make a comparison in terms of sensitivity to antiviral drug effectiveness is more difficult. If there is a large asymmetry in the effectiveness of the two drugs, the treatment and prophylaxis strategy can yield a poorer result compared to the random allocation strategy.
Consider a situation where the 10% stockpile (drug 2) turns out to be ineffective (e t → 1, e i2 → 1, e t2 → 1). For the random allocation strategy, this is of little consequence. The 90% stockpile (drug 1) is used for both treatment and prophylaxis and the result is essentially that given by a single-drug strategy. For the treatment and prophylaxis strategy however, we essentially have a prophylaxis only single-drug strategy. Our previous work (McCaw and McVernon, 2007) has shown that a "synergistic" effect is observed when effective treatment (low e t ) is introduced in a population receiving prophylaxis. Thereby, for efficacious drug 1 (low e i , e i1 and e t1 ) and poor drug 2 (high e t , e i2 and e t2 ) in Figure S16b we find that a random allocation strategy provides a longer time to median infection than does a treatment and prophylaxis strategy.  Figure S16: For a 90-10 stockpile split, we compare the random allocation (blue) and treatment and prophylaxis (red) models. The two surfaces are those already presented in Figures S10 and S14. For the blue surface, e i1 and e t1 vary on the drug 1 axis, and e i2 and e t2 on the drug 2 axis. For the red surface, e i varies on the drug 1 axis and e t varies on the drug 2 axis. The impact on susceptibility of each drug is kept fixed throughout at e s = 0.5 for both strategies. If we also tie e s to e i (for drug 1 in the treatment and prophylaxis model, and for both drugs in the random allocation model) we find a qualitatively similar result. a. Cumulative proportion resistant at stockpile expiry. The random allocation surface lies above the treatment and prophylaxis surface for most parameter values. b. Time of median infection. For the most part, the treatment and prophylaxis strategy provides a longer time to median infection. Only when drug 2 (the 10% drug) is ineffective and drug 1 is effective does the random allocation model provide a longer time to median infection.

A Multi-drug model construction
This technical appendix provides the diagrams and R 0 calculations for the multi-drug models.

A.1 Random allocation model
The possible flows are shown in the following figures. Figure S17: Tree diagram for wild strain exposures (prophylaxis part) (1 − ρ t2 )

A.2 Treatment and Prophylaxis model
The possible flows are shown in the following figures.
We define the force of infection row vectors (1 × 18) (note the transpose): (1 − K rtp ) Figure S27: Tree diagram for multi-resistant strain exposures We construct column vectors F w , F rp , F rt and F rtp as follows: and