The Timing and Targeting of Treatment in Influenza Pandemics Influences the Emergence of Resistance in Structured Populations

Antiviral resistance in influenza is rampant and has the possibility of causing major morbidity and mortality. Previous models have identified treatment regimes to minimize total infections and keep resistance low. However, the bulk of these studies have ignored stochasticity and heterogeneous contact structures. Here we develop a network model of influenza transmission with treatment and resistance, and present both standard mean-field approximations as well as simulated dynamics. We find differences in the final epidemic sizes for identical transmission parameters (bistability) leading to different optimal treatment timing depending on the number initially infected. We also find, contrary to previous results, that treatment targeted by number of contacts per individual (node degree) gives rise to more resistance at lower levels of treatment than non-targeted treatment. Finally we highlight important differences between the two methods of analysis (mean-field versus stochastic simulations), and show where traditional mean-field approximations fail. Our results have important implications not only for the timing and distribution of influenza chemotherapy, but also for mathematical epidemiological modeling in general. Antiviral resistance in influenza may carry large consequences for pandemic mitigation efforts, and models ignoring contact heterogeneity and stochasticity may provide misleading policy recommendations.

S = −(β u I u + β t I t )S − β r I r Ṡ I u = (β u I u + β t I t )(1 − ρ)S − γ u I u I t = (β u I u + β t I t )ρ(1 − c)S − γ t I ṫ I r = (β u I u + β t I t )ρc S + β r I r S − γ r I ṙ R = γ(I u + I t + I r ) (1)

Mean Field Network Model Equations
To include individual heterogeneity we employ a network model of disease transmission. Here, in contrast to the ODEs in (1), one typically needs to introduce a higher-order compartmentalization where nodes are distinguished not only by their state, but also by their degree. Hence, instead of one equation for the fraction of susceptible individuals S(t) at time t, we write an infinite number of equations for the fraction of susceptible nodes of degree k, S k (t), at time t. We obtaiṅ S k (t) = −k (β u I u + β t I t + β r I r ) S k (t) (2) where I x is the probability, that a randomly chosen link among the k belonging to a susceptible node leads to an infectious individual of type x ∈ {untreated, treated, resistant}. Similarly, Our ability to evaluate this mean-field quantity directly depends on our ability to correctly follow the links stemming from susceptible nodes. To this end, we write a set of five ODEs for the density of the possible links (where a link between nodes of state X and Y is denoted as [XY ]):  Figure 1: Degree distributions used in the main text. The heavy-tail distribution (black line) corresponds to an initial binomial regime leading into a power-law tail with exponential cutoff (for mean degree k ∼ 4.1). The binomial network is given by the red line and has the same mean degree.
[SI r ] = − (β u I u +β t I t +β r I r ) k s +β r +γ r [SI r ] + 2 (β u I u ρc +β t I t ρc + β r I r ) k s [SS] (10) where k s is the average excess degree of susceptible nodes. Equations (2)-(11) are minimal to describe the system in the sense that they are sufficient to calculate all the mean-field quantities on which they depend. Through a simple averaging procedure, one obtains: We evaluate (2)-(11) numerically to explore the transmission dynamics of this system.

Derivation of the Critical Manifold
If β r < β eff = (1 − ρ)β u + ρ(1 − c)β t , and treatment effects only transmissibility, we assume that it is always better to treat than not to treat. We also assume that we are not under the epidemic threshold of either the resistant or untreated strains, as the disease would then die out even without treatment. In opposition, the case where β r > β eff , and both the resistant and untreated wild strains are above the threshold, can lead to an important dilemma: is treatment more likely to be effective at reducing total incidence or at creating resistance? This is problematic because the resistant strain has the opportunity to lead to a bigger epidemic than the untreated disease through treatment failure, so that treatment will actually have worsened the situation.
Here, we calculate the critical initial conditions of treatment I(0) (i.e. the number of individuals in the I u state at the beginning of treatment) for which treatment has approximately a 50/50 chance of being efficient or dangerous.
Above the epidemic threshold of the resistant strain, one individual infected with the resistant strain could spark an epidemic. Thus we simply calculate the probability of getting at least one mutated strain given that we begin treatment when I(0) random individuals are currently infectious. We assume that we are below the epidemic threshold of the treated disease, as an epidemic of the resistant strain is always the more likely outcome.
Below this threshold, each of the I(0) infectious creates a microscopic epidemic of size 1+ T k , where k and k are the average degree and excess degree of the network, respectively, and T is the total probability of transmission. Removing the original untreated infectious individuals leads to I(0) T k 1−T k new infections. From these, the probability of getting at least one mutation is directly obtained by The critical initial conditions of treatment I c (0) is that I(0), for which expected epidemic sizes are equal with and without treatment. The critical probability of mutation P c := P muta (I c (0)) must thus satisfy where R u , R t and R r are the expected epidemic size for the untreated, treated without mutation and treated with mutation cases, respectively. Solving for I c (0) gives Finally, note that the total probability of transmission T is easily approximated from the effective transmission rate of the treated disease, i.e. β eff = (1 − ρ)β u + ρ(1 − c)β t ,

Details of Network Simulations
To perform MC simulations of the model, we have generated networks of size N with the degree distribution {p k } in Figure 3 via the following numerical algorithm: i. draw a sample {k i } of size N from distribution {p k } under the condition that N i=1 k i is even; ii. for each i, produce k i stubs tagged as i; iii. randomly link all stubs in pairs (i, j), thus linking nodes i and j.
This is the so-called Configuration Model with allowed loops and self-links [3]. Each and every network generated by this procedure is accepted and kept in the results, as they are part of the canonical ensemble considered by the mean-field approach of the formalism. For every generated network, a fraction I(0) of individuals are randomly chosen to be initially infectious and the dynamics are then simulated in a discrete time propagation simulation valid for a time step ∆t → 0 (we choose ∆t such that β x ∆t and γ x ∆t are less than 10 −3 for all x): i. at each ∆t, every susceptible neighbor S of every infectious individual I x is infected with probability β x ∆t; ii. wild strain infections are treated with probability ρ leading to mutation with probability c; iii. at each ∆t every infectious individual I x recovers with probability γ x ∆t.

Effects of Network Structure
The main results of the paper are robust to changes in network structure. Supplemental  Figure 1). Treatment when the system is above the critical manifold results in widespread resistance, while treatment below the critical manifold reduces epidemic size. When treatment is initiated early, resistance is minimized (Supplemental Figure 3). g.

Late Treatment
Supplemental Figure 2: Effects of Treatment Timing on Binomial Network. Effects of treatment when initial conditions are above (panels b, c, d) and below (panels e, f, g) the critical manifold. Panel a shows the critical manifold (dashed grey line) and the total number infected (red line) from the mean-field approximations, with each large dot corresponding to the panels at right. Solid lines correspond to mean-field approximations, and points correspond to means of 100,000 simulations on networks of size 250,000. Horizontal black line corresponds to a mean of 1 infected individual in a network of 250,000 over 100,000 simulations. Above the critical manifold, treatment induces large epidemics of the resistant strain (panels c and d). In the effective treatment regime (below the critical manifold) early treatment reduced the number infected by 35-fold (panel f ) while late treatment reduces cases similarly but with much higher resistance (panel g). Parameters: β u = 4 · 10 −4 , β r /β u = 1.2, β t /β u = 0.3, γ u = γ t = γ r = 10 −3 , ρ = 0.6, and c = 1/500.

Proportion Treated
Supplemental Figure 3: Comparison of Random and Targeted Treatment on the Binomial Network. Panel a shows the final size for wild-type, resistant and both infections as a function of percentage treated, ρ, for targeted (dashed lines) and random (solid lines) treatment regimes. We see a transition from wild-type to resistant infections at a lower treatment percentage in the targeted treatment regime. Panel b shows the percent of total infection that is the resistant strain for the targeted (dashed line) and random (solid line) treatment. Parameters: β u = 6 · 10 −4 , β t = 1.8 · 10 −4 , β r = 3 · 10 −4 , γ u = γ t = γ r = 10 −3 , and c = 1/500.