HIV protease, an aspartyl protease crucial to the life cycle of HIV, is the target of many drug development programs. Though many protease inhibitors are on the market, protease eventually evades these drugs by mutating at a rapid pace and building drug resistance. The drug resistance mutations, called primary mutations, are often destabilizing to the enzyme and this loss of stability has to be compensated for. Using a coarse-grained biophysical energy model together with statistical inference methods, we observe that accessory mutations of charged residues increase protein stability, playing a key role in compensating for destabilizing primary drug resistance mutations. Increased stability is intimately related to correlations between electrostatic mutations – uncorrelated mutations would strongly destabilize the enzyme. Additionally, statistical modeling indicates that the network of correlated electrostatic mutations has a simple topology and has evolved to minimize frustrated interactions. The model' s statistical coupling parameters reflect this lack of frustration and strongly distinguish like-charge electrostatic interactions from unlike-charge interactions for of the most significantly correlated double mutants. Finally, we demonstrate that our model has considerable predictive power and can be used to predict complex mutation patterns, that have not yet been observed due to finite sample size effects, and which are likely to exist within the larger patient population whose virus has not yet been sequenced.
Proteins evolve through random mutagenesis and their evolutionary selection is constrained by structural, functional and environmental factors [1]. Thermodynamic stability is by far the most important structural factor, as most proteins need to be folded in order to function. The stability range for each protein, however, is narrow and is estimated experimentally to be around 10 kcal/mol, which is of the order of three hydrogen bonds [2]. As a result of this marginal stability, proteins operate “on a knife' s edge” [3], whereby a single highly deleterious mutation could potentially lead to decreased stability and loss of activity [4]. By the same token, a single stabilizing mutation could be advantageous from an evolutionary point of view. For example, more stable forms of cytochrome P450 allowed for greater exploration of mutational space in directed evolution experiments than sequences without stabilizing mutations [5]. This increased “evolvability” is not just limited to directed evolution experiments, but may be a general property of proteins evolving under selective pressure [6]. In fact, recent experimental work on HIV protease has shown that accessory mutations compensate for the loss of stability due to destabilizing primary drug resistance mutations, helping the virus evade drugs [7]. This stabilizing effect can have an external source as well: Hsp90, a molecular chaperone, buffers deleterious mutations, allowing for polymorphisms to appear and new traits to evolve [8]. As a result of this work and prior research by other groups, it is now widely recognized that thermodynamic stability is intimately linked with the evolvability of a protein [9]–[11]. Even though the process of mutagenesis is random, the genetic and structural constraints mentioned above, coupled with functional selection, ensure that certain mutations in evolving proteins are associated with each other in a highly non-random fashion [12]. These correlated mutations are an inherent property of evolving amino acid sequences, and an evolutionary signature of viable proteins. A multitude of methods have been developed to identify such pairs and groups of mutations [13], some of which have been applied to HIV protease sequences to locate pairs or groups of coevolving residues [14]–[16]. Our previous work on higher-order correlations showed that for HIV-1 protease, including at least pair correlations is essential for reproducing statistical patterns of primary and accessory mutations observed in protease sequences from patients undergoing anti-retroviral therapy [17]. It is tempting to attribute sequence correlations to effects arising from protein stability constraints [18], and several groups have tried to connect sequence correlations with protein energetics on a detailed atomic level. For example, Zhang et al. applied Bayesian networks to infer therapeutically relevant and conditionally dependent sets of resistance mutations in HIV protease and reverse transcriptase and then used molecular simulations to model the specific interactions that cause resistance [19]. Ranganathan et al. have attempted to explain mutational coevolution by connecting statistical free energies from multiple sequence alignments to differences in experimental folding free energies [20]. Some of these results have been difficult to replicate [21], and are still a topic of active debate in the community [22], [23]. Thus while studies that link mutational correlations to thermodynamic constraints have made great progress [12], [19], [24]–[26], a consensus linking protein energetics and mutational correlation patterns has not yet emerged. These observations have motivated us to explore how correlated mutations in HIV protease are connected to energetics via their impact on protein stability. Since current methods for predicting stability changes upon mutation based on detailed atomic models are not sufficiently accurate [27], we have chosen to focus instead on the electrostatic part of the total energy for which a coarse-grained model of electrostatics is appropriate. We find that this model captures many important effects of mutations on energetics and stability of HIV protease. We show that the average electrostatic stabilization of HIV protease increases with the number of electrostatic mutations (an electrostatic mutation changes the charge of that mutating residue relative to the wild-type residue), consistent with the hypothesis that accessory electrostatic mutations buffer the destabilizing effects of primary drug-resistance mutations, most of which are non-electrostatic mutations and are therefore not modelled here. We demonstrate that correlations among electrostatic mutations are critical for stabilization; uncorrelated mutations would strongly destabilize the protein. We show that our method, which employs both electrostatic calculations and sequence analysis based on statistical inference techniques, can be used as a predictive tool for novel mutational patterns that have not yet been observed. Finally, we comment on the structure of the electrostatic mutation network of HIV protease. Energy landscape theory, which provided the framework for understanding protein folding through funnels, introduced the concept of a smooth, minimally frustrated landscape for foldable, natural proteins [28], [29]. Our results indicate that the electrostatic interaction network is minimally frustrated as is evident in the derived statistical coupling parameters which strongly predict the underlying charge patterns, providing additional evidence that proteins have evolved to minimize frustrated interactions.
Our analysis of electrostatic mutation patterns is based on the alignment of HIV protease sequences from Christopher Lee' s HIV Positive Selection Mutation Database (http: //bioinfo. mbi. ucla. edu/HIV) [30]. Each amino acid sequence in the Lee database is converted into a charge signature, which is a three letter alphabet representation of that sequence (+, −, n) corresponding to positively-charged, negatively-charged, and neutral residues. These charge signatures are compared to the wild-type charge signature to determine electrostatic mutations. We examined all primary, accessory and polymorphic drug resistance mutation positions (as designated by the Stanford HIV database [31]) and limited our analysis to a subset of 18 positions whose charged state mutates above a threshold frequency of 0. 01%. Our model therefore includes more than 380 million states or unique charge signatures involving these 18 positions (Figure 1). Of the 18,9 are sites which have been characterized as primary or accessory drug resistance mutations while the rest are sites labeled as polymorphic mutations. Mutations are labelled “polymorphic” if they are observed to mutate in the absence of drugs and whose compensatory effect has not yet been experimentally verified, even though drugs may have a significant affect on their correlations with other mutating residues [31]. If we divide this database of charge signatures into subsets with 1,2, 3… electrostatic mutations and calculate the electrostatic contribution to the average folding free energy, for each subset, we find that on average the stability of the folded state increases by kcal/mol from 1 to 3 mutations and maintains this level of stabilization beyond 3 mutations (Figure 2, black curve). Since selective pressure in the presence of inhibitors often leads to destabilizing primary drug resistance mutations [32], [33], the observed increase in electrostatic stability is due to energetic compensation: destabilizing mutations occur due to selective pressure and electrostatically active residues provide a “reservoir of stability”. The observed stabilization requires not only the correct frequencies of occurrence of each of the three possible charge states at each position, but also the presence of correlations. Generating random sequences with equal mutation frequencies for the three charge states results in a substantial destabilization of the protein (Figure 2, red curve). Introducing observed frequencies of occurrence of each charge at every position improves the stabilization relative to the previous model with equal mutation frequencies, but still results in substantial destabilization (Figure 2, blue curve). We refer to this latter model as the independent model as it generates an alignment in which mutations at each position occur independently with correct frequencies. If pair correlations are introduced by preserving the observed joint mutation frequencies (see Methods), substantial protein stabilization occurs, and the energies predicted by this pair correlation model (Figure 2, green curve) become comparable to the energies of the observed sequences. The magnitude of the difference between the observed and pair correlation model average energies is less than 2 kcal/mol for sequences with mutations, suggesting that introducing pair correlations is sufficient for explaining the observed energetic stabilization trends. Overall, the difference between the independent model and the pair correlation model is statistically significant given a sample size corresponding to the numbers of sequences observed in the database with those number of mutations (e. g. for sequences with 4 or fewer mutations). We find that the observed electrostatic stabilization can be attributed to a relatively small number of low-energy signatures which are highly unlikely under the independent model but become very probable once pair correlations are introduced (Figure S1). For example, the well-studied pair of primary drug resistance mutations, D30N-N88D [15], [32], which occurs 2220 times in the Lee database, contributes to the kcal/mol stabilization of the pair correlation model relative to the independent model shown in Figure 2. Together, these top 10 double mutants account for of the stabilization of the pair correlation model relative to the independent model. Table S3 lists the most statistically deviated pairs and Figure S7 depicts the distance between the pairs on the structure of HIV protease. It is interesting to note that 4 out of the top 10 most correlated pairs are greater than 10 Angstroms apart. With increasing numbers of mutations however, the stabilization spreads among multiple patterns (Figure S1). For 3 electrostatic mutations, the top contributor D30N-N37D-N88D is responsible for 20%, while the top 10 signatures account for 64% of the kcal/mol stabilization of the pair correlation model relative to the independent model. For 4 mutations, K20I-D30N-N37D-N88D accounts for 10%, while the top 10 signatures account for 33% of the stabilization. For 5 mutations, K20I-D30N-E35Q-N37D-N88D accounts for 17% and the top 10 signatures account for 36% of the stabilization. These stabilizing charge patterns are also strongly associated with protease inhibitor therapies, as determined by our drug association analysis (see Text S1). Most protease drug association studies focus on point mutations or pairs of mutations [30], [34]. Our drug association analysis allows us to examine the significance of drug association for patterns of more than two mutations. Table S1 lists the most significant associations between drugs and charge patterns of 2,3 and 4 electrostatic mutations with the highest probabilities as predicted by the pair correlation model. Most of the patterns listed are strongly associated with at least one drug and several are associated with many drugs. For example the D30N-N88D double mutant and the D30N-N37D-N88D triple mutant are both strongly associated with Nelfinavir monotherapy and Indinavir-Nelfinavir combination therapy with. We also find strong association between drugs and patterns predicted by the pair correlation model with more than three mutations. For example K20I-D30N-N37D-N88D and K20I-D30N-E35Q-N88D are both associated with Indinavir-Nelfinavir combination therapy while K20I-D30N-H69Q-N88D is associated with Ritonavir-Nelfinavir therapy with. The ability to predict drug resistant mutation patterns is of great therapeutic relevance. Approaches to predicting drug resistance mutations based on biophysical modeling have recently been proposed [35], [36]. In contrast, our statistical-inference based approach which includes pair correlations among mutations, allows us to predict the probabilities of arbitrary charge signatures, many of which have not yet been experimentally observed. Figure S2 shows that most of the sequences with less than 5 mutations, whose probabilities are significantly enhanced by pair correlations, are observed in the Lee database, indicating that these mutational patterns are routinely utilized by the virus. However, for 6 mutations the most stabilizing pattern, K20I-D30N-E35Q-N37D-Q58E-N88D, was not observed (Figure S1). The probability of this pattern under the pair correlation model is, too small to appear frequently in a database of sequences due to finite sample size effects. However, if the size of the database were to increase five-fold, the probability of observing at least one copy of this pattern would be. The proportion of sequences not observed in the Lee database with significantly enhanced pair correlation model probabilities increases greatly with the number of mutations, of which it is likely that many are not observed because of finite sample size effects (Figure S2). In order to test our ability to predict novel patterns of favorable electrostatic mutations unobserved in the Lee database due to finite sample size effects, we examined the contents of a separate database, the drug-annotated Stanford database which contains HIV protease subtype B sequences from various sources [31]. Figure 3 plots the probabilities of sequences using the pair correlation model, , as a function of the observed probabilities in the Lee database. Sequences are also shaded according to a gradient that represents how often the sequence occurs in the Stanford database, relative to the Lee database. The plot shows that sequences with the highest predicted probabilities that are unobserved in the Lee database are largely shaded red, indicating that most are observed in the Stanford database. In fact, of the top 25 most probable sequences predicted by the pair correlation model that are not found in the Lee database, 19 are present in the Stanford database (Table S2). Most of these sequences are also significantly associated with drug therapies (e. g.). As the predicted probability decreases, the shading of the dots gradually changes to green, indicating that sequences with the lowest predicted probabilities are unobserved in both the Lee and the Stanford database. Thus our approach exhibits considerable predictive power. If we examine other sequences in the tail of the Lee probability distribution that are observed once, twice, three times etc in the Lee database, we notice a similar trend; sequences with higher predicted probabilities are present in the Stanford database at much higher frequencies than sequences with lower predicted probabilities, even though the pair correlation model was parameterized on the Lee database. This result strongly suggests that the pair correlation model is a much better predictor of actual sequence probabilities than using the sequence counts from the databases themselves, because of finite sample size effects. In other words, the tail of the distribution is very well represented by the pair correlation model. Determining the statistical field and coupling parameters (written as and for simplicity) that best fit the pair correlation model given a set of observed univariate and bivariate marginals (and), is known in the literature as the inverse Ising problem. As described in the Methods, we iteratively determine these parameters using a graph-theoretic inference algorithm, called belief propagation (BP) [37]–[39], a method which has recently been applied by other research groups to study protein conformational entropy and protein-protein interactions [40], [41]. Within the BP framework, we apply a mean field model which includes pair correlations in the Bethe approximation to estimate the bivariate marginals, [38], [42], [43]. In the statistical physics community “mean field” in often used to refer to a class of approximations whereby the free energy of the system is written in terms of the marginals up to a given order, the corresponding mean field model at the pair correlation level is the Bethe approximation. However, it is well known that while is exact and converges to on acyclic networks, only approximates and can become unstable on cyclic networks [44]. For the electrostatic correlation network of HIV protease, we observe that the belief propagation algorithm converges quickly to the observed bivariate marginals (Figure S3). The convergence towards the observed probabilities for triplets and larger multiplets is also well approximated, a result that is non-trivial since the Bethe approximation is a pair-level approximation and does not guarantee the convergence for marginals beyond pairs [45]. For trivariate marginals, the correlation coefficient between and is 0. 98 while for four mutations, the correlation coefficient between and is 0. 90 (Figure S3). This close correlation between observed and predicted marginals argues for a simple network structure and suggests that for this system, the Bethe approximation is a good approximation and by implication the electrostatic mutation network is minimally frustrated. Another strong indicator of the lack of frustration in the electrostatic mutation network of HIV protease is that the statistical coupling parameters of the (Bethe) mean field model are able to distinguish like-charge patterns from unlike-charge patterns with high accuracy (Figure 4). We find that the sign of, a quantity derived from the sequence analysis alone, is able to correctly predict the charge patterns for of the top 35 most significantly correlated charge pairs (.), reflecting the evolutionary optimization of the protein electrostatic interaction network (Figure 4). Moreover, Figure 4 also indicates that the magnitude of correlates with the spatial distance between residues. Of the top ten pairs of residues with the largest statistical coupling parameters, nine are situated close to one another (Å) in the folded structure of the HIV protease homodimer. In this context, we note that Morcos et al. [46] used a similar approach to infer spatial contacts between residues of many proteins, through an analysis of the coupling parameters of a corresponding mean field model.
Our results suggest that electrostatic interactions play an important role in the coevolution of mutations in HIV protease. The extent to which electrostatics influences protein stability has been the subject of debate in the literature [47]–[53]. One experimental study of protease has minimized the role of electrostatics in favor of the impact of compensatory mutations on protein flexibility [54]. Others have suggested that buried charges play a more important role in protein function than stability [55], [56]. We recognize that electrostatics is only part of the total energy, and that contributions to stability from van der Waals interactions, hydrogen bonding and hydrophobic effects are significant. Nonetheless, long-range electrostatics is likely to have a substantial effect on protein stability [50], [52], [57]–[59]. Our results support the proposal that the presence of correlations among electrostatic mutations arises from the contraints imposed by the need to maintain the stability of the folded protein. HIV protease is under strong selection pressure from drugs. As a result of the initial build up of drug resistance, protease becomes less stable [32]. We hypothesize that electrostatic mutations not only bring the protein back to a more viable state, but may give the protein more “breathing room” on the evolutionary fitness landscape. Manipulating the charge distribution of HIV protease is complex and we find that uncorrelated mutations would tend to strongly destabilize the enzyme, contrary to the stability gain observed in the database. Therefore, we propose that sets of electrostatic mutations occur together, increasing the “evolvability” of a protein by providing a “reservoir of stability” which allows it to escape epistatic traps along evolutionary pathways towards higher fitness [18]. The absence of frustration could reflect evolutionary optimization of the electrostatic interaction network in HIV protease under selection pressure from drugs, or it could be a general property of protein electrostatic interaction networks. Indeed, natural proteins tend not to be frustrated systems [29], [60] – they are fine-tuned biological machines with restricted evolutionary pathways [18]. Within these pathways, proteins are highly robust and the physics underlying their folding display a kind of simplicity [29], [61]–[63]. Our conclusions based on a coarse-grained electrostatics model combined with statistical inference techniques reflect this lack of frustration. In future work we will study our algorithm on interaction networks with a larger alphabet size and different network topologies. In this context, we note that Balakrishnan et al. used an alternative learning algorithm that a solved similar problem but determined the optimal graph topology for the network [64]. Our statistical analysis of HIV sequences captures biophysical constraints in the form of a statistical network of correlated mutations. Even though the model is based on pairwise correlations, it captures the higher-order effects and correctly predicts the probabilities of sequences found in the tail of the distribution. The fact that many of these patterns are also strongly associated with protease inhibitors from patients undergoing antiretroviral therapy, highlights the clinical relevance of our method. Other mutation patterns that we predict are likely to exist within patients whose virus has not yet been sequenced. Having knowledge of these unique, but as yet unobserved patterns, can be important for the design of future inhibitors to combat drug resistance. In this work, we go beyond a purely statistical approach to modeling patterns of electrostatic mutations, and show that the statistical results are entirely consistent when viewed in the context of a structure based energy model. Though electrostatics is only part of the total energy, our work has highlighted its importance and provided support for the proposal that correlated electrostatic mutations provide a reservoir of stability for HIV protease as it builds resistance to drugs.
45,161 aligned HIV-1 DNA nucleotide sequences were downloaded from Christopher Lee' s HIV Positive Selection Mutation Database (http: //bioinfo. mbi. ucla. edu/HIV) on March 4th, 2008 [30]. This database of sequences, which we call the Lee database, consists primarily of HIV-1 subtype B samples sequenced by Specialty Laboratories Inc from 1999 to mid-2002. These sequences are not annotated. [30]. The amino acid sequences were converted into strings of characters “n”, “−” and “+”, indicating whether a given residue is neutral, negatively or positively charged at pH = 6 (i. e. His, Arg, and Lys are positively charged, while Asp and Glu are negatively charged). A second database of subtype B sequences, which we call the Stanford database, was downloaded from the Stanford HIV database on April 7th, 2010 [31]. This database consists of drug annotated sequences collected primarily from more than 900 literature and GenBank references. This drug-annotated dataset was used to associate correlated mutation patterns with specific anti-retroviral therapies (see SI). The univariate and bivariate marginals extracted from the Lee database and from the Stanford database are effectively the same (correlation coefficient of 0. 999), indicating that our results would be unchanged if we used the Stanford database to parameterize the model. Moreover, at a 20 letter amino acid alphabet level, there is little redundancy between the databases as only 7. 63% of the sequences are present in both databases. To locate electrostatic mutations, the resulting charge signatures were compared to the HIV-1 subtype B consensus sequence from the Los Alamos National Laboratory HIV sequence database. This consensus sequence was used to define the wild-type charge signature. We define an electrostatic mutation as an amino acid mutation which changes the charge at a certain position along the protein sequence, relative to the wild-type amino acid at that position (e. g. D30N and N88D). In contrast, the L90M and R8K are not considered to be electrostatic mutations. We examined all the primary, accessory and polymorphic drug resistance mutations positions (as designated by the Stanford database [31]) and included all electrostatic mutations above a threshold frequency of 0. 01%. The 18 positions included are the primary drug resistance mutatation sites D30 and N88, accessory mutation sites K20, E34, E35, K43, Q58, L63, and Q92, and polymorphic mutation sites Q7, T12, G16, Q18, N37, Q61, H69, K70, and I72. The electrostatic energy of protein folding was estimated as (1) where and are electrostatic free energies of the folded and unfolded states, computed using an Analytical Generalized Born (AGB) model [65]. The folded state electrostatic free energy was calculated by placing unit charges corresponding to a particular charge signature onto the most-distal sidechain carbon atom of the corresponding wild-type amino acid within a dimer crystal structure (PDB ID 1NH0). All other sidechain atoms remain neutral, although a partial charge dipole of is placed on every backbone amide and carbonyl group to retain the helix dipole effects [66]. Our approximation of the denatured state is a maximally extended structural representation of chain A from 1NH0, with backbone dihedral angles set to (except for prolines) and sidechain rotamer states set to all-trans. Similarly to the folded state, charges on the unfolded state are placed on the most-distal sidechain carbon atom and backbone dipoles are switched on. See SI Methods for further implementation details. As in our previous work [17], we make use of a Potts model to capture the effects of pair interactions between residues. Since our electrostatic data consists of sequences with three possible charge states at each site, we use a 3-letter alphabet (+, −, n), for positively-charged, negatively charged, and neutral residues. Including all three charge states in our study leads to possible charge signatures for positions. For each signature we calculate, the independent model probability, and, the pair correlation model probability. Specifically, we fit the frequencies of charge states at each position and the joint frequencies of charge states at pairs of positions to the 3-state Potts model: (2) where is a sequence of' s,' s or' s of length, are position indices, and are the fitting parameters for the fields and couplings, and is the partition function. The independent model, obtained by setting all, corresponds to. For the equal frequency model in Figure 2, we set. The joint probability distribution given by the Potts model has the largest entropy constrained by the univariate (independent model) or both univariate () and bivariate () marginals from the data [67]. To solve the inverse Ising problem, we implemented an efficient graph-theoretic inference algorithm called belief propagation (BP) [37]–[39]. Our algorithm employs a two-step procedure: first, all the univariate and bivariate marginals are determined for a given set of and in the Bethe approximation within BP [38], [42], [43]. Second, the predicted marginals are compared to the observed marginals to determine updated and via gradient descent [39]. See the supporting information for further information about the algorithm.
HIV is incurable because its enzymes evolve rapidly by developing resistance mutations to retroviral inhibitors. Most of these mutations work synergistically, but the biophysical basis behind their cooperation is not well understood. Our work addresses these important issues by bridging the gap between the statistical modeling of HIV protease subtype B sequences with the energetics of mutations involving charged amino acids by showing that electrostatic stability is intimately related to correlations. Moreover, we demonstrate that our statistical model has considerable predictive power and can be used to predict complex mutation patterns that have not yet been observed due to the finite sizes of the current sequence databases. In other words, as the database size increases, our model has the ability to predict the identities of the high probability mutations patterns, which are more likely to be observed. Knowing which currently unobserved mutations are more likely to be observed can be very advantageous in combating the disease.
medicine
protein interactions
protein folding
protein structure
sequence analysis
infectious diseases
proteins
hiv
biology
biophysics
physics
biochemistry
protein chemistry
viral diseases
computational biology
macromolecular structure analysis
Correlated Electrostatic Mutations Provide a Reservoir of Stability in HIV Protease