Interchain hydrophobic clustering promotes rigidity in HIV-1 protease flap dynamics: new insights from Molecular Dynamics

The dynamics of HIV-1 protease (HIV-pr), a drug target for HIV infection, has been studied extensively by both computational and experimental methods. The flap dynamics of HIV-pr is considered to be more important for better ligand binding and enzymatic actions. Moreover, it has been demonstrated that the drug-induced mutations can change the flap dynamics of HIV-pr affecting the binding affinity of the ligands. Therefore, detailed understanding of flap dynamics is essential for designing better inhibitors. Previous computational investigations observed significant variation in the flap opening in nanosecond time scale indicating that the dynamics is highly sensitive to the simulation protocols. To understand the sensitivity of the flap dynamics on the force field and simulation protocol, molecular dynamics simulations of HIV-pr have been performed with two different AMBER force fields, ff99 and ff02. Two different trajectories (20 ns each) were obtained using the ff99 and ff02 force field. The results showed polarizable force field (ff02) make the flap tighter than the nonpolarizable force field (ff99). Some polar interactions and hydrogen bonds involving flap residues were found to be stronger with ff02 force field. The formation of interchain hydrophobic cluster (between flap tip of one chain and active site wall of another chain) was found to be dominant in the semi-open structures obtained from the simulations irrespective of the force field. It is proposed that an inhibitor, which will promote this interchain hydrophobic clustering, may make the flaps more rigid, and presumably the effect of mutation would be small on ligand binding.


Introducion
Nowadays, there are many drugs available against human immunodeficiency virus (HIV). These drugs are aimed at targeting different HIV proteins at different stages of the viral life cycle. Mainly three HIV proteins, HIV-Reverse Transcriptase (HIV-RT), HIV-protease (HIV-pr), and HIV-gp41 are primarily targeted. Of these three targets, HIV-pr acts at the late stage of infection by cleaving the Gag and Gag-Pol polyproteins to produce mature infectious virions for the continuation of the viral life cycle. Although there are several protease inhibitors approved by the Food and Drug Administration (FDA), the current major problem is drug resistance caused by drug-induced mutations on HIV-pr. Most of the FDA approved drugs, including the recent one (Darunavir), are becoming resistant against HIV-pr. So designing drugs having less susceptibility to mutations considered to be a major problem. To combat this problem, knowledge about the molecular mechanism of drug resistance is essential.
HIV-pr belongs to the aspartyl protease family which is structurally homodimeric and symmetric. Active site of this protein is covered by two glycine-rich flexible flaps that control access to the active site ( Figure 1). The active site of the protease is accompanied by the presence of two catalytic triads (Asp25-Thr26-Gly27 in chain-A and Asp25′-Thr26′-Gly27′ in chain-B) having the functional Asp residues located at the dimer interface. HIV-pr contains 198 residues with 99 residues in each chain. The residues are marked as 1-99 and 100-198 (or 1′-99′) for chain-A and B, respectively. Flap (residues 43-58 and 142-157), flap elbow (residues 35-42 and 134-141), fulcrum (residues 11-22 and 110-121), cantilever (residues 59-75 and 158-174), and active site regions for both the chains-A and B are represented in Figure 1. In this paper, both types of residue numbering for chain-B (100-198 or 1′-99′) are used. X-ray diffraction studies (Lapatto et al., 1989;Spinelli, Liu, Alzari, Hirel, & Poljak, 1991;Wlodawer et al., 1989) apparently shown that in the ligand-bound form, the flaps are on top of each other resulting in a closed conformation. On the other hand, in the apo-form of HIV-pr, the flaps are widely away from each other resulting in semi-open conformation (Louis, Dyda, Nashed, Kimmel, & Davies, 1998;Prabu-Jeyabalan, Nalivaika, & Schiffer, 2000). Stable structural differences exist between the bound and unbound states of the protein ( Figure 2). However, it has been suggested that crystal-packing effects may be responsible for the specific semi-open conformations (Layten, Hornak, & Simmerling, 2006) seen in crystal structures. Semi-open conformations may have a wide variety of structures in the absence of crystal-packing effects.
The drug-induced mutations in HIV-pr structure has been found to occur near the active site and also away from the active site. (Muzammil, Ross, & Freire, 2003) The mutations near the active site influence the direct interaction between the protein and the inhibitor to a greater extent, and thus the free energy of binding can become unfavorable. However, the molecular mechanism of nonactive site mutations is more intriguing. How can a mutation away from the active site (ligand-binding site) affect the ligand binding at the active site? The current consensus from NMR and computational studies is that certain mutations change the conformational dynamics, especially the flap dynamics of HIV-pr. This in turn affects significantly the binding affinity of the ligands. Some mutations may confer both direct and indirect effects (Bandyopadhyay & Meher, 2006;Meher & Wang, 2012). There were also a number of studies carried out to reveal the molecular mechanisms of drug resistance in various HIV-pr mutants (Chen, Weber, & Harrison, 2004;Chen, Zhang, Liu, & Zhang, 2010;Dirauf, Meiselbach, & Sticht, 2010;Hou, McLaughlin, & Wang, 2008;Hu, Zhu, Zhang, Wang, & Zhang, 2010;Kar & Knecht, 2012;Muzammil et al., 2007;Prabu-Jeyabalan, Nalivaika, & Schiffer, 2002;Soares et al., 2010). But the complete understanding of drug-resistance mutations of HIV-pr is still lacking. (Ali et al., 2010) Hence, a molecular-level understanding of drug resistance needs the knowledge of both direct and indirect effects of mutations.
Torchia group from NMR studies showed that there are two time scales of flap dynamics of HIV-pr (Ishima, Freedberg, Wang, Louis, & Torchia, 1999). One is flap tip dynamics, involving only the tip of the flaps, in the nanosecond range. The other one is the opening and closing of the flaps, which is much slower, in the range of micro-to milli-second time scale. The effect of mutation on flap dynamics has been investigated by   (Maschera et al., 1996). Their findings inferred that, the mutants G48V, L90M, and L90M/G48V has lower drug-binding affinity. This is due to the increase or decrease in flap opening or closing rates, respectively. To probe the dynamics at atomistic level, molecular dynamics (MD) simulation has been a major tool to employ these days. Therefore, we used MD to investigate and understand flap dynamics. We find large number of computational studies on HIV-pr dynamics investigating both the wild type and mutants. Most of the studies have used MD simulation. These works looked into flap dynamics using simulations with different lengths and starting from structures of both apo (Meher, Satish Kumar, Sharma, & Bandyopadhyay, 2012) and liganded form of the protein.
For instance, Piana et al. performed ab initio MD and classical MD simulations of the HIV-pr complex and found that the flap dynamics affects and control the distance between the substrate and the catalytic aspartates  and hence the activation free energy barrier of the enzymatic reaction (Piana, Carloni, & Parrinello, 2002;Piana, Carloni, & Rothlisberger, 2002). Perryman, Lin, and McCammon (2004) suggested that the flap dynamics of V82F/I84V mutant is significantly different from that of the wild type. Schiffer and coworkers examined the flap opening to a greater extent (Foulkes et al., 2006;Scott & Schiffer, 2000). Ode et al. investigated the effect of nonactive site mutation L90M and found that the mutation changes the interaction between the side chain of 90 and 25 (Ode, Neya, Hata, Sugiura, & Hoshino, 2006). Our previous work looked at the I47V mutations and we found that along with subtle change in dynamics, change in a specific interaction between the side chain of residue 47 in chain-B and the ligand is likely the main reason for drug resistance (Bandyopadhyay & Meher, 2006).
One important point is that most of the previous MD simulations were done with explicit waters in the nanosecond time scale, much shorter than the actual time scale of flap opening. Interestingly, not all simulations reported the same extent of flap opening in nanosecond time scale. Later, it has been suggested that this differential flap opening is probably due to differences in the force field used and/or of the simulation protocol. In an interesting work, Simmerling group used implicit solvent with a low viscosity to model the flap opening (Hornak, Okur, Rizzo, & Simmerling, 2006). Since low viscosity allows the faster folding rates of proteins, transitions from close to open and back was seen in the simulation. In another work, Toth et al. also used implicit solvation model and found very large flap opening in only few ns of simulation (Toth & Borics, 2006). However, direct correlation of the actual time of this conformational transition to the time taken in the implicit solvent simulation is nontrivial.
While there are a number of issues related to the flap opening and closing in HIV-pr for both wild type and mutant and their effects on drug resistance, a somewhat different approach is taken to investigate the flap motion of HIV-pr in the present work. In a nice work, Carlson group showed that initial solvent distribution around HIV-pr and improper equilibration can influence the flap dynamics. (Meagher & Carlson, 2005). Moreover, biasness of a single trajectory, which has been used by most of the authors, cannot be ruled out. To investigate the effect of force field on the flap dynamics of HIV-pr, we have performed simulations on apo HIV-pr using two force fields, one AMBER ff99 (a nonpolarizable force field) (Wang, Cieplak, & Kollman, 2000) and the other AMBER ff02 (a polarizable force field) (Cieplak, Caldwell, & Kollman, 2001) as implemented in the AMBER suite of programs. Two different explicit water models, TIP3P (Jorgensen, Chandrasekhar, Madura, Roger, & Klein, 1983) and POL3 (Caldwell & Kollman, 1995), were used with ff99 and ff02, respectively. Two trajectories with different initial velocity distributions were run for each of the force field to minimize the potential bias of a single trajectory. Henceforth, the nonpolarizble and polarizable simulations will be described by nopol and pol, respectively. An equivalent simulation with the ff99 force field was performed for the complexed HIV-pr/JE-2147 system.
Our simulations were started from a ligand-bound structure of HIV-pr after deleting the ligand. It can be argued that the simulation time scale is too short for the apo-protein to reach equilibrium. However, this is done intentionally with the objective of getting close to at least partly open structure in a short time scale like nanosecond. Taking the structure of apo-protein (which is already in the semi-open state), it is difficult to get further flap opening in nanosecond range. The same protocol of removing ligand from a closed structure has been carried out in number of studies in the past. (Bandyopadhyay & Meher, 2006;Hornak et al., 2006;Perryman et al., 2004).
A detail investigation of the flap dynamics in the five trajectories gave several new insights both on the difference between the two force fields used and on the mechanism of flap movement. Especially, a number of structural interactions responsible for flap movement came out from this study by comparing the five trajectories. Generally, pol simulations make the flap dynamics tighter than the nopol ones as observed in one of our previously published article (Meher, Kumar, & Bandyopadhyay, 2009). There are number of reasons for the tightness of the pol trajectories. First, it was found that the interflap and several intraflap H-bonds are stronger in the case of pol simulations. Also, the polar interactions among the flap residues, far apart, are stronger in case of pol trajectories. Hydrogen bond with water molecules also shows slightly different patterns for the two simulations. While there are several differences between the two trajectories in each set of pol or no-pol, the qualitative conclusions are similar. While running two simulations of 20 ns, each is not enough to say the final word on the difference of these two force-ields, ff02 samples better conformations than ff99 force field. It does also show that subtle differences between two force fields can have an impact on flap dynamics, which depends on several factors such as the energy function used in the simulation (force field), the primary integrator for time propagation, initial velocity distribution, initial solvent distribution, and convergence of the trajectories.
One important observation came from this study is the importance of hydrophobic cluster formation between active site flap and active site wall of the other chain (i.e. 50-81′ or 50′-81). This is found to happen irrespective of the force field used. Flaps move in such a way that this hydrophobic interaction is maximized. Though both sideways movement and at a direction perpendicular to it are possible for the flaps, in our simulation two flaps mostly move in such a way that they are still on top of each other. In this orientation, if one flap bends, the hydrophobic cluster can still be formed even if the flaps move as far as 11.5 Å. This is a major finding of the present work that in semi-open structures, flaps can have various forms but the major driving force is the maximization of hydrophobic cluster between two chains. This result gives rise to a proposal of designing inhibitors which will facilitate the hydrophobic interactions between the two chains. The current inhibitors facilitate the hydrophobic clustering within a chain.

Methods
All five MD simulations were started using a highresolution (1.09 Å) crystal structure of HIV-pr with JE-2147 ligand-bound form (pdb ID: 1KZK) (Reiling, Endres, Dauber, Craik, & Stroud, 2002). Ligand from the protein was removed to create the apo-protein starting structure. The leap module of AMBER 9 program package (Case et al., 2006) was used to prepare the system for simulation. The ff99 and ff02 force fields were used for the nonpolarizable and polarizable simulations, respectively. Each structure was solvated in a water box such that there were 10 Å of water surrounding the protein on all sides (box volume of 83.511 Â 62.546 Â 68.898 Å 3 containing more than 8000 water molecules). The plot variation of simulation box dimension with time is given in the supplementary material section ( Figure S1), of which consistent axes shows the stability of the box. A default cutoff of 8.0 Å was used for Lennard-Jones interactions, and the long-range electrostatic interactions were calculated with the particle mesh ewald method (Essmann et al., 1995). The TIP3P and POL3 models were used to represent the water molecules with ff99 and ff02 force fields, respectively. The net positive charge on the system was neutralized through the addition of chloride ions.
For the simulation of the apo-protein, charged state of both the catalytic Asp-25 and Asp-124 were kept deprotonated. However, for the JE-2147 (JE2) complexed proteins, one of the catalytic Aspartates in chain-B (Asp124) was protonated by the addition of one proton to the OD2 atom of Asp124. Constant temperature and pressure conditions in the simulation were achieved by coupling the system to a Berendsen's thermostat and barostat (Berendsen, Postma, Van Gunsteren, DiNola, & Haak, 1984). Bonds involving the hydrogen atoms were constrained to their equilibrium position with the SHAKE algorithm. The trajectories having two different initial velocities were generated by giving different seeds for the random numbers used to generate the velocity distribution.
The system was minimized for two-phase minimization. In the first phase, the system was minimized giving restraints (30 kcal/mol/Å 2 ) to protein and crystallographic waters for 5000 steps with subsequent second phase minimization of the whole system. Then, the system was heated to 300 K over 25 ps with a 1 fs time step. The protein atoms were restrained with force constant of 30 kcal/mol/Å 2 at NVT ensemble. After that the force constant was reduced by 10 kcal/mol/Å 2 in each step to reach the unrestrained structure in three steps of 10 ps each with NVT ensemble. Then, the system was switched over to NPT ensemble and equilibrated without any restraints for 20 ps with a 1 fs time step and subsequent 160 ps with a 2 fs time step. The system was equilibrated in total of 235 ps. The time step for both nopol and pol production trajectories were 2 fs. All five trajectories were run for 20 ns. Henceforth, the five trajectories will be denoted by NOPOL-1, NOPOL-2, NOPOL-3 and POL-1, POL-2. Where, for instance, NOPOL-1 and NOPOL-2 means first and second trajectories of the nopol simulations and POL-1 and POL-2 are the first and second trajectories of the pol simulations. NOPOL-3 is the only simulation trajectory with JE-2147 complexed protein. Due to the unavailability of a suitable polarizable potential for the ligand molecule, that is compatible with AMBER ff02, the ligand-bound simulation was ignored for pol simulation. Although it is possible to develop ff02 potential for the ligand, a fine balance among ff02 force field for protein, pol3 force field for solvent and ff02 force field for ligand is required. We plan to investigate the polarization effect of the ligand in a future work. The convergence of energies, temperature, pressure, and global RMSD was used to verify the stability of the systems. The hydrogen bonds (H-bonds) were analyzed using the ptraj module of AMBER program. Formation of the H-bonds depends on the distance and angle cutoff as follows: (a) distance between proton donor and acceptor atoms were 63.5 to 64.0 Å, and (b) the angle between donor-H … acceptor was P120°. Graphic visualization and presentation of protein structures were done using Chimera (www.cgl.ucsf.edu/chimera).

Results
The results are presented in the following manner. First a general outcome of the force field quality by evaluating different parameters like RMSD, B-factors, S2 order parameters and stereo-chemical calculations for amino acids in the proteins. Secondly, the dynamics of the flaps using different structural parameters is presented with emphasis on the difference between the nopol and pol trajectories (significant difference between the two trajectories within pol or nopol is also discussed). Next, the interactions likely to be responsible for these differences are presented. Finally, some general observations of flap dynamics, present in all the trajectories irrespective of nopol or pol, is discussed.

Root mean square deviation (RMSD) of the Cα atoms
Before examining the flap movement, stability of all the trajectories were monitored by plotting the RMSD values of the C α atoms for the whole protein ( Figure 3). The values are in the range of 1.0-2.3 Å in all the trajectories, similar to other simulations on HIV-pr, ensuring stable trajectories. However, the RMSD values for the JE2 bound trajectory (NOPOL-3) are found to be lesser than the unbound trajectories (NOPOL-1 and NOPOL-2).

B-factors analysis
The B-factors (or isotropic temperature factors) can be taken as a measure to check the relative thermal/vibrational motion or flexibility of different parts of the protein structure. In order to check the flexibility in our protein, we have calculated and compared the B-factors. Figure 4 shows a comparison between X-ray structure B-factor and calculated B-factors for all the ff99 and ff02 trajectories. It can be seen that the calculated B-factors are much higher than that of X-ray structure B-factors and ff99 has higher values as compared to ff02, which apparently confirms the higher flexibility of protein in nonpolarizable simulation and rigidity in the polarizable one. While the B-factors for the two pol simulations are quite similar, it is far different within the nopol simulations. Hence, it may be possible to say that the description of protein dynamics through ff99 force field is poor.

S2 (N-H) order parameter analysis
Generalized order parameter (S2) for the N-H bond vector is a measure to evaluate the flexibility in proteins. S2 values calculated from the MD simulation trajectory using time-dependent correlation motion function typically agrees well with NMR S2 values (Chandrasekhar, Clore, Szabo, Gronenborn, & Brooks, 1992). This is often useful to validate the conformational dynamics of proteins. The S2 order parameter is an indicator of  the X-ray crystal structure (1KZK) of the protein HIV-pr. It showed that, calculated B-factors are much higher than the X-ray structure B-factor, although the general pattern is similar with only few exceptions. B-factors for NOPOL-1 and NOPOL-2 are much different than the POL-1 and POL-2 simulations, while NOPOL-3 B-factor is reduced may be due to the ligand binding.
protein backbone motions in computationally feasible time scales. The flexible regions or highly mobile regions in the protein can be inferred from the S2 values (lower the S2 value, higher is the flexibility and vice versa). In the present study, the S2 values were predicted by the S2 prediction server (Zhang & Bruschweiler, 2002) utilizing 1000 protein models obtained in an equal interval of 20 ps from the long 20 ns trajectories. The S2 values were averaged over the two chains of HIV-pr from both the nonpolarizable (ff99) and polarizable (ff02) force fields. Figure 5 shows the comparison of the predicted S2 values for both POL (S2 POL-1 & -2 ) and all NOPOL (S2 NOPOL-1, -2 & -3 ) simulations. It is markedly clear that the S2 NOPOL values are much lower than the S2 POL values, at least in five different regions of the protein, which are indicated as black solid arrows in the figure. As the S2 NOPOL values are lower, it makes the HIV-pr highly flexible as compared to the models sampled from the polarizable (POL-1 & 2) force field, ff02. On the other hand, the description for S2 NOPOL-3 is more or less different from other simulations showing the rigid behavior with lower S2 values in certain regions of the protein. In fact, this is as expected, because binding of the ligand makes the flap closure and partial reduction of protein dynamics.

Analysis of the Ramachandran plot
The stereochemical geometry of each residue in the protein models were diagnosed and analyzed by plotting the Ramachandran plot for the averaged structure over 1000 protein structures sampled through ff99 and ff02 force fields. Dihedral angles on the plot for the calculated structures shows many differences, according to which residues in the most favored regions for NOPOL-1, -2, and -3 are 89.1, 88.5, and 86.5%, respectively, whereas for POL-1 and POL-2, it is 90.4 and 92.3%, respectively ( Figure 6 and Table 1). Ramachandran plot for the JE2-bound protein (NOPOL-3) is given in the supplementary material section ( Figure S2). The percentage ensures the geometrically acceptable quality of proteins sampled from polarizable force field is greater than the nonpolarizable one.
From our earlier study on HIV-1 protease conformational dynamics, we found that the ff99 force field describes the protein dynamics poorly as compare to the ff99SB and ff03 force fields . However, that study was limited to the use of nonpolarizable force field and unliganded form of the protein only. In the current study, we explored and evaluated the protein conformations in both liganded and unliganded form which were simulated in both polarizable (ff02) and nonpolarizable (ff99) force fields. We observed that, the ff99 force field makes the protein more flexible as compared to ff02 which is evident from the comparison of the Cα RMSD values, B-factors and the S2 order parameters. The stereochemical geometry of protein models from both the force fields were evaluated by comparing the percentage of occupancy of residues in the most favored regions in Ramachandran plot. The percentage (Table 1) was found to be higher in the case of ff02 than ff99 in both ligand bound (NOPOL-3) and unbound (NOPOL-1 & -2) conformations. Hence, it may be possible to propose that ff02 samples better protein models than ff99.

Flap movement
3.3.1. Root mean square deviation (RMSD) of the flap Figure 7 shows the RMSD of the flap residues (43-58 and 43′-58′) for both chains for all five trajectories. The maximum range of RMSD is slightly more than 2.5 Å with values mostly within 0.5 and 2.5 Å with considerable overlap among the four JE2-unbound (except NOPOL-3) trajectories. RMSD for the NOPOL-3 simulation shows reduced flap dynamics behavior as expected. According to Hornak et al. (2006)

Flap tip movement
However, the flap tip movement is expected to occur in a short time scale. To check the flap dynamics behavior and if there is any difference among all the simulated trajectories, we have examined the angles involving the C α atoms of the following sets of three residues, Gly48-Gly49-Ile50, Gly49-Ile50-Gly51, Ile50-Gly51-Gly52, and Gly51-Gly52-Phe53 in both the chains. The results showed that angle Gly48-Gly49-Ile50 has the maximum variation among the five trajectories followed by the 49-50-51 angle. The mean and standard deviation (SD) of the TriCa angle Gly48-Gly49-Ile50 are 139.87°and 13.5°for NOPOL-1, 130.47°and 14.0°for NOPOL-2, 144.84°and 4.32°for NOPOL-3, 145.28°and 4.3°for POL-1, and 143.32°and 8.3°for POL-2. A list of mean and SD values of the other angles calculated are given in the supplementary section material (Table S1). It indicates that the fluctuations of flap tip region in nopol simulations are quantitatively higher than pol simulations. The angles vary a lot between the two trajectories within the nopol simulations. Figure 8(a) shows the time series of 48-49-50 angle for all the nopol trajectories. It is clearly seen that the trajectories sample angles from around 100°to little less than 160°. More interesting is the difference between the two trajectories. The first nopol trajectory (NOPOL-1) jumps between the small angle region and large angle region several times in the first 8 ns, with residence time is more in the small angles. On the other hand, the second nopol trajectory (NOPOL-2) spends more time in the smaller angle region in the initial 15 ns and with occasional jumps to the larger angle region. These findings support the highly irregular structure of the flap tip in the ns time scale. Also, it shows the difference of flap tip angle with two different MD trajectories. However, the angle (Gly48-Gly49-Ile50) remains to be stable in the upper angle    region after 15 ns time in both the NOPOL-1 and -2 cases. Additionally, the TriCa angle for JE2-bound (NOPOL-3) simulation has a consistent angular value at around 140°to 155°showing its rigid behavior upon binding of the ligand. Figure 8(b) shows the time series plot of 48-49-50 angle for the two pol trajectories. It is clearly seen that both the trajectories almost exclusively sample the high angle structures, with the second pol trajectory (POL-2) spends little time in the small angle structures around 3 ns, 12 ns, and 16 ns region and then came back to the high angle region. Here, the difference between the two trajectories reduced though clearly some differences exist. Although there are differences in two trajectories within a set of pol or nopol, the qualitative conclusions are valid if any one trajectory in one set is considered.
Flap tip curling is suggested to be related to flap opening (Perryman et al., 2004;Scott & Schiffer, 2000). Perryman et al. found that smaller flap tip angle generally indicate larger flap opening. However, Hornak et al. (2006) suggested that since the time scales of flap tip motion and full flap opening are orders of magnitude different, the direct correlation between the flap curling and flap opening may not always exist. We can conclude from the examination of the TriCa angle 48-49-50 that there are distinct differences in nopol and pol trajectories (and also within the two nopol trajectories). Nopol trajectories have a larger variation of TriCa angles going from small (low) to large (high) region and coming back, which can be observed from the superimposed images from representative structures of NOPOL-1 (Figure 9).
In chain-B, the difference between the five trajectories for angle (147-148-149) is less ( Figure S3). One pol trajectory (POL-1) is almost always sampled the high angle region, while the other four including NOPOL-3 sampled mostly the larger angles with occasional small angles region. This suggests the differential motion of each of the chains in HIV-pr. Figure 10 shows the histograms of the I50C α -I149C α distance for all the five trajectories. It can be seen that both the NOPOL (1 & 2) and POL (1 & 2) histograms are relatively similar. For both the sets of simulations, there are two peaks at least. The smaller peaks have the mean around 7 and 8 Å, while the larger peaks are around 9 Å. However, NOPOL-3 simulation has a single peak whose mean is around 6 Å. It clearly indicates that, upon binding of the ligand JE2 to the active site, the movement of the flaps is reduced and the NOPOL-3 simulation samples mostly closed structures. This reduced flap dynamics is remarkably different than the other simulations. However, for both the sets of nopol and pol simulations, the flap tip-flap tip distances span up to 11.5 Å for some snapshots. The histograms indicate that the flap-flap distance is mostly less than 7 Å for the NOPOL-3 trajectory, while for the JE2-unboud trajectories there are two important conformations (at around 7 and 9 Å) around which the flap-flap distance fluctuates. Examination of the time series plot showed that in most of the trajectories, the switching between the two conformations has happened repeatedly. According to the definition used by Toth et al. (2006), open structures should have flap-flap distance more than 10 Å and semi-open is between 6 and 10 Å. As par that definition, semi-open structures are predominant in our simulation.

Flap-active site
Flap-active site distance has been examined by several authors to get an idea of the size of the active site. Figure 11(a) shows the histograms of the I50C α -D25C β distance for all the trajectories. It is seen that the mean values of the two nopol trajectories (NOPOL-1 & -2) are slightly more than that of the two pol (POL-1 & -2) trajectories. The maximum values the two nopol trajectories sample are in the range of 15-18 Å, while the maximum  in pol simulation is 16 Å. However, the difference between the nopol and pol is more pronounced here as compared to the flap-tip-flap-tip distance. For the chain-B, all of the trajectories are mostly overlapping (Figure 11(b)) with mean around 14 and 15 Å. However, the NOPOL-3 trajectory samples the flap-active site distance on an average much less than other JE2-unbound simulations, with a mean of $14 Å indicating the structures are in closed conformations as far as the flap-active site distance is concerned.

Movement of other regions
We have also examined several other regions of HIV-pr, namely flap elbow, fulcrum, cantilever, and active site wall regions (see Figure 1) and found that the flap elbow has the maximum fluctuation.  Figure S4(a) and (b)) similar feature as the RMSD of the flaps. It was found that for chain-A, the two nopol trajectories have smaller values of RMSD than the two pol trajectories. The RMSD mean for NOPOL-1 (NOPOL-2) are 0.48 (0.52) Å as compared to the POL-1 (POL-2) having 0.59 (0.60) Å, respectively. It means that the two nopol trajectories have higher RMSD values in the flap region, as discussed earlier. This may be due to the fact that both the flaps and flap elbows have anticorrelations motion. For chain-B, the two nopol trajectories are similar and have smaller RMSD values than the POL-1 trajectory. Interestingly, POL-2 trajectory shows a much higher RMSD among the five trajectories, similar to large change in the flapactive site distance. The RMSD mean for NOPOL-1 (NOPOL-2) are 0.73 (0.65) Å as compared to the POL-1 (POL-2) having 0.79 (1.30) Å, respectively. Moreover, RMSDs for the JE2 bound protein simulation (NOPOL-3) are almost always higher for chain-A, but for chain-B, it fluctuates in between 0.3 and 1.5 Å. Mean values for NOPOL-3 chain-A and B are 0.80 and 0.83 Å, respectively.

Interactions responsible for flap opening
In the previous section, various structural parameters related to flap dynamics are discussed for all the five trajectories. In general, the pol trajectories show a tighter flap than the nopol trajectories as far as the flap-flap distance and the variation of Gly48-Gly49-Ile50 angle are concerned. In this section, different interactions in the flap regions are investigated. The flap is consisted of mainly hydrophobic residues and several Glycines which contribute to its extreme floppy character. Possibilities of hydrogen bonding exist between interflap residues 50 and 51′ (and 50′ and 51). Intraflap hydrogen bonds (H-bond) involving residues from two sides of the β-hairpin loop are also possible. Toth et al. examined the weak polar interactions involving CH and O=C between residues coming from same flap and two different flaps and proposed to be responsible for the stability of the semi-open structures. Also, H-bond between water molecules and N and/or O atoms of the residues in flap is possible. Moreover, formation of a hydrophobic cluster involving residues 50, 81 has been investigated by Scott and Schiffer (2000) and Toth and Borics (2006). All these interactions, which may be playing important roles in flap movement, are investigated in the present work. Table 2 shows the details of H-bond occupancy and average distance between the donor and acceptor atoms. The interflap (51-50′ and 50-51′) H-bonds are present for a very small timemaximum is approximately 18% in all four JE2-unbound trajectories. However, interestingly, for the pol trajectories, the interflap H-bond distances are slightly less (3.4-3.6 Å as compared to  3.8-4.0 Å) than the two nopol trajectories. For intraflap case, H-bonds are similar for both pol and nopol except for the 52-49 H-bond. For pol, it is essentially always present with a shorter distance between 52 O and 49 N, whereas for nopol it is present 55 and 70% in the two trajectories with longer distances (2.9 Å in pol as compared to 3.5 and 3.3 Å in the two nopols). This implies that the H-bonding between 52 O and 49 N is stronger in the case of pol. It is interesting to note that the difference in 49-52 H-bonding pattern is consistent with the larger flap tip movement for nopol trajectories. This indicates a model in which the flap tip of β-hairpin loop is fluctuating most, whereas other parts have much less fluctuation.

Weak polar interactions in the two flaps
There are several weak polar interactions present between flap residues. We have examined (a) G49A Oc=o and I50B Hc α (and G49B Oc=o and I50A Hc α ) and (b) G51A(B) Oc=o to F53B(A) H δ (H ɛ )) for all the five trajectories. There are distinct patterns for G49A-I50B interactions for the different trajectories. Figure 12 shows that for the initial 7 and 8 ns, in the two nopol (NOPOL-1 & -2) trajectories, the G49A Oc=o -I50B Hc α distance is clearly several angstrom more than that observed in the two pol trajectories. In the pol trajectories, the distance is around 8 Å most of the time, while in the nopol (NOPOL-1 & -2) trajectories the distance is mostly between 9 and 14 Å. For the other combination (i.e. G49B and I50A), the differences are reduced. The interaction between 51 and 53 is found to have similar for all large four JE2-unbound trajectories presumably because of the rotation of the phenyl ring of Phenylalanine. These results indicate that it is likely that the pol force field is causing a stronger polar interaction between the residues 49 and 50 present in two different flaps.
However, the interaction seems to be very stable for NOPOL-3, with a range around 7 Å, throughout the simulation period.

Hydrophobic cluster
There exists a hydrophobic cluster involving the sidechains of I50B-I50A-P81B (and I50A-I50B-P81A) in the crystal structure. It is of interest to check whether this cluster exists during the simulations and what happens when the flaps separate from each other. Importance of hydrophobic cluster has been suggested by previous workers. Scott and Schiffer (2000) suggested that the formation of hydrophobic cluster within each chain of HIV-pr opens up the active site for entrance of the ligand. Toth and Borics (2006) also discussed hydrophobic cluster formation in their implicit solvent simulations. They classified the structures based on the hydrophobic cluster formed with different residues. Recently, it was found that the rearrangement of the hydrophobic core is associated with the drug resistance of HIV-pr and may modulate the core flexibility (Mittal, Cai, Nalam, Bolon, & Schiffer, 2012). We found that hydrophobic cluster is preserved almost always in all four JE2-unbound trajectories, even when the flap-flap distance is large. A closer examination revealed that there are a variety of ways flap can move away from each other, and most of the times it tries to preserve the hydrophobic cluster involving 50B-50A-81B or at least 50A-81B. That is hydrophobic cluster between the two chains is found to be preserved for almost all the semi-open structures obtained from the simulations. The flaps move not just sideways but at a direction perpendicular to it, so that the flap-flap distances increases but on bending flap tip can still form hydrophobic cluster with the Proline of the other chain. Figure 13(a) shows the histogram of the distance between I50AC β -P81BC γ . It is clearly seen that the distributions are partially overlapped for all trajectories. For POL-1 and NOPOL-3 simulations, there is a single peak at around 4-6 Å as compared to the double peaks (around 4, 5 and 9 Å) for NOPOL-1 and -2 and POL-2. It indicates that, more number of closely interacting hydrophobic clusters are formed in the protein models sampled from pol simulations as compared to nopol simulations. Although, a certain percentage of the POL-2 sampled structures forms cluster in a larger range with average distance of 9 Å. Despite this, the clustering is closer and larger in case of pol simulations as compared to nopol as seen from the Figure 11(b). Figure 13 shows the histogram of the distance between I50BC β -P81AC γ . It was observed that the clustering event is more prominent here, except the NOPOL-1. There is little percentage of structures that form clusters in larger range at around 8 and 9 Å. As discussed before, the flap-flap distance does not differ significantly for pol and nopol  Table 3 shows the relative proportion of these structures in the four JE2-unbound simulations. Figure 14(a) and (b) shows the L-L and L-S structures. It is clearly seen that the hydrophobic cluster is preserved even when flap-flap distance is large for the L-S structure. Interestingly, the distance between I50-P81 in the same chain varies widely (8-17 Å in the four JE2-unbound trajectories). Several other interesting pictures came out from the hydrophobic cluster formation even for large values (as long as 11.5 Å) of flap-flap distance. This supports the findings of Scott and Schiffer(2000). First, since the flap has to bend to form the hydrophobic contact with the Pro of the other chain, the size of the cavity can be reduced as compared to that of the structures with no hydrophobic cluster formation. In other words, flap-flap distance alone may not always be a correct indicator of access to the active site. The mechanism by which the flap opens can be thought like this. Going from close to semi-open conformations the flaps try to preserve the hydrophobic cluster giving rise to a variety of structures. With more and more opening, ultimately, the hydrophobic cluster breaks resulting in open structures. While it might be difficult to find this mechanism experimentally, mutation at position 81 is likely to make the flapopening faster.

Effect of water
Effect of water molecules on the flap movement is examined. It was found that there is one water molecule forms a strong H-bond with oxygen atom of residues I50 (50′), G48 (48′), G49(G49′), G51(G51′) in about 30-70% the time depending on the residues in all four JE2-unbound trajectories. Though the general feature of this is similar for both nopol and pol, on average the H-bond distance is slightly shorter for pol trajectories.

Discussion
Several interesting features, both involving difference between the two force fields and general features independent of the force field were obtained from the present work. Outcome of the simulations show that in general the pol simulations make the flap movement more restricted than the nopol trajectories. This was found by calculating the flap-flap distance and the Gly48-Gly49-Ile50 curling angle. A closer examination revealed that the electrostatic interaction and H-bond pattern involving the flap residues differ for pol and nopol simulations. Weak polar interaction between carbonyl oxygen of G49 and alpha hydrogen of I149 is more prominent for pol simulations. Also, one intra-flap H-bond (between G52 and G49) is stronger in pol. Effect of water molecules around the active site is also examined. Water molecules in general are forming a slightly shorter H-bond in pol simulations.
Whether this difference between ff99 + TIP3P and ff02 + POL3 is generally true is difficult to say just by investigating one protein, it is likely to say that polarizable ff02 force filed samples better than non-polarizable ff99. However, the difference between fixed charge model and ff02 force field for K + -Benzene system has been found before (Ponder & Case, 2003). Also, the difference between TIP3P and POL3 has been found for water dimer geometry and interaction energy (Ponder & Case, 2003).
Of course, each of the above-mentioned effects alone is not enough to change the flap dynamics significantly. However, a combination of these is likely to have some effect of the flap movement. Since, the flap movement also depends on the other region of the protein in a cooperative way; there must be other factors as well with the factors mentioned above.
Some general conclusions on flap dynamics of HIV-pr also came out from the present study. The most important is the preserved hydrophobic cluster involving flap and active site wall residues even when flap-flap distance is considerably long. The hydrophobic cluster involving 50A-50B-81A (and the other combination) is present in the closed structure. When the flap starts to open up, it tries to preserve at least the 50B-81A even when the flap tip-flap tip distance is as long as 11.5 Å. This becomes possible when flaps do not move sideways rather at a direction perpendicular to it. Based on the degree of flap opening and the possibility of hydrophobic cluster formation, the semi-open structures can be classified. The 50B-81A hydrophobic cluster is also formed, but its population is less than the one formed with 50A-81B.
The conclusion regarding the hydrophobic cluster is more likely similar to that of Schiffer's work. Our more extensive simulation clearly indicates the importance of hydrophobic cluster formation in flap opening. It is likely that mutation of P81 may change the flap dynamics, which can be verified experimentally. Based on this hydrophobic cluster formation, inhibitors can be designed which will facilitate the hydrophobic cluster formation between the two chains. Our work do support Toth's proposal that the semi-open state can have a large variety of conformations. However, it is not just weak polar interactions which influence the flap movement. Hbonding within a flap, H-bond with water and finally hydrophobic cluster involving flap tip and active site wall residue are contributing to the extremely complex dynamics of flaps. It is certainly possible to have nonlocal interactions influencing the flap movement, which is not examined in the present work.
Regarding the observations from different force fields, it can be said that it is very likely that subtle differences between force field can influence flap dynamics. Two trajectories with same force field show that generally the qualitative conclusions are similar for both, with few exceptions.
One interesting observation is that unlike the implicit solvent simulation in Hornak et al. (2006), we found that the Gly49-Gly52 have H-bonding for substantial amount time (almost always for pol trajectories) for nopol trajectories.
There is one pertinent question regarding the similarity of the four JE2-unbound trajectoriesis it possible that the trajectories move to the different regions of the conformation space of HIV-pr and hence giving different results. In other words, the differences we have seen for pol and nopol are because of the force-fields or because of sampling of very different structures. Though a full clustering of structure is not done in the present work, RMSD of the whole protein, RMSD of the flap and other regions indicate that the global structural features are not very different among the four JE2-unbound trajectories. Only the local structural features differ. We have also examined the distances of six Cl À ions used in simulation from the protein structure. It was found that on average all Cl À ions are 25 Å or more apart from the residues present at the protein surface. It is unlikely that the Cl À ions can cause the difference seen between pol and no-pol.

Conclusion
An extensive MD simulation study is done to investigate the flap dynamics of HIV-pr. The motivation of the present work comes from the fact that though there are a number of works on HIV-pr dynamics, there has been considerable difference in the flap opening in ns time scale. It is likely that force field used in the simulations and simulation protocols may influence the extent of flap opening. In the present work, flap dynamics of HIV-pr is investigated by comparing two different force fields, AMBER ff99 (a nonpolarizable force field) and AMBER ff02 (a polarizable force field). To minimize biasness from a single trajectory, two different trajectories (with different initial velocity distributions) for each force field were run. All the simulations were started after an extensive equilibration procedure and continued for 20 ns. It was found that structures obtained from nopol and pol simulations have some important differences. Within nopol and pol trajectories, there are several differences; however, qualitatively, ff02 (pol) simulation of descriptions for protein dynamics is better than ff99 (nopol). It has been found that in general flap movement is tighter in the case of pol simulations. Both flap-flap distance and flap angle curling has more fluctuation for nopol. Examination of the possible interactions responsible for this difference revealed that for pol simulations some H-bond and polar interactions between atoms in the flaps have shorter distance as compared to that in nopol. In other words, the polar interaction and some H-bond are slightly stronger in the case of pol. Interaction between water and flap residues were also examined. That also showed the H-bonds with water molecules have slightly shorter distance for pol trajectories. In general for other force fields, subtle effects can lead to different ways of flap opening. Especially the flap tip is very sensitive to subtle effects of simulation protocol and force field. Since the number of simulations on HIV-pr will continue to be a major research area; the results should be interpreted carefully. One of the most important findings of the present work is the retention of hydrophobic cluster between the flap tip and active site wall residues even when the flap-flap distance is as long as 11.5 Å. The flaps move in such a way that on bending the tip can come close to the active site wall residues. Knowledge of this interchain hydrophobic cluster can be used to design inhibitors facilitating the interchain hydrophobic clustering. This may make the two chains more rigid with respect to each other. The effect of mutation is likely to be little in such a condition.

Supplementary material
The supplementary material for this paper is available online at http://dx.doi.10.1080/07391102.2013.795873.