Effect of Flap Mutation I54L/M in Inhibition of

Human immunodeficiency virus type 1 (HIV-1) encodes a 22 kDa homodimeric aspartic protease that cleaves the gag and pol viral polyproteins at nine specific sites and thus permits viral maturation (Kohl et al., 1988). HIV-1 protease is one of the most extensively studied enzymes, both experimentally and computationally. Hence, the invention of drugs that restrict proteolytic processing by the protease is a prime goal in the treatment of HIV-1 infection. Nine commercially available drugs are based on the inhibition of this enzyme to treat acquired immune deficiency syndrome (AIDS) (Wensing et al., 2010). The clinical success of these drugs is hindered by the frequent occurrence of drug-resistant mutations of the HIV-1 protease. The resistant strains encode mutated proteases with reduced binding affinity for the inhibitors but with sufficient enzymatic activity to process the viral polyproteins. However, the emergence of active, drug-resistant mutants of the protease has limited the long-term effectiveness of these inhibitors. Thus, the design of HIV-1 protease inhibitors with novel mechanisms of action and reduced resistance liabilities is still desirable.


Introduction
Human immunodeficiency virus type 1 (HIV-1) encodes a 22 kDa homodimeric aspartic protease that cleaves the gag and pol viral polyproteins at nine specific sites and thus permits viral maturation (Kohl et al., 1988). HIV-1 protease is one of the most extensively studied enzymes, both experimentally and computationally. Hence, the invention of drugs that restrict proteolytic processing by the protease is a prime goal in the treatment of HIV-1 infection. Nine commercially available drugs are based on the inhibition of this enzyme to treat acquired immune deficiency syndrome (AIDS) (Wensing et al., 2010). The clinical success of these drugs is hindered by the frequent occurrence of drug-resistant mutations of the HIV-1 protease. The resistant strains encode mutated proteases with reduced binding affinity for the inhibitors but with sufficient enzymatic activity to process the viral polyproteins. However, the emergence of active, drug-resistant mutants of the protease has limited the long-term effectiveness of these inhibitors. Thus, the design of HIV-1 protease inhibitors with novel mechanisms of action and reduced resistance liabilities is still desirable.
An important structure feature of the aspartic proteinases is the short antiparallel β-sheet with a turn known as flap. The homodimer HIV-1 protease molecule possesses two flaps. The flaps of the protease (residues 33-62) must open for segments of the polyprotein to access the active site. Once the appropriate region of the polyprotein is in the active site, the flaps must close over the substrate sequence for cleavage to occur. The catalytic residues of HIV-1 protease are located at the bottom of a large open cavity, the ceiling of which is composed of 2 flaps, one from each monomer. Although the structure and function of HIV-1 PR have been studied for over 20 years, questions remain regarding the conformations and dynamics of the flaps that cover the active site cavity. The flaps of the protease must open for segments of the polyprotein to access the active site. Once the appropriate region of the polyprotein is in the active site, the flaps must close over the substrate sequence for cleavage to occur, it is shown in Figure 1. Fascinatingly, the flaps are reversed into a closed conformation in the ligand bound protease compared to the ligand free enzyme, which suggests that an extensive conformational change takes place upon ligand binding (Wlodawer et al., 1989;Miller et al., 1989;Rick et al., 1998). Understanding the issues that govern HIV-1 protease flap mutation and effect to the binding residues has profound implications for elucidating the detailed mechanism of this enzyme and in the design of new therapeutic agents, such as allosteric inhibitors intended to interfere with flap opening and thereby with enzymatic function. Thus, numerous prior computational studies have aimed at understanding flap opening dynamics. Reported results from activated molecular dynamics (MD) simulations in the gas phase that involved forcing the atomic coordinates for non flap regions of a closed structure to the semi open state (Collins et al., 1995). Darunavir, formerly TMC-114, is the latest approved protease inhibitor (PI) for the treatment of HIV infection. It was originally designed to be active against HIV strains resistant to other currently available PIs. In the latest International AIDS Society-USA panel list (www.iasusa.org, last update in December 2008), a total of 11 mutations were defined as specifically associated with Darunavir resistance. They were segregated as major (I50V, I54M, L76V and I84V) or minor (V11I, V32I, L33F, I47V, I54L, G73S and L89V) resistance mutations (Johnson et al., 2008). The enzyme shows the resistance, if the mutation is able to change the 3D space of binding cleft.
The mechanism behind Darunavir resistance inhibition is not well *Corresponding author: Rao Sethumadhavan, Bio-informatics Division ,School of Bio Sciences and Technology (SBST), Vellore Institute of Technology University, Vellore, Tamil Nadu, India, Tel: +91-416-2202522, 2202523; Fax: +91-416-2243092; E-mail: rsethumadhavan@vit.ac.in understood, but the variants do have mutations in both the active site and nonactive site regions. The inhibition can be accomplished by disturbing favorable ligand protein interactions in the binding pocket. However, accumulating evidence suggests that the flap region also plays an important role in influencing both substrate and inhibitor binding/stability (Scott and Schiffer, 2000;Clemente et al., 2004). Another possibility for resistance is that the active site specificity itself is reduced because of a mutation(s) in the flap region preventing flap closure, inducing an expansion of the active site cavity.
We carry out computer simulations in the hope of understanding the properties of assemblies of molecules in terms of their structure and the microscopic interactions between them. Molecular dynamics simulations can be used to predict the effects of mutations on the protease structure. Molecular dynamics simulations of HIV protease have been used to evaluate the catalytic mechanism and the structural flexibility (Harrison and Weber, 1994;Liu et al., 1996;Trylska, 2004). Here, molecular dynamics simulations of solvated native and two mutants protease were run in order to better model the solution behavior and to evaluate the physical basis for drug resistance.
In this study we carry out docking and binding energetics analysis of HIV-1 native and flap mutant structures with respect to Darunavir. To understand the above mention analysis result we performed explicit solvent molecular dynamics (MD) simulations starting from a native HIV-1 protease crystal structure and two resistance flap mutant (I54L and I54M). Overall, the native and mutants proteases ultimately deviate from their respective crystal structures to a similar extent. However, the native exhibits a more localized structural change than the mutants that begins around 2 ns after initiation of the simulation. The mutants alter its overall configuration from the crystal structure very quickly relative to the native. The fast configuration change in structure of mutants, alters 3D conformation of binding residues might be the cause of resistance toward Darunavir.

Data set
We selected native HIV-1 protease (PDB ID 1ODW) (Kervinen et al., 1996) structu res from Brookhaven Protein Data Bank (Berman et al., 2000), and one small molecule/inhibitor, Darunavir for our investigation. We induced two point mutation at 54 th position isoleucine to leucine and methionin in native structure to build mutant structures and the energy was minimized by normad server (Lindahl et al., 2006) with default parameters. The native structure also subject to the same in silico mutation methods, where the native residue is mutated to itself, to ensure that this mutation procedure is not introducing an artifact. We called I54I as mutant 1 and I54M as mutant 2 in our paper. Isoleucin locating at 54 th position at flap region and above mutant at this position produce the resistance against Darunavir (King et al., 2004;Ghosh et al., 2007). The SMILES string of Darunavir were collected from PubChem, a database maintained in NCBI (Feldman et al., 2006), and submitted to CORINA (www.molecularnetworks.com/online demos/corina demo.html) for constructing the 3D structure of small molecule.

Determination of binding pocket and sequence comparison
An accurate description of protein shape derived from protein structure is necessary to establish an understanding of proteinligand interactions, which in turn will lead to improved methods for protein-ligand docking and binding site analysis. Shape descriptors representing protein structure, such as depth (Coleman et al., 2005;Nayal and Honig, 2006), surface curvature (Coleman and Sharp, 2006), extreme elevation (Agarwal et al., 2004), surface area and volume (Liang et al., 1998), have been used extensively to identify, study and compare protein ligand interactions, protein-protein interactions and the respective binding sites. We have used CASTp server (Dundas et al., 2006), to determine the binding pocket. CASTp server uses the weighted Delaunay triangulation and the alpha complex for shape measurements. It provides identification and measurements of surface accessible pockets as well as interior inaccessible cavities, for proteins and other molecules. It measures analytically the area (in Å 2 ) and volume (in Å 3 ) of each pocket and cavity, both in solvent accessible surface (SA, Richards' surface) and molecular surface (MS, Connolly's surface). We set probe radius of 1.4 Å as default value. The detected pockets from these algorithms ranked with their volumes and areas. For ligand molecules internal protein cavities appear to be a favored binding sites (Eckenhoff and Johansson, 1997), and the cavity volume may play an important role in the strength of the guest molecule-host cavity interaction. In CASTp analysis, we have chosen the first pocket for our analysis, based on volume of the pocket. Sequence comparison is one of the parameter to understand the residual arrangement of protein. We performed multiple sequence alignment by ClustalW (Thompson et al., 1994) to check the position of active site residues in native and flap mutants.

Computation of docking score and interaction energy between the inhibitor and HIV-1 protease enzyme
We used the web server PatchDock (Duhovny et al., 2005) to compute the scores of docked complexes. 3D coordinates of all the individual HIV-1 protease and the inhibitor was submitted in PDB format with default parameters. In docking Asp29, Asp30, Asp128 and Asp129 residues acting as binding residues and participate in ligand interaction (Ghosh et al., 2007;Feldman et al., 2006) . The neighboring residues were given to address the active site inside receptor and better access of ligand molecule. The underlying principle of this server is based on molecular shape representation, surface patch matching plus filtering and scoring (Duhovny et al., 2002). It is aimed at finding docking transformations that yield good molecular shape complementarity. Such transformations, when applied, induce both wide interface areas and small amounts of steric clashes. A wide interface is ensured to include several matched local features of the docked molecules that have complementary characteristics. The Patch-Dock algorithm divides the Connolly dot surface representation (Connolly 1983a;Connolly 1983b) of the molecules into concave, convex and flat patches. Then complementary patches are matched in order to generate candidate transformations. Each candidate transformation is further evaluated by a scoring function that considers both geometric fit and atomic desolvation energy (Zhang et al., 1997). Finally, an RMSD (root mean square deviation) clustering is applied to the candidate solutions to discard redundant solutions. The main reason behind PatchDock's high efficiency is its fast transformational search, which is driven by local feature matching rather than brute force searching of the six-dimensional transformation space. It further speeds up the computational processing time by utilizing advanced data structures and spatial pattern detection techniques, such as geometric hashing and pose clustering. The interaction energy of the complex was calculated by PEARL web server (Han et al., 2006). The server compute total ligand-receptor interaction energy and its components are computed by an atomic level molecular mechanicsbased force field involving intermolecular van der Waals, electrostatic and hydrogen bond interactions between the binding molecule and its receptor. The electrostatic energy is particularly well suited for JCSB/Vol.3 Issue 4 analyzing recognition process because it is a physically meaningful representation of how a molecule is perceived by another molecule in its vicinity. It is through their mutual electrostatic interactions that the two molecules involved in the interaction first 'see' each other. The electrostatic energy is the only long-range non-bonded interaction and its value is calculated for pairs of atoms that are separated by three or more bonds. Thus it is assumed that the well-separated atoms are interacting mainly electro-statically since no cut-off distance was used in the calculation of the electrostatic interactions (Politzer and Muray, 1991; Naray-Szabo and Ferenczy, 1995). The negative value of electrostatic energy enables better interaction and vice-versa. The proper formation of hydrogen bonds and van der Waals contacts require complementarity of the surfaces involved. Such surfaces must be able to pack closely together, creating many contact points, and charged atoms must be properly positioned to make electrostatic bonds. Thus van der Waals and polar interactions contribute to the dynamic stability of the ligand-receptor complex (Chothia and Janin, 1975;Getzoff et al., 1992;Jacob, 1992) Individual receptor and the ligand molecule were given as input for performing the docking experiments. The default RMSD value (4.00 Å) was used and we given the receptor binding site residues information docking experiments. The complex structure file was given as input to PEARLS web server to perform the interaction analysis.

Solvent accessibility analysis
Solvent accessibility is the ratio between the solvent accessible surface area of a residue in three-dimensional structures and that in an extended tripeptide conformation (Fraczkiewicz and Braun, 1998). It indicates the packing arrangement of residues. The solvent accessible surface area is defined as the locus of the centre of the solvent molecule as it rolls over the van der Waals surface of the protein. It is typically calculated using the 'rolling ball' algorithm (Shrake and Rupley, 1973). The accessible molecular surface is determined for a probe with a radius of 1.4 Ả. It is the major intermediate step to understand the structure-function relationships of proteins. The location of amino acid residues in protein structures paves a way to determine their influence on interaction with other residues/surrounding medium and hence the folding and stability of proteins.
We obtained the solvent accessibility information by using ASAView (Ahmad et al., 2004). Solvent accessibility was divided predominantly into buried and exposed region indicating, the least accessibility and high accessibility of the amino acid residues to the solvent, respectively (Gilis and Rooman, 1996; Gilis and Rooman, 1997). The term buried surface area has two related meanings first the surface buried away from solvent when the protein folds, and second the surface buried away from solvent when two proteins or subunits associate to form a complex. We have analyzed native and mutants and observed the solvent accessibility of binding residues, before and after docking analysis. Crystallographic waters were not included. The protein was solvated in a cubic 0.9 nm of 15960 SPC (Berendsen et al., 1981) water molecules. At physiological pH, the protein is positively charged, thus in order to make the simulation system electrically neutral, we added four chloride ions (Cl−) to the simulation box using the "genion" tool that accompanies Gromacs. The simulation system was set up as an NPT ensemble, i.e., constant number of particles (N), constant pressure (P) and constant temperature (T). All protein atoms were at a distance equal to 1.0 nm from the box edges. The system was subjected to energy minimization for 1000 steps by steepest descent. The minimized system was equilibrated for 100 ps each at 300 K by position restrained molecular dynamics simulation in order to permit waters to diffuse into the cavity to prevent the collapse (Meagher and Carlson, 2005). The equilibrated systems were then subjected to molecular dynamics simulations for 1500 ps each at 300 K. In all simulations, the temperature was kept constant at 300 K with a Berendsen thermostat (Berendsen et al., 1984). The particle mesh Ewald method (Essmann et al., 1995) was used to treat longrange Coulombic interactions and the simulations performed using the SANDER module (Case et al., 2002). The ionization states of the residues were set appropriate to pH 7 with all histidines assumed neutral. The SHAKE algorithm was used to constrain bond lengths involving hydrogens, permitting a time step of 2 fs. Van der Waals and coulomb interactions were truncated at 1.0 nm. The non-bonded pair list was updated every 10 steps and conformations were stored every 0.5 ps. Other analyses were performed using scripts included with the Gromacs (Hess et al., 2008) distribution. The visual analysis of protein structures was carried out using (Rasmol Sayle and Millner-White, 1995).

Molecular dynamics simulation
The major focus of this study is to compare the dynamic behaviors of native and mutant protein at 300 K. To that end we compare root mean square fluctuation (RMSF) of carbon alpha and root mean square deviation (RMSD) of backbone structure of protein between the trajectories generated at 300 K to explore the flexible nature of mutant.
In order to generate the three-dimensional backbone root mean square deviation (RMSD) and RMSF of carbon alpha carbon of the system were plotted for all three simulations using GRACE program (version 5.1.22).

Result
We report the pocket volume of native and mutant HIV-1 protease enzymes. Appropriate cavity volume gives better accommodation of ligand molecule inside the pocket. We found considerable difference in pocket volume of native, mutant 1 and mutant 2 as 299.8, 940.7 and 959.3 (in Å 3 ) respectively. We evaluated solvent accessibility and performed the molecular dynamics simulation of native and mutants structure to explore the conformation changes with respect to time. We highlighted the RMSF of Asp29, Asp30, Asp128 and Asp129 act as binding residue for Darunavir in HIV-1 protease enzyme (King et al., 2004;Ghosh et al., 2007). Our present study also endorses this experimentally observed fact. During sequence comparison analysis binding residues position was found identical in native and mutants, it is shown in Figure 2. Our investigation showed the behavior of protein-ligand complex of HIV-1 protease enzyme with Darunavir. Docked native HIV-1 protease complex is shown in Figure 1. The docking score of complex of Darunavir with native, mutant 1 and mutant 2 are found to be 7208, 6800 and 6532 respectively (Table 1).
Van der waals and electrostatic energy contributed major part with -0.18 and -5.13 kcal/mol in total ligand receptor interaction energy of native HIV-protease complex that is -7.28 kcal/mol. Van der waals and electrostatic energy confirms significant amount of complementarities between native receptor and Darunavir. While in mutant 1 has less van der waals and less electrostatics energy as -1.28 and -3.17 kcal/mol to total interaction energy -6.34 kcal/mol. Further reduction in both energies in mutant 2 given total interaction energy of -4.80 kcal/mol (Table 1). It proves that deviation from actual position of binding residues (Asp29, Asp30, Asp128' and Asp129') may be the reason of reduction in interaction energy and gradual decrease in docking score.
We found that binding residues in all the mutants is located in buried regions before docking and exposed region after docking (Table 1). During the process of docking the buried binding residues are being exposed and the cavity is occupied with inhibitor molecule.
We have performed molecular dynamics simulation of each native and mutant to observe the deviation of each atom of mutant as well as native structure in water with respect to time to explore the flexible nature of binding residues. In this analysis we evaluated RMSD of backbone and RMSF of C-alpha carbon.
In the trajectory of native RMSD change significantly up to 200 ps after which it remains the same till 13 ns and increases slightly and attains a value ~0.15 nm at 1500 ps. In the trajectory of mutant 1 structure backbone RMSD increases slightly from the starting conformation, which fluctuates, between 0.1 nm and ~0.18 nm during the simulation. Mutant 2 RMSD does not change significantly up to 200 ps after which it change from 0.13 nm to 0.17 nm at 220 ps and attain a value 0.22 nm at 1500 ps. These results are summarized in Figure 3.
RMS fluctuations have been calculated for 1500 ps. We computed the RMS fluctuations per residue over the production run of native and mutant structures, It is shown in Figure 4. We found that the  Table 1: Total cavity volume, Docking score, ligand-receptor interaction energy and solvent accessibility analysis at binding residues of wild and flap mutants.   In native structure the RMS fluctuation at binding residues in chain A and chain B was observed minimum and maximum in mutant 2 structure. Mutant 1 made intermediate RMS fluctuation (Table 2).

Discussion
Pocket analysis shown the expansion of cavity form native to mutant structures. This expanded cavity reduces the fit of inhibitors and lowers the binding affinity by decreasing van der Waals contacts and hydrogen bonding. Higher the docking score shows better complementarity between receptor and drug molecule. Docking score of native type structure was highest as compare to mutant type, it indicates good agreement between receptor and Darunavir. Lowest docking score of mutant 2 indicate least complementarity. We computed the total interaction energy of the complex significantly contributed by van der waals and electrostatic interaction energy between HIV-1 protease (native and mutant) and Darunavir. Lower energy of the complex shows the better interaction and good agreement between receptor and ligand molecule. Estimates of VDW interaction energy were computed to provide a theoretical quantitative assessment for the ligand-protease non-bonded interactions. The two factors, namely, score difference in docking process and variation in binding energy of complex, probably correspond to conformational alteration of binding residues due to flap mutations. Solvent accessibility of amino acids involved in structural and functional role for native and mutants of HIV-1 protease was analyzed. The occupancy of inhibitor during docking process further indicates the flexible behavior of residues of mutant structure. It proves that the increased flexibility govern by point mutation enable resistance behavior. It further provoked us to study the dynamics to understand the flexible behavior of HIV-1 protease enzyme. Since the docking process is quite related to binding residues so we have highlighted root mean square fluctuation of C-alpha carbon of binding residues in molecular dynamics simulation.
The RMSD of the backbone atoms of the protein from the starting structure over the course of simulation may be used as a measure of the conformational stability of a protein during that simulation. The two phenomenons, variation in RMSD value of residues in 3D space and decrease level of binding energy is due to mutation at flap residual position 54 th . We confirmed the above observation by studying the RMSF of binding residues. Compare to mutants, native structure shown less deviation with respect to time. Hence native have highest docking score and good interaction energy with Darunavir.
Root mean-square fluctuations (RMSF) analysis helps to understand the molecular flexibility nature of C-alpha carbon atom of binding residues which could influence the binding affinity toward the ligand molecule. A more detailed picture of differences in residue mobility within and between simulations can be obtained from graphs of the RMSF of Ca atoms relative to the average structure. Native has highest docking score and lowest fluctuation at alpha carbon atom of binding residues as compared to other two of mutants studied namely, mutant 1 and mutant 2. Mutant 2 showed the least docking score, thereby, expressing high fluctuation in the C-alpha carbon atom.
During whole simulation, native structure showed less fluctuation. This integral behavior of native HIV-1 protease gives good accommodation to inhibitor inside the pocket and better interaction with it. Hence native has significant susceptibility toward Darunavir. But flap mutations at 54 th position to leucine/methionine are able to alter the structural position of binding residues from normal. Mutants have shown less interaction energy, docking score and more fluctuation of binding residues then native. Thus in-vitro, mutants have considerable amount of Darunavir resistance.

Conclusion
In this study, mechanism of Darunavir resistance flap mutants and dynamics of HIV-1 protease binding residues conformations was investigated using molecular dynamics simulations. Our analysis revealed docking score with Darunavir varied from 7556 to 6880. Decrease in docking score suggested the alteration of binding residues in 3D space. RMSD and RMS fluctuation values with respect to native type supported the docking score order. On the basis of docking calculation, binding energy and molecular dynamics simulation we conclude the resistance behavior of mutants was due to alteration in 3D structure at binding site residues.
The mutations at 54 th position in flap region as a result of conformational change in binding residues often contribute to Darunavir resistance. However, several of these residues are not part of the active site cavity, and their essential role in causing drug resistance could possibly be rationalized if this conformational change actually occurs. Trapping HIV-1 protease in this inactive conformation would provide a unique opportunity for future drug design The structures insight obtained with these simulations could be considered as a starting point for more refined all atom models and simulations aimed at designing drugs that target different sites far from the active one, suggesting novel therapeutic strategies.