Hedging against Antiviral Resistance during the Next Influenza Pandemic Using Small Stockpiles of an Alternative Chemotherapy

Mathematically simulating an influenza pandemic, Joseph Wu and colleagues predict that using a secondary antiviral drug early in local epidemics would reduce global emergence of resistance to the primary stockpiled drug.


Discrete-event stochastic simulation for a single population
We consider three drug regimens: otherwise.
where ε is the reduction in infectiousness provided by antiviral treatment (see Figure B). β is calibrated to yield the desired 0 R using the formula where Presymptomatic D , Asymptomatic D , and Symptomatic D are the durations of the presymptomatic, asymptomatic, and symptomatic stages .
To illustrate how we simulate emergence of resistance (see Figure B for the resistance emergence model), suppose the drug regimen is Mono throughout the pandemic. Upon onset of symptoms, an infected individual is treated with drug A with probability T p . If the individual is treated, antiviral resistance emerges during his symptomatic stage with probability A p . We assume that the rate of emergence within the individual A r is constant throughout the symptomatic stage. As such, the time to emergence within an individual A T is an exponential random variable with mean 1 A r . We compute A r by solving the equation where both sides of the equation represent the probability that resistance does not emerge. When an individual becomes symptomatic and is treated with drug A, Symptomatic D and A T are drawn from their distributions and compared.
If the latter is smaller, then resistance emerges A T time units after treatment is initiated for the individual. Otherwise, the individual does not develop resistance. This model implicitly assumes that the hazard of emergence of resistance is constant over the infectious course of an individual once treatment begins, but results are weakly sensitive to this assumption since the overall probability of emergence A p is fixed by assumption; hence the assumption of constant hazard affects only when, not whether, emergence of resistance occurs in a particular host. Although later emergence of resistance in individuals should reduce the probability of transmission, it is unlikely that our conclusions are sensitive to our more simple model structure: a slight reduction in transmission from the initial resistant case due to the timing of its appearance must be dominated by the range of values over several orders of magnitude we use for de novo rates of emergence ( A p and B p ) in our sensitivity analysis.
Emergence of dual resistance within a single individual (i.e. ) under combination chemotherapy is simulated using the same mechanism.

Discrete-time stochastic meta-population simulation for global network
The discrete-event simulation above is an exact algorithm for simulating the stochastic epidemic dynamics in a single closed population. However, the computational requirement of this method becomes prohibitive when applied to the many connected populations. We use a discrete-time simulation as an alternative, which is computationally more efficient at this scale. A time-step of 0.25 days is chosen for the discrete-time simulation. With this time-step, the results given by the discrete-time simulation and its discrete-event counterpart are statistically indistinguishable when simulating the pandemic in a single population.

The algorithm
To simulate the global spread of resistance, we first specify the initial conditions at time 0 t = . For example, if Hong Kong (with population index H ) is the first city infected (the "source") with y initial infections with wild-type virus, then is generated from the latent duration distribution specified in the natural history model (see Figure A).
After the initial conditions have been specified, Steps 1 to 4 below are then executed for each time k t k t = Δ , 1, 2, k = … until 365 k t = days: Step 1: System states from the previous time step The state variables at time k t are initiated with the state variables from the previous time step: 1 ( , , , ) ( , , , Step 2: Inter-population travel Data from the International Air Transport Association for 2002 is used to construct an air travel volume matrix , where ij w represents the rate at which a resident from population i travels to population j . An individual who travels to a foreign population returns to his home population after R T days. We assume that R T is exponentially distributed with mean 7 days. We illustrate below how inter-population travel is implemented for the susceptibles. The mechanism is exactly the same for all individuals. Then the ( ) ij k S t 's are updated as follows: Step 3: Emergence of resistance Let ( ) α ⋅ be the function that maps emergence probabilities to emergence rates (e.g., ( ) ( ) are updated accordingly after these emergence events are generated.
Step 4: Infection and disease progression Infection. The instantaneous force of infection for strain s in population i at ⎦ be the number of strain s infections among residents of population i who are in population j at time k t . These infection events are generated with the averaged force of infection ( ) with status ( ,, , ,,) Presymptomatic s None t i j Δ progress to the next disease stage. Upon exiting the presymptomatic stage, these infected individual can become asymptomatic, symptomatic and untreated, or symptomatic and treated. Let Then the following updates are performed:

Reduction in infectiousness under combination chemotherapy
In the main text, we assumed that if a symptomatic individual infected with strain s was treated with antiviral regimen d , his infectiousness was reduced by a proportion 0.67 ε = if strain s was sensitive to at least one drug in regimen d [1]. However, combination chemotherapy is likely to be more efficacious than monotherapy in reducing infectiousness [2]. If we set 1 ε = when a wild-type infection was treated with combination chemotherapy, results (not shown) were indistinguishable from those presented in the main text. Thus, our conclusions are not sensitive to higher values of ε under combination chemotherapy.

Household-based antiviral prophylaxis
In the main text, we assumed that large-scale antiviral intervention involved only treatment of symptomatic individuals. Recent studies on mitigation of influenza pandemics suggested that adding targeted antiviral prophylaxis could further reduce the attack rate [3,4] and many countries are now considering this option (http://www.hhs.gov/pandemicflu/plan/sup7.html). Here, we show that our conclusions in the main text hold when household-based antiviral prophylaxis is implemented in addition to treatment of symptomatic individuals. For illustration, we only present the results for ECC. Results for SMC are similar to that for ECC with high synergy.
To model household-based antiviral prophylaxis, we used an individual-based network simulation that we have previously built to study mitigation of influenza pandemics [4]. See Ref 4 for the details on model specifications. To make direct comparison with the results presented in the main text, we made the following modifications to the simulation: 1. The multi-strain structure and the process of resistance emergence were modeled as described in Discrete-event stochastic simulation for a single population above. Rates of resistance emergence during the presymptomatic and asymptomatic stages were the same as that during the symptomatic stage. Rates of resistance emergence were the same for both treatment and prophylaxis (a worst-case assumption). 2. The effect of hospitalization was ignored. 3. Household-based quarantine was not modeled. 4. Household-based antiviral prophylaxis was assumed to take place as follows. Upon onset of symptoms, an individual who was not already in the intervention program reported with probability 0.4 T p = and received antiviral treatment. In addition, all household members of this individual were immediately recruited into the program and received antivirals: non-symptomatic members received prophylactic regimen for 5 days while symptomatic members (who did not report symptoms earlier) received treatment regimen. The "prophylactic clock" of 5 days was renewed whenever a previously non-symptomatic household member became symptomatic. 5. Antiviral prophylaxis reduced susceptibility by 85% [1]. 6. A 1% stockpile coverage of drug B was available for ECC (as in the main text).
Figure D-A shows that this network model (with household-and workplacestructures) and the simpler homogeneous mixing model in the main text gave essentially the same predictions on AR and RAR when antiviral intervention involved only treatment of symptomatic individuals. Therefore, household-and workplace-structures had little effect on the emergence and spread of resistance. When household-based antiviral prophylaxis was added, AR actually increased unless A p < 10 -3 . This is because the addition of prophylaxis required a higher level of antiviral use (approximately tripling the stockpile requirement) which speeded up the emergence of antiviral resistance at the population level.  Figure 1D. This was because when prophylaxis was implemented, antivirals were depleted faster. As a result, the wild-type prevalence was smaller (compared to the treatment-only scenario) when the 1 % stockpile coverage of drug B was used up, and therefore the subsequent emergence and spread of resistance was more severe. To restore the same level of effectiveness of ECC (as in Figure 1D), the stockpile of drug B needed to be increased by the same proportion of increase in drug A required by prophylaxis (compared to treatment only). Since the stockpile of drug A needed to be approximately tripled when household-based prophylaxis was added, the reductions in AR and RAR were similar to that presented in the main text when the stockpile coverage of drug B was boosted to 3% ( Figure  D-C).
These results suggested that our conclusion in the main text remains valid when household-based antiviral prophylaxis is implemented in addition to treatment.

The number of populations that implement antiviral intervention and SMC/ECC
In the main text, we assumed that 28 out of 105 cities in the global network implement antiviral intervention with Hong Kong as the source. To show that the conclusions made in the main text are robust against this assumption, we varied the number of populations that implement antiviral intervention, M, from 8 to 64. In Figure G, we simulated the scenarios in Figure 3A-D for each value of M. For scenario C and D (as in Figure 3), we also varied the proportion of the M populations that implement SMC (which we denoted by q ) from 0.25 to 0.75. The policies adopted by Hong Kong, London, New York and Geneva in each scenario are the same as in Figure 3. In each realization, the source population was chosen randomly from the network (subject to the exception that (i) the source cannot be London in scenario c because London implements monotherapy, and (ii) the source cannot be Hong Kong, New York and Geneva in scenario d because these populations implement SMC) Figure G shows that the conclusions drawn from Figure 3 are applicable in all scenarios regardless of the values of M and q .

When will SMC/ECC fail in downstream populations if the source does implement SMC/ECC?
In the main text, we claimed that if the source population did implement SMC, SMC became much less effective in downstream populations only if most of the neighbors (i.e. those populations that were immediately downstream) of the source failed to control resistance, which was likely only if most of these neighbors implemented only monotherapy and A p was high. Here, we provide results to support this claim. We considered the same scenarios as in Figure  3A-C but with Kinshasa as the source of infection in the network with 30 wildtype seeds on day 0. Kinshasa was chosen as the source because it has the fewest neighbors in our network model (six neighbors, which are Addis Ababa, Brussels, Johannesburg, Lagos, Nairobi, and Paris). In each scenario, we increased the number of populations that implement monotherapy by n , where 0, 2, 4,6 n = . We conducted two sets of simulations: (a) the n additional monotherapy populations were chosen randomly from the network; (b) the n additional monotherapy populations were chosen randomly from the neighbors of Kinshasa. Figure H shows that for a given A p , the spread of resistance was essentially insensitive to the value of n in Set (a) (which is consistent with the results in Figure G) but increased considerably with n in Set (b). In Set (b), the effectiveness of SMC diminished significantly only when n was large (relatively to the total number of neighbors) and A p was larger than 0.05. Similar results were obtained for a different source population. Taken together, these results give support to our claim that SMC becomes much less effective in downstream populations only if most of the neighbors of the source fail to control resistance, which is unlikely unless most these neighbors implement only monotherapy and A p is very high.

Zanamivir is an effective candidate for the secondary antiviral though it does not reduce infectiousness and is not licensed for use in small children
In the main text, we claimed that although zanamivir is not licensed for treatment in children less than 7 years of age and may not be able to reduce infectiousness, these drawbacks have little impact on the effectiveness of SMC and ECC. We justify this claim in this section.
Empirical evidence suggested that zanamivir provides no reduction in infectiousness [5]. However, our results have shown that reduction in infectiousness by the secondary antiviral is not required for SMC or ECC to be an effective hedge: if the secondary antiviral provides no reduction in infectiousness, then SMC is the same as delaying large-scale antiviral intervention, which we have shown to be an effective hedge ( Figure 1C). For ECC, we have assumed that the reduction in infectiousness under combination chemotherapy is the same as that provided by the primary antiviral alone ( 66% ε = ), which would be the case if the secondary antiviral provides no reduction in infectiousness and does not interfere with the action of the primary antiviral. Therefore, zanamivir should not be ruled out as a suitable secondary antiviral on the basis of a potential lack of reduction in infectiousness.
Another potential problem with using zanamivir as the secondary antiviral is that zanamivir is not licensed for treatment in children less than 7 years of age. However, the need for treating this special group with the primary antiviral (oseltamivir) has little impact on the overall effectiveness of SMC and ECC. For instance, this group accounts for ~10% of the population in the US and treating this group with the primary drug under SMC is similar to reducing A p by 10 times at the population level during the early phase of large-scale antiviral intervention. Further, SMC and ECC will still reduce the spread of Aresistance during the early phase by lowering the infectiousness of treated individuals who do not belong to the special group and are infected with an Aresistant strain that is not dually resistant. These effects together ensure that the overall effectiveness of SMC and ECC will only be slightly attenuated by the need for treating the special group with the primary antiviral.

Resistance emergence due to reassortment
Recent establishment of the oseltamivir-resistant H1N1 strain [6,7,8,9,10,11,12] has elevated the concern that an oseltamivir-sensitive pandemic strain may gain resistance via reassortment with circulating oseltamivir-resistant strains. We have shown in Figure 1C that the spread of resistance is limited if the cumulative number of infections at the time resistance emerges (in a population of 6.8 million) is 100,000 or above. Therefore, as long as the probability of resistance emergence via any means other than drug pressure is less than 0.00001 for every infected case, our conclusions in this study will hold. To put this threshold in the context of reassortment with seasonal strains, consider a region where the epidemic curve of seasonal influenza is relatively flat with an annual attack rate of 15% (e.g. Singapore). For resistance to emerge via reassortment between the pandemic strain and a circulating seasonal strain, an individual must be coinfected by both strains. Assume that the pandemic has no effect on the transmission of seasonal influenza, which is a conservative assumption favoring reassortment. In this scenario, if an individual is infected with the pandemic strain, the probability that he will be coinfected by a seasonal strain during his pandemic infection is on the order of 0.001 (0.15/365*duration of pandemic infection in days). If we multiply this by the probability of the occurrence of a reassortment event that gives rise to a resistant and fully fit strain in a coinfected individual (say, 1 in 100), the contribution of reassortment to resistance emergence will likely satisfy the 0.00001-threshold and is small compared to drug-induced resistance in our base case.

Deterministic Version of the Stochastic Model
The modeling equations for the deterministic model are shown below:  All parameters and symbols in the above system have been defined earlier in the section "Algorithms" (except that the mean duration for disease stage i is denoted here by i D instead of E[D i ]). Table A. Assumptions for key unknown transmission parameters. Ranges presented here were used for sensitivity analyses. Parameters were translated into symptomatic and asymptomatic durations and relative infectiousness for different disease stages (see Figure A). Broad sensitivity analysis. Lower bound must be greater than proportion of cases asymptomatic.   Table 1 shows that the deterministic model is not necessarily an accurate estimate of the mean behavior of the stochastic model.   ECC outperformed SMC (i.e. AR(ECC)−AR(SMC) < 0%) in 53% of the scenarios in part C, which was much higher than the 22% in Figure 2C. The reason is that most of the values of s generated were very high (72% and 86% of the s values generated were greater than 0.99 and 0.9. respectively). See Figure F for a more detailed examination of the effect of synergy on the relative performance of ECC and SMC. D-F This column is the same as Figure 2 except 1 T p = . ECC outperformed SMC (i.e. AR(ECC)−AR(SMC) < 0%) in 59% of the scenarios, which was much higher than the 22% in Figure  2C. The reason is that higher treatment probability tended to increase the proportion of scenarios where containment was likely under ECC but not under SMC."  Figure  E). Panels in the left column are analogous to Figure 2C, and parts C and F of Figure E. Panels in the right column show the sensitivity of AR(ECC)-AR(SMC) to the synergy parameter (note that on the x-axis, 1 s − is plotted on log-scale in d while s is plotted on linear scale for the others). By comparing the proportion of scenarios for which ECC outperformed SMC in each row, we conclude that SMC is likely to outperform ECC unless synergy is very high (> 0.95).

Figure G. Sensitivity of the effectiveness of SMC against the number of populations that implemented antiviral intervention (M) and monotherapy (qM).
Vertical bars indicate 95% prediction intervals with triangles indicating the upper and lower bounds and horizontal dashes indicating the means (from 500 realizations). In each panel, the four colors correspond to the scenarios in Figure 3 (red -scenario A, green -scenario B, blue -scenario C, cyan -scenario D). Columns A-C correspond to different values of M. Each row corresponds to a different population (the same as in Figure 3).  ECC0 and ECC1 refer to ECC with s = 0 and s = 1, respectively. A 1% drug B coverage, as in the main text. B 5% drug B coverage. C 20% drug B coverage, which was sufficient to give combination chemotherapy to all treated cases throughout the entire epidemic.  Figure 1C, we argued that in a population of 6.8 million, the spread of resistance can be substantially reduced if the emergence of resistance can be delayed until the cumulative (wild-type) attack rate had reached 1.5% or above. This figure shows that the same conclusion applied for population sizes up to 20 million as long as the probability of emergence of resistance was within reasonable range (<0.05). Therefore, if 1% drug B coverage was sufficient for ECC and SMC to be effective in a population of 6.8 million (as in the main text), the same level of coverage was sufficient for all population sizes of interest in the global model. Consequently, in the global model, we assumed 1% drug B coverage for all populations that implemented SMC or ECC.     Figure 1 in the main text. A This chart is the same as Figure 1A with the outcomes in the deterministic model shown by the red curves. B The singlepopulation epidemic in the deterministic model with 0.01 A p = (brown shaded area corresponds to resistant incidence; blue shaded area corresponds to wild-type incidence; AR, attack rate; RAR, resistance attack rate). C AR and RAR as functions of the cumulative number of wild-type infections at the time resistance emerged (colors correspond to the value of A p as per diamonds in part a; in the absence of interventions AR=73% and RAR=0%; in the absence of resistance, AR=56% and RAR=0%). D Efficacy of early combination chemotherapy (ECC, antiviral synergies of 0 s = and 1 s = ) and sequential multi-drug chemotherapy (SMC) in reducing the attack rate.