Superinfection and cure of infected cells as mechanisms for hepatitis C virus adaptation and persistence.

RNA viruses exist as a genetically diverse quasispecies with extraordinary ability to adapt to abrupt changes in the host environment. However, the molecular mechanisms that contribute to their rapid adaptation and persistence in vivo are not well studied. Here, we probe hepatitis C virus (HCV) persistence by analyzing clinical samples taken from subjects who were treated with a second-generation HCV protease inhibitor. Frequent longitudinal viral load determinations and large-scale single-genome sequence analyses revealed rapid antiviral resistance development, and surprisingly, dynamic turnover of dominant drug-resistant mutant populations long after treatment cessation. We fitted mathematical models to both the viral load and the viral sequencing data, and the results provided strong support for the critical roles that superinfection and cure of infected cells play in facilitating the rapid turnover and persistence of viral populations. More broadly, our results highlight the importance of considering viral dynamics and competition at the intracellular level in understanding rapid viral adaptation. Thus, we propose a theoretical framework integrating viral and molecular mechanisms to explain rapid viral evolution, resistance, and persistence despite antiviral treatment and host immune responses.

Many RNA viruses, including global pathogens of major medical importance such as hepatitis C virus (HCV), HIV, and ebolavirus, exist within infected hosts as large populations of genetically related viral variants commonly referred to as a quasispecies (1). The diverse nature of the viral quasispecies allows viral populations to evolve rapidly to adapt to abrupt changes in the host environment (2); examples include escape from adaptive immune pressures (3, 4) and development of resistance to antivirals (5, 6). Despite the importance of viral adaptation, our understanding of how the viral quasispecies responds to selection pressure and the underlying molecular mechanisms supporting viral adaptation in vivo is limited by a lack of frequent longitudinal viral sequence data generated by methods that retain linkage across genes and genomes (7).

Previously, much attention has been devoted to estimating the effectiveness of antivirals against viral variants in the quasispecies (8–10); however, the mechanisms that drive an adaptive/fitter variant (e.g., a drug-resistant mutant) at an initially low frequency to a high frequency are less well studied. Theoretical work suggests that fitter variants at a low frequency cannot expand unless there are sufficient numbers of target cells, that is, “replication space” (5, 11–14). This is analogous to the concept of competitive release in ecological “niche” theory, that is, a species is not able to expand unless there is a niche vacated by other species (15). Previous models have in general assumed that, once a variant infects a cell, the cell is occupied by the variant for its remaining lifetime, and thus the replication space is provided by generation of new target cells (5, 8, 13, 16–18). Under this framework, viral variants compete for infecting newly generated target cells, and thus, the expansion of a low-frequency fitter mutant is determined by how quickly infected cells die and are replaced by new target cells. This framework has been successfully applied to understand HIV drug resistance (6, 13, 16, 18). However, the source of replication space for other rapidly evolving viruses, and in particular HCV, is unclear.

HCV is a positive-strand RNA virus belonging to the flaviviridae, a family of medically important single-strand RNA viruses. HCV infects ∼80–170 million people worldwide (19). Chronic infection can lead to cirrhosis and hepatocellular carcinoma. Current treatments for HCV infection involve combinations of direct acting antivirals (DAAs), and they have achieved remarkably high cure rates (20–22). Ironically, because of this clinical effectiveness, combination treatments cannot generally be used to probe the molecular mechanisms that contribute to virus persistence. However, this information may be critical for designing effective vaccination strategies for HCV and therapeutic regimens and vaccines for other members of the flaviviridae family. Here, we analyzed clinical samples collected from a phase I/IIa monotherapy trial (23) of a second-generation HCV protease inhibitor, MK-5172 (also called grazoprevir) (22, 24–26), that is extremely potent but susceptible to drug resistance development through multiple viral genetic pathways of escape. Measuring viral loads and analyzing viral populations longitudinally by single-genome sequencing (7) allowed us to follow the population dynamics and evolutionary responses of HCV to drug pressure and its removal in vivo at high resolution. We analyzed this unique dataset by combining mathematical modeling and phylogenetic analysis [termed “phylodynamic” analysis (27)]. We provide strong evidence that cure and superinfection of infected cells contribute essentially to the replication space, allowing resistance mutant expansion and evolution.

Results

Study Subjects and Single-Genome Sequence Analysis.

We studied eight participants from a dose-ranging monotherapy trial of MK-5172 (23). Five were treated with MK-5172 (50 or 800 mg) once daily for 7 d (subjects 1–5, genotypes 1a, n = 4; 1b, n = 1) and three were treated with placebo (subjects 6–8, genotypes 1a, n = 2; 1b, n = 1) (). All of the participants were followed for 7–8 wk after treatment cessation. We measured viral loads longitudinally and analyzed viral populations using single-genome sequencing (7).

Samples taken immediately before treatment initiation, designated as day 0, were analyzed in all eight subjects for genotypic complexity and baseline drug resistance mutations (DRMs). At each sampling timepoint, ∼100 viral genomes were sequenced in each subject (Fig. 1 and ). Day 0 sequences of each subject showed patterns of virus diversity typical of chronic infection with broad genotypic heterogeneity (). Maximum pairwise diversity and mean pairwise diversity ranged from 1.17 to 5.64% and 0.53 to 3.45%, respectively. No drug resistance-associated mutations were identified in the pretreatment sequences, although a Q80K substitution was found in all sequences in subjects 2 and 3. Q80K is a common polymorphism that confers no resistance to MK-5172.

Fig. 1.

Sequential plasma virus load and sequences from subject 1. (A) Time course of treatment with MK-5172 (shaded area, days 1–7), viral load determinations (blue solid dots), and viral sequence analyses (open circles at days 0, 14, 27, and 56). (B) A maximum-likelihood (ML) phylogenetic tree of all viral sequences sampled from subject 1 from all time points. Tree tips are color coded according to the known resistant mutations they bear. (C) ML phylogenetic trees of viral sequences sampled from subject 1 at each time point.

In the five subjects treated with MK-5172, viral load decreased over 5 logs and went below (or close to) the limit of quantification (=25 IU/mL) during the 7-d treatment period (Fig. 1 and ). At the end of this period, viral load was too low to permit sequencing analyses. By day 27, 20 d after treatment was stopped, the viral load rebounded to levels (106 to 107 IU/mL) similar to the baseline level before treatment in all treated subjects and could be quantified as early as day 14 (7 d after stopping treatment) in subjects 1 and 4. Follow-up single-genome sequencing analyses were conducted on samples from days 14 to 62.

After therapy was stopped, the viral populations that emerged carried drug resistance mutations. Surprisingly, the dynamics of mutant virus populations exhibited stringent population bottlenecks and rapid turnover of the dominant resistant mutants during the follow-up period, even in the presence of relatively stable viral load (e.g., subject 1 between days 27 and 56; Fig. 1). Extremely rapid and continuous turnover of the resistant mutants, after therapy cessation, was observed in subjects 1 and 4 (Fig. 1 and ). In contrast to these dramatic changes in plasma viral RNA (vRNA) load in subjects treated with MK-5172, the three subjects who received placebo had no changes in plasma vRNA levels and no sequences containing DRMs. Viral sequences from these individuals were obtained on days 0 and 27 or 34 (). Below, we describe the evolutionary patterns in subject 1 in detail. The patterns in other treated subjects are described in .

In subject 1, none of 103 single genome sequences on day 0 contained drug resistance-associated mutations (Fig. 1). Conversely, on day 14 (i.e., 7 d after therapy was stopped), all 111 single-genome sequences analyzed contained one or more drug resistance-associated mutations. Fifty-six of the 111 sequences contained single A156T or V mutations listed in color on the right-hand side of Fig. 1 and two sequences contained a single R155G mutation. These 58 sequences containing single drug resistance mutations were widely dispersed throughout the phylogenetic tree of pretreatment sequences (Fig. 1). This indicates that nearly all of these drug-resistant sequences represented distinct viral lineages that traced back 14 d earlier to infected hepatocytes that carried A156T or A156V mutations in widely divergent genetic backgrounds. In marked contrast, 51 of 111 sequences (46%) contained a R155W-plus-A156G double mutation, in a nearly homogeneous genetic background where sequence diversity ranged from 0 to 1 nucleotide substitution over a span of 2,210 bp (<0.1% diversity). This result suggests that ∼46% of hepatocytes, releasing HCV virions into the circulation at 14 d after the initiation of MK-5172 treatment (and 7 d after its discontinuation), were infected by HCV genomes that had evolved from a single R155W-plus-A156G drug-resistant viral sequence lineage. In turn, this lineage most likely emanated from a single or a few productively infected hepatocytes that existed before treatment initiation (Fig. 1). This inference is drawn based of the identity or near identity (1 mutation in 2,210 nt) of the genetic background in which the R155W-plus-A156G drug resistance mutation expanded. These findings are consistent with model predictions of the expected frequency of a single mutation or any set of two mutations in an HCV vRNA quasispecies in the absence of selection (5, 28). In addition, one sequence with Y56H-plus-D168V and one with Y56H-plus-D168N double mutations were also detected. At day 27 when the plasma viral load had returned to baseline (>106 vRNA/mL), there was a striking contraction in both the R155W-plus-A156G and A156T or V populations and a remarkable expansion of the discrete homogeneous Y56H-plus-D168V and Y56H-plus-D168N double-mutant populations (already detected at day 14). A single dominant Y56H-plus-D168N double-mutant population and three dominant Y56H-plus-D168V double-mutant populations together represented 70% of the sequences. This finding is all of the more remarkable given that, by this time, steady-state plasma vRNA load had returned to its pretreatment setpoint of >106 vRNA copies per mL and drug therapy had been stopped 20 d previously. Two additional monophyletic lineages, each containing the double-mutant D168E-plus-F169I, were present on day 27. By day 56, the D168E-plus-F169I mutant lineages (first detected at day 27) expanded to comprise 79% of the sequences, whereas the previously dominant Y56H-plus-D168N and Y56H-plus-D168V lineages contracted and now comprised only 12% of the sequences. Four wild-type sequences, that is, sequences with no known resistance mutations, were detected for the first time posttreatment at day 56.

Fitting a Baseline Viral Dynamic Model to the Viral Load and Sequence Data.

The rapid viral load decline during treatment and continuous turnover of resistant mutants after treatment cessation observed in the treated subjects raise intriguing questions, such as how do resistance clones expand so rapidly and how can resistant mutants dominate the viral population for such a long period after treatment cessation without being replaced by the nonresistant virus? To address these questions, we first constructed a “baseline” multistrain HCV model, similar to a previous standard viral dynamic model that incorporates competition for target cells (5) (). The strains in the model are grouped according to the shared drug resistance mutations in the clinical samples as described in the section above and as color-coded in Fig. 1 and (see , for the mutants modeled). In this model, we assume that infected cells are lost at a per-capita rate δ and that target cells proliferate following a logistic growth term to replace them. Cure and superinfection of infected cells are not included in this model.

We fitted this model to the data from all five treated subjects. This model describes both the viral load and the sequence data well (). However, the mean loss rate of infected cells, δ, estimated across all five treated subjects was 0.68 d−1 with a SD of 0.09 d−1 (see , for the best-fit values of δ). In our model, the loss of infected cells corresponds to death of these cells. A large death rate is needed to explain the rapid second-phase viral load decline during treatment and rapid turnover of dominant mutants after treatment. However, such a large death rate of infected cells (0.68 d−1) does not seem physiological, as it implies a half-life for infected hepatocytes of 1 d. In stark contrast, the death rate of infected cells has been estimated as 0.14 d−1 from IFN therapy (29, 30) and previously used by Guedj et al. (31) and Lau et al. (21) to explain viral declines seen with other DAAs. To confirm the necessity of a high loss rate in this model, we refitted the model fixing δ at 0.14 d−1 and found this alternative parameterization explains the data poorly (). These results strongly suggest that other mechanisms (in addition to death of infected cells) must play essential roles in driving the patterns seen in the data.

Potential Roles for Infected Cell Cure and Intracellular Viral Competition Through Superinfection.

We hypothesize that cure of infected cells is the primary cause for the loss of infected cells during the rapid second-phase decline observed in our data and other studies (5, 8, 32). In addition, we hypothesize that the continuous turnover of dominant mutants in the absence of treatment is mostly driven by intracellular competition among viral strains, as a result of superinfection (rather than competition for newly generated target cells). Previously, cure of infected cells has been demonstrated in in vitro studies (33–37). In vivo, the rapid second-phase viral declines under treatment were suggested to be attributable to cure of infected cells (32, 38). For superinfection, in vitro experiments showed that multiple viruses can enter the same cell, although limiting host resources necessary for viral replication and/or translation may restrict active production of multiple viruses (39, 40). Interestingly, a recent study demonstrated that a fitter HCV strain can enter already-infected cells and outcompete the resident strain (41), suggesting that multiple HCVs can enter a cell and compete for intracellular resources.

To test these hypotheses, we extended the baseline model by incorporating infected cell cure and superinfection, and fixing the rate of infected cell death at 0.14 d−1 as estimated under IFN therapy (29, 30) (). Here, we assumed that infected cells are cured and become target cells again under MK-5172 treatment at per-capita rate , where is a rate constant and is the drug efficacy against the ith strain (see for detail). To model the intracellular viral competition due to superinfection, we assume, for simplicity, that once an infected cell is superinfected, the fitter strain can outcompete the less fit strain intracellularly, and the cell becomes a cell infected by the fitter strain at rate ksuper.

We used three variations of the extended model to test the importance of cure and superinfection in explaining the kinetic patterns in the data: (i) we allowed cure but no superinfection (denoted as the “cure” model), or (ii) superinfection but no cure (denoted as the “superinfection” model), or (iii) both processes (denoted as the “full” model) (see and Fig. 2, respectively, for model fits and for the best-fit parameter values). To statistically test the importance of cure or/and superinfection in explaining the clinical data, we performed model selection using the corrected Akaike information criterion (AICc) (Table 1). The results suggest that cure of infected cells is needed to explain the rapid second-phase viral load decline during treatment seen in all five subjects (). It also helps to explain the rapid selection of the resistant mutants seen at day 14 in subjects 1 and 4 (compare the fits between solid lines and dashed lines to data in ). Based on the best-fit parameters in the best model for each subject (), we estimated the mean rate of cure of cells infected by baseline viruses under MK-5172 treatment [], where is the drug efficacy against the baseline viruses) to be 0.56 d−1. On average, it takes about 1.9 d to cure a cell under MK-5172 treatment () and the loss of infected cells is primarily through cure rather than death of infected cells (assumed to occur at rate of 0.14 d−1).

Fig. 2.

Best fit of the “full” model (lines) to the clinical data (circles) from five subjects treated with MK-5172. (A–E) Schematic diagrams of the evolutionary dynamics in subjects 1–5, respectively. A–E, the data and simulation results for viral loads are shown in open circles and black lines, respectively, on the Left; the data and simulation results for mutant frequencies are shown in colored open circles and lines, respectively, on the Right. The color coding for each mutant considered is shown in .

Table 1.

Summary of the model characteristics and the fitting results (i.e., AICc scores) of each model for each subject

Model	Model characteristics		Fitting results (AICc)
	kcure	ksuper	Subject 1	Subject 2	Subject 3	Subject 4	Subject 5	Total
Baseline model with δ = 0.14 d−1	0.0	0.0	−16.2	−22.1	−8.1	5.0	−24.9	−66.3
Cure model	Fitted	0.0	−32.8	−53.2	−47.6	−42.6	−86.5	−262.7
Superinfection model	0.0	Fitted	−12.0	−17.8	−9.5	−5.1	−21.2	−65.6
Full model	Fitted	Fitted	−63.8	−51.8	−44.5	−90.6	−83.6	−334.3

Bolded AICc scores denote the best model fit among all models for the five subjects.

The AICc results suggest that intracellular competition through superinfection is needed to explain the rapid and continuous turnover of dominant resistant strains after treatment stops in subjects 1 and 4 (). This is because after treatment cessation, drug is eliminated and drug-induced cure is no longer possible. As a result, new target cells become available at a much slower rate than the rate during treatment. However, superinfection allows fitter strains to enter already-infected cells and compete intracellularly, greatly increasing the rate at which a fitter strain rises to a high frequency. In contrast, without superinfection and intracellular competition, a fitter strain is predicted to increase in frequency at a much lower rate (see the changes in the frequencies of the mutant in green during the period of days 30–60 in ).

We performed uncertainty analysis in the best model for each subject using likelihood profiling (42) (). These analyses indicated that the data allow us to estimate accurately the effectiveness of MK-5172 against the baseline virus, ε1, the rate of viral clearance, c, and the cure rate kcure in all five subjects, and the value of superinfection constant ksuper in subjects 1 and 4. In addition, we explored whether models with alternative assumptions can explain the data. This included a model assuming a constant rate of target cell generation and models assuming a DAA-independent cure of infected cells (). Fitting results show that the full model presented in the main text is the best model, and the conclusions about the role of superinfection and cure of infected cells are robust against these model variations.

Compensatory Mutations, Clonal Interference of HCV-Resistant Mutants.

To understand better the evolutionary dynamics of the viral strains considered in our model, we analyzed how the strains at different time points are related to each other. Here, we define strains phenotypically, that is, each “strain” consists of a group of related viral sequences in the clinical samples sharing drug resistance mutations. We calculated the genetic distances between sequences belonging to different groups at the different sampling time points (see and ). We then combined the mutational pattern derived from the sequence data with the results from fitting the dynamical models to generate a diagram that summarizes the evolutionary dynamics of the drug-resistant strains over the period of the study (Fig. 3).

Fig. 3.

Evolutionary dynamics of HCV before and after treatment with MK-5172. (A–E) Schematic diagrams of the evolutionary dynamics in subjects 1–5, respectively. Each mutant is denoted as a colored circle, and the color coding is shown at the left-hand side of the y axis in each panel. The numbers within the circles denote the fitness of the mutant relative to the baseline strain in the absence of treatment according to the best-fit parameter values in the best model for each subject. The x axis shows the time when sequence data are taken. The size of the circle is scaled according to the frequency of the mutant in a given sample. Solid arrows show the strain (where the arrow starts; say, strain a) to which the majority of the sequences in a mutant strain of interest (where the arrow ends; say, strain b) are most closely related, and the numbers on each arrow show the number of sequences in strain b that are mostly closely related to strain a and the total number of sequences in strain b.

The diagram reveals how the HCV population responds to MK-5172 treatment. First, viral load rebounded rapidly in subjects 1 and 4 to high levels that permitted sequencing by day 14, presumably because of selection of the resistant mutants bearing A156T/V or R155G/W in subject 1 and A156T/V in subject 4. We estimated that these mutants are highly resistant, yet have high fitness costs (Fig. 3), which is consistent with results shown in previous in vitro studies (24–26). Because of their short genetic distance to the wild-type viruses (), they are likely present at low frequencies before treatment. During and shortly after treatment, they rise to high frequencies transiently before being replaced by other resistant mutants with higher fitness while drug is cleared. Second, mutants with the Y56H mutation are observed transiently in subjects 1, 3, and 4 between days 27 and 34. Our analysis suggests that the transient appearance of these mutants is due to their intermediate resistance and fitness values. These mutants gave way to fitter strains with resistant mutations at position 168 (e.g., D168E, D168A, or D168Y; see refs. 24 and 25) that dominated the viral population in subjects 1–4 at the end of follow-up period, that is, 7–8 wk after treatment cessation. These dynamics resemble the classical pattern, termed “clonal interference,” as suggested by population genetic models for the evolution of asexual organisms (43, 44).

The fitnesses of the resistant mutants observed at late time points are predicted to be close to or higher than those of the baseline viruses, which remained at low frequencies in subjects 1–4 (Fig. 3). A likely explanation for the late appearance of these high-fitness variants is that mutation at position 168 on the baseline virus background resulted in low fitness (and thus these mutants are at low frequency before treatment). However, growth of these low-fitness resistant mutants generated compensatory mutations that increased their fitness, and allowed their frequency to increase and stabilize at a high level in the population.

Last, we observe that drug-sensitive viruses, that is, viruses with no known drug-resistant mutations, reappeared in the sequences at the last time point of follow-up in subjects 1, 4, and 5. Statistical analysis shows that these viruses are likely to be descendants of drug-sensitive viruses present before treatment rather than reversions from resistant variants (). This suggests that the drug-sensitive viruses were suppressed to low frequencies under treatment and only grow to high frequencies long after treatment cessation.

Discussion

Here, we have analyzed the evolutionary dynamics of the HCV population in response to 7-d treatment with the protease inhibitor grazoprevir (MK-5172) in five subjects using single-genome sequencing and phylodynamic modeling. HCV-resistant variants underwent rapid expansion after drug treatment with surprisingly fast sequential turnover of dominant resistant mutants in the 7- to 8-wk follow-up period, a distinct feature of viral persistence. Using mathematical models, we showed that superinfection and cure of infected cells are likely to be crucial mechanisms driving the extremely rapid expansion and turnover of resistant variants observed in the clinical data. Although superinfection and cure of infected cells have been demonstrated in vitro (33–37, 41), our work suggests that these processes occur and play important roles in driving viral adaptation and persistence in vivo.

We estimated the rate of cure of cells infected by wild-type viruses under grazoprevir to be 0.56 d−1 on average. Further, it takes about 1.9 d on average to cure a cell under MK-5172 treatment. This high rate of infected cell cure under treatment is consistent with the clinical data in this study and from previous studies (8, 30, 32, 45), where rapid second-phase viral load declines have been observed. Previous mathematical models suggested that the second-phase decline reflects both the death rate of infected cells and the rate of intracellular vRNA degradation (31, 32, 38). Thus, this rapid cure of infected cells is likely to be a result of the rapid loss of intracellular HCV RNAs under treatment. Our model shows that, due to the rapid cure, drug-sensitive viruses were at a much lower frequency than resistant mutants when treatment stopped. That coupled with the rapid expansion of resistant mutants and generation of compensatory mutations led to the disappearance or a delay in the reappearance of drug-sensitive virus.

Superinfection has been demonstrated in vivo indirectly by the identification of recombinant forms of HCV in patients (46–49). Although recombination between HCV genomes is apparently rare (50, 51), this does not mean that superinfection of hepatocytes is infrequent, because the formation of distinct intracellular HCV replication complexes and a processive RNA polymerase may prevent recombination. In line with a recent in vitro study showing that a fitter HCV strain can enter already-infected cells and outcompete the resident strain (41), our results suggest that superinfection and intracellular competition occur frequently in vivo, especially when different strains exhibit a difference in intracellular replication fitness.

We argue that superinfection and intracellular competition are critical for the establishment and expansion of adaptive mutants and thus for the persistence of HCV under selection pressures. For viruses that largely transmit by cell-to-cell infection, such as HCV (52–55), transmission is often limited from an infected cell to a few neighboring cells. When most of the cells are infected by resident viral variants, without superinfection, the probability that an adaptive mutant finds an uninfected cell is low. Thus, adaptive mutants would likely be lost, and for those mutants that do establish infection, the rate at which they expand would be low. In contrast, with superinfection and intracellular competition, an adaptive mutant can enter all neighboring cells to compete with any resident virus, which substantially increases the probability of establishment and the rate at which the adaptive mutant expands. While previous theoretical work discussed the role of superinfection in viral evolution at the host level (56, 57), our work highlights the important role that superinfection plays in facilitating viral adaptation at the cellular level.

By combining phylogenetic analysis and mathematical modeling, we provided an integrated understanding of the HCV evolutionary and population dynamics. First, the HCV within-host evolution of HCV resembles the classical pattern for the evolution of asexual organisms as suggested by population genetic models, termed clonal interference (43, 44). This emphasizes the usefulness of population genetic models in understanding viral evolutionary dynamics. Second, our results suggest that compensatory mutations are likely to be generated after treatment, allowing for stabilization and persistence of resistant mutants in the viral population long after treatment cessation. This highlights the need to maintain a high level of adherence to combination therapies of DAAs. Treatment interruption would allow rapid expansion of viral mutants resistant to one or more DAAs and possible generation of compensatory mutations that stabilize the resistant mutations at a high frequency, and thus creating opportunity for transmission of resistant mutants. In situations where treatment is interrupted frequently, multiple mutations may accumulate such that the evolved virus may become resistant to all DAAs in the combination therapy (12).

Although beyond the scope of this work, a couple of predictions of our study can be further tested experimentally. First, our models suggest that superinfection and intracellular competition are important mechanisms facilitating the rapid expansion and turnover of resistant mutants in vivo after treatment cessation. Although testing this prediction in vivo is not possible, intracellular viral competition experiments as in ref. 41 can be performed to compete pairs of laboratory viruses with the resistant mutations identified in this study. The results would reconfirm the occurrence of superinfection and determine the relative fitness of those mutants at the intracellular scale. Second, our analysis suggests that compensatory mutations are involved in stabilizing the resistant mutations at amino acid 168 seen in subjects 1–4. This can be tested by first performing whole-genome sequencing on the viral isolates and then identifying the shared mutations of the viruses with a mutation at position 168 in each subject. These shared mutations would be candidate compensatory mutations for the mutation at amino acid 168. Introducing one or a combination of these mutations into viral replicons, and then testing and comparing the replication fitness of the replicons (as shown in refs. 24–26), would be a means to identify compensatory mutations.

To summarize, the results in our study elucidate important, yet previously unrecognized, mechanisms that we suggest operate in vivo, drive within-host viral evolution, and allow the adaptation and persistence of HCV in the face of drug pressure and immune response. A previous framework frequently used in viral dynamic models (5, 8, 13, 16–18) assumed that viruses compete for the replication space needed for adaptation to changing host environments only through infection of newly generated target cells (Fig. 4). Consequently, the rate of adaptation was limited by the rate of target cell generation. In contrast, our framework emphasizes that both superinfection and cure of infected cells (in addition to infection of new target cells) can be important mechanisms contributing to the replication space needed for mutant virus expansion (Fig. 4). This mechanism can greatly accelerate the rate of within-host viral adaptation beyond the rate set by the generation of target cells. This framework may be important in elucidating the contributions and mechanisms of HCV containment and/or elimination mediated by the host immune system that occur naturally or following vaccination (58). More broadly, our work points toward the important role molecular mechanisms play in facilitating viral evolution and persistence. Frequent cell-to-cell infection and high multiplicity of infection, that is, mechanisms that provide the means for intracellular competition, have been reported recently for HIV, HCV, influenza, and other viruses (54, 59–61). We argue that our framework and the consideration of viral competition at both the intracellular and host level in particular, is critical for the understanding of the evolutionary and adaptive dynamics of viral population at the within-host level (62, 63). Thus, it opens the door to the next generation of viral dynamic models and quantitative frameworks to understand viral evolution, interpret clinical datasets, and predict treatment outcomes.

Fig. 4.

Conceptual frameworks for virus evolution under antiviral pressure. Resistant mutants expand through occupying/competing for available replication space. (A) Previous modeling mostly assumed that, once a cell is infected by a virus, the cell remains infected until death. The replication space arises through generation of new target cells. Under this framework, resistant mutants expand through infection of newly generated target cells, and the rate of the increase of mutant frequency is mostly set by the rate at which infected cells die and are replaced by newly generated cells. (B) Our results suggest a conceptual framework where the replication space arises from multiple sources. In the presence of potent antivirals, the level of intracellular HCV RNAs decreases, leading to cure of infection in some cells. Replication space thus arises from both newly generated cells and cured cells. In addition, superinfection makes replication space available by allowing resistant viruses to enter an already-infected cell and compete for intracellular resources. Thus, cure and superinfection of cells allow resistant mutant expansion and turnover to occur at a much faster rate than the rate set by the death and replacement of infected cells only.

Methods

Experimental Model and Subject Details.

For this study, Merck Sharp and Dohme Corporation provided preexisting, deidentified human blood plasma specimens collected from a multiinstitutional experimental drug treatment protocol (MK-5172-004) designed to test the antiviral activity of the HCV protease inhibitor grazoprevir/MK-5172 (ClinicalTrials.gov identifier NCT00998985). The specimens were collected, stored, and inventoried with study subjects’ informed consent and approval for future use in the analysis of the effects of grazoprevir/MK-5172 on HCV replication, although the specimens were not collected specifically for the research project described here. Based on these stipulations, the research described in the current report was deemed by the Institutional Review Board of the University of Pennsylvania not to constitute human subjects in research.

Single-Genome Sequencing.

The partial NS2 and complete NS3 and NS4A gene sequences from all eight subjects were generated using the single-genome sequencing method previously described (50). The primers used for the nested PCR included the following: (i) genotype 1a: first-round sense primer 1aNS2.F2, 5′-ACCCGRCTTTGGTATTTGACATCACC-3′ (nucleotides 2983–3008, H77); first-round antisense primer 1aNS4B.R3, 5′-TATTGTATCCCACTGATGAAGTTCCACAT-3′ (nucleotides 5634–5662, H77); second-round sense primer 1aNS2.F3, 5′-AAAGTGCCCTACTTYGTGCGCGT-3′ (nucleotides 3063–3085, H77); and second-round antisense primer 1aNS4B.R4, 5′-AGGGCCTTCTGCTTGAACTGCTC-3′ (nucleotides 5517–5539, H77); (ii) genotype 1b: first-round sense primer 1bNS2.F2, 5′-GCCCGTCGTCTTYTCTGACATGGA-3′ (nucleotides 3257–3280, H77); first-round antisense primer 1bNS4B.R3, 5′-TTCCACATGTGCTTCGCCCA-3′ (nucleotides 5622–5641, H77); second-round sense primer 1bNS2.F3, 5′-TCATCACCTGGGGGGCAGACA-3′ (nucleotides 3289–3309, H77); and second-round antisense primer 1bNS4B.R4, 5′-CGAGCGCCTTCTGCTTGAATTG-3′ (nucleotides 5520–5541, H77).

Viral Load and Sequence Diversity Analysis.

The vRNA measurement was done by Merck and Co., Inc., using the Roche Cobas TaqMAN 2.0 assay with a lower limit of quantification of 25 IU/mL and a limit of detection of 10 IU/mL. A total of 2,755 NS2–NS3–NS4 genomes were generated and analyzed from the eight subjects. The median number of sequences analyzed per time point was 107 (mean, 106; range, 43–144) (). Sequences alignments were initially made with clustalW and then hand-checked using Geneious for correct codon alignments. The maximum, minimum, median, and mean diversity for each sequence set was calculated based on Hamming distance (). For each subject, phylogenetic trees were generated by maximum-likelihood methods using PhyML. The combined phylogenetic tree and the tree from each time point were rooted with consensus sequences from day 0.

Construction of a Baseline Multistrain HCV Model.

We construct a viral dynamic model for within-host infection by multiple strains of HCV based on previous published models (5, 29, 64). In this model, we keep track of n different HCV strains. The ordinary differential equations (ODEs) for the model are as follows:

Target cells (T), presumably a subset of hepatocytes, die at a per-capita rate, d. Existing target cells can proliferate, and this proliferation is modeled using a logistic term, , as in Rong et al. (5), where is a proliferation rate constant, I are cells infected with virus strain i, N is the concentration of hepatocytes that are not target cells, and Tmax is the liver carrying capacity. Target cells are infected by the ith strain of HCV () at rate . Infected cells die at per-capita rate, . Drug-sensitive viruses, represented by strain i = 1, are produced from infected cells at rate p per cell in the absence of treatment. We assume that the ith strain has a fitness of r relative to the drug-sensitive viruses (r1 = 1). This difference in fitness is reflected only in the differences in viral production in our model. Thus, the production rate of the ith strain is in the absence of treatment. Under treatment, we assume that viral production of the ith strain is reduced by a factor , where corresponds to a 100% effective drug. The drug effectiveness ε is modeled as a function of the drug concentration and the EC50 value for each strain (EC50,) according to an Emax model (65). We assume that the drug concentration, D, stays constant during the first 7 d of treatment, since drug concentration reaches its maximal concentration within 2–4 h (23). After treatment stops on day 7, we assume the drug concentration declines exponentially at rate, w, obtained previously (23). All viruses are cleared at per-capita rate c. The values of the parameters are shown in .

Construction of Models with Cure and Superinfection of Infected Cells.

We constructed a model incorporating cure and superinfection of infected cells by extending the basic model. The ODEs for this model are as follows:

To model cure of infected cells, we assumed that the rate of cure of infected cells is linearly dependent on the log10 of the efficacy of the drug, , and the rate constant is kcure. This assumption is based on a previous finding with the protease inhibitor telaprevir that the rate of second phase of viral load decline increased linearly with − (32).

The superinfection we track here includes the infection of a cell already infected with strain i with a virus of strain j, and the subsequent conversion of the cell into a cell that produces virus of strain j. We assume, for simplicity, that, if a fitter virus, strain i, enters a cell infected with a less fit strain, strain j (modeled using the term ), it can outcompete strain j, and then the cell is converted to a cell productively infected by strain i. We further assume the rate of this conversion is linearly dependent on the fitness difference () of the two strains, , and ksuper is a constant for the efficiency of this conversion.