The impact of antiretroviral therapy on population-level virulence evolution of HIV-1

In HIV-infected patients, an individual's set point viral load (SPVL) strongly predicts disease progression. Some think that SPVL is evolving, indicating that the virulence of the virus may be changing, but the data are not consistent. In addition, the widespread use of antiretroviral therapy (ART) has the potential to drive virulence evolution. We develop a simple deterministic model designed to answer the following questions: what are the expected patterns of virulence change in the initial decades of an epidemic? Could administration of ART drive changes in virulence evolution and, what is the potential size and direction of this effect? We find that even without ART we would not expect monotonic changes in average virulence. Transient decreases in virulence following the peak of an epidemic are not necessarily indicative of eventual evolution to avirulence. In the short term, we would expect widespread ART to cause limited downward pressure on virulence. In the long term, the direction of the effect is determined by a threshold condition, which we define. We conclude that, given the surpassing benefits of ART to the individual and in reducing onward transmission, virulence evolution considerations need have little bearing on how we treat.

In this section we cover in more detail the calculations using the PDE model that help us to understand when ART is expected to favour the less virulent strain. We then relate the PDE model to the ODE model at equilibrium and solve the full ODE model at equilibrium. This is done in order to calculate 'total number of AIDS-related deaths per untreated infected person per year' at equilibrium and its relation to the relative importance of low CD4 transmissions.

ODE model equations
Recall the equations for the ODE model with ART as presented in (2) in the main text. S represents the population of susceptible (uninfected) individuals, I are the various infected classes, where the superscripts H and L indicate high and low (resp.) SPVL and the subscripts S1 and S2 index strain type, with u and t distinguishing between individuals who will remain untreated or will be treated with ART. To shorten notation we denote I j i = I j i,t + I j i,u and define λ i = (β H cI H i + β L cI L i )/N . T is the population of individuals receiving ART, assumed to be noninfectious.

PDE model definition and solution at epidemic equilibrium
To recap, in the PDE model both treatment and AIDS-related death are dependent on age of infection. Therefore each infected class is no longer split into 'treated' and 'untreated'. Let τ denote time since infection, and define y H 1 (t, τ ) to be the number of high SPVL virulent strain (S1) infections of infection-age τ at time t, so that ∞ 0 y H S1 (t, τ )dτ = I H S1 (t) and similarly for y L S1 , y H S2 , y L S2 . Let an average person with a high SPVL infection become classified as 'low CD4' at τ = T H and an average person with a low SPVL infection become classified as 'low CD4' at τ = T L . Let d H and d L be the AIDS-related death rates of patients with high and low SPVL infections respectively. Then the full set of PDEs describing the spread of the two strains is: if τ > T H y H S1 (t, 0) = λ S1 S ∂y L S1 ∂t + ∂y L S1 ∂τ = −(µ + m)y L S1 (t, τ ) + my L S2 (t, τ ) if τ ≤ T L −(µ + d L + m)y L S1 (t, τ ) + my L S2 (t, τ ) if τ > T L y L S1 (t, 0) = 0 where λ S1 = (β H cI S1 H +β L cI S1 L )/N and λ S2 = (β H cI S2 H +β L cI S2 L )/N . There is also a discontinuity at τ = T i where a proportion k of infected patients reaching that infected age are put onto treatment (and hence removed from the infectious class): At equilibrium, when ∂y j i ∂t = 0, the above reduce to a set of ODEs in τ : where S * is the equilibrium value of S and so λ S1 S * , λ S2 S * are now constants (independent of τ ). We may solve these for y j S1 + y j S2 and y j S1 − y j S2 : Recall that the purpose of introducing age-of-infection dependent death rates and treatment was to determine for which parameter sets the following inequality is true Pr(infected with strain S1 | low CD4) > Pr(infected with strain S1) Here the sample space is the population of currently infected individuals and we are concerned with how these probabilities compare in the absence of ART (k = 0). This can be written as: We may now use the expressions in (3) for y H S1 (τ ) etc., with k = 0, to calculate both sides of (4) at equilibrium. First the left hand side: and similarly to before the integral of all four populations is: Substituting these all into equation 3, cross multiplying and comparing the λ S1 λ S1 , λ S1 λ S2 and λ S2 λ S2 terms on both sides we obtain To simplify further we set m = 0, which is reasonable when m is very small compared to µ. So, with this assumption we are able to obtain

Relating the ODE and PDE models at equilibrium
As covered briefly in the main text, we wish to relate the parameters a j and p j in the ODE model to the parameters T j and d j in the PDE model. We do this by equating the two models at equilibrium, i.e. ensuring that at equilibrium the sizes of both the high SPVL and low SPVL infected classes in the ODE model match those in the PDE model. So, for j = H, L: Substituting in k = 0 and k = 1 shows that in order to obtain a condition that is independent of k we must have µ + d j = p j (µ + a j )/(p j − a j ). Substituting this back in: So we must choose a i and p j to satisfy

Equilibrium population size calculations for the ODE model
To find the equilibrium population class sizes we first re-write the model equations in (4) in matrix notation. Let dy be the vector of derivatives, ordered as in (4). Similarly let y be the vector of population sizes. Then define the matrices Λ and M such that At equilibrium dy = 0 and so we wish to solve ΛyS * /N * = −My, where S * and N * are the equilibrium values of S and N . The solution for y that we want is the eigenvector of M −1 Λ corresponding to the eigenvalue with the largest absolute value. Now let δ 1 = µ + a H + m, δ 2 = µ + p H + m, δ 3 = µ + a L + m, δ 4 = µ + p L + m so that and we can calculate the inverse wherek = (1 − k). We can calculate the eigenvalues and eigenvectors of this matrix computationally. If y * is the eigenvector corresponding to the largest eigenvalue then the equilibrium value of total AIDS-related deaths per untreated infected person per year is , as plotted in figure 5 (main text).
2 Low heritability of SPVL

Methods
In this new model an infection initiated with strain S1 has high SPVL with frequency f ∈ (0.5, 1) (whereas before f = 1) and remains at the same SPVL regardless of within host strain switching (as before). The same is true for infections initiated by strain S2 and low SPVL. The lower f is, the lower the heritability of SPVL in the model. The adjusted model equations are: Note that a low strain switching rate, m, is still required to avoid one strain necessarily becoming extinct.
The R 0 values are now calculated as (approximating using m = 0 again):

infections generated by high SPVL infection) + Pr(new infection has low SPVL) × (avg. # infections generated by low SPVL infection)
Hence: And the equilibrium matrix becomes:

Results: Magnitude and direction of the ART effect
We have not re-done the PDE analysis for this model since we expect the direction of the ARTeffect at equilibrium to remain the same as in the original model. This is because, even though high SPVL infections don't correspond exactly to strain S1 infections, it is still the case that if 'low CD4' transmissions occur proportionally more in high SPVL infections (compared to low SPVL infections) then they are of greater importance to the onward transmission of S1. Therefore, so long as f > 0.5, the direction of the long-term ART-effect will remain the same.
In supplementary figures 1-3 we have redone figures 3, 5 and 6 (resp.) of the main text using this model and f = 0.8. We see that, in both the short-and the long-term, the magnitude of the ART-effect on virulence is lessened. This is because reducing the heritability brings the number of high and low SPVL infections closer together, and so there is less scope for change in overall virulence. Aside from this the results are qualitatively the same; the figures show that the patterns of virulence change in the absence of ART and the direction of the ART-effect in the short-and long-term are all very similar to in the original model. Parameters and initial conditions are as in figure 3 (main text) but with f = 0.8. A, C: In the absence of ART R S1 0 = 2.4, R S2 0 = 2.2. B, D: In the absence of ART R S1 0 = 2.2, R S2 0 = 2.4. Solid lines show the output of the model without ART while dotted lines plot trajectories for the model when ART is introduced at 30% coverage 40 years into the epidemic. A, B: The number of strain S1 infections through time is shown in red and the number of strain S2 infections in blue. With reduced heritability, the numbers of S1 and S2 infections are similar to in the case where f = 1 (figure 3, main text). C, D: However since each strain now gives rise to a mixture of high and low SPVL infections, the numbers of high and low SPVL infections are closer to each other and this results in the smaller magnitude variations in virulence. The shape of virulence change over time is similar to when f = 1. The effect of ART is still to increase virulence compared to the model without ART, but this effect is also reduced in magnitude when heritability is reduced.
Total number of AIDS-related deaths per untreated infected person per year is plotted on the y axis. Since each strain now gives rise to a mixture of high and low virulence infections, the overall virulence values that occur when one strain dominates and the other is maintained at mutation-selection balance are less extreme than in figure 5 (main text) when f = 1. Hence the potential effect of ART is reduced in magnitude while the direction of the effect remains the same.
Total number of AIDS-related deaths per untreated infected person per year is plotted on the y axis. Note that points are not plotted for k = 1 because when all patients are treated there are no AIDS-related deaths and so overall virulence is not well defined at this point. Since each strain now gives rise to a mixture of high and low virulence infections, changes in the strain balance have a less severe effect on overall virulence and so the potential effect of ART is reduced compared to in figure 6 (main text) when f = 1. The direction of the ART-effect is to decrease overall virulence at this point in time for all parameter sets, as before.

Methods
In this model we add in an 'acute' stage of infection, A, which lasts for on average 1/α years. Individuals in the acute classes have γ-fold higher infectiousness. When these individuals transition into the standard infected classes they are at that point subdivided unto those who will or won't later receive treatment. The ODE model equations are shown below.
The R 0 values for this model are calculated as follows (in the absence of ART, and using the approximation m = 0 as before): To examine how changing the transmission profile alters the relative importance of low CD4 transmissions we switch to using the PDE form of the model. The partial differential equations for the infected populations are no different to in the original PDE model (equations 4, main text), except that λ i changes. Let acute infection last for A years. In the simple m = 0 case we are now interested in when the following inequality holds: where, as before (with m = 0) Note that the end of the acute period (t = A) is before the start of the treatment period for both high and low virulence infections (t = T H , T L ). Hence Writing (6) in terms of the parameters of the model: After multiplying out, the λ S1 λ S1 terms all cancel and we can then divide by λ S1 λ S2 S 2 to get: With these parameters, the condition for ART to benefit the low virulence strain is: (where as when we previously assumed constant probability of transmission throughout infection the RHS was (e µT H − 1)/(e µT L − 1)). Since e µT H and e µT L are both a little greater than 1, the fraction on the RHS is larger than previously. Hence, when we account for increased transmission in acute infection, then d L and d H have to be closer in size than they did before in order to see a benefit to the low virulence strain. Or, in other words, there are fewer cases in which the effect of ART is to drive a decrease in virulence (as seen in supplementary figure 5 compared with figure 5 (main text)).  figure 5 (main text) plus α = 2, γ = 13. Recall B H A H is kept constant while B L A L is varied from low (red line) to high (yellow line). A: R S1 0 = 2.0, R S2 0 = 1.8. B: R S1 0 = 1.8, R S2 0 = 2.0. Total number of AIDS-related deaths per untreated infected person per year is plotted on the y axis. Since transmissions late in infection are now of lesser relative importance to both strains, ART very rarely has a large effect on virulence. There are fewer cases in which the effect of ART is to drive a decrease in virulence. is kept constant while B L A L is varied from low (red line) to high (yellow line). Initial conditions are as given in figure 3 (main text). A: R S1 0 = 2.0, R S2 0 = 1.8. B: R S1 0 = 1.8, R S2 0 = 2.0. Total number of AIDS-related deaths per untreated infected person per year is plotted on the y axis. As seen in the equilibrium calculations, the effect of ART on virulence in this model is smaller in magnitude compared to the model with constant transmission profile (figure 6, main text) and there are fewer cases in which the effect of ART is to drive a decrease in virulence.

Modelling increased transmission amongst a subsection of the population 4.1 Methods
For this extension of the model, we split the population into a 'core' and 'periphery', where there will be much higher contact between individuals in the core and the core will be much smaller than the periphery. We assume no movement between the core and periphery compartments. So the basic structure of the model is: where, now specifying the viral strain, i, and we now have three contact rates: c c is the contact rate within the core, c p is the contact rate within the periphery, and c b is the contact rate between the core and the periphery. Hence c c > c p , c b . The spread of the infection within each compartment is described by the equations in (1), with λ i replaced by λ i,c or λ i,p as appropriate. We introduce ART to the core only.
Let c b = c p , so it is only within the core that there are higher contact rates. We wish to use this model to simulate the dynamics of generalised and concentrated epidemics and compare the two. For the generalised epidemic we let c p = 1 2 c c and then for the concentrated epidemic c p = 1 100 c c .

Comparison of generalised and concentrated epidemics
The results in figure 7 show that in the generalised epidemic the dynamics of the core and the periphery are almost identical. This situation is not much different from our original model. Therefore it is not surprising that virulence change over time is very similar to in figure 3 (main text), both in the absence and presence of ART. In contrast, in figure 8 we see that in the concentrated epidemic the dynamics of the periphery lag behind those of the core and the fraction of individuals infected in the periphery is more than 10 times smaller than in the core. Note that in the concentrated model the forces of infection are dominated by the following terms: Hence the spread of infection in the core is influenced very little by the periphery, whereas in the periphery the incidence of infection is almost entirely determined by the prevalence of infection in the core and the force of infection is low. These observations explain the strain dynamics seen in figure 8.
Despite the differences in dynamics between the core and periphery, we see in figure 8 that patterns of virulence change in the first 200 years of the epidemic are very similar to in the generalised epidemic, both in the absence and presence of ART. We wish to gain an understanding of what we would expect the ART-effect in the long term (at equilibrium) to be in the concentrated model. First note that in general we would expect the magnitude of the effect to be slightly greater because we are targeting therapy at the individuals who contribute most to onward transmission (the core). We would expect the direction of the ART effect to be similar to that in the original model (since the relative importance of 'low CD4' transmissions to each strain hasn't changed). so long as the ratio of high:low SPVL infections in the core (where ART is applied) is similar to that in the periphery. We can use the approximated forces of infection above to do a straightforward calculation of the equilibrium position of the model (in the absence of ART): Similar calculations for S2 in the core and S1 and S2 in the periphery yield Combining the core and periphery pairs of equations and rearranging gives (µ + a H + m)I H S1,c − m 2 µ+a H +m I H S1,c Initially there are no infections in the periphery, but in the core I H S1,u (0) = I L S2,u (0) = 1 and all other infected populations are 0 initially. A, C: In the absence of ART, the R 0 s of the core (where c = 1) are set at R S1 0 = 2.4, R S2 0 = 2.2. In other words, we choose β H , β L such that 2.4 = β H /(µ + a H ), 2.2 = β L /(µ + a L ). B, D: In the absence of ART, the R 0 s of the core (c = 1) are set at R S1 0 = 2.2, R S2 0 = 2.4. Solid lines show the output of the model without ART while dotted lines plot trajectories for the model when ART is introduced 40 years into the epidemic. Only the core is treated (at coverage of 30%). A, B: The number of strain S1 infections through time is shown in red and the number of strain S2 infections in blue. In this simulation of a generalised epidemic, though the dynamics are driven by transmission in the core, the dynamics in the periphery closely match those in the core. C, D: Patterns of virulence change over time are qualitatively similar to the original model, though the peak is smaller and the magnitude of the ART-effect is lesser since only the core is treated.  Figure 8: Behaviour of the core-periphery model under the assumptions of a concentrated epidemic, with and without ART. Parameters are as in table 1 but with the differing contact rates c c = 1, c p = c b = 0.01. The core is 20 times smaller than the periphery, so S c (0) = 2000, S p (0) = 40000 and also B c = 50 and B p = 1000. Initially there are no infections in the periphery, but in the core I H S1,u (0) = I L S2,u (0) = 1 and all other infected populations are 0 initially. A, C: In the absence of ART, the R 0 s of the core (c = 1) are set at R S1 0 = 2.4, R S2 0 = 2.2. B, D: In the absence of ART, the R 0 s of the core (c = 1) are set at R S1 0 = 2.2, R S2 0 = 2.4. Solid lines show the output of the model without ART while dotted lines plot trajectories for the model when ART is introduced 40 years into the epidemic. Only the core is treated (at coverage of 30%). A, B: The number of strain S1 infections through time is shown in red and the number of strain S2 infections in blue. Compared to the generalised epidemic the peak of infections is later due to reduced transmission in the periphery, and the number of infections is much lower in the periphery. The force of infection in the periphery is very dependent on the proportion of infected individuals in the core and the shape of the epidemic peak reflects this. C, D: Patterns of virulence change are still very similar to those seen in the original model, except that the low virulence strain (S2) starts to dominate earlier in (D)

Results
Under this treatment protocol, the low virulence strain is always affected most when ART is introduced because a greater proportion of its onward transmissions are effectively removed (compared to the onward transmissions of the high virulence strain). Therefore in the long term (at equilibrium), ART always drives an increase in virulence, whatever the value of the parameter p j . The magnitude of the effect is smaller when p j is large, i.e. when patients are treated particularly early in infection, because in that case both strains are (more equally) heavily affected by ART.
In the short term, the effect of ART is either to benefit the high virulence strain (S1), in which case the magnitude is small because S1 dominates in the absence of ART anyway, or to benefit the low virulence strain. However, the benefit to the low virulence strain occurs in the (unrealistic) cases where ART is causing the eventual extinction of both strains -since the turnover of the S2infected population is lower it decays at a slower rate and so the average virulence of an infection rises in these cases. : The relationship between ART coverage and virulence at equilibrium when patients are treated as soon as possible after infection, regardless of disease progression. A: R S1 0 = 2.4, R S2 0 = 2.2. B: R S1 0 = 2.2, R S2 0 = 2.4. Total number of AIDS-related deaths per untreated infected person per year is plotted on the y axis. Treating independent of disease progression translates to setting p L = p H in the ODE model. Different coloured lines here represent different rates at which individuals in the 'treated' class are treated, from p j = 2 i.e. an average of 6 months of infection before initiating ART (blue) to p j = 1/4 i.e. an average of 4 years of infection before initiating ART (red). The remaining parameters are the same as in Table 1 (main text). Note that when everyone is treated almost immediately, both strains take a big hit to onward transmission, where as when individuals are treated after four years the low virulence strain, S2, takes a much bigger hit, hence overall virulence is higher the lower p j is. are shown, with total number of AIDS-related deaths per untreated infected person per year once again used as a proxy for virulence A: R S1 0 = 2.4, R S2 0 = 2.2. B: R S1 0 = 2.2, R S2 0 = 2.4. As in supplementary figure 9, different coloured lines represent different rates at which individuals in the 'treated' class are treated, from p j = 2 i.e. an average of 6 months of infection before initiating ART (blue) to p j = 1/4 i.e. an average of 4 years of infection before initiating ART (red). The remaining parameters and initial conditions are the same as in figure 3 and table 1 (main text). When p j is low then even at this stage in the epidemic the effect of ART is to benefit the high virulence strain (S1). For higher values of p j and ART coverage, k, ART results in the eventual extinction of both strains. Since the death rate of the S2-infected population is lower than that of the S1-infected population this can lead to a rise in the proportion of S2 compared to S1 -infected individuals in the period leading to extinction, and hence a decrease in virulence by this measure. The points representing parameter sets where this extinction has already occurred at 90 years (I S1 (90) + I S2 (90) < 1) are not plotted.

Additional sensitivity analyses
In this section we conduct sensitivity analyses in which we vary the R 0 values (figures 11-13), the timing of ART initiation relative to the peak of the epidemic (figures 14 & 15) and the initial conditions ( figure 16). These sensitivity analyses are based on the original model as presented in the main text. Table 1 summarises the results of all these sensitivity analyses and the analyses in sections 2-4 in terms of the effect of the various modifications on the key results regarding virulence evolution and the ART-effect. A, C: In the absence of ART, the R 0 s are set at R S1 0 = 1.6, R S2 0 = 1.4. B, D: In the absence of ART, the R 0 s are set at R S1 0 = 1.4, R S2 0 = 1.6. Solid lines show the output of the model without ART while dotted lines plot trajectories for the model when ART is introduced at 30% coverage 55 years into the epidemic. The number of strain S1 infections through time is shown in red and the number of strain S2 infections in blue. A, B: The peak of the epidemic now occurs much later and the progress of the epidemic is in general much slower. In the case where R S1 0 < R S2 0 , the high virulence strain is unable to spread through the population as quickly (B) and its domination in the early stages is significantly reduced compared to in simulations with higher R 0 s. C, D: Despite this, patterns of virulence change are still quantitatively similar -a first-phase rise in virulence is followed by a second phase decline (though this is modest in (C)), and the third phase is determined by the balance of the R 0 s. There is a decrease in virulence in the second half-century of the epidemic in all cases, this is most severe when R S2 0 is much greater than R S1 0 and very slight when R S1 0 is greater than R S2 0 . ART always results in a decrease in virulence (compared to the model with no ART) and this ART effect is largest when R S2 0 is a little larger than R S1 0 . Magnitude smaller. ART more often benefits high virulence strain.
Large effects rarer. Fewer cases of decrease in virulence in particular.
Low heritability (4-6) Clearer decrease in virulence after peak. Smaller magnitude of changes in virulence in general.
Magnitude similar in range for this range of parameters. Direction: Sometimes an increase in virulence when ART initiation is later.
Larger magnitude possible at low coverage. Direction: Always an increase in virulence.
High R 0 s (11) Similar patterns seen but twice as fast.
Magnitude slightly larger. Direction still negative.
Magnitude and direction largely unchanged (still determined by equation 8, main text).
Low R 0 s (12) Similar patterns but epidemic progresses much slower.
Magnitude and direction seem similar.
Magnitude and direction unchanged.
Negligible change in magnitude and direction.
Negligible change in magnitude and direction.
Varying initial conditions (16) First phase (before peak) varies depending on initial conditions. Second phase decrease sometimes shorter in duration.
Magnitude can be smaller. Direction unchanged.