Epistasis and entrenchment of drug resistance in HIV-1 subtype B

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: October 8, 2019 (This version)
Accepted: September 9, 2019
Received: July 25, 2019

1. Of interest
Shared enhancer gene regulatory networks between wound and oncogenic programs

Swann Floc'hlay, Ramya Balaji ... Stein Aerts

Research Article Updated Jun 2, 2023
Further reading

Article
Figures and data
Abstract
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

The development of drug resistance in HIV is the result of primary mutations whose effects on viral fitness depend on the entire genetic background, a phenomenon called ‘epistasis’. Based on protein sequences derived from drug-experienced patients in the Stanford HIV database, we use a co-evolutionary (Potts) Hamiltonian model to provide direct confirmation of epistasis involving many simultaneous mutations. Building on earlier work, we show that primary mutations leading to drug resistance can become highly favored (or entrenched) by the complex mutation patterns arising in response to drug therapy despite being disfavored in the wild-type background, and provide the first confirmation of entrenchment for all three drug-target proteins: protease, reverse transcriptase, and integrase; a comparative analysis reveals that NNRTI-induced mutations behave differently from the others. We further show that the likelihood of resistance mutations can vary widely in patient populations, and from the population average compared to specific molecular clones.

https://doi.org/10.7554/eLife.50524.001

Introduction

HIV mutates rapidly as it jumps from host to host, acquiring resistance to each host’s distinct immune response and applied drug regimen. Drug-resistance mutations (DRMs) arise when the virus evolves under selective pressure due to antiretroviral therapy (ART). Primary DRMs often incur a fitness penalty which is then compensated for by accompanying associated mutations (Heeney et al., 2006; Shafer and Schapiro, 2008). With the use of current robust inhibitors in drug therapy, the drug-resistance mutation patterns in HIV have become increasingly more complex (Richman et al., 2004a; Iyidogan and Anderson, 2014) often leading to ART failure in patients. Resistance is estimated to develop in up to 50% of patients undergoing monotherapy (Richman et al., 2004b) and up to 30% of patients receiving current combination antiretroviral therapy (c-ART) (Gupta et al., 2008). The primary drug targets in treatment of HIV are the enzymes coded by the pol gene, reverse transcriptase (RT), protease (PR), and integrase (IN). A large number of sequences of HIV are available for RT, PR, and IN for patients who have been treated during the past nearly 30 years, and this information permits critical sequence-based informatic analysis of drug resistance.

The selective pressure of drug therapy modulates patterns of correlated mutations at residue positions which are both near and distal from the active site (Chang and Torbett, 2011; Haq et al., 2012; Flynn et al., 2015; Yilmaz and Schiffer, 2017). A mutation’s impact on the stability or fitness of a protein however is dependent on the entire genetic background in which it occurs: a phenomenon known as ’epistasis’. Drug resistance develops as these mutations accumulate, providing the virus a fitness benefit in the presence of drug pressure, with a complex interplay in the roles of primary and secondary mutations (Yilmaz and Schiffer, 2017; Ragland et al., 2017). When a primary resistance mutation is incurred in the context of a wild-type background, there is usually a fitness penalty associated with it. In backgrounds with more (accessory) mutations however, the fitness penalty decreases and on average, the primary mutation can become more likely than the wild-type residue. Because the beneficial effects of the associated mutations depend on the primary mutation, with the accumulation of (accessory) mutations, the reversion of the primary mutation can become increasingly deleterious, leading to a type of evolutionary ‘entrenchment’ of the primary mutation (Pollock et al., 2012; Shah et al., 2015; McCandlish et al., 2016). The entrenchment effect on a primary mutation can be very strong on average, and is in fact, modulated by the collective effect of the entire sequence background.

The effective modeling of epistasis is then critical to the identification and understanding of the drug and immune pressure mediated mutational combinations that give rise to drug-resistant, stable viruses. Experimental techniques to assess the effect of multiple mutations on phenotype have proven effective (Troyer et al., 2009; da Silva et al., 2010; Liu et al., 2013), but functional assays to test all possible combinations are impossible because of the vast size of the mutational space. Co-evolutionary information derived from multiple sequence alignments (MSAs) of related protein sequences have also served as a basis for building models for protein structure and fitness (Göbel et al., 1994; Lockless and Ranganathan, 1999; Morcos et al., 2011; Hinkley et al., 2011; Haq et al., 2012; Ferguson et al., 2013; Mann et al., 2014; Jacquin et al., 2016; Hopf et al., 2017; Tubiana et al., 2019). A subset of such models, called Potts statistical models (Levy et al., 2017) (a generalization of Ising spin glass models), have been used successfully to predict the tertiary and quaternary structure of proteins (Morcos et al., 2011; Marks et al., 2012; Sułkowska et al., 2012; Sutto et al., 2015; Haldane et al., 2016; Levy et al., 2017; Sjodt et al., 2018) as well as viral protein stability and fitness (Shekhar et al., 2013; Barton et al., 2016a; Hopf et al., 2017; Flynn et al., 2017; Levy et al., 2017; Louie et al., 2018).

Recent studies using Potts models inferred the epistatic interactions involved in the evolution of drug resistance in HIV-1 protease (Flynn et al., 2017). The authors identified sequence backgrounds with differing patterns of associated mutations which were predicted to either strongly entrench or strongly disfavor particular primary mutations in HIV-1 protease. Of crucial importance is the observation that whether a mutation favors or disfavors a primary resistance mutation with which it is associated, depends on the entire sequence background. Here we confirm the predictions in Flynn et al. (2017) providing the first direct confirmations of ‘entrenchment’ by the entire sequence background and show how such models capture the epistatic interactions that lead to drug resistance in each of the three predominant target enzymes in HIV subjected to selective pressure of ART: protease, reverse transcriptase, and integrase. We show that it is possible to predict which mutations a particular background will support in the drug-experienced population and how conducive that background is towards that mutation. Understanding the relationship between the likelihood of a drug-resistance mutation and the specific sequence background is important for the interpretation of mutagenesis experiments, which are carried out in the background of specific reference sequences such as subtype B molecular clones NL4-3, HXB2, and others (Martinez-Picado et al., 1999; An et al., 1999). Our results demonstrate that the background-dependent likelihood of a primary mutation and the strength of compensatory behavior can vary widely in patient populations, and from population averages to predominant subtype B molecular clones such as NL4-3 or HXB2.

It is useful to review some basic properties of Potts models in this context. A Potts model is a maximum-entropy model of the likelihood of a dataset multiple sequence alignment (MSA), constrained to predict the bivariate (pairwise) residue frequencies between all pairs of positions in the alignment, in other words it is designed to capture the observed pairwise mutational correlations caused by epistasis. A central quantity known as the ‘statistical’ energy of a sequence $S$ (see Equation 1, Materials and methods) is commonly interpreted to be proportional to fitness; the model predicts that sequences will appear in the dataset with probability $P (S) \propto e^{- E (S)}$ , such that sequences with favorable statistical energies are more prevalent in the MSA. A key feature of the Potts model is that the effect of a mutation on fitness and $E (S)$ is ‘background-dependent’, as a single change at one position causes a difference in $L - 1$ of the coupling values in the sum for $E (S)$ , and a mutation at one position will affect mutations at all other positions both directly and indirectly through chains of interactions involving one or many intermediate residues. Because of this, the Potts model predicts complex patterns of mutational correlations, even though the effective energy model includes only pairwise terms.

Results

Epistasis: the effect of the background on primary resistance mutations in HIV

We begin by demonstrating the importance of the genetic background in determining the frequencies of primary mutations in HIV and verifying this by comparison to the empirical background dependence of primary mutations observed in the Stanford HIV database. We do this using a background-dependent computation of residue biases. If we mutate a sequence $S$ at position $i$ to a residue $α$ , calling the mutated sequence $S_{α}^{i}$ , we can also calculate the mutant sequence’s probability $P (S_{α}^{i})$ according to the Potts model. Further, we can calculate the relative probability of any mutation at position $i$ of sequence $S$ , as $ℙ (S, i, α) = P (S_{α}^{i}) / \sum_{β} P (S_{β}^{i})$ . Notably, $ℙ (S, i, α)$ depends on the background comprised of all positions except $i$ . We will call this background $S ∖ i$ representing the sequence $S$ with an unspecified residue at position $i$ . Given a very large sample of sequences such that the background $S ∖ i$ is sampled multiple times, a residue $α$ will appear at position $i$ in that background with frequency $ℙ (S, i, α)$ . In practice, obtaining such large samples is difficult. But equivalently, if a background $S ∖ i$ is sampled from nature without knowledge of the residue at position $i$ , the Potts model can be used to predict that the residue $α$ will be at the missing position with a likelihood $ℙ (S, i, α)$ . We emphasize that the model was not fit using this information. While the Potts Hamiltonian model only contains fields and pairwise coupling terms; $ℙ (S, i, α)$ involves combinations of couplings of residues at $i$ with all other positions and cannot be expressed as a simple function of the correlations between mutations at pairs of positions.

Using $ℙ (S, i, α)$ we can classify sequences by how likely a residue $α$ is to appear at a position $i$ (Figure 1A). We divide the dataset MSA into 10 groups of sequences ordered by $ℙ (S, i, α)$ , and compare the Potts predicted frequency of the residue in each group to its observed frequency. Figure 1B shows good agreement for the likelihood of M at position 90 in HIV-1 protease, an important drug-resistance mutation. Figure 1D shows that the average absolute error between Potts predicted and observed frequencies is very small for all major drug-resistance mutations in HIV PR, RT and IN, establishing that the Potts model is a good predictor of epistasis in viral sequence backgrounds. These results provide one of the most direct verifications of the complex patterns (involving the effect of many mutations) of epistatic interactions involving drug-resistance mutations in HIV. This result also shows that the background dependence can be very strong, as many sequences can have likelihoods of $ℙ (S, i, α) > 0.9$ while others have $ℙ (S, i, α) < 0.1$ as is seen for the DRM L90M in HIV PR (Figure 1B).

Figure 1 with 1 supplement see all

Download asset Open asset

The Potts model predicts residue frequencies.

(A) Schematic showing that the Potts model can be used to classify sequences by how likely a residue $α$ is to appear at a position $i$ in a sequence $S$ using the background-dependent probability, $ℙ (S, i, α)$ . (B) The observed frequency of the resistance mutation L90M in HIV-1 drug-experienced proteases matches the Potts model predicted frequencies in sequence clusters binned according to the Potts-predicted frequencies in steps of 0.1 (blue circles) with statistical error (green). Diameters of the circles represent the number of sequences. Inset shows the significant overlap in Hamming distance for sequences with low predicted mutant frequencies between 0.2 and 0.3 (blue) and high, 0.7 and 0.8 (pink) depicting the difficulty of such a classification based on Hamming distance. (C) The receiver operator characteristic (ROC) curve comparing the Potts model and Hamming distances as classifiers of mutational probabilities for L90M in HIV protease. (D) The average absolute error between the observed mutational frequencies and the Potts-predicted frequencies for the major drug-resistance mutations in HIV-1 in three drug targets. The average absolute error is calculated by binning the sequences in ascending order of their predicted frequencies such that there are roughly equal number of sequences in each bin as shown in Figure 1—figure supplement 1, and averaging over the absolute error in each bin.

https://doi.org/10.7554/eLife.50524.002

The epistatic predictions of the Potts model cannot be simply explained by the number of background mutations (or Hamming distance) relative to the wild-type (consensus wild-type subtype B amino acid sequence). It is not simply the number of mutations that bias residues, but particular patterns of mutations. The inset figure of Figure 1B shows that there is significant overlap in the distributions of the number of mutations between two groups of sequences in which the probability of observing the mutation L90M is ~3x smaller in one than the other (‘blue’ compared to ‘red’) illustrating that a prediction of residue frequencies based on Hamming distance alone is insufficient. We compare the residue-prediction ability of the Potts model using $ℙ (S, i, α)$ to that of a ‘Hamming’ model which predicts residue frequencies purely based on the number of mutations. This is plotted in Figure 1C for the L90M mutation in HIV protease. There is a significant performance gain for the Potts model with a true positive rate (TPR) of 0.81 at a false positive rate (FPR) of 0.1, compared to a TPR of 0.46 at an FPR of 0.1 for the Hamming distance model, again showing that the latter model is much less predictive.

‘Entrenchment’ of a primary resistance mutation and its verification

The log likelihood of a mutation occurring at a position $i$ in a given sequence background relative to the wild-type residue (consensus) at that position can be calculated from the change in the Potts statistical energy on acquiring that mutation from the wild-type residue at that site as $Δ E^{i} = E_{w i l d - t y p e}^{i} - E_{m u t a t i o n}^{i}$ . If $Δ E^{i} > 0$ , this implies that the mutation is more favorable in the given sequence background and that there is a fitness penalty associated for reversion back to the wild-type residue; while $Δ E^{i} < 0$ implies that the wild-type residue is more favorable and the mutation is likely to revert to the wild-type residue. A single resistance mutation in the wild-type background is almost always disfavored relative to the wild-type residue, and is generally disfavored ( $Δ E < 0$ ) in backgrounds which are close to the wild-type (i.e carry a few accessory mutations). However, as the number of mutations in the background increases, the preference for the wild-type residue at position $i$ decreases on average, until eventually the mutant is preferred at that position. We describe a drug-resistance mutation at position $i$ as being ‘entrenched’ in a particular sequence S if $Δ E^{i} > 0$ at that position. When a drug-resistance mutation is entrenched at $i$ in sequence S, this means that the drug-resistance mutation is more probable than the wild-type residue at $i$ in sequence S.

Figure 2A illustrates entrenchment for the drug-resistance mutation L90M in HIV-1 proteases. The preference to revert as a function of the number of mutations at PI-associated sites is shown as box plots. Each box plot shows the median and mean values of the preference to revert for sequences, conditional on the indicated number of mutations at PI-associated sites in the corresponding column, as well as the region containing half of the sequences (shown in gray bounded by blue/red boxes); the whiskers of the box plot indicate the range of $Δ E$ values for sequences in the tails of the distributions conditional on the number of mutations at PI-associated sites. For drug-experienced viral PR sequence backgrounds that contain eight or fewer PI-associated mutations, (Figure 2A left area in blue), the L90M mutation is disfavored on average and likely to revert, with more than a 10 fold preference to revert in backgrounds with no inhibitor-associated mutations. However, the box plots show that the sequences, conditional on a fixed number of PI-associated mutations, typically span a large range of $Δ E$ s. For example, the mean and median values of $Δ E$ for L90M in those sequences which contain exactly eight PI-associated mutations is slightly smaller than zero (meaning a probability of reverting slightly greater than one), the range of $Δ E$ values spans from about −6 to +5.4 (the preference to revert spans the range ∼217x to 0.003x). Figure 2A shows that for each box plot corresponding to nine or more PI-associated mutations, the mean and median value of $Δ E > 0$ , again for each box plot $Δ E$ or equivalently, the preference to revert, assumes a large range. It is noteworthy that in the consensus background (one of the sequences containing no PI-associated mutations), the predicted probability of observing the consensus residue L90 is ~40x greater than observing the drug-resistance mutation 90M in the Stanford database of drug-experienced sequences. This prediction for the relative likelihood of L90M in the consensus background is qualitatively consistent with what is actually observed in the drug-experienced Stanford dataset (40 sequences are observed to contain L90 and 2 sequences are observed to contain 90M). Figure 2 does not however, represent a time ordering for the acquisition of the primary DRM, instead it follows the general favorability and likelihood of the mutation in sequence backgrounds conditional on the number of inhibitor-associated mutations; and the mutation L90M on average, is favored over the wild-type residue and is entrenched (unlikely to revert) only in sequence backgrounds with 9 or more PI-associated mutations (seen in Figure 2A for sequences in red).

Figure 2

Download asset Open asset

Entrenchment and the effect of epistasis on the favorability of a primary resistance mutation.

(A) $Δ E^{90} = E_{w i l d - t y p e}^{90} - E_{m u t a t i o n}^{90}$ change for sequences, conditional on the number of PI-associated mutations is shown as boxplots annotating the first, second and third interquartile range. Whiskers extend to 1.5 times the interquartile range with outliers marked as ’x’s and the mean values are marked as squares. The left ordinate scale shows the relative probability of reversion $(e^{- Δ E})$ , and the right shows $Δ E$ . Sequences whose energy difference fall above $Δ E = 0$ (dashed line) are entrenching backgrounds favoring the mutation. Sequence backgrounds where the mutation is favored on average are shown in red, the others in blue. The mutation L90M becomes favorable on average when there are about 9 PI-associated mutations, but there is a wide range of favorability and ‘which’ PI-associated mutations are present play an important role in determining if the primary mutation is favored/disfavored. The highlighted regions (white with dark border) show there are many sequence backgrounds with between 7 and 14 mutations in which L90M is either ‘highly entrenched or favored’ (top) or ‘highly disfavored’ (bottom). (B) Distributions of number of PI-associated mutations in sequences in the ‘highly entrenching’ and ‘highly disfavoring’ regions from panel A have a large overlap, again showing that entrenchment is not primarily determined by number of associated mutations. Prediction of the likelihood of M based on the number of PI-associated mutations alone for these sequence backgrounds where M is highly entrenched (red) or highly disfavored (blue) would be especially difficult due to the large overlap between them.

https://doi.org/10.7554/eLife.50524.004

For a quantitative comparison of the Potts model predictions of ‘entrenchment’ with the observed statistics concerning the relative likelihood of L90M in the Stanford drug-experienced dataset, it is necessary to aggregate the sequence statistics since the dataset size is too small to compare predictions with observations on a sequence by sequence basis. We compare the observed and predicted residue frequencies for different subsets of sequences as shown in Table 1 for L90M in HIV protease. The agreement between the predictions of the Potts model and the observed sequence statistics is very good. For the set of highly entrenching sequences ( $Δ E^{i} > 1 σ$ with $σ$ being the standard deviation in $Δ E$ ) the L90M mutation is predicted to be present in this set with frequency 97.6% whereas it is observed to be present with frequency 97.3% (see Table 1 second row). In the backgrounds classified by the Potts model as ‘highly disfavoring’, the mutation is predicted to occur with frequency 3.5% whereas the observed frequency is 1.9%. The L90M mutation is observed to be enriched in the sequences classified by the Potts model as ‘highly entrenching’ over those classified as ‘highly disfavoring’ the mutation by a factor of ~1907.

Table 1

Confirmation of the Potts model entrenchment predictions for the mutation L90M in HIV-1 protease.

Analytical predictions of likelihoods of the mutation (using the Potts model) are shown along with the corresponding observed frequencies in different subsets of our dataset classified as entrenching ( $Δ E > 0$ ) or disfavoring ( $Δ E < 0$ ) for the mutation. The agreement between the predicted and observed frequencies is remarkable and serves as a confirmation for the Potts model entrenchment predictions.

https://doi.org/10.7554/eLife.50524.005

Sequences	Classification	Number of sequences	Predicted frequency	Observed frequency
All	Entrenched	1475	82.7%	83.9%
All	Disfavored	3283	12.7%	11.9%
All with $\| Δ E \| \geq 1 σ$ of $Δ E = 0$	Entrenched	560	97.6%	97.3%
All with $\| Δ E \| \geq 1 σ$ of $Δ E = 0$	Disfavored	1444	3.5%	1.9%
With between 7 and 14 mutations and $\| Δ E \| \geq 1 σ$ of $Δ E = 0$	Entrenched	470	97.5%	97.0%
	Disfavored	239	2.8%	2.1%

Entrenchment cannot be predicted based simply on the number of mutations, particularly for sequences with an intermediate number of associated mutations in the background. To illustrate this we consider sequences for L90M in HIV-1 PR, with between 7 and 14 inhibitor-associated mutations which correspond to the two highlighted regions (white with dark border) in the center of Figure 2A which contain roughly equal numbers of sequences (see Table 1 third row agreement between predictions and observations is excellent). There is significant overlap of the number of PI-associated mutations in the two distributions as seen in Figure 2B, which makes clear why it is difficult to classify them based on the number of associated mutations alone. This also serves to emphasize that the probability of observing a drug-resistance mutation depends on the specific pattern of mutations that appears in the background and not just the total number of associated mutations in the background; this is also implied from the variance of the Potts $Δ E$ (lengths of the box and whiskers in Figure 2A) for sequences, for instance, with 13 to 17 PI-associated mutations, where some sequence backgrounds with even such high numbers of associated mutations disfavor the primary DRM while others entrench it.

Comparative entrenchment of resistance mutations in protease, reverse transcriptase, and integrase

Entrenchment of a resistance mutation depends crucially on the specific background in which it is acquired. In this section, we compare the effect of ‘background dependence’ of primary DRMs in each of the three major HIV drug target proteins, reverse transcriptase, protease and integrase for the four drug classes: nucleoside analog reverse transcriptase inhibitors (NRTIs), non-nucleoside analog reverse transcriptase inhibitors (NNRTIs), protease inhibitors (PIs) and integrase strand transfer inhibitors (INSTIs) (Table 2). Due to the large number of primary resistance mutations, we restrict our analysis to only those occurring at ~1% mutation frequencies or more. A mutation is defined to be ‘entrenched in the population’ if at least 50% of the sequences which contain that mutation have a change in Potts statistical energy score $Δ E > 0$ .

Table 2

Entrenchment in the population (of sequences carrying the mutation) for four major classes of resistance mutations in HIV-1 subtype B: A DRM is entrenched in the population if at least ~50% or more of the sequences containing the mutation entrench it (i.e $Δ E = E_{w i l d} - E m u t > 0$ ).

We catalog entrenchment in the population (of sequences carrying the mutation) for all primary DRMs appearing at ~1% frequency or more, and our study reveals mutations in response to NNRTIs are much less entrenched in the population (of sequences carrying the mutation) than others.

https://doi.org/10.7554/eLife.50524.006

	NRTIs	NNRTIs	PIs	INSTIs
Years in therapy	32	23	24	12
Number of primary DRMs	18	15	13	15
DRMs entrenched in the population (of sequences carrying the mutation)	11 (61.1%)	1 ( 6.7%)	10 (76.9%)	13 (86.7 %)
Number of DRMs conferring high-level resistance	7	11	7	6
High-level resistance DRMs entrenched in the population (of sequences carrying the mutation)	3 (42.8%)	1 (9.1 %)	5 (71.4%)	5 (83.3 %)

Source: Results shown in this table are based on calculations for ‘entrenchment in the population’ (of sequences carrying the mutation) for each primary DRM occurring at ~1% frequency or more in HIV-1 RT, PR, and IN, respectively as shown in Table 2—source data 1, Table 2—source data 2, Table 2—source data 3 and Table 2—source data 4. Resistance levels are determined according to Stanford HIVDB (Rhee et al., 2003; Shafer, 2006) mutation scores for PIs, NRTIs, NNRTIs, and INSTIs.

Table 2—source data 1 Table showing entrenchment in the population (of sequences carrying the mutation) for primary resistance mutations against NRTIs.: https://doi.org/10.7554/eLife.50524.007
Download elife-50524-table2-data1-v1.docx
Table 2—source data 2 Table showing entrenchment in the population (of sequences carrying the mutation) for primary resistance mutations against NNRTIs.: https://doi.org/10.7554/eLife.50524.008
Download elife-50524-table2-data2-v1.docx
Table 2—source data 3 Table showing entrenchment in the population (of sequences carrying the mutation) for primary resistance mutations against PIs.: https://doi.org/10.7554/eLife.50524.009
Download elife-50524-table2-data3-v1.docx
Table 2—source data 4 Table showing entrenchment in the population (of sequences carrying the mutation) for primary resistance mutations against INSTIs.: https://doi.org/10.7554/eLife.50524.010
Download elife-50524-table2-data4-v1.docx

If a resistance mutation is ‘entrenched in the population’, it means that the majority of sequences containing that mutation are unlikely to revert. These sequences would then support the development of drug resistance in the population as it evolves under drug selection pressure. Most drug-resistance mutations against all four classes of drugs are usually present only at low frequencies (much lower than 50% with the one exception of M184V in RT) in the drug-experienced dataset (Table 2—source datas 1, 2, 3, 4); but many of these mutations (with the exception of those conferring resistance to the NNRTIs) are highly ‘entrenched in the population’ of sequences where they do occur (Table 2), meaning that the mutation is more likely than the wild-type residue in these particular backgrounds, with low probability of reversion. The INSTI-selected resistance mutation Q148R in integrase, for example, is present in only about 5% of the patient population (in our dataset), but the mutation is entrenched in more than 85% of these sequences which have the mutation. These entrenching sequences, even though they are only a small fraction (~4.3%) of the total population, can then support the persistence of drug resistance associated with that mutation as the population evolves under drug selection pressure.

Over the last decade, viral IN has emerged as the major new focus of antiretroviral therapy. IN has been subjected to drug pressure for a considerably shorter period of time than other HIV target proteins, yet IN drug-resistance mutations are strongly ‘entrenched in the population’ of sequences containing the mutation (Table 2) similar to DRMs against PIs, and NRTIs which have been used in therapy for much longer. Figure 3 illustrates the degree (strength or value of $Δ E$ , and frequency or number of times $Δ E > 0$ ) of entrenchment for some key resistance mutations occurring in the catalytic core domain (CCD) of IN. The mutations Q148R and Q148H for example, reduce susceptibilities to the IN strand transfer inhibitors (INSTIs) raltegravir (RAL) and elvitegravir (EVG) by 40–100 fold and 5–10 fold, respectively (Goethals et al., 2010; Fransen et al., 2009; Abram et al., 2013; Goethals et al., 2008). Together, they are present in ~20% of the treated population, but they are strongly entrenched in the sequences in which they occur, with the probability of reversion being much smaller than 1%. Figure 3 shows that the strength of entrenchment can be high for many of these mutations. For example, mutations at site 148 are entrenched in almost all sequences with this mutation, ensuring the persistence of drug resistance associated with these mutations. Comparatively, resistance mutations in RT associated with the non-nucleoside analog reverse transcriptase inhibitors (NNRTIs) are much less ‘entrenched in the population’ of sequences containing the mutation (Table 2). Only 1 (K103N/T) of the 15 mutations rendering resistance to NNRTIs shows any significant ‘entrenchment’ by the background.

Figure 3

Download asset Open asset

Degree of entrenchment of key resistance mutations occurring in the catalytic core domain (CCD) of HIV-1 IN.

The change in Potts statistical energy $Δ E$ for some of the key resistance mutations occurring in the catalytic core domain (CCD) of IN, is plotted as a function of the rank of mutation-carrying sequences, ranked in descending order of their favorability towards the mutant. Plot shows the degree of entrenchment for these mutations. For example, Q148R is highly entrenched in almost all sequences carrying it, whereas, G163K is entrenched ( $Δ E > 0$ ) in only about half of the sequences carrying G163K.

https://doi.org/10.7554/eLife.50524.011

The striking observation that NNRTI-resistance mutations are less entrenched than those for NRTIs, PIs, and INSTIs is noteworthy (Table 2). One clear difference is that NNRTIs are the only allosteric inhibitors among the four target classes being examined here. NRTIs, PIs, and INSTIs exert their action by binding to the catalytic active sites of the target enzymes. One characteristic of active sites is that they contain invariant (or at least highly conserved) residues, and it is possible that nearby mutations interact strongly with those and their potential entrenchment may be influenced by the presence of the conserved residues. Another shared characteristic of the PI, RT, and IN active site is their catalytic residues are all acidic (PR: 2 Asp; RT 3 Asp; IN 1 Glu, 2 Asp). The NNRTI-binding site has considerably more sequence variation, with only Trp229 being invariant. Perhaps the overall variability of sequence in this region makes it more difficult for residues to become entrenched. Another key difference is that the variety of conformational states of the NNRTI-binding region is broader than for the active sites. The pocket can be closed, open, hyperextended open, and open with a shift of the primer grip, with varying degrees of solvent exposure of numerous internal hydrophobic residues (Sarafianos et al., 2009).

Entrenchment of a DRM makes it hard to escape drug resistance once the mutation is highly favored in its background. Table 2 shows the ‘entrenchment in the population’ of individual DRMs in sequences that carry the DRM. On aggregating over all possible DRMS, the fraction of sequences containing at least one primary DRM is greater than 60% for all four classes of DRMS (see Table 3 second row). When considering the fraction of sequences carrying at least one entrenched DRM, again the fraction is close to 50% or greater for three of the four DRM classes; the prominent exception is NNRTIs for which only 28.1% carry at least one entrenched primary DRM. This implies that the persistence of drug resistance against NNRTIs as a class (compared to the others) is less likely in the population than the other three classes of HIV drugs.

Table 3

Entrenchment of at least one primary DRM for each of the four drug classes in HIV-1: Table shows the percentage of drug-experienced sequences which contain at least one primary DRM, and the percentage of drug-experienced sequences which contain at least one primary DRM such that the DRM is entrenched by its respective background.

This is shown separately for resistance mutations occurring in response to each of the four HIV drug classes. A significantly lower percentage (~28%) of the patient sequences carry at least one entrenched resistance mutation conferring resistance to the NNRTIs when compared to other drugs (~50%–65%).

https://doi.org/10.7554/eLife.50524.012

	NRTIs	NNRTIs	PI	Insti
Number of sequences in MSA	19194	19194	4758	1220
% of total sequences containing at least one primary DRM	78.5%	62.4%	64.8%	61.7%
% of total sequences with at least one entrenched primary DRM	64.7%	28.1%	50.8%	47.2 %

Fitness and degree of entrenchment of observed resistance mutations: why some mutations are seen and others are not

The evolutionary trajectory of a resistance mutation depends on the balance between the fitness cost of evolving that mutation in its particular background, and the advantage it provides in evading drug pressure. It has been shown for HIV protease that a majority of the possible amino acid substitutions are observed rarely or not at all in isolates, indicating that the protein function is under strong purifying selection (Boucher et al., 2019). Observed resistance mutations arising with lower fitness costs which allow the virus to escape drug pressure with the least effect on its ability to propagate infection are therefore expected to become evolutionarily ‘entrenched’ by their respective sequence-backgrounds. Figure 4 shows the distribution of likelihoods as measured by the Potts $Δ E = E_{w i l d} - E_{m u t}$ scores at drug-resistance sites in sequences where the corresponding drug-resistance mutations are present in HIV integrase. The distribution of likelihoods for observed DRMs in sequences where they are present is seen in ‘green’, while the ‘blue’ distribution shows the distribution of likelihoods of mutations to other possible residue types in these sequences (containing the DRMs) at the sites where the DRM is observed. The distribution of the likelihoods of observed IN DRMs (green in Figure 4) shows that these drug-resistance mutations are generally more favorable than the wild-type residue as they have a $Δ E > 0$ in their respective backgrounds where the mutations are present. The relative displacement of the green from the blue distribution shows that in sequences which contain the drug-resistance mutation, the observed mutations are predicted on average to be more likely than the consensus residue type and that mutations to other residue types at those sites in the same sequence backgrounds are even less likely than the consensus residue type, on average. Hence, the observed resistance mutations arise in these drug-experienced backgrounds and help enable the virus to evade drug pressure. The more prohibitive fitness costs associated with the appearance of rare/unobserved mutations (blue distribution in Figure 4) at the same sites support the suggestion that drug-resistance mutations generally confer resistance at the least fitness cost. Each distribution also has a wide range of $Δ E$ scores illustrating the wide variation in the favorabilities of different mutations in different sequence backgrounds.

Figure 4

Download asset Open asset

Distribution of Potts $Δ E$ scores for key residues associated with drug resistance in HIV-1 IN.

The distribution of the Potts $Δ E (E_{w i l d} - E_{m u t})$ scores for sequences carrying the particular resistance mutation are shown in ’green’ for the most frequently observed INSTI selected resistance mutations in HIV IN, and in ’blue’ for all other possible mutations at the same sites. Other possible mutations include rarely observed or unobserved mutations. The histograms show the differential distribution of $Δ E$ scores for observed vs. unobserved/rare mutations at 15 primary mutation sites associated with evolving drug resistance in HIV-1 IN. The green (observed) and blue (unobserved/rare) distributions are normalized to the total number of primary DRMs in IN in the Stanford HIVDB and the total number of other possible mutations at the same sites, respectively. The mean $Δ E$ scores for observed vs. unobserved mutations are +2.11 and −5.58, respectively ( $p < 0.001$ ). The wide distribution of $Δ E$ scores also illustrates the role of the background in which the resistance mutation occurs.

https://doi.org/10.7554/eLife.50524.013

Molecular clones and the effects of specific backgrounds

In vitro fitness and drug susceptibility assays for viral proteins of HIV are often based on mutagenesis experiments, which are performed with specific molecular clones of the virus such as NL4-3, HXB2, LAI IIIB, etc (Hu and Kuritzkes, 2010; Abram et al., 2013). Mutagenesis experiments performed with specific clones, however, can only explore a limited region of the fitness and mutational landscapes available to the viral proteins accessible through the particular clones. These measurements are not necessarily representative of the fitness effects and the favorabilities of primary resistance mutations in the presence of in vivo selective pressures and entrenchment in viral background(s) of diverse patient populations as contained in the Stanford HIV drug resistance database (HIVDB, https://hivdb.stanford.edu/).

Figure 5 compares the predicted log likelihoods (Potts $Δ E$ s also associated with the degree of entrenchment) of primary drug-resistance mutations in integrase and reverse transcriptase averaged over the drug-experienced patient sequences in the Stanford HIVDB (denoted as $< Δ E_{p a t i e n t} >$ ) with the corresponding values in the specific backgrounds of molecular clones NL4-3 (panels A and D), HXB2 (panels B and E) and the subtype B consensus sequence (panels C and F), respectively. The likelihood and degree of entrenchment of a primary resistance mutation depends strongly on the specific background in which it is accrued and as observed in Figure 5, these effects studied in the genetic background of a specific molecular clone does not accurately (correlates only moderately) represent the effects averaged over drug-experienced patient populations. Among the three specific sequences analyzed, NL4-3 seems to be somewhat more representative of drug-experienced patient populations than the other two. That the consensus B sequence is not the most representative serves to highlight that epistatic effects can vary significantly from the average to patient populations. In fact, the molecular clones likely represent patient backgrounds where the mutation is disfavored and likely absent, rather than where the mutation is entrenched.

Figure 5

Download asset Open asset

Entrenchment and favorability of key resistance mutations in specific backgrounds.

The $Δ E$ change in Potts energy of a sequence is used as the measure of ‘entrenchment’ and favorability of key resistance mutations in HIV-1 NL4-3, HXB2, and the subtype-B consensus sequences, respectively, shown as a function of the average $Δ E$ in drug-experienced HIV-1 subtype B patient populations ( $< Δ E_{p a t i e n t} >$ ) in the Stanford HIVDB for viral integrase (**A,B,C**) and reverse transcriptase (**D,E,F**). In each case, the Pearson correlation coefficients are indicated. The protein sequences for the molecular clones NL4-3 and HXB2 are obtained from GenBank with accession number AF324493.2 and K03455.1, respectively with protein ids for the pol polyprotein as AAK08484.2 and AAB50259.1, respectively. The subtype B consensus sequence is obtained from the Stanford HIVDB. The degree of ‘entrenchment’ in these subtype B strains is often not representative of the average entrenchment effects in a patient population or even the most representative background from a patient population.

https://doi.org/10.7554/eLife.50524.014

Figure 6 shows the distribution of the change in Potts energy, $Δ E$ (log likelihood) values for the DRMS N155H (panel A) and G140S (panel B) in HIV IN in different drug-experienced sequence backgrounds in the Stanford HIVDB along with the $Δ E$ change in the specific backgrounds of NL4-3, HXB-2 and the subtype B consensus sequences. In terms of the change in Potts energy, which represents the degree of entrenchment of a mutation and its overall favorability in a given background, there is a clear distinction between sequence backgrounds where the mutation is present and mostly favored (or ‘entrenched’) from those where the mutation is absent. Key resistance mutations are typically disfavored in molecular clones in line with the fitness effects reported in Hu and Kuritzkes (2010) and in line with sequence backgrounds that do not carry the mutation (Figure 6). An exception to the general rule is observed for the IN mutation N155H in NL4-3 (Figure 6A), where the virus, despite not carrying the mutation, is favorable towards it.

Figure 6

Download asset Open asset

The particular sequence background in which a resistance mutation occurs affects the degree of entrenchment, with often a clear distinction between sequences where the mutation is present (green) versus absent (blue) shown here for the mutations N155H (A) and G140S (B) in integrase.

The degree of entrenchment for subtype B consensus, NL4-3 and HXB2 are shown as ’black dashed’, ’red’ and ’magenta’ lines, respectively. The Potts entrenchment score manifests as a clear distinction between backgrounds where a particular mutation is observed from ones where the mutation is absent, with the former more likely to present a distinct fitness advantage towards the mutation.

https://doi.org/10.7554/eLife.50524.015

Effective epistasis in the presence of drug selection pressure

HIV evolves rapidly with studies indicating that in the absence of drug pressure, the virus in a single patient explores the majority of all point mutations many times daily (Coffin, 1995; Perelson et al., 1996). The presence of drug selection pressure can however bias the virus to explore regions of the mutational landscape that allow the virus to evade the drug, which are not accessible otherwise. In this section, we discuss how drug treatment leads to a mutational landscape that reflects both intrinsic epistatic effects and the effects of drug selection pressure; we use the term ‘effective epistasis’ to refer to the combined effects of both.

As HIV jumps from host to host, it adapts to the applied selection pressure from each host’s immune response and drug regimen if the host is currently undergoing retroviral treatment. In the case of immune response selection pressure, (Shekhar et al., 2013; Ferguson et al., 2013) suggest that due to the diversity of host immune responses among the HIV population, selective effects of immune pressure are averaged so that a Potts model fit to sequences from many different hosts can effectively capture the ‘intrinsic’ fitness landscape of the virus. In the case of drug selection pressure, the correlated mutation patterns which are contained in MSAs built on drug-experienced datasets reflect a combination of intrinsic epistatic effects and possibly correlations induced by drug selection pressure. To distinguish between intrinsic epistasis and correlated drug selection pressure, we have built a second Potts statistical energy model on an MSA constructed from the Stanford HIV drug-naive dataset and compared the two models.

Figure 7 shows that the effects of point mutations captured by the changes in Potts statistical energies, $Δ E = E_{w i l d - t y p e} - E_{m u t a t i o n}$ when scored using the two models (the drug-naive model fit to an MSA constructed using HIV sequences in the Stanford database without prior drug exposure and the drug-experienced model fit to an MSA constructed with drug-experienced sequences) are highly correlated. The changes in the ‘field’ terms are responsible for the shift in $Δ E$ s, whereas the changes in the ‘couplings’ determine the correlation coefficient. Overall, the high correlation (Pearson correlation coefficient > 0.8) between the probabilities of observing the drug-resistance mutation relative to the wild-type in a given background (measured as a function of $Δ E$ ), when the Potts model is parameterized on the drug-naive dataset as compared with the drug-experienced dataset shows that ‘intrinsic’ epistatic effects have a large influence on the virus evolving under drug selection pressure. A detailed comparison of the two models will be presented elsewhere.

Figure 7

Download asset Open asset

Comparison of Potts models parameterized on drug-naive and drug-experienced HIV protein sequences.

The comparison of the effects of point mutations is shown in terms of Potts $Δ E$ scores (which forms the basis of our study) using two different Potts models, parameterized on drug-naive vs drug-experienced sequences for PR drug-resistance mutations M46I (A), I54V (B), V82A (C), and L90M (D), all of which appear with at least a frequency of 0.25% in both datasets. Bin shading for the 2D histogram scatter plots shown here scales logarithmically with the number of sequences whose scores fall into each bin. To obtain $Δ E$ scores, the sequences are scored using a drug-naive model vs. a drug-experienced model. The $Δ E$ scores are highly correlated with a Pearson correlation coefficient of 0.82 ( $p < 0.001$ ), 0.93 ( $p < 0.001$ ), 0.85 ( $p < 0.001$ ), and 0.82 ( $p < 0.001$ ), respectively.

https://doi.org/10.7554/eLife.50524.016

As a consequence of fitting a Potts model on drug-experienced sequences, there may be apparent mutational correlations induced due to the varying selective pressures of different drug regimens. However, the large overlap in resistance mutation sites among different drugs of the same class (Wensing et al., 2017) as seen in Figure 8 as well as the large number of possible resistance mutations that can arise in response to any one drug, and the multi-drug resistance mutation patterns induced by many drugs (Condra et al., 1995; Hertogs et al., 2000; Delaugerre et al., 2001), all tend to minimize the effects of such ‘spurious’ correlations in the drug-experienced dataset. Furthermore, if such ‘spurious’ correlations picked up from a drug-experienced dataset strongly affected the drug-experienced Potts model, we would not expect to see the ‘effective’ epistatic effects on $Δ E$ inferred by the drug-experienced model to be highly correlated with the ‘intrinsic’ epistatic effects of the drug-naive model as observed in Figure 7.

Figure 8

Download asset Open asset

Drug-pressure associated mutations are largely common between drugs of the same class.

The mutations (both primary and associated) arising in drug-treated HIV proteases in response to inhibitor treatment are shown corresponding to each protease drug. The diameters of the circles represent the number of mutations at that site that occurred in sequences treated with the particular drug. Most mutations that occur in response to one drug are seen to have occurred when treated with another drug of the same class, showing that the ‘spurious correlations’ (that could be picked up by a Potts model built on a mixture of patient sequences treated with different drugs of the same class, if the mutations occurring in response to one drug are not at all observed in response to another), are minimal.

https://doi.org/10.7554/eLife.50524.017

Discussion

The evolution of viruses like HIV under drug and immune selective pressures induces correlated mutations due to constraints on the structural stability and fitness (ability to assemble, replicate, and propagate infection) of the virus (Theys et al., 2018). This is a manifestation of the epistatic interactions in the viral genome. In fact, it has been shown that long-range epistasis can shift a protein’s mutational tolerance during HIV evolution (Haddox et al., 2018). Epistasis has recently become a major focus in structural biology and genomics, and co-evolutionary information encoded in collections of protein sequences have been used to infer epistatic couplings between them (Hinkley et al., 2011; Ferguson et al., 2013; Mann et al., 2014; Hopf et al., 2017; Figliuzzi et al., 2016; Butler et al., 2016; Haldane et al., 2018). In the current study, we have used the correlated mutations encoded in an MSA of drug-experienced HIV-1 sequences to parameterize a Potts Hamiltonian model of sequence statistical energies, and used it to infer the epistatic interactions leading to drug resistance in HIV-1 Subtype B. We first confirmed the Potts model’s ability to accurately predict the likelihood of a mutation based on a Potts statistical energy analysis. The predictive power of the model is established by verifying predictions of mutation prevalence based on its background using aggregate sequence statistics from the MSA, illustrating the crucial role of the background in giving rise to a mutation. This also provides the most direct verification of the complex patterns of epistasis which lead to drug resistance in HIV.

Entrenchment of functionally important mutations has been shown to be a general feature in systems exhibiting epistasis (Pollock et al., 2012). The entrenchment of primary resistance mutations described in this work suggests that epistasis plays a major role in the evolution of HIV under drug selection pressure. Both primary and accessory drug-resistance mutations exhibit strong epistatic interactions and therefore, entrenchment is a likely mechanism by which drug-resistance mutations accumulate within the population and in the persistence of drug resistance. Building on previous work by Flynn et al. (2017) identifying sequence backgrounds which were predicted to strongly entrench or strongly disfavor particular primary resistance mutations in HIV-1 PR, we first confirmed the Potts model entrenchment predictions using aggregate statistics from the MSA, and then extended it to the two other major HIV drug target enzymes, RT and IN.

It has been shown by Barton et al. (2016b) and Chen and Kardar (2019) that large differences in escape times for HIV patients targeting the same epitope can be explained by differences in the background mutations contained in the sequence of the virus strains that infected these patients. A similar result can also be expected for the virus subjected to antiretroviral therapy. If the background mutations have strong disfavorable couplings with the primary escape mutation, the escape mutation will likely take much longer to establish itself in the population. Due to the nature of entrenchment, we expect this to be reflected in the entrenchment score of a primary mutation in a given background. Resistance mutations, particularly ones associated with high levels of drug resistance, can be expected to occur with a significant degree of entrenchment in most drug-experienced backgrounds containing the mutation. Therefore, in the context of patient populations, a study of the entrenchment of primary resistance mutations is very relevant. The most entrenched mutations are the ones which are at some local maxima in the fitness landscape. Accumulating correlated mutations as observed in Figure 2A, can unlock pathways to these local fitness maxima (Gupta and Adami, 2016). But of crucial importance is the background in which these mutations occur. The local maxima can be 100 times more favorable in particular backgrounds, and these highly resistant sequences pose a significant risk for transmission of drug resistance to new hosts. Entrenchment of resistance mutations can thus play a central role in the persistence of drug-resistant viruses, and many of the most strongly entrenched mutations have been shown to revert slowly in drug-naive patients with transmitted resistance or in drug-experienced patients after withdrawal of ART (Izopet et al., 2000; Yang et al., 2015; Gandhi et al., 2003; Borman et al., 1996).

In vitro studies of the acquisition of drug resistance are based on mutagenesis experiments usually performed on molecular clones such as NL4-3, HXB2, IIIB, etc. We show that the epistatic effects on major drug-resistance mutations in molecular clones can be very different from what is observed in patient populations. Thus, fitness studies in particular clones may not be representative of fitness in drug-experienced patient populations or even the consensus background of a patient population.

Our comparison of the statistical energy differences $Δ E$ obtained from two Potts models of HIV protease sequence co-variation, one constructed from an MSA of drug-naive sequences in the Stanford HIV database, and the other constructed from an MSA of drug-experienced sequences, shows that the effect of the sequence background on the likelihood of observing a drug-resistance mutation compared with a consensus residue is highly correlated between the drug-experienced and drug-naive datasets (Figure 7). This supports our conclusion that ‘intrinsic’ epistatic effects have a strong influence on the virus evolving under drug selection pressure. While beyond the scope of the present work, this observation invites further analysis and opens up new avenues of research designed to interrogate in a quantitative way how the application of drug selection pressure changes the intrinsic mutational landscape. In order to identify HIV sequences which are most responsible for conferring drug resistance, we need to know the probabilities of observing highly entrenching sequences in the drug-naive population as well as in the drug-experienced population. Multi-canonical reweighting techniques have been developed for this purpose (Tan et al., 2016). While the present work does not deal with the time development of drug resistance or identify evolutionary paths, our construction of Potts models from both the drug-naive and drug-experienced populations, sets the stage for kinetic Monte-Carlo studies of pathways by which drug resistance is acquired.

Overall, the analysis presented here provides a framework for further studies of the fitness of HIV proteins evolving under drug selection pressure. The Potts model parameterized on the drug-experienced dataset is a powerful classifier that can identify complex sequence patterns that highly favor (entrench) or disfavor each drug-resistance mutation. Elucidating the epistatic effects of key resistance mutations has the potential to impact future HIV therapies and their entrenchment may be an important factor reinforcing the emergence of drug-resistant strains. The degree of entrenchment of key resistance mutations can provide insights about the fitness of these mutations in the context of a patient HIV reservoir; and their differential degree of entrenchment in different viral sequence backgrounds has the potential to impact future drug design strategies for the next generation of c-ART based on the patient viral reservoir, and suggests that the future of these strategies may very well involve personalized medicine based on an analysis of each patient’s individual viral genomic background.

Materials and methods

In this section, we give a brief introduction to the Potts Hamiltonian model and the motivation behind the model. The Potts model is a probabilistic model of sequence co-variation built on the single and pairwise site amino-acid frequencies of a protein MSA, and aimed to describe the observation probabilities of the various states of the system (sequences in the MSA).

To build the inference model, the goal is to approximate the unknown empirical probability distribution $P (S)$ which best describes HIV-1 sequences $S$ of length $L$ , where each residue is encoded in an alphabet $Q$ , by a model probability distribution $P^{m} (S)$ as in Mora and Bialek (2011). We choose the ‘least biased’ or maximum entropy distribution as the model distribution. Such distributions maximizing the entropy have been previously derived in Mézard and Mora (2009); Weigt et al. (2009); Morcos et al. (2011); Ferguson et al. (2013); Barton et al. (2016a) with the constraint that the empirical univariate and bivariate marginal distributions are effectively preserved. We follow a derivation of the maximum entropy model using Lagrange multipliers as in Mora and Bialek (2011) and Ferguson et al. (2013). Our maximum entropy model takes the form of an exponential distribution given by:

E (S) = \sum_{i}^{L} h_{S_{i}}^{i} + \sum_{i < j}^{L (L - 1) / 2} J_{S_{i} S_{j}}^{i j}

P^{m} (S) \propto e x p (- E (S))

where the quantity $E (S)$ is the Potts Hamiltonian determining the statistical energy of a sequence $S$ of length $L$ , the model parameters $h_{S_{i}}^{i}$ called ‘fields’ represent the statistical energy of residue $S_{i}$ at position i in sequence $S$ , and $J_{S_{i} S_{j}}^{i j}$ are ‘couplings’ representing the energy contribution of a position pair $i, j$ . In this form, the Potts Hamiltonian consists of $L Q$ field parameters $h_{S_{i}}^{i}$ and $(\begin{matrix} L \\ 2 \end{matrix}) Q^{2}$ coupling parameters $J_{S_{i} S_{j}}^{i j}$ , and for the exponential distribution $P^{m} \propto e x p (- E)$ , negative fields and couplings signify favored amino acids. The change in Potts energy due to mutating a residue $α$ at position $i$ in $S$ to $β$ is then given by:

Δ E (S_{α \to β}^{i}) = E (S_{α}^{i}) - E (S_{β}^{i}) = h_{α}^{i} - h_{β}^{i} + \sum_{j \neq i}^{L} J_{α S_{j}}^{i j} - J_{β S_{j}}^{i j}

In this form, a $Δ E (S_{α \to β}^{i}) > 0$ implies that the residue $β$ is more favorable than residue $α$ at position $i$ for the given sequence $S$ . If $α$ represents the wild-type residue at $i$ and $β$ the mutant, then the mutant is favorable over the wild-type if $Δ E > 0$ for the change, and vice versa.

The methodology followed in this analysis is similar to the one followed in Flynn et al. (2017) for HIV-1 PR (for further details on derivation and description of the model parameters see their Materials and methods section as well as the SI). The sample size of the MSA plays a critical role in determining the quality and effectiveness of the model (Haldane and Levy, 2019) and we confirm that the models are fit using sufficient data with minimal overfitting.

Data processing and mutation classification

Request a detailed protocol

The method begins with the collection of sequence and patient as well as reference information from the Stanford University HIV Drug resistance database (HIVDB, https://hivdb.stanford.edu) (Rhee et al., 2003; Shafer, 2006). In this work, we have used the Stanford HIVDB genotype-rx (https: //hivdb.stanford.edu/pages/genotype-rx.html) to obtain protein sequences for HIV target proteins. Alternatively, downloadable sequence datasets are also available on the Stanford HIVDB ( https://hivdb.stanford.edu/pages/geno-rx-datasets.html). The filtering criteria we used are: HIV-1, subtype B and nonCRF, and drug-experienced (# of PI = 1–9 for PR, # of NRTI = 1–9 and # of NNRTI = 1–4 for RT, and # of INST = 1–3 for IN), removal of mixtures, and unambiguous amino acid sequences (amino acids are ‘-ACDEFGHIKLMNPQRSTVWY’). For RT, sequences with exposure to both NRTIs and NNRTIs were selected. An alternate search for RT sequences exposed to only NRTIs or NNRTIs would return a vastly smaller number of isolates (5398 and 80 isolates, respectively for sequences exposed to only NRTIs or only NNRTIs as compared to 22446 isolates for sequences exposed to both). Sequences with insertions (‘#’) and deletions (‘~’) are removed. MSA columns and rows with more than 1% gaps (‘.’) are removed. This resulted in a final MSA size of $N = 4758$ sequences of length $L = 93$ for PR, $N = 19194$ sequences of length $L = 188$ for RT, and $N = 1220$ sequences of length $L = 263$ for IN. (Sequence datasets now available in the Stanford HIVDB are further updated, last on February, 2019). To retain enough sequence coverage in the MSA, we removed residues: residues 1–38 and residue 227 onwards for RT (see Appendix 1—figure 1 of Appendix), and 264 onwards for IN. For this reason, some interesting DRMs like F227I/L/V/C, L234I, P236L or N348I (NNRTI affected) for RT are not amenable for our analysis. For drug-naive PR model, sequences are obtained from the Stanford HIVDB using the criterion # of PI = 0, and after filtering of inserts, deletes, rows or columns with too many ‘gap’ characters, etc results in an MSA of $N = 14969$ sequences of length $L = 93$ . The protein sequences for the molecular clones NL4-3 and HXB2 are obtained from GenBank (Clark et al., 2016) with accession number AF324493.2 and K03455.1, respectively. The protein ids for the pol polyprotein are AAK08484.2 and AAB50259.1, respectively. The subtype B consensus sequence is obtained from the Los Alamos HIV sequence database (Foley et al., 2018) consensus and ancestral sequence alignments (https://www.hiv.lanl.gov/content/sequence/HIV/CONSENSUS/Consensus.html, last updated August 2004), and is also available in Appendix 1 of the Stanford HIVDB release notes page (https://hivdb.stanford.edu/page/release-notes/). Sequences are given weights reciprocal to the number of sequences contributed by each patient such that if multiple sequences are obtained from a single patient, the effective number of sequences obtained from any single patient is still 1. Sequences from different patients are assumed to be independent. It has been shown that phylogenetic trees of drug-naive and drug-experienced HIV-1 patients exhibit star-like phylogenies (Keele et al., 2008; Gupta and Adami, 2016), and thus, phylogenetic corrections are not required. Potts models of different HIV-1 protein sequences under immune pressure have also previously been parameterized without phylogenetic corrections (Ferguson et al., 2013; Mann et al., 2014; Butler et al., 2016).

In this work, mutations in the HIV genome have been classified into three main classes: primary or signature drug-resistance mutations, accessory or secondary mutations, and polymorphic mutations. The subtype B consensus sequence, which is derived from an alignment of subtype B sequences maintained at the Los Alamos HIV Sequence Database (https://www.hiv.lanl.gov), and is a commonly used reference sequence to which new sequences are compared, is used as our reference sequence and is referred to as the ‘consensus wild-type’ throughout the text. In accordance with current literature Wensing et al. (2017) and the Stanford HIV Drug Resistance Database (Rhee et al., 2003; Shafer, 2006), any mutation that (i) affects in vitro drug-susceptibility, (ii) occurs commonly in persons experiencing virological failure, and (iii) occurs with fairly low extent of polymorphism (mutations occurring commonly as natural variants) among untreated individuals, is classified as a major or primary drug-resistance mutation. In contrast, mutations with little or no effect on drug susceptibility are classified as accessory or secondary. Such mutations may reduce drug susceptibility or increase replication fitness only in combination with a primary mutation. Polymorphic mutations are mutations occurring as natural variants typically in antiretroviral drug (ARV) naive patients. Polymorphic mutations may occur in the absence of drug pressure and usually have little effect on ARV susceptibility when they occur without other DRMs. In this regard, polymorphic mutations affecting ARV susceptibility in combination with other DRMs can be classified as accessory.

Alphabet reduction

Request a detailed protocol

It has been shown that a reduced grouping of alphabets based on statistical properties can still capture most of the information in the full 20 letter alphabet while decreasing the number of degrees of freedom and thereby, leading to a more efficient model inference (Barton et al., 2016a; Haldane et al., 2016). All the possible alphabet reductions from 21 amino acid characters (20+1 gap) to 20 amino acid characters at any site $i$ are enumerated for all pair of positions $i j (j \neq i)$ , by summing the bivariate marginals for each of the $(\begin{matrix} 21 \\ 2 \end{matrix})$ possible combinations of amino acid characters, and choosing the alphabet grouping which minimizes the root mean square difference (RMSD) in mutual information (MI):

M I_{R M S D} = \sqrt{\frac{1}{N} \sum_{i j} {(M I_{Q = 21}^{i j} - M I_{Q = 20}^{i j})}^{2}}

between all pairs of positions $i j (j \neq i)$ . The process is then iteratively carried out from 20 amino acid characters to 19, and so on until the desired reduction in amino acid characters is reached. Due to residue conservation at many sites in HIV-1, several previous studies have used a binary alphabet to extract meaningful information from HIV sequences (Wu et al., 2003; Ferguson et al., 2013; Flynn et al., 2015). A binary alphabet however, marginalizes the information at a site to only the wild-type and mutant residues, with the loss of some informative distinctions between residues particularly at sites acquiring multiple mutations. To strike a balance between loss of information due alphabet reduction and reduction of the number of degrees of freedom for effective model inference, we use a reduced alphabet of 4 letters in line with Flynn et al. (2017). Our four letter alphabet reduction gives a Pearson’s $R^{2}$ coefficient of 0.995, 0.984, and 0.980 for protease, reverse transcriptase, and integrase, respectively between the $M I$ of bivariate marginal distributions with the full 21 letter alphabet and the reduced four letter alphabet, representing a minimal loss of information due to the reduction.

Using the reduced alphabet, the original MSA is then re-encoded and the bivariate marginals are recalculated. Small pseudocounts are then added to the bivariate marginals, as described by Haldane et al. (2016) to make up for sampling bias or limit divergence in the following inference procedure.

Model inference

Request a detailed protocol

The goal of the model inference is to find a suitable set of Potts parameters ${h, J}$ that fully determines the Potts Hamiltonian $E (S)$ and the total probability distribution $P^{m} (S)$ given in Equation 1 and Equation 2, respectively. This is done by obtaining the set of fields and couplings ${h, J}$ , which yield bivariate marginal estimates that best reproduce the empirical bivariate marginals.

A number of techniques to this effect have been developed previously (Mézard and Mora, 2009; Weigt et al., 2009; Balakrishnan et al., 2011; Cocco and Monasson, 2011; Morcos et al., 2011; Haq et al., 2012; Jones et al., 2012; Ekeberg et al., 2013; Ferguson et al., 2013; Barton et al., 2016a). We follow the methodology given in Ferguson et al. (2013). Given a set of fields and couplings, the bivariate marginals are estimated by generating sequences through a Markov Chain Monte Carlo (MCMC) sampling with the Metropolis criterion for a generated sequence proportional to the exponentiated Potts Hamiltonian. A multidimensional Newton search algorithm is then used to find the optimal set of Potts parameters ${h, J}$ . The descent step in the Newton search is determined after comparing the bivariate marginal estimates generated from the MCMC sample with the empirical bivariate marginal distribution. Although approximations are made in the computation of the Newton steps, the advantage of this method is that it avoids making explicit approximations to the model probability distribution. The method is limited by the sampling error of the input empirical marginal distributions and can also be computationally quite intensive. Our GPU implementation of the MCMC method makes it computationally tractable without resorting to more approximate inverse inference methods. The MCMC algorithm implemented on GPUs has been used to infer Potts models of sequence covariation which are sufficiently accurate to infer higher order marginals as shown by Flynn et al. (2017); Haldane et al. (2018). For a full description of the inference technique, we refer the reader to the supplemental information of Haldane et al. (2016).

The scheme for choosing the Newton update step is important. Ferguson et al. (2013) developed a quasi-Newton parameter update approach which determines the updates to $J^{i j}$ and $h^{i}$ by inverting the system’s Jacobian. In order to simplify and speed up the computation, we take advantage of the gauge invariance of the Potts Hamiltonian. We use a fieldless gauge in which $h^{i} = 0$ for all $i$ , and we compute the expected change in the model bivariate marginals $Δ f_{m}^{i j}$ (hereafter dropping the m subscript) due to a change in $J^{i j}$ to the first order by:

Δ f_{S_{i} S_{j}}^{i j} = \sum_{k l, S_{k} S_{l}} \frac{\partial f_{S_{i} S_{j}}^{i j}}{\partial J_{S_{k} S_{l}}^{k l}} Δ J_{S_{k} S_{l}}^{k l} + \sum_{k, S_{k}} \frac{\partial f_{S_{i}, S_{j}}^{i j}}{\partial h_{S_{k}}^{k}} Δ h_{S_{k}}^{k}

Computing the Jacobian $\frac{\partial f_{S_{i} S_{j}}^{i j}}{\partial J_{S_{k} S_{l}}^{k l}}$ and inverting the linear system in Equation 5 to solve for the changes in $Δ J^{i j}$ and $Δ h^{i}$ given the $Δ f^{i j}$ , is the challenging part of the computation. We choose the $Δ f^{i j}$ as:

Δ f^{i j} = γ (f_{e m p i r i c a l}^{i j} - f^{i j})

with a small enough damping parameter $γ$ such that the linear (and other) approximations are valid.

The computational cost of fitting $(\begin{matrix} L \\ 2 \end{matrix}) \times (4 - 1)^{2} + 93 \times (4 - 1) = 372, 816$ model parameters for the smallest protein in our analysis, PR, on 2 NVIDIA K80 or 4 NVIDIA TitanX GPUs is $\approx 20 â € ‰ h$ . The methodology followed in this analysis is almost the same as done in Flynn et al. (2017) (see Materials and methods) for HIV-1 PR. For a more detailed description of data preprocessing, model inference, and comparison with other methods, we refer the reader to the SI and text of Flynn et al. (2017); Haldane et al. (2016); Haldane et al. (2018); Haldane and Levy (2019).

A repository containing the inference methodology code is available at https://github.com/ComputationalBiophysicsCollaborative/IvoGPU and the final MSAs are available at https://github.com/ComputationalBiophysicsCollaborative/elife_data (copy archived at https://github.com/elifesciences-publications/elife_data). Appendix 1 contains a figure showing the sequence coverage for RT and tables of most entrenched and least entrenched sequences with same Hamming distances, for some important DRMs in PR, RT and IN.

Statistical robustness of HIV potts models

Request a detailed protocol

Finite sampling error and overfitting can play an important role in all inference problems, and inverse Ising inference is no exception. Overfitting of the Potts model cannot be simply estimated by comparing the number of model parameters to the number of sequences being fit. Our Potts model is not directly fit to the sequences but to the bivariate marginals (pairwise residue frequencies) of the MSA, and there are equal numbers of model parameters as there are bivariate marginals. The inference problem is neither over nor under constrained. The major source of error that affects our model inference is the fact that the inference relies on pair statistics (bivariate marginals) found in a finite collection of sequences. Overfitting comes from this finite-sampling error in the bivariate marginals estimated from a finite-sized multiple sequence alignment (MSA) which then serves as the input for model inference. It has been shown by Haldane and Levy (2019) that it is possible to fit accurate Potts models to MSAs when the number of sequences in the MSA is much smaller than the number of parameters of the Potts model.

The degree of Potts model overfitting can be quantified using the ‘signal to noise ratio’, or SNR, which is a function of the sequence length $L$ , alphabet size $q$ , number of sequences in the MSA $N$ , and a measure of the degree of conservation of the MSA, $χ^{2}$ . Using an SNR analysis of the kind described in Haldane and Levy (2019), we conclude that the Potts models for HIV-1 PR and RT are clearly not overfit (SNR values are 21.6 for PR, and 43.7 for RT). Since the number of sequences in the IN MSA is considerably smaller than the others, the IN model can be more affected by overfitting and may be less reliable than others (SNR is 0.14). However, the different types of predictions based on Potts models are differently affected by finite sampling error; predictions of the effect of point mutations to a sequence, $Δ E$ , which forms the basis of this study are the most robust and least affected by finite sampling. Using the in silico tests suggested in Haldane and Levy (2019), we find that the effects of point mutations are accurately captured even in the IN model, which is the most susceptible to finite sampling errors among our three models (for a more detailed analysis of the effects of finite sampling on the predictions of the Potts model, we refer the reader to Haldane and Levy, 2019). Thus, we conclude that the MSA sample sizes for PR, RT, and IN used in this study are sufficiently large to construct Potts models for these HIV proteins that adequately reflect the effects of the sequence background on point mutations which are the central focus of this work.

Appendix 1

Appendix 1—figure 1

Download asset Open asset

Sequence coverage for RT.

Figure shows the sequence coverage (# of sequences vs the # of residues) for RT drug-experienced (both NRTI and NNRTI) sequences derived from the Stanford HIVDB (22,444 isolates from 20422 patients). For RT, sequences with exposure to both NRTIs and NNRTIs were selected as an alternate search for RT sequences exposed to only NRTIs or NNRTIs would return a vastly smaller number of isolates (5398 and 80, respectively). Sequences with insertions (‘#‘) and deletions (‘ ~‘) are removed. MSA columns and rows with more than 1% gaps (‘.’) are removed. This resulted in a final MSA size of N = 19194 sequences from 17130 persons each with length L = 188 for RT. To retain enough sequence coverage in the MSA, we removed residues: residues 1–38, and residue 227 onwards for RT. For this reason, some interesting DRMs like F227I/L/V/C, L234I, P236L or N348I (NNRTI affected) for RT are not amenable for our analysis.

https://doi.org/10.7554/eLife.50524.020

Appendix 1—table 1

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against the INSTIs

https://doi.org/10.7554/eLife.50524.021

		E92Q		L74M		N155H		Q148H		G140S
Position	Consensus	ME	LE	ME	LE	ME	LE	ME	LE	ME	LE
6	D	--	--	--	E	--	--	--	--	--	--
11	E	--	--	--	--	--	--	--	D	--	D
17	S	--	--	N	N	N	N	--	T	--	--
20	R	--	--	--	--	--	--	--	--	K	--
22	M	--	--	I	--	--	--	--	--	--	--
23	A	--	--	V	--	--	--	--	--	--	--
28	L	--	--	I	--	--	--	--	--	--	--
31	V	--	--	--	--	I	--	I	--	I	--
32	V	--	--	--	--	--	--	--	--	--	I
37	V	--	--	--	--	--	--	--	--	--	I
39	S	C	C	--	--	--	--	--	--	--	--
45	L	Q	--	V	--	--	--	--	--	--	--
50	M	--	--	M	I	--	--	I	--	L	--
63	L	I	--	--	--	--	--	--	--	--	--
72	I	L	--	--	V	--	V	--	--	--	--
74	L	M	--	M	M	--	--	--	--	--	--
92	E	Q	Q	--	--	--	--	--	--	--	--
97	T	--	--	A	A	--	--	--	--	--	--
101	L	--	I	I	--	I	--	--	I	I	I
119	S	--	--	T	R	--	--	--	P	--	--
124	T	--	N	--	--	--	A	--	N	--	--
135	I	--	V	--	--	--	--	--	V	--	--
136	K	--	--	--	--	--	--	--	--	Q	--
138	E	--	--	D	--	--	--	K	--	--	--
140	G	--	--	--	--	--	S	S	--	S	S
143	Y	--	--	R	--	--	--	--	--	--	--
148	Q	--	--	--	--	--	H	H	H	H	R
151	V	--	--	--	I	I	--	--	--	--	--
154	M	--	--	--	--	--	I	--	--	--	--
155	N	--	H	--	H	H	H	--	--	--	--
156	K	--	--	--	N	--	--	--	--	--	--
170	E	--	--	--	--	--	A	--	--	--	--
181	F	--	L	--	--	--	--	--	L	--	--
188	K	--	--	--	--	--	--	R	--	--	--
196	A	--	--	--	--	P	--	--	--	--	--
201	V	--	I	--	--	I	I	I	I	--	I
206	T	--	--	--	--	--	--	S	--	--	--
208	I	--	--	--	M	--	--	--	--	L	--
212	E	--		--	A	--	--	--	--	--	--
215	K	--	--	--	--	--	--	S	--	--	--
216	Q	N	--	--	--	--	--	--	--	--	--
220	I	--	--	--	--	--	--	--	L	--	--
230	S	--	--	--	--	--	--	G	--	--	--
232	D	--	--	--	--	N	--	--	--	--	E
234	L	I	--	--	--	--	--	--	--	--	--
256	D	--	--	E	--	E	--	--	--	E	E
HD		8	8	12	12	9	9	10	10	9	9

Residues same as consensus are shown as ‘--', mutations are shown as the one letter abbreviated alphabet encoding the mutant residue (in bold are primary drug-resistance mutations appearing at more than 1% frequency).

Appendix 1—table 2

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against NNRTIs

https://doi.org/10.7554/eLife.50524.022

		P225H		Y188L		L100I		K103N/E		Y181C/G		E138K/R
	Consensus	ME	LE	ME	LE	ME	LE	ME	LE	ME	LE	ME	LE
39	T	--	--	---	--	--	S	--	--	--	--	--	--
41	M	--	--	L	L	L	L	--	--	--	--	--	L
43	K	--	--	--	--	--	R	--	--	--	--	--	--
44	E	--	--	--	D	--	--	--	--	--	--	--	--
48	S	--	--	--	--	--	--	T	--	--	--	--	--
49	K	--	--	--	--	--	--	--	R	--	--	--	--
50	I	--	--	--	--	--	--	--	--	--	--	--	--
58	T	--	--	N	--	--	--	--	--	--	--	--	--
60	V	--	--	--	--	--	--	--	--	--	I	--	--
62	A	--	--	V	--	--	--	--	--	V	--	--	--
65	K	--	--	--	--	--	--	--	--	R	--	--	--
67	D	--	--	--	G	N	--	--	--	--	--	--	--
68	S	--	--	--	G	--	--	--	--	G	--	--	--
69	T	--	--	N	--	N	--	--	--	I	N	--	--
70	K	--	--	--	T	--	--	--	--	--	R	--	--
74	L	--	V	--	V	V	--	V	--	--	--	--	V
75	V	--	--	I	--	T	--	--	--	A	--	--	--
77	F	--	--	L	--	--	--	--	--	--	--	--	--
82	K	--	--	--	--	--	--	--	--	--	--	R	--
90	V	--	--	--	--	--	--	I	--	I	--	--	--
98	A	S	G	--	--	--	S	--	S	--	--	--	--
100	L	--	--	--	--	I	I	I	--	--	--	--	--
101	K	Q	--	--	--	--	--	--	E	E	--	E	--
103	K	N	--	--	N	N	--	N	N	--	N	R	--
104	K	--	--	--	--	--	--	--	--	N	--	--	--
106	V	--	--	I	--	--	--	--	A	--	--	--	--
108	V	--	--	--	--	--	--	--	--	I	--	--	--
116	F	--	--	W	--	--	--	--	--	--	--	--	--
118	V	--	--	--	I	--	--	--	--	--	--	--	--
122	K	--	--	--	E	E	--	E	--	--	--	E	E
123	D	--	--	--	--	--	E	--	--	--	--	--	--
126	K	--	--	--	--	--	--	--	--	--	--	R	--
135	I	L	--	--	K	M	V	--	--	--	T	--	--
138	E	--	--	--	--	--	--	--	K	--	--	K	K
139	T	--	--	--	--	--	--	--	--	--	--	E	--
142	I	--	--	--	--	--	Q	--	--	--	V	T	--
151	Q	--	--	M	--	--	--	--	--	--	--	--	--
158	A	--	S	--	--	--	--	--	--	--	--	--	S
162	S	--	--	--	--	--	C	--	--	--	--	--	A
165	T	--	I	--	--	--	K	--	--	--	--	--	--
166	K	--	--	--	--	--	--	--	Q	--	--	--	--
169	E	--	--	--	A	--	--	--	--	--	--	--	--
172	R	--	--	K	--	--	--	--	--	--	--	--	--
173	K	--	--	N	--	--	--	--	--	--	--	--	--
176	P	Q	--	--	--	--	--	--	--	--	--	--	--
177	D	--	--	--	--	--	--	N	--	--	--	--	--
178	I	--	--	--	--	--	--	M	--	--	--	M	--
179	V	--	D	--	I	--	D	--	--	--	--	--	--
180	I	--	--	--	--	--	--	--	--	--	M	--	--
181	Y	--	C	--	C	--	--	--	--	C	C	--	--
184	M	V	--	--	--	V	I	I	--	--	V	--	--
188	Y	--	--	L	L	--	--	--	--	--	L	--	--
189	V	--	--	I	--	--	--	--	--	--	--	--	--
190	G	--	A	--	S	--	--	--	A	A	--	--	--
196	G	E	--	--	--	--	--	--	--	--	--	--	--
197	Q	--	--	--	--	--	--	--	K	--	--	--	--
200	T	A	V	A	--	A	A	A	A	--	A	--	A
202	I	--	--	V	--	V	--	--	--	--	--	--	--
203	E	---	--	--	--	D	--	--	--	--	--	--	--
207	Q	--	--	--	--	--	A	--	--	--	--	--	--
210	L	--	--	--	W	W	--	--	--	--	--	--	--
211	R	--	--	--	K	--	K	--	--	K	K	--	K
214	F	--	--	L	--	--	--	--	--	--	L	--	--
215	T	--	--	Y	Y	Y	D	--	--	--	--	--	Y
219	K	--	--	--	--	E	--	--	--	Q	Q	--	--
221	H	--	--	--	--	--	--	--	--	Y	--	--	--
223	K	--	--	T	--	--	--	--	--	--	--	--	--
225	P	H	H	--	--	--	--	H	--	--	--	--	--
HD		9	9	18	18	16	16	11	11	14	14	9	9

Residues same as consensus are shown as ‘--', mutations are shown as the one letter abbreviated alphabet encoding the mutant residue (in bold are primary drug-resistance mutations appearing at more than 1% frequency).

Appendix 1—table 3

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against PIs

https://doi.org/10.7554/eLife.50524.023

		I50V		I84V		L90M		V82A
Position	Consensus	ME	LE	ME	LE	ME	LE	ME	LE
10	L	I	--	F	--	I	--	I	I
13	I	--	V	--	--	V	--	--	V
15	I	--	V	--	V	--	V	--	--
18	Q	--	--	--	--	H	--	--	H
19	L	--	I	--	--	--	--	--	--
20	K	--	--	--	--	I	R	--	--
24	L	--	--	--	--	--	--	I	--
30	D	--	--	--	N	--	--	--	--
32	V	--	--	--	--	--	I	--	--
33	L	--	F	--	F	--	--	F	--
34	E	Q	--	--	--	--	--	--	--
35	E	--	D	--	D	--	D	--	--
36	M	--	--	--	L	--	I	--	--
37	N	--	S	--	E	--	D	--	--
41	R	--	K	--	--	--	K	--	--
43	K	--	--	--	T	--	--	--	--
46	M	I	--	I	--	--	L	L	I
47	I	--	--	--	--	--	A	--	--
48	G	V	--	--	--	--	--	--	--
50	I	V	V	--	--	--	--	--	--
53	F	Y	--	--	--	--	--	L	--
54	I	S	--	V	--	V	--	V	--
55	K	--	--	R	--	--	--	--	--
57	R	--	--	--	--	--	K	--	--
58	Q	--	E	--	--	--	--	--	--
60	D	--	--	--	E	--	--	--	--
61	Q	--	--	--	--	--	Y	--	--
62	I	V	--	--	V	V	--	--	--
63	L	Q	P	P	P	P	--	P	P
64	I	--	--	V	--	--	--	--	--
66	I	--	--	V	--	--	--	--	--
67	C	--	--	--	--	F	--	--	--
71	A	V	T	--	--	V	--	V	--
72	I	V	V	--	--	K	--	--	--
73	G	--	--	T	--	S	--	--	--
74	T	S	--	--	--	--	--	--	--
76	L	--	--	--	--	--	--	--	V
77	V	I	--	--	--	--	I	--	--
79	P	--	--	A	--	A	--	--	--
82	V	A	I	--	--	--	I	A	A
84	I	--	--	V	V	V	--	--	V
88	N	--	--	--	D	--	--	--	--
89	L	--	V	--	--	--	--	--	I
90	L	--	--	M	M	M	M	--	--
93	I	L	L	L	--	--	--	--	--
95	C	--	--	F	--	--	--	--	--
HD		15	15	13	13	14	14	9	9

Residues same as consensus are shown as ‘--', mutations are shown as the one letter abbreviated alphabet encoding the mutant residue (in bold are primary drug-resistance mutations appearing at more than 1% frequency).

https://doi.org/10.7554/eLife.50524.019

Data availability

Sequence data analyzed in this study is obtained from the Stanford University HIV drug resistance database (https://hivdb.stanford.edu/), Los Alamos HIV sequence database (https://www.hiv.lanl.gov/content/sequence/HIV/mainpage.html). Source data tables are provided for Table 2.

The following previously published data sets were used

1. Rhee S-Y
2. Gonzales MJ
3. Kantor R
4. Betts BJ
5. Ravela J
6. Shafer RW
(2003) Stanford HIV drug resistance database
Stanford University HIV drug resistance database: Genotype-Treatment Correlations.

https://hivdb.stanford.edu/pages/geno-rx-datasets.html
1. Foley B
2. Leitner T
3. Apetrei C
4. Hahn B
5. Mizrachi I
6. Mullins J
7. Rambaut A
8. Wolinsky S
9. Korber B
(2004) Los Alamos HIV sequence database
Consensus and Ancestral Sequence Alignments, Select 'Alignment type:Consensus/Ancestral', 'organism: HIV-1/SIVcpz', 'Pre-defined region of the genome: POL', Subtype:All', 'DNA/PRotein: Protein'.

https://hiv.lanl.gov/content/sequence/NEWALIGN/align.html

References

1. Abram ME
2. Hluhanich RM
3. Goodman DD
4. Andreatta KN
5. Margot NA
6. Ye L
7. Niedziela-Majka A
8. Barnes TL
9. Novikov N
10. Chen X
11. Svarovskaia ES
12. McColl DJ
13. White KL
14. Miller MD
(2013) Impact of primary elvitegravir Resistance-Associated mutations in HIV-1 integrase on drug susceptibility and viral replication fitness
Antimicrobial Agents and Chemotherapy 57:2654–2663.
https://doi.org/10.1128/AAC.02568-12
- Google Scholar
1. An DS
2. Morizono K
3. Li QX
4. Mao SH
5. Lu S
6. Chen IS
(1999)
An inducible human immunodeficiency virus type 1 (hiv-1) vector which effectively suppresses hiv-1 replication

Journal of Virology 73:7671–7677.
- PubMed
- Google Scholar
(2011) Learning generative models for protein fold families
Proteins: Structure, Function, and Bioinformatics 79:1061–1078.
https://doi.org/10.1002/prot.22934
- PubMed
- Google Scholar
(2016a) ACE: adaptive cluster expansion for maximum entropy graphical model inference
Bioinformatics 32:3089–3097.
https://doi.org/10.1093/bioinformatics/btw328
- PubMed
- Google Scholar
(2016b) Relative rate and location of intra-host HIV evolution to evade cellular immunity are predictable
Nature Communications 7:11660.
https://doi.org/10.1038/ncomms11660
- PubMed
- Google Scholar
(1996) Resistance of human immunodeficiency virus type 1 to protease inhibitors: selection of resistance mutations in the presence and absence of the drug
Journal of General Virology 77:419–426.
https://doi.org/10.1099/0022-1317-77-3-419
- PubMed
- Google Scholar
1. Boucher JI
2. Whitfield TW
3. Dauphin A
4. Nachum G
5. Hollins C
6. Zeldovich KB
7. Swanstrom R
8. Schiffer CA
9. Luban J
10. Bolon DNA
(2019) Constrained mutational sampling of amino acids in HIV-1 protease evolution
Molecular Biology and Evolution 36:798–810.
https://doi.org/10.1093/molbev/msz022
- PubMed
- Google Scholar
(2016) Identification of drug resistance mutations in HIV from constraints on natural evolution
Physical Review E 93:022412.
https://doi.org/10.1103/PhysRevE.93.022412
- PubMed
- Google Scholar
1. Chang MW
2. Torbett BE
(2011) Accessory mutations maintain stability in drug-resistant HIV-1 protease
Journal of Molecular Biology 410:756–760.
https://doi.org/10.1016/j.jmb.2011.03.038
- PubMed
- Google Scholar
Preprint
1. Chen H
2. Kardar M
(2019) Mean-field computational approach to hiv dynamics on a fitness landscape
bioRxiv.
https://doi.org/10.1101/518704
- Google Scholar
(2016) GenBank
Nucleic Acids Research 44:D67–D72.
https://doi.org/10.1093/nar/gkv1276
- PubMed
- Google Scholar
1. Cocco S
2. Monasson R
(2011) Adaptive cluster expansion for inferring boltzmann machines with noisy data
Physical Review Letters 106:090601.
https://doi.org/10.1103/PhysRevLett.106.090601
- PubMed
- Google Scholar
1. Coffin JM
(1995) HIV population dynamics in vivo: implications for genetic variation, pathogenesis, and therapy
Science 267:483–489.
https://doi.org/10.1126/science.7824947
- PubMed
- Google Scholar
1. Condra JH
2. Schleif WA
3. Blahy OM
4. Gabryelski LJ
5. Graham DJ
6. Quintero JC
7. Rhodes A
8. Robbins HL
9. Roth E
10. Shivaprakash M
(1995) In vivo emergence of HIV-1 variants resistant to multiple protease inhibitors
Nature 374:569–571.
https://doi.org/10.1038/374569a0
- PubMed
- Google Scholar
(2010) Fitness epistasis and constraints on adaptation in a human immunodeficiency virus type 1 protein region
Genetics 185:293–303.
https://doi.org/10.1534/genetics.109.112458
- PubMed
- Google Scholar
1. Delaugerre C
2. Rohban R
3. Simon A
4. Mouroux M
5. Tricot C
6. Agher R
7. Huraux JM
8. Katlama C
9. Calvez V
(2001) Resistance profile and cross-resistance of HIV-1 among patients failing a non-nucleoside reverse transcriptase inhibitor-containing regimen
Journal of Medical Virology 65:445–448.
https://doi.org/10.1002/jmv.2055
- PubMed
- Google Scholar
1. Ekeberg M
2. Lövkvist C
3. Lan Y
4. Weigt M
5. Aurell E
(2013) Improved contact prediction in proteins: using pseudolikelihoods to infer potts models
Physical Review E 87:012707.
https://doi.org/10.1103/PhysRevE.87.012707
- Google Scholar
(2013) Translating HIV sequences into quantitative fitness landscapes predicts viral vulnerabilities for rational immunogen design
Immunity 38:606–617.
https://doi.org/10.1016/j.immuni.2012.11.022
- Google Scholar
(2016) Coevolutionary landscape inference and the Context-Dependence of mutations in Beta-Lactamase TEM-1
Molecular Biology and Evolution 33:268–280.
https://doi.org/10.1093/molbev/msv211
- PubMed
- Google Scholar
1. Flynn WF
2. Chang MW
3. Tan Z
4. Oliveira G
5. Yuan J
6. Okulicz JF
7. Torbett BE
8. Levy RM
(2015) Deep sequencing of protease inhibitor resistant HIV patient isolates reveals patterns of correlated mutations in gag and protease
PLOS Computational Biology 11:e1004249.
https://doi.org/10.1371/journal.pcbi.1004249
- PubMed
- Google Scholar
1. Flynn WF
2. Haldane A
3. Torbett BE
4. Levy RM
(2017) Inference of epistatic effects leading to entrenchment and drug resistance in HIV-1 protease
Molecular Biology and Evolution 34:1291–1306.
https://doi.org/10.1093/molbev/msx095
- PubMed
- Google Scholar
Report
1. Foley B
2. Leitner T
3. Apetrei C
4. Hahn B
5. Mizrachi I
6. Mullins J
7. Rambaut A
8. Wolinsky S
9. Korber B
(2018) Hiv Sequence Compendium 2018
US Department of Energy.
https://doi.org/10.2172/1458915
- Google Scholar
1. Fransen S
2. Gupta S
3. Danovich R
4. Hazuda D
5. Miller M
6. Witmer M
7. Petropoulos CJ
8. Huang W
(2009) Loss of raltegravir susceptibility by human immunodeficiency virus type 1 is conferred via multiple nonoverlapping genetic pathways
Journal of Virology 83:11440–11446.
https://doi.org/10.1128/JVI.01168-09
- PubMed
- Google Scholar
1. Gandhi RT
2. Wurcel A
3. Rosenberg ES
4. Johnston MN
5. Hellmann N
6. Bates M
7. Hirsch MS
8. Walker BD
(2003) Progressive reversion of human immunodeficiency virus type 1 resistance mutations in vivo after transmission of a multiply drug-resistant virus
Clinical Infectious Diseases 37:1693–1698.
https://doi.org/10.1086/379773
- PubMed
- Google Scholar
(1994) Correlated mutations and residue contacts in proteins
Proteins: Structure, Function, and Genetics 18:309–317.
https://doi.org/10.1002/prot.340180402
- PubMed
- Google Scholar
1. Goethals O
2. Clayton R
3. Van Ginderen M
4. Vereycken I
5. Wagemans E
6. Geluykens P
7. Dockx K
8. Strijbos R
9. Smits V
10. Vos A
11. Meersseman G
12. Jochmans D
13. Vermeire K
14. Schols D
15. Hallenberger S
16. Hertogs K
(2008) Resistance mutations in human immunodeficiency virus type 1 integrase selected with elvitegravir confer reduced susceptibility to a wide range of integrase inhibitors
Journal of Virology 82:10366–10374.
https://doi.org/10.1128/JVI.00470-08
- PubMed
- Google Scholar
1. Goethals O
2. Vos A
3. Van Ginderen M
4. Geluykens P
5. Smits V
6. Schols D
7. Hertogs K
8. Clayton R
(2010) Primary mutations selected in vitro with raltegravir confer large fold changes in susceptibility to first-generation integrase inhibitors, but minor fold changes to inhibitors with second-generation resistance profiles
Virology 402:338–346.
https://doi.org/10.1016/j.virol.2010.03.034
- PubMed
- Google Scholar
1. Gupta R
2. Hill A
3. Sawyer AW
4. Pillay D
(2008) Emergence of drug resistance in HIV type 1-infected patients after receipt of first-line highly active antiretroviral therapy: a systematic review of clinical trials
Clinical Infectious Diseases 47:712–722.
https://doi.org/10.1086/590943
- PubMed
- Google Scholar
1. Gupta A
2. Adami C
(2016) Strong selection significantly increases epistatic interactions in the Long-Term evolution of a protein
PLOS Genetics 12:e1005960.
https://doi.org/10.1371/journal.pgen.1005960
- PubMed
- Google Scholar
(2018) Mapping mutational effects along the evolutionary landscape of HIV envelope
eLife 7:e34420.
https://doi.org/10.7554/eLife.34420
- PubMed
- Google Scholar
1. Haldane A
2. Flynn WF
3. He P
4. Vijayan RS
5. Levy RM
(2016) Structural propensities of kinase family proteins from a potts model of residue co-variation
Protein Science 25:1378–1384.
https://doi.org/10.1002/pro.2954
- PubMed
- Google Scholar
1. Haldane A
2. Flynn WF
3. He P
4. Levy RM
(2018) Coevolutionary landscape of kinase family proteins: sequence probabilities and functional motifs
Biophysical Journal 114:21–31.
https://doi.org/10.1016/j.bpj.2017.10.028
- PubMed
- Google Scholar
1. Haldane A
2. Levy RM
(2019) Influence of multiple-sequence-alignment depth on potts statistical models of protein covariation
Physical Review E 99:032405.
https://doi.org/10.1103/PhysRevE.99.032405
- PubMed
- Google Scholar
1. Haq O
2. Andrec M
3. Morozov AV
4. Levy RM
(2012) Correlated electrostatic mutations provide a reservoir of stability in HIV protease
PLOS Computational Biology 8:e1002675.
https://doi.org/10.1371/journal.pcbi.1002675
- PubMed
- Google Scholar
(2006) Origins of HIV and the evolution of resistance to AIDS
Science 313:462–466.
https://doi.org/10.1126/science.1123016
- PubMed
- Google Scholar
1. Hertogs K
2. Bloor S
3. Kemp SD
4. Van den Eynde C
5. Alcorn TM
6. Pauwels R
7. Van Houtte M
8. Staszewski S
9. Miller V
10. Larder BA
(2000) Phenotypic and genotypic analysis of clinical HIV-1 isolates reveals extensive protease inhibitor cross-resistance: a survey of over 6000 samples
Aids 14:1203–1210.
https://doi.org/10.1097/00002030-200006160-00018
- PubMed
- Google Scholar
(2011) A systems analysis of mutational effects in HIV-1 protease and reverse transcriptase
Nature Genetics 43:487–489.
https://doi.org/10.1038/ng.795
- PubMed
- Google Scholar
1. Hopf TA
2. Ingraham JB
3. Poelwijk FJ
4. Schärfe CP
5. Springer M
6. Sander C
7. Marks DS
(2017) Mutation effects predicted from sequence co-variation
Nature Biotechnology 35:128–135.
https://doi.org/10.1038/nbt.3769
- PubMed
- Google Scholar
1. Hu Z
2. Kuritzkes DR
(2010) Effect of raltegravir resistance mutations in HIV-1 integrase on viral fitness
JAIDS Journal of Acquired Immune Deficiency Syndromes 55:148–155.
https://doi.org/10.1097/QAI.0b013e3181e9a87a
- PubMed
- Google Scholar
1. Iyidogan P
2. Anderson KS
(2014) Current perspectives on HIV-1 antiretroviral drug resistance
Viruses 6:4095–4139.
https://doi.org/10.3390/v6104095
- PubMed
- Google Scholar
1. Izopet J
2. Massip P
3. Souyris C
4. Sandres K
5. Puissant B
6. Obadia M
7. Pasquier C
8. Bonnet E
9. Marchou B
10. Puel J
(2000) Shift in HIV resistance genotype after treatment interruption and short-term antiviral effect following a new salvage regimen
Aids 14:2247–2255.
https://doi.org/10.1097/00002030-200010200-00005
- PubMed
- Google Scholar
(2016) Benchmarking inverse statistical approaches for protein structure and design with exactly solvable models
PLOS Computational Biology 12:e1004889.
https://doi.org/10.1371/journal.pcbi.1004889
- PubMed
- Google Scholar
(2012) PSICOV: precise structural contact prediction using sparse inverse covariance estimation on large multiple sequence alignments
Bioinformatics 28:184–190.
https://doi.org/10.1093/bioinformatics/btr638
- PubMed
- Google Scholar
1. Keele BF
2. Giorgi EE
3. Salazar-Gonzalez JF
4. Decker JM
5. Pham KT
6. Salazar MG
7. Sun C
8. Grayson T
9. Wang S
10. Li H
11. Wei X
12. Jiang C
13. Kirchherr JL
14. Gao F
15. Anderson JA
16. Ping LH
17. Swanstrom R
18. Tomaras GD
19. Blattner WA
20. Goepfert PA
21. Kilby JM
22. Saag MS
23. Delwart EL
24. Busch MP
25. Cohen MS
26. Montefiori DC
27. Haynes BF
28. Gaschen B
29. Athreya GS
30. Lee HY
31. Wood N
32. Seoighe C
33. Perelson AS
34. Bhattacharya T
35. Korber BT
36. Hahn BH
37. Shaw GM
(2008) Identification and characterization of transmitted and early founder virus envelopes in primary HIV-1 infection
PNAS 105:7552–7557.
https://doi.org/10.1073/pnas.0802203105
- PubMed
- Google Scholar
(2017) Potts hamiltonian models of protein co-variation, free energy landscapes, and evolutionary fitness
Current Opinion in Structural Biology 43:55–62.
https://doi.org/10.1016/j.sbi.2016.11.004
- PubMed
- Google Scholar
1. Liu Y
2. Holte S
3. Rao U
4. McClure J
5. Konopa P
6. Swain JV
7. Lanxon-Cookson E
8. Kim M
9. Chen L
10. Mullins JI
(2013) A sensitive real-time PCR based assay to estimate the impact of amino acid substitutions on the competitive replication fitness of human immunodeficiency virus type 1 in cell culture
Journal of Virological Methods 189:157–166.
https://doi.org/10.1016/j.jviromet.2012.10.016
- PubMed
- Google Scholar
1. Lockless SW
2. Ranganathan R
(1999) Evolutionarily conserved pathways of energetic connectivity in protein families
Science 286:295–299.
https://doi.org/10.1126/science.286.5438.295
- PubMed
- Google Scholar
(2018) Fitness landscape of the human immunodeficiency virus envelope protein that is targeted by antibodies
PNAS 115:E564–E573.
https://doi.org/10.1073/pnas.1717765115
- PubMed
- Google Scholar
1. Mann JK
2. Barton JP
3. Ferguson AL
4. Omarjee S
5. Walker BD
6. Chakraborty A
7. Ndung'u T
(2014) The fitness landscape of HIV-1 gag: advanced modeling approaches and validation of model predictions by in vitro testing
PLOS Computational Biology 10:e1003776.
https://doi.org/10.1371/journal.pcbi.1003776
- PubMed
- Google Scholar
(2012) Protein structure prediction from sequence variation
Nature Biotechnology 30:1072–1080.
https://doi.org/10.1038/nbt.2419
- PubMed
- Google Scholar
(1999)
Human immunodeficiency virus type 1 cloning vectors for antiretroviral resistance testing

Journal of Clinical Microbiology 37:2943–2951.
- PubMed
- Google Scholar
(2016) Epistasis and the dynamics of reversion in molecular evolution
Genetics 203:1335–1351.
https://doi.org/10.1534/genetics.116.188961
- PubMed
- Google Scholar
1. Mézard M
2. Mora T
(2009) Constraint satisfaction problems and neural networks: a statistical physics perspective
Journal of Physiology-Paris 103:107–113.
https://doi.org/10.1016/j.jphysparis.2009.05.013
- PubMed
- Google Scholar
1. Mora T
2. Bialek W
(2011) Are biological systems poised at Criticality?
Journal of Statistical Physics 144:268–302.
https://doi.org/10.1007/s10955-011-0229-4
- Google Scholar
1. Morcos F
2. Pagnani A
3. Lunt B
4. Bertolino A
5. Marks DS
6. Sander C
7. Zecchina R
8. Onuchic JN
9. Hwa T
10. Weigt M
(2011) Direct-coupling analysis of residue coevolution captures native contacts across many protein families
PNAS 108:E1293–E1301.
https://doi.org/10.1073/pnas.1111471108
- PubMed
- Google Scholar
(1996) HIV-1 dynamics in vivo: virion clearance rate, infected cell life-span, and viral generation time
Science 271:1582–1586.
https://doi.org/10.1126/science.271.5255.1582
- PubMed
- Google Scholar
(2012) Amino acid coevolution induces an evolutionary stokes shift
PNAS 109:E1352–E1359.
https://doi.org/10.1073/pnas.1120084109
- PubMed
- Google Scholar
(2017) Elucidating the interdependence of drug resistance from combinations of mutations
Journal of Chemical Theory and Computation 13:5671–5682.
https://doi.org/10.1021/acs.jctc.7b00601
- PubMed
- Google Scholar
1. Rhee SY
2. Gonzales MJ
3. Kantor R
4. Betts BJ
5. Ravela J
6. Shafer RW
(2003) Human immunodeficiency virus reverse transcriptase and protease sequence database
Nucleic Acids Research 31:298–303.
https://doi.org/10.1093/nar/gkg100
- PubMed
- Google Scholar
1. Richman DD
2. Little SJ
3. Smith DM
4. Wrin T
5. Petropoulos C
6. Wong JK
(2004a)
Hiv evolution and escape

Transactions of the American Clinical and Climatological Association 115:289.
- PubMed
- Google Scholar
1. Richman DD
2. Morton SC
3. Wrin T
4. Hellmann N
5. Berry S
6. Shapiro MF
7. Bozzette SA
(2004b) The prevalence of antiretroviral drug resistance in the united states
Aids 18:1393–1401.
https://doi.org/10.1097/01.aids.0000131310.52526.c7
- PubMed
- Google Scholar
1. Sarafianos SG
2. Marchand B
3. Das K
4. Himmel DM
5. Parniak MA
6. Hughes SH
7. Arnold E
(2009) Structure and function of HIV-1 reverse transcriptase: molecular mechanisms of polymerization and inhibition
Journal of Molecular Biology 385:693–713.
https://doi.org/10.1016/j.jmb.2008.10.071
- PubMed
- Google Scholar
1. Shafer RW
(2006) Rationale and uses of a public HIV drug-resistance database
The Journal of Infectious Diseases 194:S51–S58.
https://doi.org/10.1086/505356
- PubMed
- Google Scholar
1. Shafer RW
2. Schapiro JM
(2008)
Hiv-1 drug resistance mutations: an updated framework for the second decade of haart

AIDS Reviews 10:67.
- PubMed
- Google Scholar
(2015) Contingency and entrenchment in protein evolution under purifying selection
PNAS 112:E3226–E3235.
https://doi.org/10.1073/pnas.1412933112
- PubMed
- Google Scholar
(2013) Spin models inferred from patient-derived viral sequence data faithfully describe HIV fitness landscapes
Physical Review E 88:062705.
https://doi.org/10.1103/PhysRevE.88.062705
- Google Scholar
1. Sjodt M
2. Brock K
3. Dobihal G
4. Rohs PDA
5. Green AG
6. Hopf TA
7. Meeske AJ
8. Srisuknimit V
9. Kahne D
10. Walker S
11. Marks DS
12. Bernhardt TG
13. Rudner DZ
14. Kruse AC
(2018) Structure of the peptidoglycan polymerase RodA resolved by evolutionary coupling analysis
Nature 556:118–121.
https://doi.org/10.1038/nature25985
- PubMed
- Google Scholar
1. Sułkowska JI
2. Morcos F
3. Weigt M
4. Hwa T
5. Onuchic JN
(2012) Genomics-aided structure prediction
PNAS 109:10340–10345.
https://doi.org/10.1073/pnas.1207864109
- PubMed
- Google Scholar
(2015) From residue coevolution to protein conformational ensembles and functional dynamics
PNAS 112:13567–13572.
https://doi.org/10.1073/pnas.1508584112
- PubMed
- Google Scholar
1. Tan Z
2. Xia J
3. Zhang BW
4. Levy RM
(2016) Locally weighted histogram analysis and stochastic solution for large-scale multi-state free energy estimation
The Journal of Chemical Physics 144:034107.
https://doi.org/10.1063/1.4939768
- PubMed
- Google Scholar
(2018) The impact of HIV-1 within-host evolution on transmission dynamics
Current Opinion in Virology 28:92–101.
https://doi.org/10.1016/j.coviro.2017.12.001
- PubMed
- Google Scholar
1. Troyer RM
2. McNevin J
3. Liu Y
4. Zhang SC
5. Krizan RW
6. Abraha A
7. Tebit DM
8. Zhao H
9. Avila S
10. Lobritz MA
11. McElrath MJ
12. Le Gall S
13. Mullins JI
14. Arts EJ
(2009) Variable fitness impact of HIV-1 escape mutations to cytotoxic T lymphocyte (CTL) response
PLOS Pathogens 5:e1000365.
https://doi.org/10.1371/journal.ppat.1000365
- PubMed
- Google Scholar
(2019) Learning protein constitutive motifs from sequence data
eLife 8:e39397.
https://doi.org/10.7554/eLife.39397
- PubMed
- Google Scholar
1. Weigt M
2. White RA
3. Szurmant H
4. Hoch JA
5. Hwa T
(2009) Identification of direct residue contacts in protein-protein interaction by message passing
PNAS 106:67–72.
https://doi.org/10.1073/pnas.0805923106
- PubMed
- Google Scholar
1. Wensing AM
2. Calvez V
3. Günthard HF
4. Johnson VA
5. Paredes R
6. Pillay D
7. Shafer RW
8. Richman DD
(2017)
2017 update of the drug resistance mutations in hiv-1

Topics in Antiviral Medicine 24:132–133.
- PubMed
- Google Scholar
1. Wu TD
2. Schiffer CA
3. Gonzales MJ
4. Taylor J
5. Kantor R
6. Chou S
7. Israelski D
8. Zolopa AR
9. Fessel WJ
10. Shafer RW
(2003) Mutation patterns and structural correlates in human immunodeficiency virus type 1 protease following different protease inhibitor treatments
Journal of Virology 77:4836–4847.
https://doi.org/10.1128/jvi.77.8.4836-4847.2003
- PubMed
- Google Scholar
(2015) Persistence of transmitted HIV-1 drug resistance mutations associated with fitness costs and viral genetic backgrounds
PLOS Pathogens 11:e1004722.
https://doi.org/10.1371/journal.ppat.1004722
- PubMed
- Google Scholar
1. Yilmaz NK
2. Schiffer CA
(2017)
Antimicrobial Drug Resistance

535–544, Drug resistance to hiv-1 protease inhibitors: molecular mechanisms and substrate coevolution, Antimicrobial Drug Resistance, Springer, 10.1007/978-3-319-46718-4_35.
- Google Scholar

Decision letter

Patricia J Wittkopp

Senior and Reviewing Editor; University of Michigan, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The first decision letter after peer review is shown below.]

Thank you for submitting your work entitled "Epistasis and entrenchment of drug resistance in HIV-1 subtype B" for consideration by eLife. Your article has been reviewed by a Senior Editor, a Reviewing Editor, and two reviewers. The reviewers have opted to remain anonymous.

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

The topic of your paper is very important, and indeed epistasis and the effects of coupled mutations in HIV proteins on the fitness of HIV is documented. However, as the reviews below note, because your study is based on individuals undergoing treatment, it is impossible to separate correlated selection pressures from epistatic effects.

Reviewer #1:

Here Biswas and collaborators apply coevolutionary analysis methods to study the entrenchment of HIV drug resistance mutations in the protease, integrase, and reverse transcriptase proteins. Entrenchment refers to an increasing preference for the mutant (i.e., resistant) residue as epistatic accessory mutations are accumulated. The authors demonstrate an impressive ability to predict the frequency of resistance mutations in different genetic backgrounds. They then characterize how entrenchment depends on the number of accumulated resistance-associated mutations, and which specific mutations are most strongly associated with entrenchment.

Specific comments:

1) One might reasonably hypothesize that the reason that drug resistance mutations are more likely to appear together in the Potts model is because the model is capturing correlated selection pressures rather than true epistasis. In other words, in patients who are treated by one drug, there is selection for the accumulation of resistance mutations for that drug. No such selection is present in patients not treated by that drug. Thus, whether the resistance/accessory mutations are truly epistatic or not, they will tend to appear either clustered together in sequences or not at all. This could then be reflected as a network of positive epistatic couplings in the Potts model. Butler et al., cited by the authors, used sequence data obtained before the widespread use of protease inhibitors to infer an Ising model that contained positive epistatic couplings between resistance sites. In principle, this approach should limit the effect of correlated selection pressures due to differences in drug treatment among different individuals. Butler and collaborators further found some external evidence for true epistasis between a resistance mutation at site 84 in protease and other accessory mutations. In the present work, given that the MSA contains sequences from individuals who are drug-experienced, are there independent lines of evidence that support the hypothesis of true epistatic interactions between associated resistance mutations? And more broadly, how would one deal with the confounding influence of correlated selection pressures?

2) As a possible suggestion, one might expect that inferring Potts parameters using the pseudolikelihood method would provide an even better fit of the probability of different mutations at a given site as a function of the sequence background, since this is precisely what the pseudolikelihood algorithm is optimizing. In contrast, the canonical maximum entropy approach must simultaneously fit all low-order correlations.

3) Some of the figures and tables are challenging to interpret. For example, in Supplementary file 1 both the most (ME) and least (LE) residues at position 17 for the L74M mutation are N. This result also appears for other mutations and positions in other tables. This needs further explanation. The plot in the inset of Figure 1B is also a bit challenging to interpret at first, and the small size of the numbers/labels in the inset plot makes it difficult to read at standard resolution. The ordering of resistance mutations in Figure 1D is confusing. The linear order would suggest order along the sequence, but mutations within each sequence are presented in reverse order along the protein sequence. Integrase is also placed before reverse transcriptase. This figure would likely be more clear if the same color scheme was maintained, but the mutations were reordered according to their positions in Pol.

Reviewer #2:

The study analyzes molecular co-evolution in sequences of HIV from drug-exposed patients, applying a Potts model to infer fitness effects of mutations and their dependence on the genetic background. The authors use the inferred model to analyze the notion of entrenchment - that a mutation tends to become more difficult to revert, over time, due the epistatic interactions with subsequent substitutions in the genetic background on which it first arose. The authors report this trend in three different HIV proteins that are targets of various drug treatments.

The intellectual basis of this study is sound and interesting - that epistasis can shape and channel molecular evolution - and it builds on a broader theoretical literature that has demonstrated such phenomenon for protein evolution in general (e.g. Pollack et al., 2012; Shah et al., 2015). Demonstrating this effect for real-world sampled protein sequences, in the setting of HIV protein evolution in response to drugs, is an especially nice project - as it would provide real-world relevance, in a clinical setting, for a broad body of theoretical work.

Despite the enthusiasm for the topics, I am disappointed by a large number of inconsistencies in the reported results - which do not, in the end, seem to demonstrate the type of entrenchment the authors start out to explore. And I'm also disappointed by the presentation of the model and analysis, which is ambiguous or unclear in both main text and Materials and methods section. I am left without a clear understanding of what the authors claim to have demonstrated, or what specific model inferences they are reporting in most of their figures.

These two problems - whether their analysis actually addresses the established concept of entrenchment, and whether the writing/notation is clear enough to explain what analyses they actually performed – occur throughout the manuscript, and so I will address instances of these two problems in a linear fashion throughout the work:

1) The Potts model appears to be massively overfit. Even after artificially reducing the alphabet size from 20 acids to 4, the authors seem to be fitting roughly 380,000 parameters to account for all pairwise interactions between sites. And yet, as far as a I can tell, the authors have only ~20,000 sequences. So, they are fitting an order of magnitude more parameters than they have data. Even so, they make a big deal of the fact the model makes good "predictions" for the frequency of sequences observed (Figure 1B)- but it is not clear if these prediction were made with cross-validation (e.g. fitting on one part of the data, and prediction for another part), and little discussion of overfitting is provided, if any. In any case, the notion that pairwise maxent Potts-like models can do a decent job of explaining sequence variability is not a new result and has been shown in dozens of papers.

1b) The authors discuss their Potts model as providing a general probability for the observing a mutation alpha at site i in a genetic background of sequence S. But, in fact, the model depends only on pairwise interactions, the full sequence S. The authors should clarify this in the main text.

1c) It is not clear what reference sequence is used as the "wildtype" S throughout the analysis. Is this a consensus sequence from the database? This should be explained in the main text.

2) The main results on entrenchment (Figure 2) are not described precisely and seem to contradict the notion of entrenchment to which the authors refer in their introduction. Entrenchment according to Pollack et al., occurs when a substitution is initially neutral or slightly beneficial (so that its reversion is neutral or only slightly deleterious), but when subsequent substitutions in the protein render the focal mutation highly deleterious to revert. And yet the mutations discussed here are drug resistance mutations arising in the presence of drugs - and so the primary mutation must be net highly beneficial, even if it partly disrupts HIV protein function. (Indeed Figure 4 shows that the primary DRMS are favored.) And so, prima facie, this situation is quite different to Pollack at all - because the reversion is immediately disfavored, even without changes to the genetic background. Presumably the authors intend to show that the DRMs become *even more deleterious* to revert as subsequent substitutions accrue, although they do not say this in the main text and repeatedly refer to the DRMs as carrying fitness costs as the time of their emergence.

2b) Figure 2A, showing the main entrenchment results according to the Potts model, seems to contradict the notion the primary mutation is a resistance site that is beneficial at the time of its substitution (as in Figure 4). The figure reports a nearly 100-fold preference to revert the primary DRM L90M substitution in the background in which it arose - i.e. the mutation is inferred to carry a huge fitness cost (and its reversion a benefit) at the time of its origin. This does not make sense to me, as a reader, as the DRM (even if it contains some costs for protein function) must be net adaptive at the time of its substitution in a patient. Why does Figure 2A show that the drug resistance site is strongly disfavored in the genetic background in which it initially arose?

2c) The authors should define in terms of the ΔE notation when, precisely, is being plotted on the y-axis of Figure 4A.

2d) The authors state "entrenchment of key resistance mutations is likely to play an important role in emergence of drug-resistance viral strains". I cannot understand what the authors mean by this. Certainly, the initial *emergence* of a drug resistance strain depends on the adaptive benefit of DRM on the background in which it arises - and has nothing to do with future changes to the genetic background. Perhaps the authors mean that the "persistence" of the DRM over subsequent evolution depends on entrenchment?

3) [Presentation issue] The authors discuss the stabilization or destabilization of mutations in their Introduction and throughout (Table 1). I believe by "stabilization" they mean epistasis with subsequent substitutions that render the primary mutation costly to revert. But this is a terrible choice of words, because the large literature on entrenchment (Pollack, Shah etc.) has been developed for models of thermostability - and, indeed, entrenchment has been observed to arise from interactions between sites for protein stability. The authors do not seem to be using the word "stability" here to mean thermostability, which is tremendously confusing when referring to a literature that is all about the effects of mutations, and their interactions, on thermostability.

4) The authors make a big deal of their analysis of "entrenchment of DRMs" in the population - meaning that drug resistance mutations are beneficial in over 50% of the sequences sampled/fit. I do not find this surprising or even attributable to entrenchment at all - because, after all, the sequences being used to fit the model are sequences from drug-experienced HIV populations, where the drug resistance mutation will be net beneficial *at the time of its original occurrence* and even without any further substitutions. This is shown directly in Figure 4, for example.

5) When the authors report that they have "verification" of entrenchment I assumed they meant they had direct data on fitness measurements for a mutation in one genetic background, and in a subsequent background. But in fact, the "verification" is simply relative to the fitted Potts model. By contrast, when the authors do compare the model to actual measurements of fitness effects (subsection “Molecular clones and the effects of specific backgrounds”), the authors seem to be saying that the predictions of the Potts model do not match the empirical measurements, and they attribute this to limitation of the mutagenesis experiments, as opposed to a deficiency of the Potts model. Comparison to data would have been the only thing that really convinces me that any statistical effects of their model are "verified".

5b) Also, I had a hard time understanding what conclusions the authors draw from Figure 5, or even what the axes really denote in Figure 5. The notation Δ E_patient does not seem to be defined in the main text. The authors seem to be showing effects of mutations predicted by the Potts model fitted to patient data versus their predicted effects in the specific genetic backgrounds used for mutagenesis experiments - but it's not clear if that's what they're plotting, in fact, or what we are supposed to conclude from Figure 5.

https://doi.org/10.7554/eLife.50524.029

Author response

[Editors' note: the author responses to the re-review follow.]

Reviewer #1:

Here Biswas and collaborators apply coevolutionary analysis methods to study the entrenchment of HIV drug resistance mutations in the protease, integrase, and reverse transcriptase proteins. Entrenchment refers to an increasing preference for the mutant (i.e., resistant) residue as epistatic accessory mutations are accumulated. The authors demonstrate an impressive ability to predict the frequency of resistance mutations in different genetic backgrounds. They then characterize how entrenchment depends on the number of accumulated resistance-associated mutations, and which specific mutations are most strongly associated with entrenchment.

Specific comments:

1) One might reasonably hypothesize that the reason that drug resistance mutations are more likely to appear together in the Potts model is because the model is capturing correlated selection pressures rather than true epistasis. In other words, in patients who are treated by one drug, there is selection for the accumulation of resistance mutations for that drug. No such selection is present in patients not treated by that drug. Thus, whether the resistance/accessory mutations are truly epistatic or not, they will tend to appear either clustered together in sequences or not at all. This could then be reflected as a network of positive epistatic couplings in the Potts model. Butler et al., cited by the authors, used sequence data obtained before the widespread use of protease inhibitors to infer an Ising model that contained positive epistatic couplings between resistance sites. In principle, this approach should limit the effect of correlated selection pressures due to differences in drug treatment among different individuals. Butler and collaborators further found some external evidence for true epistasis between a resistance mutation at site 84 in protease and other accessory mutations. In the present work, given that the MSA contains sequences from individuals who are drug-experienced, are there independent lines of evidence that support the hypothesis of true epistatic interactions between associated resistance mutations? And more broadly, how would one deal with the confounding influence of correlated selection pressures?

The reviewer’s comments suggest that we need to address three important questions: (1) Can we distinguish between correlations induced by “true epistastic interactions” and those due to correlated selection pressure induced by drugs; (2) Do our conclusions about the entrenchment of drug resistance mutations depend on whether the origins of the mutational correlation patterns can be separated into “intrinsic epistasis” and effects due to correlated selection pressure induced by drugs; (3) Is the Potts model we build from the drug experienced dataset contaminated by artifacts associated with different patients being given different drug regimens.

1) Can we distinguish between correlations induced by “true” or “intrinsic” epistasis and those induced by correlated selection pressure induced by drugs?

To begin with, we note that while reviewer #1 refers to epistatic effects which are independent of drug selection pressure as “true epistasis”, we prefer to describe them as “intrinsic epistasis”, since correlated mutations induced by drug selection pressure may also be considered to be a form of epistasis. We use the term “effective epistasis” to refer to mutational correlations induced by the combined effects of “intrinsic epistasis” and correlated selection pressure due to drugs. We have added a new subsection “Effective epistasis in the presence of drug selection pressure” where we discuss the distinction between “intrinsic epistasis” and correlated selection pressure.

In the case of immune response selection pressure Shekhar et al., 2013 suggest that due to the diversity of host immune responses among the HIV population, a Potts model fit to sequences from many different hosts can effectively capture the “intrinsic” fitness landscape of the virus because the selective effects of immune pressure are averaged in a way that can be modeled by adjusting the fields of the Potts model, leaving the couplings unchanged.

We have available to us from the Stanford HIV database, a drug experienced sequence data set and a drug naïve dataset. We have fit a Potts model for HIV-1 Subtype B Protease to each dataset individually. Comparisons of the Potts statistical energy values ΔE predicted by the drug naïve and drug experienced models for mutating from the wild type to the drug resistance mutations for four HIV PR DRMS as a function of the background sequence are presented in subsection “Effective epistasis in the presence of drug selection pressure” of the revised manuscript. There is a high correlation (r > 0.8) between the probabilities of observing the drug resistance mutation relative to the wild type in a given background, when the Potts model is parameterized on the drug naïve dataset as compared with the drug experienced dataset. This means that “intrinsic” epistatic effects have a large influence on the virus evolving under drug selection pressure.

It is possible to carry out a more quantitative analysis of how the application of drug selection pressure changes the intrinsic mutational correlation patterns and probabilities of observing DRMs in specific sequence backgrounds using multi-canonical reweighting techniques familiar to researchers who work with Ising or Potts models. While such an analysis is beyond the scope of the present work, in a future communication we will provide a detailed comparison of the drug naïve and drug experienced Potts models and the relationship between them, as such an analysis can provide insights into sequence patterns which are most strongly associated with the transmission of drug resistance.

2) Do our conclusions about the entrenchment of drug resistance mutations depend on whether the origins of the mutational correlation patterns can be separated into “intrinsic epistasis” and effects due to correlated selection pressure induced by drugs?

As far as the goal of the current work is concerned, to demonstrate that the Potts model can accurately predict the effect of the sequence background on the probability of observing a resistance mutation in the drug experienced population, and therefore identify those sequences which are highly favoring (entrenching) or highly disfavoring the mutation, it is not necessary to separate the effects of intrinsic epistasis from those due to correlated selection pressure induced by drugs. The data we are comparing the predictions of our Potts models to are the observed frequencies (prevalences) in the Stanford HIV drug experienced sequence database of drug resistance mutations, in sets of sequence backgrounds classified by our model as either favoring or disfavoring the DRM. These are non-trivial predictions; the observations confirm the predictions. Put another way, the Potts model parameterized on the drug experienced dataset is a powerful classifier that can identify complex sequence patterns that highly favor (entrench) or disfavor each drug resistant mutant. It should be noted that the frequency of any individual drug resistance mutation (DRM) appearing in the drug experienced Stanford HIV database is typically small. See our response to reviewer #2 comment 4) for additional discussion of the distinction between the frequency of a DRM appearing in the drug experienced population and the “entrenchment” of that DRM.

3) Is the Potts model we build from the drug experienced dataset contaminated by an artifact associated with different patients being given different drug regimens? This could result in spurious apparent correlations predicted by the Potts model.

Reviewer #1 also expresses a concern that if the selection pressures induced by different drugs are significantly different, then this could lead to artifacts when parameterizing a Potts model on a drug experienced dataset, as the model would potentially be sensitive to the details of how many patients in the dataset are exposed to each drug. While we cannot definitively rule out such an artifact, two lines of evidence suggest that if such effects are present, they are small.

First, we note that HIV drugs are frequently given in combination, and that many of the drugs have similar resistance patterns. We include a comparison of these patterns for different Protease inhibitors in Figure 8 of the revised manuscript. Secondly, as noted above there is a strong correlation between the ΔE values predicted using the Potts model parameterized on the drug naïve dataset and the one parameterized on the drug experienced set. If the details concerning the number of individuals which were given specific drugs and combination therapies strongly influenced the Potts model parameterized on the drug experienced dataset, we would not expect the “effective epistatic” effects on ΔE inferred by the drug experienced model to be highly correlated with the “intrinsic epistatic” effects of the drug naïve model. These observations mitigate possible artifacts associated with the fact that individual HIV patients represented in the drug experienced dataset may be undergoing different drug therapies. We added a statement about in subsection “Effective epistasis in the presence of drug selection pressure” of the revised manuscript.

2) As a possible suggestion, one might expect that inferring Potts parameters using the pseudolikelihood method would provide an even better fit of the probability of different mutations at a given site as a function of the sequence background, since this is precisely what the pseudolikelihood algorithm is optimizing. In contrast, the canonical maximum entropy approach must simultaneously fit all low-order correlations.

We are using an approach which optimizes the maximum entropy Boltzmann likelihood function by fitting the bivariate marginals of the Potts model to the observed bivariate marginals using MCMC sampling. This is a more compute intensive approach to inferring the parameters of the Potts Hamiltonian, but it is also more accurate as the errors in the bivariate and univariate marginals are much smaller than those obtained by the pseudolikelihood method. This has been demonstrated in the literature; for instance, the supplementary information in Barton et al., 2016 shows large relative errors in the generated univariate marginals and the pairwise correlations using the pseudolikelihood method, and in Jacquin et al., 2016, Barton et al., 2016, it is shown that the pseudolikelihood method is not "generative", meaning sequences generated using it have modified bivariate marginals. The agreement we obtain using our MCMC GPU implementation of inverse inference, between predicted and observed frequencies of drug resistance mutations in different sets of backgrounds is excellent, and probably cannot be matched using more approximate inverse inference methods like the pseudo-likelihood.

3) Some of the figures and tables are challenging to interpret. For example, in Supplementary file 1 both the most (ME) and least (LE) residues at position 17 for the L74M mutation are N. This result also appears for other mutations and positions in other tables. This needs further explanation. The plot in the inset of Figure 1B is also a bit challenging to interpret at first, and the small size of the numbers/labels in the inset plot makes it difficult to read at standard resolution. The ordering of resistance mutations in Figure 1D is confusing. The linear order would suggest order along the sequence, but mutations within each sequence are presented in reverse order along the protein sequence. Integrase is also placed before reverse transcriptase. This figure would likely be more clear if the same color scheme was maintained, but the mutations were reordered according to their positions in Pol.

Supplementary file 1 in the updated supplement presents a pair of sequences for each mutation; one being the most entrenching for that mutation and the other the least entrenching (the probability of observing the mutation relative to the wildtype at this position is the smallest in the least entrenching background). There will be positions at which the same accessory mutation appears in both the most and least entrenching backgrounds; for example, the mutation 17N, which is accessory to L74M. The background sequences which are the most and least entrenching for any resistance mutation is a collective property of the set of all mutations in the corresponding background; any particular accessory mutation may appear in both backgrounds.

As per the reviewer’s suggestion, we have changed the ordering of mutations in Figure 1D according to their positions in Pol and resized Figure 1 for clarity.

Reviewer #2:

The study analyzes molecular co-evolution in sequences of HIV from drug-exposed patients, applying a Potts model to infer fitness effects of mutations and their dependence on the genetic background. The authors use the inferred model to analyze the notion of entrenchment – that a mutation tends to become more difficult to revert, over time, due the epistatic interactions with subsequent substitutions in the genetic background on which it first arose. The authors report this trend in three different HIV proteins that are targets of various drug treatments.

The intellectual basis of this study is sound and interesting – that epistasis can shape and channel molecular evolution – and it builds on a broader theoretical literature that has demonstrated such phenomenon for protein evolution in general (e.g. Pollack et al., 2012; Shah et al., 2015). Demonstrating this effect for real-world sampled protein sequences, in the setting of HIV protein evolution in response to drugs, is an especially nice project – as it would provide real-world relevance, in a clinical setting, for a broad body of theoretical work.

We appreciate the reviewer’s enthusiasm for the topic of the paper.

Despite the enthusiasm for the topics, I am disappointed by a large number of inconsistencies in the reported results – which do not, in the end, seem to demonstrate the type of entrenchment the authors start out to explore. And I'm also disappointed by the presentation of the model and analysis, which is ambiguous or unclear in both main text and Materials and methods section. I am left without a clear understanding of what the authors claim to have demonstrated, or what specific model inferences they are reporting in most of their figures.

These two problems – whether their analysis actually addresses the established concept of entrenchment, and whether the writing/notation is clear enough to explain what analyses they actually performed – occur throughout the manuscript, and so I will address instances of these two problems in a linear fashion throughout the work:

We have made revisions to the manuscript to make it clearer following the suggestions of the reviewers, and in light of what appears to be a serious misunderstanding of our work at several places by reviewer #2. We do not believe that there are inconsistencies in the reported results; it is apparent from the reviewer’s comments that we have done an inadequate job of describing our work. Several statements made by referee #2 are incorrect. We have revised the manuscript to address the issues and concerns raised by referee #2 as described below.

1) The Potts model appears to be massively overfit. Even after artificially reducing the alphabet size from 20 acids to 4, the authors seem to be fitting roughly 380,000 parameters to account for all pairwise interactions between sites. And yet, as far as a I can tell, the authors have only ~20,000 sequences. So, they are fitting an order of magnitude more parameters than they have data. Even so, they make a big deal of the fact the model makes good "predictions" for the frequency of sequences observed (Figure 1B) - but it is not clear if these prediction were made with cross-validation (e.g. fitting on one part of the data, and prediction for another part), and little discussion of overfitting is provided, if any. In any case, the notion that pairwise maxent Potts-like models can do a decent job of explaining sequence variability is not a new result and has been shown in dozens of papers.

Overfitting is an important issue and we appreciate the suggestion that we should discuss it more thoroughly. We have recently published a detailed analysis of overfitting effects in the context of inverse inference of Potts model parameters from protein multiple-sequence alignments (Haldane and Levy, 2019). We have added a discussion of overfitting effects to the main text in subsection “Statistical Robustness of HIV Potts models” which we summarize here.

The reviewer suggests that we are overfitting our model because the number of model parameters (θ) is much larger than the number of sequences in our dataset (N). However, we explain in Haldane and Levy, 2019 why this is not the correct comparison to make. First, our model is not fit directly to the sequences but rather to the bivariate marginals of the MSA, and there are an equal number of model parameters θ as there are bivariate marginals. The inference problem is therefore neither overconstrained nor underconstrained. In addition, while for (linear) regression overfitting occurs when the number of datapoints is less than the number of model parameters, this is not the case for inference methods in general.

In Haldane and Levy, 2019 we explain how overfitting ultimately arises due to finite-sampling error in the bivariate marginal frequencies used as input to the inference procedure, which are computed from the MSA of N sequences. Each bivariate marginal f is estimated from a sample of size N, and its statistical error is the multinomial mean-square error f(1-f)/N. This statistical error in the bivariate marginals leads to error in the inferred model parameters.

Based on this fact, in Haldane and Levy, 2019 we go on to derive the dependence of model quality on the number of parameters and the number of sequences in the MSA, and we show how one can fit accurate Potts models to MSAs containing only thousands of effective sequences with negligible overfitting, even though there are many more Potts model parameters than sequences. We illustrate this using a simple-to-understand toy model and demonstrate it numerically for realistic Potts models. We find that the degree of overfitting can be quantified using the "signal-to-noise ratio" (SNR), which is a function of the number of model parameters θ, the degree of conservation of the MSA as measured by a value called Χ², and the number of sequences N. While the SNR is related to the ratio of number of parameters θ to number of sequences N, it contains other terms including Χ^2.. In order to demonstrate that our Potts models for HIV proteins are not overfit, we have added a discussion of the SNR of our models to the revised manuscript showing that they have suitable SNR and satisfy the criteria described in Haldane and Levy, 2019 for estimating statistical energy differences without overfitting. Our paper (Haldane and Levy, 2019) also discusses cross-validation of Potts statistical energies and quantities which can be derived from the statistical energies such as ΔE.

We think reviewer #2’s comment “In any case, the notion that pairwise maxent Potts-like models can do a decent job of explaining sequence variability is not a new result and has been shown in dozens of papers.” is inappropriate. We have been active and publishing in this field since 2012 and are well versed in the literature which is extensively cited. Our analysis of the effects of the sequence background on the likelihood of drug resistance mutations in HIV-1, which first appeared in 2017 (Flynn et al., 2017) and which is presented in a more complete form in the current manuscript where the predictions are verified, is highly novel, there is nothing comparable in the literature. Our analysis of the entrenchment of DRMS for HIV-1 is not simply a calculation of the fitness of single point mutations, which has been reported by several groups working with Potts models of sequence variability. In order to demonstrate the entrenchment of a DRM by the background it is necessary to be able to capture the effects of many simultaneous mutations on DRM sites and how those effects change as the sequence background changes. This is related to our ability to model higher order marginals of the distribution, which is possible using our MCMC method to infer the Potts model. It would be quite difficult to obtain results concerning the entrenchment of DRMS of comparable quality using more approximate inference methods.

1b) The authors discuss their Potts model as providing a general probability for the observing a mutation alpha at site i in a genetic background of sequence S. But, in fact, the model depends only on pairwise interactions, the full sequence S. The authors should clarify this in the main text.

While the Potts Hamiltonian model contains only pairwise interaction terms, the Potts effective potential function induces higher order correlations and the probability of a given sequence cannot be decomposed into a sum of products of lower order terms. In other words, pairs of residues interact with each other both directly (through a direct coupling in the Hamiltonian) and indirectly through chains of interactions involving one or more intermediate residues. As a consequence, the higher order marginals of the distribution cannot be expressed as an explicit function of the pair correlations. Our ability to predict the effect of the background on the likelihood of a drug resistance mutation is related to the ability of our MCMC algorithm to infer a Potts model of sequence co-variation which is sufficiently accurate to infer higher order marginals of the distribution as described in Flynn et al., 2017; Haldane et al., 2018. The Potts model Hamiltonian with pairwise terms is both necessary and sufficient to capture many observable effects that involve higher order correlations, like the effect of the sequence background on the likelihood of a DRM. Whether or not Hamiltonian functions which include triplet or even higher order terms are needed to capture other properties of the HIV sequence ensembles is a question that is not strictly germane to our paper, and in any case is best studied first in toy models because of the data limitations inherent in HIV sequence databases. We have added a discussion of this point to the main text in subsection “Epistasis: the effect of the background on primary resistance mutations in HIV” and subsection “Model inference”.

1c) It is not clear what reference sequence is used as the "wildtype" S throughout the analysis. Is this a consensus sequence from the database? This should be explained in the main text.

For all three proteins (PR, RT, IN) the reference sequence referred to as “wildtype” is the consensus subtype B sequence, as obtained from an alignment of subtype B sequences maintained in the Los Alamos HIV sequence database, which is frequently used as the reference for mutation studies. Where the consensus sequence can be obtained is stated the revised Materials and methods section. Since this reference sequence is often referred to as the “consensus wild-type subtype B” sequence, we have adopted this terminology in the main text and in the Materials and methods section.

2) The main results on entrenchment (Figure 2) are not described precisely and seem to contradict the notion of entrenchment to which the authors refer in their introduction. Entrenchment according to Pollack et al. occurs when a substitution is initially neutral or slightly beneficial (so that its reversion is neutral or only slightly deleterious), but when subsequent substitutions in the protein render the focal mutation highly deleterious to revert. And yet the mutations discussed here are drug resistance mutations arising in the presence of drugs – and so the primary mutation must be net highly beneficial, even if it partly disrupts HIV protein function. (Indeed Figure 4 shows that the primary DRMS are favored.) And so, prima facie, this situation is quite different to Pollack at all – because the reversion is immediately disfavored, even without changes to the genetic background. Presumably the authors intend to show that the DRMs become *even more deleterious* to revert as subsequent substitutions accrue, although they do not say this in the main text and repeatedly refer to the DRMs as carrying fitness costs as the time of their emergence.

2b) Figure 2A, showing the main entrenchment results according to the Potts model, seems to contradict the notion the primary mutation is a resistance site that is beneficial at the time of its substitution (as in Figure 4). The figure reports a nearly 100-fold preference to revert the primary DRM L90M substitution in the background in which it arose – i.e. the mutation is inferred to carry a huge fitness cost (and its reversion a benefit) at the time of its origin. This does not make sense to me, as a reader, as the DRM (even if it contains some costs for protein function) must be net adaptive at the time of its substitution in a patient. Why does Figure 2A show that the drug resistance site is strongly disfavored in the genetic background in which it initially arose?

We use the term entrenchment to refer to the fact that the sequence background has a large effect on the relative likelihood of observing a drug resistance mutation vs. a consensus wildtype residue at a DRM site. For some sequence backgrounds with several accessory mutations the background is highly favorable for the DRM; we refer to these mutation patterns as “entrenching”, while for other backgrounds with the same total number of associated mutations, those backgrounds highly disfavor the DRM. Our Potts Hamiltonian model is able to predict these sequence specific effects (entrenching vs. disfavoring DRMS) with high fidelity, even for sets of sequences conditional on having the same total Hamming distance from the wild type consensus sequence. This is an achievement that would not be possible using more approximate inverse inference techniques to construct the Potts Hamiltonian.

Perhaps the reviewer’s questions 2 and 2b are due to a misunderstanding concerning the information contained in Figure 2? It appears that the reviewer is interpreting Figure 2 as a time ordering of the appearance of the L90M mutation in backgrounds with increasing numbers of mutations. There is no information about time ordering of the appearance of L90M in Figure 2. We are not currently using the Potts model to infer information about the kinetics by which drug resistance is acquired. In fact, the drug resistance mutation 90M appears in the Stanford drug experienced dataset only 2 times in the wild type consensus background, while the wild-type consensus residue L90 appears 40 times in the wild type consensus background in the same dataset. Even though L90M is a drug resistance mutation, it is deleterious in the consensus wildtype background because of intrinsic epistatic effects, and according to the Stanford drug experienced database, the L90M mutation in the consensus background is ~ 20x less likely to be observed than the wild-type residue (L90) in the same consensus background sequence.

While Figure 2 should not be thought of as a time ordering of the appearance of L90M in a population, we can speculate about kinetics given the information in Figure 2. A possible scenario would be that L90M first appears in backgrounds where the L90M mutation is almost neutral, there are many such sequences with a Hamming distance of ~7 to 10 mutations from the consensus wildtype. Many of these backgrounds are almost neutral for L90M. Additional mutations, sometimes accompanied by reversions, will in general increase the likelihood of L90M leading to stronger entrenchment. Such a scenario is consistent in a general way with the time dependent model of entrenchment described by Shah et al., (2015) and Pollock et al., (2012), even though our current model does not include kinetics. But there are differences with the scenario envisioned by these authors in that the initial HIV-1 DRM mutation in the wildtype consensus sequence need not be almost neutral, since the path to entrenchment of drug resistance mutations need not pass through that sequence containing 1 DRM in the consensus background (the paper by Shah et al., does describe a similar scenario). We have revised the main text (subsection ““Entrenchment” of a primary resistance mutation and its verification”, Discussion section) to make it clear that our model is not a kinetic model, and that the L90M mutation is only rarely observed in the wild type consensus background in the drug experienced Stanford database. We plan to investigate the time evolution of the appearance of drug resistance in specific backgrounds using explicitly kinetic models in upcoming work.

2c) The authors should define in terms of the ΔE notation when, precisely, is being plotted on the y-axis of Figure 4A.

It is not clear whether the reviewer is referring to Figure 2A or Figure 4? Concerning the Y-axis on the rhs of Figure 2A, D,E is the difference in the Potts statistical energy of the sequence with the consensus wild type residue L at position 90 minus the Potts statistical energy of the drug resistance mutation M at position 90. A positive value ΔE means the mutant is favored. ΔE is the log of the probability of observing the mutant M in the background relative to the wild type residue in the background. The box plots show the statistics (median, and quartiles) for ΔE.

Figure 4 shows the distributions of ΔE values for key residues associated with drug resistance in HIV-1 IN. The green histogram shows the calculated ΔE for every drug resistance position in sequences where the drug resistance mutation is present, while the blue histogram shows the calculated ΔE to all other residue types in the same sequence background. The relative displacement of the green from the blue distributions shows that in sequences which contain the drug resistance mutation, the mutation is predicted on average to be more likely than the consensus residue type in that background (green distribution), and mutations to other residue types in that same sequence background are less likely (blue distribution). The green histogram is normalized to the total drug resistance mutation count in the Stanford HIV database (the total number of drug resistance mutations in each sequence times the total number of sequences), while the blue histogram is normalized to the total number of other possible mutations excluding the drug resistance mutation. This is made clear in the main text in subsection “Fitness and degree of entrenchment of observed resistance mutations: why some mutations are seen and others are not” and the legend for Figure 4.

2d) The authors state "entrenchment of key resistance mutations is likely to play an important role in emergence of drug-resistance viral strains". I cannot understand what the authors mean by this. Certainly, the initial *emergence* of a drug resistance strain depends on the adaptive benefit of DRM on the background in which it arises – and has nothing to do with future changes to the genetic background. Perhaps the authors mean that the "persistence" of the DRM over subsequent evolution depends on entrenchment?

We can understand why a reader might be confused by the statement "entrenchment of key resistance mutations is likely to play an important role in emergence of drug-resistance viral strains". We have revised the text in the Discussion section to provide a better explanation of what we intended to say concerning the relationship between entrenchment and the emergence of drug resistance. In the new text we state clearly that in this manuscript we are not modeling the time development of drug resistance or identifying evolutionary paths by which drug resistance is acquired. However, our construction of Potts models for both the drug naïve and drug experienced populations, sets the stage for kinetic studies of the pathways by which drug resistance is acquired. In order to identify HIV sequences which are most responsible for conferring drug resistance, we need to know the probabilities of observing highly entrenching sequences in both the drug naïve and drug experienced populations. As we discuss in the Discussion section of the revised manuscript, multicanonical reweighting techniques that can be applied to Ising and Potts models have been developed for this purpose.

Concerning the reviewer’s statement “Certainly the initial *emergence* of a drug resistance strain depends on the adaptive benefit of DRM on the background in which it arises”, we emphasize again that the initial “emergence” of a drug resistance mutation does not need to appear in the consensus wild type background (see our response to reviewer #2 comment 2 and 2b above), and for example the Protease resistance mutation L90M is rarely observed in the wildtype consensus background in the Stanford drug experienced dataset.

3) [Presentation issue] The authors discuss the stabilization or destabilization of mutations in their Introduction and throughout (Table 1). I believe by "stabilization" they mean epistasis with subsequent substitutions that render the primary mutation costly to revert. But this is a terrible choice of words, because the large literature on entrenchment (Pollack, Shah etc.) has been developed for models of thermostability – and, indeed, entrenchment has been observed to arise from interactions between sites for protein stability. The authors do not seem to be using the word "stability" here to mean thermostability, which is tremendously confusing when referring to a literature that is all about the effects of mutations, and their interactions, on thermostability.

We have removed the words “stabilizing” and “destabilizing” when referring to mutations as we agree that these words are often used in the context of the effects of a mutation on “protein stability” or “thermostability”. Instead we use the terms “favoring” or “disfavoring” to describe the correlated effects of sequence backgrounds on the probability of a drug resistance mutation appearing in that background.

The Potts model infers the probabilities of sequences in a population; several papers have characterized how well Potts statistical energies track fitness. There are many different ways that mutations can affect fitness; they include mutational effects on thermostability, replicative capacity, viral infectivity, etc. All of these factors will contribute to the prevalence of the sequence in the Stanford HIV database. (Also, see our answer to reviewer #2 comment 4 next).

4) The authors make a big deal of their analysis of "entrenchment of DRMs" in the population – meaning that drug resistance mutations are beneficial in over 50% of the sequences sampled/fit. I do not find this surprising or even attributable to entrenchment at all – because, after all, the sequences being used to fit the model are sequences from drug-experienced HIV populations, where the drug resistance mutation will be net beneficial *at the time of its original occurrence* and even without any further substitutions. This is shown directly in Figure 4, for example.

Regrettably, the reviewer appears to misunderstand the points we make about “entrenchment of DRMS” and therefore misunderstand the main theme of our paper. The reviewer’s statement that “drug resistance mutations are beneficial in over 50% of the sequences sampled/fit” is incorrect. This is clear simply from an analysis of the frequencies (univariate marginals) of the four classes of drug resistance mutations in the three target proteins (RT mutations against NRTIs and NNRTIs, PR mutations against PIs, and IN mutations against INSTIs). In fact, for all four classes of DRMs, the frequency with which each DRM appears in the Stanford HIV drug experienced database is almost always under 50% (there is one exception, the NRTI resistance mutation M184V); and the frequencies of individual DRMS in the drug experienced dataset are usually much smaller than 50%.

We use the term “Entrenchment of a DRM in the population” to refer to the fact that for most DRM mutations (the exception being NNRTIs), we observe that the DRM is favored over the consensus residue type at the corresponding DRM position (i.e. the Potts energy change is ΔE >0 at that position, and this depends on the residues at every other position) for more than 50% of the sequences which contain the DRM. If a mutation is entrenched in the population it is then unlikely or difficult to revert in more than 50% of the sequences in which the DRM is present.

This information was contained in Table 2—source data 1, Table 2—source data 2, Table 2—source data 3, and Table 2—source data 4, of the original manuscript. We have revised the tables to make the point clearer and revised the main text in subsection “Comparative entrenchment of resistance mutations in protease, reverse transcriptase, and integrase” to point this out. Most sequence backgrounds are unfavorable for any specific DRM. Consider an example. For the π mutation D30N in Protease, only 8% of the sequences in the Stanford HIVDB contain this π mutation, but 66% of the sequences with the D30N mutation have sequence backgrounds which are entrenching, meaning that there is a greater than even chance of observing D30N in 66% of the sequences containing D30N (see Table 2—source data 3). Consider another example, the INSTI-selected mutation Q148H which occurs with frequency 19% in the Stanford Drug Experienced dataset, but in 98% of these sequences (containing Q148H), the Potts model predicts that the Q148H mutation is more likely than not (i.e. that it is entrenched in these sequences). We have revised the text in subsection “Comparative entrenchment of resistance mutations in protease, reverse transcriptase, and integrase” to make it very clear that there is a distinction between the frequency of observing a DRM (which is usually low), and the likelihood that the DRM is “entrenched” in the sequences which contain that mutation (which is usually high).

Entrenchment in the population is better defined with an example the second paragraph of subsection “Comparative entrenchment of resistance mutations in protease, reverse transcriptase, and integrase”. It is also discussed in the fifth paragraph of subsection “Comparative entrenchment of resistance mutations in protease, reverse transcriptase, and integrase” leading to a discussion of sequences containing at least one entrenched DRM.

5) When the authors report that they have "verification" of entrenchment I assumed they meant they had direct data on fitness measurements for a mutation in one genetic background, and in a subsequent background. But in fact, the "verification" is simply relative to the fitted Potts model. By contrast, when the authors do compare the model to actual measurements of fitness effects (subsection “Molecular clones and the effects of speci1c backgrounds”), the authors seem to be saying that the predictions of the Potts model do not match the empirical measurements, and they attribute this to limitation of the mutagenesis experiments, as opposed to a deficiency of the Potts model. Comparison to data would have been the only thing that really convinces me that any statistical effects of their model are "verified".

5b) Also, I had a hard time understanding what conclusions the authors draw from Figure 5, or even what the axes really denote in Figure 5. The notation Δ E_patient does not seem to be defined in the main text. The authors seem to be showing effects of mutations predicted by the Potts model fitted to patient data versus their predicted effects in the specific genetic backgrounds used for mutagenesis experiments – but it's not clear if that's what they're plotting, in fact, or what we are supposed to conclude from Figure 5.

“Entrenchment” is a statement about the relative probability of seeing a drug resistance mutation as compared with the consensus residue type at a specific position, given the sequence background, which is determined by the coupling of that position to all other positions in the sequence. Our predictions of entrenchment using the Potts model are amply verified using the Potts model calculations of the change in the statistical energy ΔE as a classifier to predict the sequence backgrounds in which the drug resistance mutation is likely to be present (or absent), and then comparing the predictions with the observations from the Stanford HIV drug experienced dataset which are aggregated into groups according to the predicted values of ΔE. The comparison between our predictions and observations are presented in Table 1. We cannot emphasize too strongly that we do not fit a Potts model to observations about which sequence backgrounds are likely to favor or disfavor a DRM; the Potts model is only fit to the bivariate marginals of a very high dimensional distribution (the MSA).

A classifier of sequence backgrounds as “entrenching” or “disfavoring” towards DRMs where the predicted frequencies match the observed frequencies so well as shown in Table 1 is highly non-trivial and extremely unlikely to have happened by chance (p-value close to zero).

The subsection “Molecular clones and the effects of specific backgrounds” in the revised text presents a discussion of the results shown in Figure 5 and Figure 6. In Figure 5 we compare the predicted log likelihood (ΔE) of observing a drug resistance mutation in three specific sequence backgrounds – the molecular clones NL4-3, and HXB2, and the consensus sequence background – with the corresponding values calculated as an average over the drug experienced sequences in the Stanford HIV-1 dataset <ΔE_patients>. We make the following points: (a) the log likelihood values of a drug resistance mutation appearing in the three specific sequence backgrounds NL4-3, HXB2, and consensus, correlates moderately with the average value in the drug experienced population and NL4-3 is most representative of the drug experienced population. Figure 5 shows that the log likelihood values are negative for most drug resistance mutations in the NL4-3 and HXB2 molecular clones, and indeed most drug resistance mutations are not present in these clones.

The reviewer’s claim that “By contrast, when the authors do compare the model to actual measurements of fitness effects (subsection “Molecular clones and the effects of speci1c backgrounds”), the authors seem to be saying that the predictions of the Potts model do not match the empirical measurement” is not correct, that is not what we are saying. As stated again for emphasis, our paper is focused on an analysis of how the sequence background affects the likelihood of observing a drug resistance mutation. We verify the predictions by comparison with the prevalence of drug resistance mutations observed in the Stanford database in particular backgrounds when the statistics are aggregated according to the predicted likelihood of observing the mutation in sequence backgrounds with Potts ΔE scores in different ranges. Potts statistical energy scores ΔE are sometimes used as a proxy for organism “fitness”; this can take on different meanings depending on which experiments (e.g. thermal stability, replicative capacity, etc.) are used to interrogate for “fitness” (Wargo and Kurath, 2012). We note that many key resistance mutations are predicted by the Potts model to be disfavored in the NL4-3 and HXB2 molecular clones, this is consistent with literature reports of the effects of these drug resistance mutations on experimental probes of fitness using the NL4-3 and HXB2 clones for these experiments (Hu and Kuritzkes, 2010; Abram et al., 2013).

Figure 6 shows the distributions of Potts statistical energy differences ΔE for two IN drug resistance mutations N155H and G140S calculated over the set of drug experienced sequences in the Stanford database. It can be seen that on average the statistical energies are much more favorable in sequence backgrounds where the drug mutation appears as compared with those where it does not. For comparison the statistical energy differences for the NL4-3 and HXB2 clones are also shown. The Potts statistical energy of a sequence is commonly interpreted to be proportional to fitness. Using changes in statistical energy as a proxy for changes in fitness when a DRM is substituted for a wild-type consensus residue, the results shown in Figure 6 imply that the fitness effects of the N155H and G140S IN mutations measured through mutagenesis experiments on specific molecular clones are likely to be quite different from the effects averaged over a distribution of sequence backgrounds represented in the Stanford HIVDB.

https://doi.org/10.7554/eLife.50524.030

Article and author information

Author details

Avik Biswas
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Physics, Temple University, Philadelphia, United States
Contribution
Conceptualization, Data curation, Formal analysis, Writing—original draft, Writing—review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-3519-3944
Allan Haldane
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Physics, Temple University, Philadelphia, United States
Contribution
Conceptualization, Data curation, Formal analysis, Writing—original draft, Writing—review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-8343-1994
Eddy Arnold
1. Center for Advanced Biotechnology and Medicine, Rutgers University, Piscataway, United States
2. Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, United States
Contribution
Formal analysis, Writing—original draft

Competing interests
No competing interests declared
Ronald M Levy
1. Center for Biophysics and Computational Biology, Temple University, Philadelphia, United States
2. Department of Physics, Temple University, Philadelphia, United States
3. Department of Chemistry, Temple University, Philadelphia, United States
Contribution
Conceptualization, Formal analysis, Writing—original draft, Writing—review and editing

For correspondence
ronlevy@temple.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8696-5177

Funding

National Institutes of Health (U54-GM103368)

Ronald M Levy
Eddy Arnold

National Institutes of Health (R01-GM030580)

Ronald M Levy

National Institutes of Health (S10OD020095)

Ronald M Levy

National Institutes of Health (1R35GM132090)

Ronald M Levy

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the very supportive collaborative environment provided by the HIV Interaction and Viral Evolution (HIVE) Center at the Scripps Research Institute (http://hivescripps.edu, last accessed January, 2019). We thank Alan Engelman for his comments and suggestions for revisions of a draft of this manuscript.

Senior and Reviewing Editor

Patricia J Wittkopp, University of Michigan, United States

Publication history

Received: July 25, 2019
Accepted: September 9, 2019
Version of Record published: October 8, 2019 (version 1)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

1,010

Page views
128

Downloads
15

Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Avik Biswas
Allan Haldane
Eddy Arnold
Ronald M Levy

(2019)

Epistasis and entrenchment of drug resistance in HIV-1 subtype B

eLife 8:e50524.

https://doi.org/10.7554/eLife.50524

Categories and tags

Research organism

Virus

Share this article

Cite this article

The Potts model predicts residue frequencies.

Entrenchment and the effect of epistasis on the favorability of a primary resistance mutation.

Confirmation of the Potts model entrenchment predictions for the mutation L90M in HIV-1 protease.

Entrenchment in the population (of sequences carrying the mutation) for four major classes of resistance mutations in HIV-1 subtype B: A DRM is entrenched in the population if at least ~50% or more of the sequences containing the mutation entrench it (i.e Δ⁢E=Ew⁢i⁢l⁢d-E⁢m⁢u⁢t>0).

Table 2—source data 1

Table 2—source data 2

Table 2—source data 3

Table 2—source data 4

Degree of entrenchment of key resistance mutations occurring in the catalytic core domain (CCD) of HIV-1 IN.

Distribution of Potts Δ⁢E scores for key residues associated with drug resistance in HIV-1 IN.

Entrenchment and favorability of key resistance mutations in specific backgrounds.

The particular sequence background in which a resistance mutation occurs affects the degree of entrenchment, with often a clear distinction between sequences where the mutation is present (green) versus absent (blue) shown here for the mutations N155H (A) and G140S (B) in integrase.

Comparison of Potts models parameterized on drug-naive and drug-experienced HIV protein sequences.

Drug-pressure associated mutations are largely common between drugs of the same class.

Sequence coverage for RT.

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against the INSTIs

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against NNRTIs

Most entrenched (ME) and least entrenched (LE) or most disfavoring sequence pairs with the same Hamming distances (HD) for primary resistance mutations against PIs

Author details

Avik Biswas

Contribution

Competing interests

Allan Haldane

Contribution

Competing interests

Eddy Arnold

Contribution

Competing interests

Ronald M Levy

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Research organism

Further reading

Entrenchment in the population (of sequences carrying the mutation) for four major classes of resistance mutations in HIV-1 subtype B: A DRM is entrenched in the population if at least ~50% or more of the sequences containing the mutation entrench it (i.e $Δ E = E_{w i l d} - E m u t > 0$ ).

Distribution of Potts $Δ E$ scores for key residues associated with drug resistance in HIV-1 IN.