Dynamical basis for drug resistance of HIV-1 protease.

BACKGROUND: Protease inhibitors designed to bind to protease have become major anti-AIDS drugs. Unfortunately, the emergence of viral mutations severely limits the long-term efficiency of the inhibitors. The resistance mechanism of these diversely located mutations remains unclear.
RESULTS: Here I use an elastic network model to probe the connection between the global dynamics of HIV-1 protease and the structural distribution of drug-resistance mutations. The models for study are the crystal structures of unbounded and bound (with the substrate and nine FDA approved inhibitors) forms of HIV-1 protease. Coarse-grained modeling uncovers two groups that couple either with the active site or the flap. These two groups constitute a majority of the drug-resistance residues. In addition, the significance of residues is found to be correlated with their dynamical changes in binding and the results agree well with the complete mutagenesis experiment of HIV-1 protease.
CONCLUSIONS: The dynamic study of HIV-1 protease elucidates the functional importance of common drug-resistance mutations and suggests a unifying mechanism for drug-resistance residues based on their dynamical properties. The results support the robustness of the elastic network model as a potential predictive tool for drug resistance.

Background

HIV-1 protease (human immunodeficiency virus type 1 protease) is an enzyme that plays a critical role in the virus replication cycle. It cleaves the gag and pol viral polyproteins at the active site to process viral maturation [1-3], and without HIV-1 protease the virus was found to be noninfectious [4]. Thus HIV-1 protease is widely considered the major target for AIDS treatment [5,6]. One of the most severe obstacles to protease-inhibiting drugs is the rapid emergence of protease variants. Variants are able to evolve resistance by developing a chain of mutations, and as a result limit the long-term efficiency of these drugs [7,8].

HIV-1 protease is a dimer of C2 symmetry with each monomer consisting of 99 amino acid residues. Each monomer has one α helix and two antiparallel β sheets in the secondary structure. The enzyme active site is a catalytic triad composed of Asp25-Thr26-Gly27 from each monomer. It is gated by two extended β hairpin loops (residues 46−56) known as flaps [9]. At the molecular level, resistance to protease inhibition predominantly takes the form of mutations within the protein that preferentially lower the affinity of protease inhibitors with respect to protease substrates, while still maintaining a viable catalytic activity [10]. Mutations associated with drug resistance occur within the active site as well as non-active distal sites [11].

During the past two decades, researchers and clinicians from different disciplines have made enormous efforts to investigate resistance against HIV-1 protease targeted drugs. To elucidate the molecular mechanisms of drug resistance, biochemists and molecular biologists have characterized the structure, energetics and catalytic efficiency of a large number of HIV-1 protease mutants to unravel the resistance mechanism in combination with extensive computational studies [12-15]. Moreover, drug resistance data collected from AIDS patients treated with HIV-1 protease inhibitor drugs [16-19] provide opportunities for researchers to identify resistance-related mutation patterns [20-22]. Recently there have been efforts to link protein physical and functional stability with its evolutionary dynamics [23,24].

At the heart of understanding the molecular basis of drug-resistant behaviors of HIV-1 protease is the structural distribution of resistance mutations. Presumably these mutations are not randomly located throughout the protein structure. Although different HIV-1 protease inhibitors elicit different combinations of mutation types to generate distinctive resistance levels, there are 21 most common mutations associated with resistance against all inhibitors [19]. Prediction of resistance mutations of proteins is based on either sequence or structure information [25]. Sequence-based methods predict resistance mutations by analyzing large datasets of sequences with known resistance properties. Thus the availability of those datasets is a prerequisite for such methods [22,26-28]. On the other hand, predicting mutations using protein structure has largely relied on the characterization of binding thermodynamics [29-32], as the mutations with resistance against inhibitors lower the binding affinity of inhibitors far more than that of natural substrates. The accuracy of the prediction is directly related to the accuracy of the potential function used in the calculations and the adequacy of the sampling of the protein conformational space. It is also sensitive to the error/noise in the free energy calculations [32].

Conformational dynamics play an essential role in regulating protein function [33,34]. In the past few years a deepening understanding of the relationship of protein dynamics and function has emerged [35]. Relevant to the study here is the utilization of protein dynamics to identify the sequence regions of functional importance even though their locations may be remote from the active site. Computationally there have been rapid methodological developments in relating protein dynamics to function by probing the long range communications between residues: perturbation method [36,37], clustering analysis of correlation matrix [38], network analysis [39], and energy diffusivity estimation by propagation through vibrational modes [40]. The success of these methods in reproducing experimental results as well as findings from sequence-based methods has established the validity of dynamics-based approaches [38,41].

The dynamics of HIV-1 protease, especially binding dynamics of its ligands are fundamental to the protease inhibitor design and have been a subject of intense computational study [42-49]. Because of limitations of time scale in all-atom simulations, various coarse-grained models have been used to investigate HIV-1 protease binding dynamics and kinetics, shedding light on important dynamics issues [45-49]. The main features of substrate interactions and dynamics at the active site were analyzed within the framework of the coarse-grained model [45,49]. Gaussian models were shown to describe accurately the correlated motion of HIV-1 protease residues in thermodynamic equilibrium through a series of successful comparisons with an extensive MD simulation [49,50]. There is increasing evidence relating protease's drug-resistance mutations to its dynamics. The impact of some distal mutations on catalytic function of HIV-1 protease was linked to protein flexibility [51,52]. Multi-drug resistance residues of HIV-1 protease were found to overlap the global hinge region identified from coarse-grained normal-mode analysis of the protease [53]. Nevertheless, despite extensive research efforts, a general explanation for drug-resistance mutations of HIV-1 protease is still lacking [54].

In this study, a coarse-grained elastic network model is used to investigate the dynamics of HIV-1 protease, to probe the connection between its global dynamics and the distribution of drug-resistance mutations, and to examine the potential of the dynamics-based approach as a predictive tool for drug resistance prediction, with an attempt to provide a unifying mechanistic explanation for all residues of resistance based on their dynamical properties. The crystal structures of an unbound form and bound forms with a substrate and nine FDA approved inhibitors of HIV-1 protease are used as model systems. Correlation analysis of the protease at equilibrium focuses on two functional sites of HIV-1 protease: the active site (the Asp25-Thr26-Gly27 triad [45]) and the flap (residues 45-55 [47]). The protease dynamic changes upon ligand binding are examined as well. The implications for resistance mechanisms and protein evolution are discussed.

Results

Here HIV-1 protease is represented by a coarse-grained network model, and its dynamics is examined in several X-ray crystallographic structures. The linkage between global dynamics and the distribution of drug-resistance mutations is examined first in individual unbound and bound forms, then in the dynamical differences between the unbound and bound forms. The former is a measure of the residual fluctuations in different structures, and the latter is an estimate of dynamical change caused by ligand binding.

Dynamically coupled regions identified by equilibrium correlations Unbound form

The correlation matrix, consisting of correlations of all residue pairs, captures the essence of protein dynamics. A 198 × 198 correlation matrix (Figure 1) is generated from the elastic network modeling (see Methods section) based on the unbound HIV-1 protease (PDB id 1HHP). The most conspicuous features in the figure are the beta-sheets, represented by the line across the diagonal (residues 19-24, 43-66, 69-78, 118-123, 142-165 and 168-177). Extraction of further information from the matrix requires the application of data analysis tools such as clustering algorithms. The Markov cluster (MCL) algorithm (see Methods section) is the chosen method, and the results from the MCL program for the unbound form consist of five clusters (Table 1). Cluster 1 is the largest cluster with 48 residues. It covers the core domain, which excludes both termini and the flap, and forms the scaffold surrounding the active site of the enzyme. Clusters 2 and 3 are exclusively composed of the residues located at the beginning of the N-terminus and the end of the C-terminus, respectively. The remaining two clusters concentrate on two functional sites (Figure 2). Cluster 4 contains the active site of HIV-1 protease (Asp25-Thr26-Gly27 catalytic triad), one residue at the dimerization region (Leu10), and a large segment of the C-terminal (Ile84-Cys95). Cluster 5 is mainly made up of the residues at the flap region (Pro44-Arg57) with two additional residues from the C-terminal domain (Leu76 and Gly78). These two functional sites are involved in distinct aspects of the protease function: the active site is where the enzymatic catalysis takes place, and the flap controls access to the active site. The clustering results indicate that these two functional sites belong to separate interaction networks of the protease.

Figure 1

Correlation matrix. The correlation matrix for an unbound form of HIV-1 protease (PDB id 1HHP).

Table 1

Clustering result of the unbound form

Scaffold	11 12 13 14 15 16 17 18 19 20 21 22 32 33 34 35 36 37 38 39 40 41 42 43 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 77 79 80 81 82 83
N-terminal	1 2 3 4 5 6 7 8 9
C-terminal	96 97 98 99
Coupled with the active site	10 23 24 25 26 27 28 29 30 31 84 85 86 87 88 89 90 91 92 93 94 95
Coupled with the flap region	44 45 46 47 48 49 50 51 52 53 54 55 56 57 76 78

Residues of the active site and of the flap are underlined and double-underlined, respectively.

Figure 2

Two clusters containing functional sites. Residues in Clusters 4 (red) and 5 (yellow) are displayed in an unbound form of HIV-1 protease (PDB id 1HHP). The active site is highlighted by a ball representation.

Bound forms

The ten structures of HIV-1 protease in complex with the natural substrate and nine FDA approved inhibitors are chosen to represent HIV-1 protease in bound forms. Despite their varying degrees of similarity with the unbound form, these ten structures generate very similar clustering results (Table 2). The general clustering pattern of HIV-1 protease unbound form is preserved in all the bound forms: the scaffold, the N- and C- termini, the active site and the flap. Nevertheless, inside the clusters there are reorganizations and splits of the clusters due to ligand binding. Upon binding, the tips of the flaps (residues 48-55) close and cover the active site, and the physical proximity facilitates stronger interactions between the flap region and part of the C-terminal. As a result, the cluster containing the flap grows in the bound form. The cluster containing the flap is enlarged by including residues 79-83, while part of the C-terminal (residues 87-95) is disengaged from the cluster with the active site and forms a new cluster of its own. The diversity of the ligand types of the ten structures does not induce a dramatic impact on the clusters. The number of residues in the cluster of the active site and of the flap ranges from 12 to 18 and from 19 to 24 respectively, but the number of separate sequence segments in these two clusters remains constant.

Table 2

Clustering results of the ten bound forms

PDB id	Mutation	RMSD (Å)	Coupled with the active site	Coupled with the flap region
2FNS 7(substrate)	Q7K, D25N,L63P, I64V	0.52	1023 24 25 26 27 28 29 30 31 8485 86 87	44 45 464748 49 50 51 52 5354 55 56 78 79 80 81 8283
1HPV 6(Amprenavir)	None	0.54	1023 24 25 26 27 28 29 30 31 3285 86	44 45 464748 49 50 51 52 5354 55 56 78 79 80 81 8283 84
1HXB 6(Saquinavir)	None	0.31	1023 24 25 26 27 28 29 30 31 3285 86	44 45 464748 49 50 51 52 5354 55 56 78 79 80 81 8283 84
1MUI 5(Lopinavir)	N37S	0.53	1023 24 25 26 27 28 29 30 31 328485 86	44 45 464748 49 50 51 52 5354 55 56 78 79 80 81 8283
2O4K 8(Atazanavir)	Q7K	0.35	23 24 25 26 27 28 29 30 31 32 83 8485 86 87 88 89 90	43 44 45 464748 49 50 51 52 5354 55 56 57 58 7677 78 79 80 81 82
3JVY 6(Darunavir)	Q7K, L33I*, L63I, C67A, G86A, C95A	0.38	23 24 25 26 27 28 29 30 31 3285 86	44 45 464748 49 50 51 52 5354 55 56 79 80 81 8283 84
2O4P 8(Tipranavir)	Q7K	0.39	23 24 25 26 27 28 29 30 31 32 83 8485 86 87 88 89 90	44 45 464748 49 50 51 52 5354 55 56 79 80 81 8283 84
1HXW 6(Ritonavir)	None	0.40	1023 2425 26 27 28 29 30 31 85 86 87	44 45 464748 49 50 51 52 5354 55 56 79 80 81 8283 84
2PYN 6(Nelfinavir)	D30N, A71V*	0.26	1023 24 25 26 27 28 29 30 31 328485 86	44 45 464748 49 50 51 52 5354 55 56 78 79 80 81 8283
2B7Z 7(Indinavir)	K20R*, V32I*, L33F*, M36I*, M46I*, L63P, A71V*,V82A*, I84V*, L90M*	0.32	10 22 23 24 25 26 27 28 29 8485 86	34 43 44 45 464748 49 50 51 52 5354 55 56 57 7677 78 79 80 81 8283

The type of the ligand is listed in brackets. The global average RMSD of backbone atoms between the unbound (PDB id 1HHP) and bound forms is given. The mutations of each protease are listed. Only the clusters containing the active site and the flap are shown here.

(1) *drug-resistance mutations

(2) Residues belonging to the 21 common resistant mutations are underlined.

(3) Residues shared by all the bound forms are in italics.

In summary, there exist two independent clusters containing two important functional sites of HIV-1 protease. These two networks are relatively robust to perturbations caused by the different types of ligand.

Dynamical behavior differences between the unbound and bound forms

Without directly engaging the active site, another way to influence a protein's function is to perturb the motions essential to its function. For enzymes, these essential motions are the conformational changes accompanied by the association and dissociation of ligands [55]. In HIV-1 protease the large scale open/close conformational change of the flap along the reaction pathway is the major structural reorganization induced by ligand binding [43-47]. Extensive investigations of 73 X-ray mutant and complex structures of HIV-1 protease revealed that a common and predominant dynamic behavior was found among the protease in complex with different ligands [56]. The focus of study here is the dynamical changes from the unbound to bound forms of HIV-1 protease at the residue level. The dynamical property of each residue is characterized by the sum of its couplings with other residues (see Methods). The binding-induced changes in the dynamical behaviors of individual residues were very similar among the ten bound structures (Figure 3). Overall, the perturbations to the residues are not isotropic, and the regions exhibiting the largest deviations are signaled by peaks. The ten bound proteases all share very similar locations and magnitudes of these peaks (residues 32-42, 44-57 and 77-82). The most pronounced difference between the unbound and bound protease is the loss of flexibility in the flap tip upon binding, as indicated by the highest peak (around residues 44-57) in Figure 3.

Figure 3

Fluctuation differences between the unbound and bound forms. Fluctuation difference between the unbound (PDB id 1HHP) and ten bound forms of HIV-1 protease plotted as a function of residues. Only residues with larger than average difference are plotted.

Discussion

Using a coarse-grained network model, functionally important residues of HIV-1 protease were identified based on correlation analysis of either equilibrium fluctuations or dynamical changes. Experimentally the residues of functional importance can be directly probed by mutagenesis experiments. The complete mutagenesis experiments carried out on HIV-1 protease revealed three mutationally sensitive sequence domains [57]: the active site region (Ala22-Leu33), the flap region (Ile47-Gly52) and the part of C-terminal (Thr74-Arg87). The active site and flap region have obvious functional significance. The influence the C-terminal residues exert on the catalytic cycle is most likely through long-range couplings with the functional sites. Correlation analysis suggests that some residues of the C-terminal (Ile84-Cys95) are functionally important due to their strong interactions with the active site, while the others (Leu76 and Gly78) may influence the dynamics of the protein through their couplings with the flap. The dynamical changes analysis also indicates that residues 77-82 are involved in the binding process. Structurally, residues 79-84 form a wall of the active site, and their motions were shown by previous simulations to correlate well with the open/closed dynamics of the flap [52].

Currently there are nine FDA approved protease inhibitors, and the most up-to-date clinical data indicate 21 most common drug-resistant mutation positions ("most common" defined as mutation position shared by at least two inhibitors among the nine currently approved protease inhibitors):10, 20, 24, 32, 33, 36, 46, 47, 48, 50, 53, 54, 62, 71, 73, 76, 77, 82, 84, 88 and 90 [19]. These mutations are located at the dimer interface (residues 88 and 90), the core domain (residues 10, 20, 24, 32, 33, 36, 62, 71, 73, 76, 77, 82, and 84), and the flap domain (residues 46, 47, 48, 50, 53, and 54), respectively. The clustering analysis of HIV-1 protease in various forms identified 10, 24, 32, 46, 47, 48, 50, 53, 54, 76, 77, 82, 84, 88 and 90 (15 out of the total 21 resistance sites) as those coupled either with the active site or the flap. The dynamical change analysis identifies two additional residues 33 and 36 as residues of importance. The undetected residues (20, 62, 71 and 73) are located in the hydrophobic core of the protease, and their mutations likely affect the protease activity through influencing the protease structure and stability [30]. Another noteworthy fact is that there are far more residues discovered by clustering analysis than by clinical studies (15 out of the total 34 residues in the two clusters containing functional sites are the residues of resistance). Drug resistance residues not only influence protease binding, but also generate differential binding affinity between the substrate and inhibitors. The cluster analysis can only detect residues that may influence the binding. Therefore the pool of residues identified by the clustering analysis is larger than those found with drug resistance affinity in clinical studies. The similar clustering results from different protease inhibitor complexes further suggest that global dynamics are preserved among different complexes of HIV-1 protease inhibitors. The convergence of results regardless of the conformation of the protein was also found by the Gaussian network modeling of HIV-1 protease [50]. Computational studies using atomistic MD simulations reached the same conclusion that correlation matrix based analysis does not differentiate the essential modes of motion of the protein native forms from those of the mutants [58]. Nevertheless, the dynamic approach proposed here locates a majority of drug-resistance mutations, and provides insight on the drug-resistance mechanism of HIV-1 protease. Drug resistance mutations of HIV-1 protease can be classified as active site or non-active site mutations, depending on their location within the protein. Active site mutations are located in the vicinity of the active site and directly affect the protease-inhibitor interactions. Thus their action on inhibitor binding affinity can be readily understood in structural terms. On the other hand, non-active site mutations influence binding from various distal locations and their mechanism of action is not immediately apparent. Although some residue-specific explanations and suggestions have been proposed, the overall mechanism by which these diversely located non-active site residues influence inhibitor-binding remains unclear [19]. Dynamic studies presented here suggest a simple yet general explanation for the distribution of the drug resistance residues in HIV-1 protease. Except for these residues of resistance influencing protease binding affinity by acting on the structural stability of HIV-1 protease, the drug-resistance residues belong to clusters that are either the "coupled with the active site" or "coupled with the flap". It is noteworthy that residues coupled with the functional sites (Clusters 4 and 5) do not all locate in the physical vicinity of each other. Long distance communications play an important role in mediating the interactions between the functional sites and distal residues. The residues clustered with the functional sites can exert an impact on protein function, even though they may not locate near the active site. Leu10 and Leu90, among the residues clustered with the active site (Cluster 3), are well-known residues whose mutations lead to drug resistant variants [17]. The finding that Leu76 and Gly78 are coupled with the flap region is also corroborated by a detailed protein flexibility study of HIV-1 protease which concluded that mutations at residues Thr74-Val77, although far from the active site, reduce protease activity because of their correlated motions with the flap region [55]. Traditionally, computer methods of predicting resistance mutations based on protein structure have been largely focused on energetic analysis [25], in which atomistic molecular mechanics and/or molecular dynamics are used to investigate ligand-protein binding affinity. The mechanism of non-active site mutations has largely remained a challenge for energetic analysis because of the minimal structural and energetic perturbations caused by those mutations [51]. These difficulties, however, open the door for dynamics-based study, especially via coarse-grained methods such as elastic network modeling, which provide an efficient way of sampling the global dynamics of proteins. The interaction network identified by the correlation analysis contains both active site and non-active site residues. These residues influence the inhibitor binding by coupling with the active site and the flap, regardless of their physical proximity. Thus the dynamic approach in this study is able to detect both active and non-active drug-resistance sites based on coarse-grained protein models.

It is believed that only a few amino acid sites are responsible for adaptive evolution in almost all proteins [59], although the nature of these positively selected residues is yet to be elucidated. The findings from this study indicate that the interaction networks of globally distributed residues involving the functional sites play a dominant role in the evolutionary pathways of HIV-1 protease, and they become the major sites that develop resistance under selective antiviral pressure. Pathogenic proteins such as HIV-1 protease escape the challenges to their survival imposed by drug inhibition through mutations at these amino acid residues. Structure-based modeling of proteins confirms the decisive role of physical interactions in the evolution of virus proteins and raises the possibility of constructing a protein fitness landscape by means of physical modeling of proteins.

Conclusions

This study examines the functional significance of common drug-resistance mutations of HIV-1 protease by characterizing its global dynamics using coarse-grained modeling. The calculations show that most residues of drug-resistance are coupled either with the active site or with the flap. These couplings are rather robust to the perturbation of ligand binding. These findings result in a unifying mechanism for all drug-resistance residues based on their dynamical properties. They also indicate that global dynamics of HIV-1 protease are intrinsically connected to the structural distribution of drug-resistance mutations, thus dynamic study provides a simple yet general and useful tool to examine the tendency of drug resistance of residues in addition to traditional energetic analysis.

Methods

Elastic network modelling

The elastic network model was applied according to the standard protocol [60]. The details of the correlation matrix can be found elsewhere [41]. In this study the cutoff distance Ris set to be 10Å, but the choice of Rwas shown not to noticeably affect the results based on the correlation matrix generated from the model [61]. The structures of the unbound and ligand-bound forms of HIV-1 protease with various ligands were used. For the ligand-bound protein, only the Catoms of the protein are represented by the network model and the ligand is not incorporated in the model. The online server http://ignmtest.ccbb.pitt.edu/cgi-bin/anm/anm1.cgi[60] was used to generate the correlation matrix, and the element in the correlation matrix is defined as

where ΔRand ΔRare the fluctuations of nodes i and j, respectively.

Clustering analysis

The correlation matrix is submitted for clustering analysis by the Markov cluster (MCL) algorithm [62]. The MCL algorithm is one of the most successful clustering procedures in identifying protein-protein interactions from genomic data [63] and has been shown to be robust and outperform other clustering algorithms [64]. Relevant to the study here is the application of MCL to clustering protein residues based on the interaction correlation matrix [38]. MCL finds cluster structure in graphs by performing a random walk through the graphs. The process computes the probabilities of random walks through the graph, and uses expansion and inflation to change the probabilities associated with the random walks departing from one particular node. It results in the separation of the graph into different segments. Cluster granularity is controlled by the inflation parameter which is the only variable in the MCL program used in this study. In order to reduce the noise, the correlation matrix has to be adjusted before being submitted to the MCL program. First the absolute values of the correlation coefficients are taken. Then a cutoff value is applied to produce a condensed version. Correlations less than the cutoff value are set to zero and the cutoff value is subtracted from the remaining correlations. The cutoff and the inflation parameter (0.08 and 1.5, 0.075 and 1.4 for the unbound and bound protease, respectively) were chosen to produce a total number of five clusters for the unbound protease, and six to seven clusters for the bound protease. The number of clusters is chosen to be five or six because the resulting clusters make the most physical sense. All the bound forms are subject to the same parameters.

Fluctuation difference between unbound and bound forms

The structural differences between the unbound (ligand-free) and bound (ligand-bound) proteins are usually significant. The structural changes are not to be identified with the dynamical behavior changes. The change in dynamical behavior caused by binding for residue i, ΔC, is calculated as the sum of the absolute values of the difference between the correlation coefficients of residue i in the unbound and bound forms,

In Equation 2, the sum is over all the residues in the protein. Cand C'denote the correlation coefficient of residues i and j in the unbound and bound forms, respectively.

Authors' contributions

YM is responsible for all the work related to the manuscript.