Exploring the Reasons for Decrease in Binding Affinity of HIV-2 Against HIV-1 Protease Complex Using Interaction Entropy Under Polarized Force Field

Introduction

Acquired immunodeficiency syndrome (AIDS) has been a major global health challenge due to its serious threat to human life. According to data of The Joint United Nations Programme on HIV/AIDS (UNAIDS), around 36.7 million people within the worldwide are infected with HIV in recent years. Currently, there are two major etiological causative agents of AIDS: human immunodeficiency virus type-1 (HIV-1) and type-2 (HIV-2). As the most common type, HIV-1 which is divided into four groups (M, N, O, and P) and several subtypes (Tie et al., 2012), is observed around the global. Another type, HIV-2, is firstly identified in West Africa (Guyader et al., 1987; Menéndez-Arias and Tözsér, 2008; Peterson et al., 2011), but now it is found which is gradually spreading to other parts of the world (Soriano et al., 2000; Barin et al., 2007; Kovalevsky et al., 2008).

HIV-1 protease (PR1) (Navia et al., 1989) is extremely effective as a drug target for AIDS treatment (Debouck, 1992). It plays an essential role in the life cycle of HIV by cleaving the Gag and Gag-pol non-functional polypeptides into mature and infectious HIV viral particles (Wlodawer and Vondrasek, 1998; Louis et al., 2000). PR1 inhibitors can combine preemptively the active-site cavity of PR1, and at the same time, it can prevent the hydrolysis of the viral Gag and Gal-Pol polypeptides, which results in immature and non-infectious viral particles. However, due to the lack of specifically targeted inhibition for HIV-2 protease (PR2), the study for those patients infected by HIV-2 faces significant challenges. PR1 and PR2 share ~50–70% sequence identity (Gustchina and Weber, 1991; Tözsér et al., 1991) and they have very similar tertiary structure, while PR1 inhibitors show lower efficiency and weaker inhibition for PR2 (Brower et al., 2008; Peterson et al., 2011). Therefore, the design of new inhibitor specially targeting PR2 is urgent for those treatment of patients infected by HIV-2.

The family of HIV protease (PR) is a kind of homodimer enzyme, which consists of two identical 99 amino acid monomers. The flap regions (flap A: residues 43–58 of A chain, flap B: 43'-58' of B chain) form two β hairpin structures at the active site, respectively. They control the substrate to enter or leave the substrate-binding cavity by closing or opening pocket of cavity (Miller et al., 1989; Hornak et al., 2006b). Furthermore, the bridging water W301 plays a significant role in opening and closing of the flap as well as the binding for HIV and inhibitor (Lu et al., 2006; Barillari et al., 2007; Duan et al., 2007). W301 forms four hydrogen bonds at the active site together, two with flap tip (Ile50 and Ile50′) and two with inhibitor (shown in Figure 1). In addition, PR has five alternative protonation states: unprotonation (unpro), mono protonation of Asp25 of A chain on OD1 (25OD1), mono protonation of Asp25 on OD2 (25OD2), mono protonation of Asp25′ of B chain on OD1 (25′OD1), and mono protonation of Asp25′ on OD2 (25′OD2). For PR1, the mono protonation state is generally assigned to OD2 of Asp25′ (Hou and Yu, 2007; Chen et al., 2014; Duan et al., 2015); however, the protonation state of PR2 is unknown up to now.

FIGURE 1

Figure 1. The native structures of molecular: (A) complexes of PR2 and two inhibitors (DRV is cyan, APV is magenta); (B) Inhibitor DRV; (C) Inhibitor APV.

Although many theorists and experimenters have studied the binding mechanism of PR1 and inhibitors, the research about PR2 is still deficient. Clavel et al. analyzed the nucleotide sequence of the human retrovirus associated with PR2 and found that their biological properties are conserved in spite of limited sequence homology in PR1 and PR2 system (Guyader et al., 1987). Weber et al. compared sequences and structures of PR1 and PR2 and found that the origin of difference in affinity was partly attributable to the change of Val32 in PR1 to Ile32 in PR2 (Gustchina and Weber, 1991). Freire et al. performed kinetic inhibition assays to measure the inhibition constants (Ki) of the HIV-1 protease inhibitors against the HIV-2 protease and found that inhibitors show Ki ratios ranging between 2 and 80 for HIV-2 and HIV-1 proteases, respectively (Brower et al., 2008). Recently, Regad et al. compared the binding pockets of PR1 and PR2 and found that pockets of bound PR2 were more hydrophobic with more oxygen atoms and fewer nitrogen atoms than pockets of PR1 (Triki et al., 2018). In this paper, our work focuses on the investigation for the differences of binding patterns between two PRs and inhibitors, besides discussing the origin of the decrease in affinity of two inhibitors against HIV-2 protease compared with HIV-1 protease by molecular dynamics (MD) simulation and the binding free energy calculation.

As we know, the accuracy of the dynamics results obtained from MD simulation mainly depends on the accuracy of the force field adopted. Current standard force fields, such as AMBER, GROMACS, CHARMM and so on, lack the electronic polarization effect and this leads to inaccurate and unreliable results. Recent years, some polarizable force fields (Halgren and Damm, 2001), such as AMOEBA (Ponder et al., 2010) in AMBER package, have been gradually developed, however, the practical application is more complicated than non-polarizable force field. Currently, several models including polarization effects, such as the fluctuating charge model (Banks et al., 1998), induced multipole (Ren and Ponder, 2003), Drude oscillator (Lamoureux et al., 2003), and quantum mechanical treatment (Gordon et al., 2013), have been proposes. The uncertainty of the accuracy and validity of the basic theoretical model used to derive polarizable force field still exists (Warshel et al., 2007). In our study, the polarized protein-specific charge (PPC) force field (Ji et al., 2008) based on quantum mechanical calculation is employed. It can provide more accurate electrostatic interaction for PR and inhibitor than traditional non-polarizable force field. Extensive studies have also demonstrated that the electronic polarization effect plays an essential role in MD simulations (Ji and Zhang, 2008; Duan et al., 2010, 2016a; Ji and Mei, 2014).

In the binding free energy calculation, it is widely believed that free energy perturbation (FEP) (Bash et al., 1987; Rao et al., 1987) and thermodynamic integration (TI) (Straatsma and Berendsen, 1988; Åqvist et al., 1994) are the most accurate and rigorous methods among all methods. However, the low computational efficiency and the lack of the ability to calculate absolute free energy are the common shortcomings in the two methods. In contrast, as another effective method for free energy calculation, the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) (Srinivasan et al., 1998; Wang and Kollman, 2000; Wang J. et al., 2001; Wang W. et al., 2001; Xu et al., 2013; Sun et al., 2014; Genheden and Ryde, 2015; Chen et al., 2016) has become a very popular method duo to its full functioning and efficiency in the absolute free energy calculation. However, the entropy change calculation in MM/PBSA method is suffering the developmental bottleneck. The standard normal mode (Nmode) (Nguyen and Case, 1985) in the MM/PBSA is time-consuming and approximate for the entropic calculation. The method considers the internal motions of protein as superposition of vibrations with different frequencies (Xu et al., 2011), and then calculates the vibrational entropy. Besides, the calculation of Hessian matrix of the energy (second derivative of energy) is extremely costly in multiple degrees of freedom (Yan et al., 2017). Ray Luo, the developer of MM/PBSA, finds that the normal mode approximation does not benefit too much to the quality of the MM/PBSA calculations (Wang et al., 2016). So, it is promising to further improve the accuracy of MM/PBSA methods by improving the accuracy of the calculation of entropy change using a highly efficient and reliable method. In this paper, a novel and effective method, namely interaction entropy (IE) (Duan et al., 2016b, 2017), is used to calculate the entropy change. The method is more theoretically rigorous and time-saving than traditional Nmode method (Cong et al., 2017; Yan et al., 2017; Liu et al., 2018; Qiu et al., 2018; Song et al., 2018).

In addition, as the main tool to search for hot-spots residues, the residue decomposition of binding free energy in MM/GBSA method has always been used to explore the binding mechanism of protein and ligand. However, it is difficult to decompose entropy change into each residue using traditional Nmode method. Therefore, the effect of entropic contribution on the calculation for predicting hot-spots residues is generally ignored. Fortunately, the IE method is able to calculate straightforwardly the entropy change of each residue. The combination of MM/GBSA and IE will further provide detailed and reliable information for predicting hot spots and understanding the interaction mechanism between protein and ligand.

Method

Polarized Protein-Specific Charge

PPC force field based on quantum mechanical calculation can provide accurate partial atomic charges of proteins to represent electrostatic polarization effect. With the help of molecular fragmentation using conjugated caps scheme (MFCC) method (Zhang and Zhang, 2003), we firstly cleave the protein into fragments at the peptide bond and add a pair of conjugate caps to achieve the electronic structure of protein by fully quantum mechanical (QM). Next the restrained electrostatic potentials (RESP) (Bayly et al., 1993) program is employed to fit atomic charges. Then the Poisson-Boltzmann (PB) Equation (Tannor et al., 1994) is used to calculate electrostatic solvation energy and generate discrete induced charges on cavity surface. The newly obtained charges of other residues and surface charges mimicking the solvent effect are regarded as background charges in the next cycle of QM calculation for each fragment. Finally, the solute and solvent polarize each other until solvation energy and induced charges converge.

MM/PBSA

As one of the most widely used methods, MM/PBSA has always played a significant role on calculation of binding free energy. In the MM/PBSA approach, the binding free energy ΔG_bind is calculated as follows:

$\begin{array}{cc}\begin{array}{r}\Delta G_{bind}=\Delta H-T\Delta S=\langle E_{pl}^{int}\rangle +\Delta G_{sol}-T\Delta S\end{array}& (1)\end{array}$

where the ΔH represents enthalpy change, and the $\langle E_{pl}^{int}\rangle$ represents the ensemble averaged protein-ligand interaction including electrostatic interaction and van der Waals (vdW) interaction. ΔG_sol and −TΔS represent the solvation free energy and contribution of entropy change, respectively. Among them, the ΔG_sol can be divided into two parts:

$\begin{array}{cc}\begin{array}{r}\Delta G_{sol}=\Delta G_{pb}+\Delta G_{np}\end{array}& (2)\end{array}$

where ΔG_pb and ΔG_np represent the polar and non-polar solvation free energy, respectively. The ΔG_pb can be calculated through the PB equation. In our calculation, the exterior and interior dielectric constants are set to 80 and 1, respectively. In the meantime, the ΔG_np can be calculated through the following equation:

$\begin{array}{cc}\begin{array}{r}\Delta G_{np}=\gamma \times SASA+\beta \end{array}& (3)\end{array}$

where SASA represents solvent-accessible surface area, and it can be calculated by MSMS program (Sanner et al., 1996). The numerical values of γ and β are the standard values of 0.00542 kcal/(mol · Ų) and 0.92 kcal/mol, respectively.

In our works, the calculation of enthalpy change is performed based on 100 snapshots from MD trajectories. Considering that the calculation of entropy change is extremely time-consuming, we only extract averagely 20 snapshots from the above 100 snapshots to finish the calculation of entropy change by Nmode method.

Interaction Entropy

In addition to Nmode, a new more rigorous and concise method, interaction entropy (IE), is employed to calculate entropy change. It can be defined as:

$\begin{array}{cc}\begin{array}{r}-T\Delta S=KTln\langle e^{\beta \Delta E_{pl}^{int}}\rangle \end{array}& (4)\end{array}$

where $\Delta E_{pl}^{int}$ represents the fluctuation of protein-ligand interaction energy ($E_{pl}^{int}$) around the average energy ($\langle E_{pl}^{int}\rangle$). It can be calculated as:

$\begin{array}{cc}\begin{array}{r}\Delta E_{pl}^{int}=E_{pl}^{int}-\langle E_{pl}^{int}\rangle \end{array}& (5)\end{array}$

the protein-ligand interaction energy ($E_{pl}^{int}$) consists of electrostatic interaction and vdW interactions. The efficiency of this approach lies in that the two averages $\langle E_{pl}^{int}\rangle$ and $\langle e^{\beta \Delta E_{pl}^{int}}\rangle$ can be calculated simply by the following equations:

$\begin{array}{cc}\begin{array}{r}\langle E_{pl}^{int}\rangle =\frac{1}{N}{\displaystyle \sum _i^N}{E}_{pl}^{int}(t_i)\end{array}& (6)\end{array}$

and

$\begin{array}{cc}\begin{array}{r}\langle e^{\beta \Delta E_{pl}^{int}}\rangle =\frac{1}{N}{\displaystyle \sum _i^N}{e}^{\beta \Delta E_{pl}^{int}(t_i)}\end{array}& (7)\end{array}$

where β is 1/KT.

In the IE method, the residue decomposition of entropy change (Wang et al., 2017; Yan et al., 2017) is performed by the following equations:

$\begin{array}{cc}\begin{array}{r}-T\Delta S_{rl}=KTln\langle e^{\beta \Delta E_{rl}^{int}}\rangle \end{array}& (8)\end{array}$

where $\Delta E_{rl}^{int}$ represents the fluctuation of residues-ligand interaction energy ($E_{rl}^{int}$) around the average energy ($\langle E_{rl}^{int}\rangle$).

MD Simulation

The initial structures of DRV-PR1, DRV-PR2, AVP-PR1, and AVP-PR2 is obtained from the Protein Databank (PDB entry: 4LL3, 3EBZ, 3EKV, and 3S45, respectively). All water molecules, but the bridging water molecule W301, are removed from initial structures. For PR1, the mono protonated state is assigned to OD2 of Asp25′; For PR2, the five alternative protonation states all are considered. The geometry of ligand is optimized at the HF/6-31G** level using Gaussian03 program (Frisch et al., 2003). Single point energy is calculated at the B3LYP/cc-PVTZ level to generate electrostatic potentials (ESP) and the atomic charges are fitted using the restrained ESP (RESP) (Cornell et al., 1993) method. The parameters of the protein and ligands are generated based on AMBER12SB force field and the general AMBER force field (GAFF) (Hornak et al., 2006a), respectively. The TIP3P water box is chosen as the solvent environment with 10 Å buffer between the complex and extremity of the water box. Counterions are added to keep the system electrically neutral. Energy minimization is carried out by the steepest descent method followed by the conjugate gradient minimization to remove bad contacts between the solute and solvent water molecules. Firstly, solvent water molecules are optimized by holding the complex fixed with force constant of 500 kcal/(mol · Ų). Secondly, the entire system is energy minimizes without any restrictions.

Before running MD simulation, the whole system is heated from 0 to 300 k with 10 kcal/(mol · Ų) harmonic constraints on all solute atoms for 300 ps. SHAKE (Ryckaert et al., 1977) algorithm is used to constrain all chemical bonds involving hydrogen atoms. Finally, MD simulation is performed for 100 ns in the explicit water using AMBER16 (Case et al., 2016). In the first 80 ns, the time of writing to the coordinated file is set to 4 ps. In the last 20 ns, the time of writing to the coordinated file is set to10 fs to obtain enough conformational sampling. The MD simulations using PPC force field is the same with AMBER force field except that the standard charges from AMBER force field are simply replaced by the PPC method.

Results and Analyses

Analysis of the Stability

In order to evaluate the stability of the system during MD simulations, the root mean square deviation (RMSD) of the protein backbone atoms relative to the corresponding crystal structure is calculated. And these results are shown on Figures S1, S2 in the Supporting Information. Obviously, the values of RMSD are fluctuated between 0.5 and 1.5 Å, which indicates that simulation has reached equilibrium. Moreover, the distances between W301 and inhibitor are calculated and plotted on Figures S3, S4 in the Supporting Information. It can be clearly observed in the six systems for DRV-PR2(unpro), DRV-PR2(25OD1), DRV-PR2(25OD2), APV-PR1(25′OD2), APV-PR2(unpro), and APV-PR2(25OD1), the bridging water W301 leaves the binding pocket after around 50, 80, 70, 20, 60, and 20 ns under AMBER force field, respectively. And there is no other water molecule taking its place in the subsequent MD simulation. Therefore, the above RMSD obtained from AMBER force field is calculated before the bridging water is far away from the binding pocket. In contrast, W301 tightly binds to its original position under PPC simulation using the fluctuation of hydrogen bond length for about 3.5 Å.

In order to explore the reason, those charges of bridging water molecules and the surrounding atoms that form hydrogen bonds with W301 are plotted on Figures 2, 3 under AMBER and PPC force field respectively. For most atoms, the absolute values of charges from PPC are larger than corresponding AMBER charges. The larger polarized charges can provide more attractive electrostatic interaction under PPC than AMBER force field, which results into the strong hydrogen bonds for the binding of protease-W301-inhibitor. Thus, the bridging water is very stable under PPC simulation. The phenomenon may be an excellent illustration that correctly including the polarization effect in MD simulation can give stable and reliable electrostatic environment and therefore it is of critical importance in PR-W301-inbihitor interaction.

FIGURE 2

Figure 2. The charge of bridging water molecules and surrounding atoms that form hydrogen bonds in DRV complex. MO2, MO4, and MO5 represent O2, O4, and O5 atoms of ligand molecule, respectively. WO, WH1, and WH2 represent bridging water molecule. (A) DRV-PR1(25′OD2) system. (B) DRV-PR2(unpro) system. (C) DRV-PR2(25OD1) system. (D) DRV-PR2(25OD2) system. (E) DRV-PR2(25′OD1) system. (F) DRV-PR2(25′OD2) system.

FIGURE 3

Figure 3. The charge of bridging water molecules and surrounding atoms that form hydrogen bonds in APV complex. (A) APV-PR1(25′OD2) system. (B) APV-PR2(unpro) system. (C) APV-PR2(25OD1) system. (D) APV-PR2(25OD2) system. (E) APV-PR2(25′OD1) system. (F) APV-PR2(25′OD2) system.

Protonation States of PR2

Considering that the two residues (Asp25 and Asp25′) containing charged groups are very close to each other, protonation prevents the two residues from repulsing each other because of strong Coulomb interaction. Different protonation states are validated according to the local environment. In our works, five possible protonated states of PR2 complex are performed to explore the effect of different types of protonation on structure and binding free energy.

Some average geometric parameters of the protonated region during MD simulation are shown in Table 1. Detailed results are presented at Tables S1, S2 in the Supporting Information. The selection of the six parameters (D₁~D₆) is consistent with Hou et al. (Hou and Yu, 2007). The unsigned mean error (UME) of the distance and angle is used to measure deviations between simulated structures and the native structures. The average (AVE) of six parameters (D₁~D₆) is used to measure average distance among the three (Asp25, Asp25′ and inhibitor). For DRV-PR2 complex, the average distances in the unprotonated state under AMBER and PPC are 4.18 and 3.40 Å, respectively, larger than the average distances (3.27 Å) of crystal structure and most other protonation states. For APV-PR2 complex, the average distances in the unprotonated state under AMBER and PPC are 4.69 and 4.47 Å, respectively, much larger than the average distances (3.00 Å) of crystal structure and other protonation states. The large conformation change comes from the consequence of the strong Asp-Asp repulsion. Therefore, the unprotonated state is not reasonable in our simulation.

TABLE 1

Table 1. Average values and unsigned mean errors of geometric parameters of the protonated region during MD simulation.

In contrast, these parameters of distance and angle in the 25′OD1 protonated state under PPC force field are the closest crystal structure among these protonated states. For DRV-PR2 complex, the UME of distance parameters in the 25′OD1 protonated state under PPC force field is 0.26, 0.18 Å less than the next closest state. The UME of angle parameters is 12.82°, which is also small enough. Related to APV-PR2 complex, the UMEs of distance and angle parameters in the 25′OD1 protonated state under PPC force field are 0.26 and 9.34°, respectively, 0.11 and 3.98° less than the next closest state. This indicates that the 25′OD1 protonated state under PPC force field provides better geometry structure in the protonated region than other protonated states.

The Difference in Binding Free Energies

In our works, four methods (AMBER-Nmode, AMBER-IE, PPC-Nmode, and PPC-IE) are used to calculate the binding free energy between PR and inhibitor. Under the condition of ensuring the stability of bridging water and the convergence of IE method, the binding free energy and detailed items are calculated and shown in Tables S3, S4 in the Supporting Information. It can be seen that there is considerable difference between binding free energies under different protonated state. This indicates that type of protonation does play an indispensable role in calculation of binding free energy.

In order to explore the influences of types of protonation in deep and the calculation method on binding free energy, we calculate the difference of binding free energy between PR1 and PR2 complex. Those results are shown on Table 2 and Figure 4. In Figure 4, red lines represent the difference of binding free energy of between the experimental values of PR1-inhibitor and PR2-inhibitor, which is 1.6 and 1.3 kcal/mol in DRV complexes and APV complexes, respectively. In DRV complexes, the difference of binding free energy between PR1-DRV and PR2-DRV calculated by PPC-IE method in the 25′OD1 protonated state in PR2 system is 0.78 kcal/mol. In APV complexes, it is 0.66 kcal/mol. Both of these two values are the closest to the experimental value among all combinations. And at the above condition, the rank of the predicted binding free energies is significantly consistent with the experimental rank. This indicates that binding free energy calculated by PPC-IE method in the 25′OD1 protonated state in PR2 system is the optimal choice in our study.

TABLE 2

Table 2. Binding free energy and delta of PR1 and PR2 complexes.

FIGURE 4

Figure 4. Difference between binding free energy of PR1 complex and binding free energy of PR2 complex. AM represents AMBER; Nm represents Nmode. Cyan lines represent experimental values. *The experimental value is calculated by Weber (Tie et al., 2012).

In view of the consistency of conformation and binding free energy between experimental measure and theoretical calculation, the 25′OD1 protonation model and PPC force field are used to further study.

According to Table 2, it can be observed that the computed binding free energies of DRV complexes are stronger than APV complexes in two PRs consistently. The binding free energies and detailed items of DRV and APV complexes are compared on Table 3. The G_ele+pol item contains the electrostatic interactions and polar solvation free energies. In PR1 complexes, the delta of binding free energies between APV-PR1 and DRV-PR1 is 14.13 kcal/mol. The most contribution of the difference comes from vdW interaction (delta is 6.58 kcal/mol). In PR2 complexes, the delta of binding free energies is 14.01 kcal/mol. Although those deltas of electrostatic interaction and polar solvation free energy is special significant, the two offset each other. The total delta of two is 1.33 kcal/mol, which is small enough. The most contribution still comes from vdW interaction (delta is 12.89 kcal/mol).

TABLE 3

Table 3. Binding free energy and delta of DRV and APV complexes under PPC force field.

Based on vdW interaction, the decomposition of residue is performed which is shown in Figure 5. According to Figure 5A, the vdW interaction of Ala28′ in DRV-PR1 complex is more favorable than in APV-PR1 complex with the delta value of 2.57 kcal/mol. Subsequently, the vdW interaction of per atom of Ala28′ and inhibitor is calculated and shown in Figure 6A. It can be observed that the difference of vdW interaction is mainly reflected on N, CA, CB, and C atoms. Figures 7A,B depict the geometrical positions of Ala28′ and two inhibitors based on the lowest potential energy structure from MD trajectory. These distances between above atoms and hydrophobic groups of DRV are smaller than APV, which result in stronger vdW interaction. The structural descriptions agree well with analysis of vdW interaction.

FIGURE 5

Figure 5. Decomposition of residue based on vdW interaction. (A) PR1 complexes. (B) PR2 complexes.

FIGURE 6

Figure 6. vdW interaction of per atom and inhibitor in Ala28′ of PR1 and Asp25/Asp25′ of PR2. (A) Ala28′ residues in PR1. (B) Asp25 residue in PR2; (C) Asp25′ residues in PR2.

FIGURE 7

Figure 7. The geometrical positions of Ala28′, Asp25/Asp25′ and two inhibitors based on the lowest potential energy structure from MD trajectory. (A) Ala28′ residue in DRV-PR1. (B) Ala28′ residue in APV-PR1. (C) Asp25residue in DRV-PR2. (D) Asp25 residue in APV-PR2. (E) Asp25′ residue in DRV-PR2. (F) Asp25′ residue in APV-PR2.

According to Figure 5B, the vdW interaction of Asp25 and Asp25′ in APV-PR2 complex is more unfavorable than DRV-PR2. Those deltas are 3.70 and 2.48 kcal/mol, respectively. The vdW interaction of per atom of Asp25/Asp25′ and inhibitor is calculated and shown in Figures 6B,C, respectively. The difference of vdW interaction is mainly reflected on OD2/OD1 atom in Asp25/Asp25′, respectively. Figures 7C–F depict the geometrical positions of Asp25/Asp25′ and two inhibitors based on the lowest potential energy structure from MD trajectory. These distances between OD2/OD1 atom and O3 atom of APV are smaller than in DRV complex. Excessively short distances result into stronger repulsion among these atoms. The vdW interaction OD2/OD1 atom and O3 atom in Asp25-DRV, Asp25-APV, Asp25′-DRV, and Asp25′-APV is 1.35, 5.00, 2.43, and 4.72 kcal/mol, respectively. The structural descriptions and analysis of vdW interaction are both consistent.

Bridging Water Molecule and Hydrogen Bond Network

Water molecules play a critical role in regulating the structure and function of proteins and maintaining the fundamental biological process of cells (Chaplin, 2006; Vukovic et al., 2016). Philip Ball called water “an active matrix of life for cell and molecular biology” (Ball, 2017). In this paper, the bridging water molecules and their hydrogen bond networks are analyzed to explore its impact on the binding of PR and inhibitors.

The above calculation on the binding free energy has included the effect of the bridging water W301 by considering it as part of PR. Next, in order to investigate the contribution of W301, we again calculate the binding free energy excluding the bridging water. These energies of individual items are compared with the previous energy calculated in the presence of bridging water and are shown in Table 4. Obviously, the electrostatic interaction and polar solvation free energy undergo significant changes before and after considering the bridging water. In particular, the electrostatic interaction makes strong and favorable contribution toward binding free energy at the presence of the W301. The reason is mainly from the four hydrogen bonds formed by bridging water. In PR1 complexes, the contribution of bridging water to total binding energy is −3.43 and −6.62 kcal/mol, respectively. This indicates that bridging water does play a powerful role in promoting the binding of the two PR1 complexes which is significantly consistent with other studies (Lu et al., 2006; Barillari et al., 2007; Duan et al., 2007). However, in PR2 complexes, the contribution of bridging water to total binding energy is 0.12 kcal/mol and 2.06 kcal/mol, respectively. This observation clearly shows that bridging water not only doesn't promote the binding of PR2 complexes, but also plays an inhibitory role.

TABLE 4

Table 4. Comparison of binding free energy with bridging water (WAT) and no bridging water (NW) when using the PPC-IE method in the 25′OD1 protonated state.

According to Table 4, considering that the contribution of polar solvation free energy to delta is cancelled by electrostatic interactions, vdW interaction plays a major contribution to delta of binding free energy. In our study, bridging water is treated as part of protease. Therefore, the vdW interaction of bridging water and per heavy atom of inhibitors is further calculated and they are shown in Figure 8. It can be clearly observed that vdW interaction between O4 atom of DRV and bridging water showed a stronger unfavorable contribution in PR2 complex than PR1 complex. Similarly, O5 atom of APV obtains the same conclusion as above, as well. Subsequently, the geometrical positions of bridging water and inhibitors in four complexes based on the lowest potential energy structure are depicted in Figure 9. Obviously, these distances between O4/O5 atom of inhibitors and bridging water are smaller in PR2 complexes than PR1 complexes. Excessively short distances lead to stronger repulsion between bridging water and inhibitors. The different geometrical positions may be from the binding environment where the bridging water molecules are located. It will further be analyzed in the analysis of the binding pocket section. To some extent, the unfavorable contribution of bridging water explains the reason for the decline in the effectiveness of PR2-inbibitor binding.

FIGURE 8

Figure 8. The vdW interaction of bridging water and per heavy atom of inhibitors. (A) DRV complexes. (B) APV complexes.

FIGURE 9

Figure 9. The geometrical positions of bridging water and inhibitors based on the lowest potential energy structure from MD trajectory. (A) DRV-PR1 complexes. (B) DRV-PR2 complexes. (C) APV-PR1 complexes. (D) APV-PR2 complexes.

Based on the above analysis, the electrostatic interaction makes strongly favorable contribution toward binding free energy. The main reason is from the four hydrogen bonds formed by bridging water molecules. In the next step, the analysis of hydrogen bond network of protease-W301-inhibitor is performed. Among the four hydrogen bonds, two are attached to the inhibitor and two are attached to the residues Ile50 and Ile50′. These two residues are located in the flap A and flap B tip, respectively. They control the substrate to enter or leave the substrate-binding cavity by closing or opening pocket of cavity. This means that the hydrogen bond network not only enhances electrostatic interactions, but also affects the binding of proteins and inhibitors by changing the structure and function of the binding pocket.

In our works, the average distance, angle and occupancy of each hydrogen bond during MD simulation are calculated and shown on Table 5. Of the 16 hydrogen bonds, the average distance of 13 hydrogen bonds is less than 3.5 Å, the average angle of all hydrogen bonds is greater than 120°, and the occupancy of 10 hydrogen bonds is greater than 80%. Moreover, each bridging water forms at least two stable hydrogen bonds connected to proteins and inhibitor, respectively. The above data shows the vast majority of hydrogen bonds are very stable during the simulation. This result explains the significant changes in electrostatic interactions before and after considering the bridging water.

TABLE 5

Table 5. Occupancy of hydrogen bonds of bridging water during MD simulation.

Comparison of Binding Pockets of PR1 and PR2

In order to explore the mechanism of binding affinity of inhibitor to PR2, it is essential to research comprehensively operating mode of binding pocket (Duan et al., 2015; Triki et al., 2018). As the active region of binding between protein and inhibitor, differences of binding pockets can explain essentially the reasons for the decline in the binding effectiveness of PR2 complex. The volume of binding pocket is calculated firstly through the POVME (Durrant et al., 2011, 2014) program during MD simulation. The frequency distributions of volume calculated are plotted in Figure 10. It is obviously observed that volumes of binding pocket of PR2 are both smaller compared with PR1. The average volumes of binding pocket calculated during MD simulation are 423.84, 394.00, 415.92, and 402.98 ų in DRV-PR1, DRV-PR2, APV-PR1, and APV-PR2 complexes, respectively. Our analysis shows that the significant difference of volume is mainly caused by the two possible factors: one is the difference of volume in native structure between PR1 and PR2; the other is the variation in the flexibility of residues during MD simulation. The volumes of native crystal structure are 385, 368, 353, and 360 ų, respectively, and the influence can be eliminated by calculating the delta of native volume and average volumes. Larger delta of volume indicates greater flexibility of residues during MD simulation. The delta calculated are 38.84, 26.00, 62.92, and 42.98 ų in DRV-PR1, DRV-PR2, APV-PR1, and APV-PR2 complexes, respectively. Obviously, the delta of volume in PR2 complex is still less than PR1 complex. The result indicates flexibility of residues in PR1 complex is greater during MD simulation than in PR2 complex.

FIGURE 10

Figure 10. The frequency distributions of volume of the binding pockets during MD simulation. The center of the binding pocket is set as the centroid of the inhibitor. (A) DRV complexes. (B) APV complexes.

To further analyze residual atomic flexibility, the isotropic temperature factor (B-factor) is utilize to measure atomic fluctuations of individual residues. B-factor reflects the mobility of each residue around its mean position. Figure 11 has plotted the B-factor of protein Cα atoms. The average B-factors of protein Cα in DRV-PR1, DRV-PR2, APV-PR1, and APV-PR2 complexes are 12.87, 7.69, 12.91, and 9.02 Ų, respectively. Obviously, according to above result, the B-factors of PR2 complexes are less than that of PR1 complex. Current research illustrates the mobility of PR1 complexes is greater, compared to PR2 complexes. The different behavior patterns of the two proteins explain the difference in the volume of the binding pocket.

FIGURE 11

Figure 11. The B-factor of protein Cα atoms in the four complexes during MD simulation. (A) DRV complexes. (B) APV complexes.

In addition, Figure 11 shows that the B-factor of flap A and flap B does not differ significantly between PR1 and PR2, relative to other regions. As we all know, flap dynamics not only affects the enzyme catalysis of PR, but also controls the substrate to enter or leave the substrate-binding cavity. Therefore, it is essential to further probe the local structure of this region. In this paper, the distance of the flap tip (Ile50 and Ile50′) is calculated to explore the extent of flap dynamics. The frequency distribution of distance between Ile50 (Cα) and Ile50′ (Cα) is plotted in Figure 12. The average distances are 5.95, 5.80, 5.80, and 5.70 Å in DRV-PR1, DRV-PR2, APV-PR1, and APV-PR2 complexes, respectively. It is clearly seen that the flap area with a tighter structure in PR2 complex than in PR1 complex makes the pocket of cavity tend to close and this results into a smaller volume, which in turn causes local differences of residues lining their surfaces and affects the binding mode of PR2 complex. At the same time, the tighter structure in PR2 complex leads to distances between bridging water and inhibitors shorter and thus produces stronger vdW repulsion than in PR1 complex. This explains further the reason of the unfavorable contribution of bridging water in PR2 complexes.

FIGURE 12

Figure 12. The frequency distribution of distance between Ile50 (Cα) and Ile50′ (Cα) during MD simulation. (A) DRV complexes. (B) APV complexes.

The Decomposition of Residue

Searching for hot-spot residues and exploring its interactions with inhibitors have long been an important step in understanding the binding mechanism of complex. In our work, combining the enthalpy change calculated by MM/GBSA with the entropy change calculated by IE, the binding free energy is decomposed to generate more detailed residue-inhibitor interaction spectrum.

Those results have been plotted in the Figure 13. In general, the residue-inhibitor interaction spectrums are extremely similar on the four complexes. Among them, regions of binding pocket (Ala28/Ala28′, Ile50/Ile50′, and Ile84/Ile84′ residues) provide major favorable contribution on binding free energy in the four systems. Those energies contributed from those hot-spot residues are greater than −2 kcal/mol. Further, the contribution of these six hot-spot residues is divided into electrostatic interaction, vdW interaction, polar solvation energy, non-polar solvation energy and entropy change. The result of four complexes is plotted in Figure 14. Obviously, the favorable contribution mainly comes from electrostatic interaction and vdW interaction. Among them, vdW interaction provides the most favorable contribution in the range of −2.3 to −5.8 kcal/mol. This result is basically consistent with the previous other research about PR1 (Meher and Wang, 2012).

FIGURE 13

Figure 13. Decomposition of the binding free energy on a per-residue basis. (A) DRV-PR1 complexes. (B) DRV-PR2 complexes. (C) APV-PR1 complexes. (D) APV-PR2 complexes.

FIGURE 14

Figure 14. Decomposition of the binding free energy on a per-basis into contributions from electrostatic interaction (ELE), van der Waals interaction (vdW), polar solvation energy (Pol), non-polar solvation energy (Nopol), and entropy change (-TΔS) in the six hot-spot residues. (A) DRV-PR1 complexes. (B) DRV-PR2 complexes. (C) APV-PR1 complexes. (D) APV-PR2 complexes.

For the residue-inhibitor interaction spectrum shown in Figure 13, it is worthwhile that Asp30 residue demonstrates great difference on PR1 and PR2 complexes. The energies of inhibitor-residue interaction are 1.60, 8.98, −0.06, and 9.01 kcal/mol in DRV-PR1, DRV-PR2, APV-PR1, and APV-PR2 complexes, respectively. The Asp30 residue contributes unfavorably to the binding of inhibitor against PR2. We further decomposed the binding energy of inhibitor and Asp30 residue for the four complexes shown in Figure 15. Although electrostatic interaction and polar solvation free energy show great differences between PR1 complexes and PR2 complexes. The two offset each other, and the sum of electrostatic interaction and polar solvation free energy is basically same in the two PR systems. The differences are mainly from the entropy change. When DRV is combined with PR1 and PR2, the difference of entropy change is 6.19 kcal/mol. For APV inhibitor, the differences are 8.41 kcal/mol. It may be another reason that causes the decline of the affinity in PR2 complex. The above information is mainly provided by the IE method, which emphasizes the importance of entropy change in the decomposition of residues.

FIGURE 15

Figure 15. Decomposition of the binding free energy on a per-basis into contributions from electrostatic interaction (ELE), van der Waals interaction (vdW), polar solvation energy (Pol), the sum of electrostatic interaction and polar solvation energy (ELE+Pol), non-polar solvation energy (Nopol), and entropy change (-TΔS) in the Asp30 residues. (A) DRV complexes. (B) APV complexes.

Conclusions

The current study probes the differences in binding patterns between PR1 and PR2 with inhibitors (DRV and APV) by molecular dynamics simulation. Two force fields, non-polarized traditional AMBER force field and polarized PPC force field, combined two methods for computing entropy change (traditional Nmode and newly developed IE method) are employed to explore the effect of electrostatic polarization on the simulation process and the effectiveness of binding. Our studies find that polarized force field is able to provide a relatively stable binding environment during MD simulation by analysis of stability and protonation states. In addition, polarized force field coupled with the highly efficient IE method shows more significantly consistent with the experimental result in the computation of binding free energy, compared to AMBER force field and Nmode method. Our result shows the following features:

Firstly, five different protonation states of PR2 have been studied in detail through some relevant geometric parameters and the binding free energy calculation. The 25′OD1 protonated state in B chain is found to have the more similar geometry with the crystal structure than other protonated states in two PR2 complexes consistently.

Secondly, the calculated binding free energy is significantly consistent with the experimental observation under the 25′OD1 protonated state in PR2 complexes. Further, the binding free energies of DRV complexes are stronger than APV complexes. The reason is mainly caused by vdW interaction of Ala28′ and Asp25/Asp25′. On the one hand, vdW interaction of Ala28′ in DRV-PR1 complex is more favorable than in APV-PR1 complex; on the other head, the vdW interaction of Asp25 and Asp25′ in APV-PR2 complex is more unfavorable than in DRV-PR2.

Thirdly, we investigate the important impact of bridging water molecule W301 on the effectiveness of binding and our study finds that the bridging water is capable of stably linking proteins and inhibitors through stable hydrogen bonding under PPC force field whether PR1 system or PR2 system. Besides, the bridging water contributes favorably to the binding of two PR1 complexes. However, it is unfavorable to two PR2 complexes. The reason is mainly from that tighter structure of PR2 complex, which results into stronger vdW repulsion between bridging water and inhibiter than that in the PR1 complex. To some extent, it explains the reason for the decline in the effectiveness of inhibitor against PR2.

Fourthly, analysis of the binding pocket finds the B-factor of Cα atoms and distance of flap tip in PR2 complexes all are less than that in PR1 complexes. The poorer flexibility and tighter structure of flap region in PR2 complex make the pocket of cavity closer and the volume smaller than in PR1 complex, which in turn causes local differences of residues lining their surfaces and thus affects the bonding mode of PR2 complex.

Fifthly, the method for the decomposition of residue in traditional MM/GBSA calculation usually neglects entropic contribution due to the difficulty using Nmode method. The IE method is proposed to efficiently compute the residue-specific entropic contribution to PR-inhibitor binding free energy. Using this method, the predicted hot-spot resides (Ala28/Ala28′, Ile50/Ile50′ and Ile84/Ile84′) are found to be nearly similar in the two types PR systems. However, Asp30 residue shows more unfavorable entropic contributions on PR2 complexes, compared with PR1 complexes. The phenomena is hard to obtain using traditional Nmode method. It may be another reason that results in the decline in potency of DRV and APV against PR2 relative to PR1.

We hope that current study will provide essential value and theoretical guidance for the future design of new dual inhibitors targeting HIV proteases.

Author Contributions

YC and YL outperformed the MD simulations, drafted the main text of the manuscript and prepared all the figures. KJ and SZ helped with data analysis. JZ, HL, and LD designed this study and revised the manuscript.

Funding

This work was supported by National Key R&D Program of China (Grant no. 2016YFA0501700), National Natural Science Foundation of China (Grant nos. 11774207, 11574184, 21433004, 91753103), the Natural Science Foundation of Shandong Province (ZR2016JL003), Primary Research & Development Plan of Shandong Province (No. 2017GSF18108), Innovation Program of Shanghai Municipal Education Commission (201701070005E00020), and NYU Global Seed Grant. We thank the Supercomputer Center of East China Normal University for providing us computer time.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer XA and handling Editor declared their shared affiliation.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2018.00380/full#supplementary-material

References

Åqvist, J., Medina, C., and Samuelsson, J. E. (1994). A new method for predicting binding affinity in computer-aided drug design. Protein Eng. 7, 385–391. doi: 10.1093/protein/7.3.385

PubMed Abstract | CrossRef Full Text | Google Scholar

Ball, P. (2017). Water is an active matrix of life for cell and molecular biology. Proc. Natl. Acad. Sci.U.S.A. 114, 13327–13335. doi: 10.1073/pnas.1703781114

PubMed Abstract | CrossRef Full Text | Google Scholar

Banks, J. L., Kaminski, G. A., Zhou, R., Mainz, D. T., Berne, B. J., and Friesner, R. A. (1998). Parametrizing a polarizable force field from ab initio data. I. The fluctuating point charge model. J. Chem. Phys. 110, 741–754. doi: 10.1063/1.478043

CrossRef Full Text | Google Scholar

Barillari, C., Taylor, J., Viner, R., and Essex, J. W. (2007). Classification of water molecules in protein binding sites. J. Am. Chem. Soc. 129, 2577–2587. doi: 10.1021/ja066980q

PubMed Abstract | CrossRef Full Text | Google Scholar

Barin, F., Cazein, F., Lot, F., Pillonel, J., Brunet, S., Thierry, D., et al. (2007). Prevalence of HIV-2 and HIV-1 group O infections among new HIV diagnoses in France: 2003–2006. Aids 21, 2351–2353. doi: 10.1097/QAD.0b013e3282f15637

PubMed Abstract | CrossRef Full Text | Google Scholar

Bash, P. A., Field, M. J., and Karplus, M. (1987). Free energy perturbation method for chemical reactions in the condensed phase: a dynamic approach based on a combined quantum and molecular mechanics potential. J. Am. Chem. Soc. 109, 8092–8094. doi: 10.1021/ja00260a028

CrossRef Full Text | Google Scholar

Bayly, C. I., Cieplak, P., Cornell, W., and Kollman, P. A. (1993). A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 97, 10269–10280. doi: 10.1021/j100142a004

CrossRef Full Text | Google Scholar

Brower, E. T., Bacha, U. M., Kawasaki, Y., and Freire, E. (2008). Inhibition of HIV-2 Protease by HIV-1 protease inhibitors in clinical use. Chem. Biol. Drug Des. 71, 298–305. doi: 10.1111/j.1747-0285.2008.00647.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Betz, R. M., Cerutti, D. S., Cheatham, T. E., Darden, T. A., Duke, R. E., et al. (2016). AMBER 16. University of California, San Francisco, CA.

Chaplin, M. (2006). Do we underestimate the importance of water in cell biology? Nat. Rev. Mol. Cell Biol. 7, 861–866. doi: 10.1038/nrm2021

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, F., Liu, H., Sun, H., Pan, P., Li, Y., Li, D., et al. (2016). Assessing the performance of the MM/PBSA and MM/GBSA methods. 6. Capability to predict protein–protein binding free energies and re-rank binding poses generated by protein–protein docking. Phys. Chem. Chem. Phys. 18, 22129–22139. doi: 10.1039/C6CP03670H

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Liang, Z., Wang, W., Yi, C., Zhang, S., and Zhang, Q. (2014). Revealing origin of decrease in potency of darunavir and amprenavir against HIV-2 relative to HIV-1 protease by molecular dynamics simulations. Sci. Rep. 4:6872. doi: 10.1038/srep06872

PubMed Abstract | CrossRef Full Text | Google Scholar

Cong, Y., Li, M., Feng, G., Li, Y., Wang, X., and Duan, L. (2017). Trypsin-Ligand binding affinities calculated using an effective interaction entropy method under polarized force field. Sci. Rep. 7:17708. doi: 10.1038/s41598-017-17868-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Cornell, W. D., Cieplak, P., Bayly, C. I., and Kollman, P. A. (1993). Application of RESP charges to calculate conformational energies, hydrogen bond energies, and free energies of solvation. J. Am. Chem. Soc. 115, 9620–9631. doi: 10.1021/ja00074a030

CrossRef Full Text | Google Scholar

Debouck, C. (1992). The HIV-1 protease as a therapeutic target for AIDS. AIDS Res. Hum. Retroviruses 8, 153–164. doi: 10.1089/aid.1992.8.153

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, L., Feng, G., Wang, X., Wang, L., and Zhang, Q. (2017). Effect of electrostatic polarization and bridging water on CDK2–ligand binding affinities calculated using a highly efficient interaction entropy method. Phys. Chem. Chem. Phys. 19, 10140–10152. doi: 10.1039/C7CP00841D

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, L., Liu, X., and Zhang, J. Z. (2016b). Interaction entropy: a new paradigm for highly efficient and reliable computation of protein–ligand binding free energy. J. Am. Chem. Soc. 138, 5722–5728. doi: 10.1021/jacs.6b02682

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, L. L., Feng, G. Q., and Zhang, Q. G. (2016a). Large-scale molecular dynamics simulation: effect of polarization on thrombin-ligand binding energy. Sci. Rep. 6:31488. doi: 10.1038/srep31488

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, L. L., Mei, Y., Zhang, D., Zhang, Q. G., and Zhang, J. Z. H. (2010). Folding of a helix at room temperature is critically aided by electrostatic polarization of intraprotein hydrogen bonds. J. Am. Chem. Soc. 132, 11159–11164. doi: 10.1021/ja102735g

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, L. L., Tong, Y., Mei, Y., Zhang, Q. G., and Zhang, J. Z. (2007). Quantum study of HIV-1 protease-bridge water interaction. J. Chem. Phys. 127:145101. doi: 10.1063/1.2770720

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, R., Lazim, R., and Zhang, D. (2015). Understanding the basis of I50V-induced affinity decrease in HIV-1 protease via molecular dynamics simulations using polarized force field. J. Comput. Chem. 36, 1885–1892. doi: 10.1002/jcc.24020

PubMed Abstract | CrossRef Full Text | Google Scholar

Durrant, J. D., de Oliveira, C. A., and McCammon, J. A. (2011). POVME: an algorithm for measuring binding-pocket volumes. J. Mol. Graph Model. 29, 773–776. doi: 10.1016/j.jmgm.2010.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Durrant, J. D., Votapka, L., Sørensen, J., and Amaro, R. E. (2014). POVME 2.0: an enhanced tool for determining pocket shape and volume characteristics. J. Chem. Theory Comput. 10, 5047–5056. doi: 10.1021/ct500381c

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2003). Gaussian 03, Revision, A.1. Pittsburgh, PA: Gaussian, Inc.

Google Scholar

Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 10, 449–461. doi: 10.1517/17460441.2015.1032936

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, M. S., Smith, Q. A., Xu, P., and Slipchenko, L. V. (2013). Accurate first principles model potentials for intermolecular interactions. Annu. Rev. Phys. Chem. 64, 553–578. doi: 10.1146/annurev-physchem-040412-110031

PubMed Abstract | CrossRef Full Text | Google Scholar

Gustchina, A., and Weber, I. T. (1991). Comparative analysis of the sequences and structures of HIV-1 and HIV-2 proteases. Proteins 10, 325–339. doi: 10.1002/prot.340100406

PubMed Abstract | CrossRef Full Text | Google Scholar

Guyader, M., Emerman, M., Sonigo, P., Clavel, F., Montagnier, L., and Alizon, M. (1987). Genome organization and transactivation of the human immunodeficiency virus type 2. Nature 326, 662–669. doi: 10.1038/326662a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Halgren, T. A., and Damm, W. (2001). Polarizable force fields. Curr. Opin. Struc. Biol. 11, 236–242. doi: 10.1016/S0959-440X(00)00196-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006a). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123

PubMed Abstract | CrossRef Full Text | Google Scholar

Hornak, V., Okur, A., Rizzo, R. C., and Simmerling, C. (2006b). HIV-1 protease flaps spontaneously open and reclose in molecular dynamics simulations. Proc. Natl. Acad. Sci. U.S.A. 103, 915–920. doi: 10.1073/pnas.0508452103

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, T., and Yu, R. (2007). Molecular dynamics and free energy studies on the wild-type and double mutant HIV-1 protease complexed with amprenavir and two amprenavir-related inhibitors: mechanism for binding and drug resistance. J. Med. Chem. 50, 1177–1188. doi: 10.1021/jm0609162

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, C., and Mei, Y. (2014). Some practical approaches to treating electrostatic polarization of proteins. Acc. Chem. Res. 47, 2795–2803. doi: 10.1021/ar500094n

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, C., Mei, Y., and Zhang, J. Z. (2008). Developing polarized protein-specific charges for protein dynamics: MD free energy calculation of pKa shifts for Asp26/Asp20 in thioredoxin. Biophys. J. 95, 1080–1088. doi: 10.1529/biophysj.108.131110

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, C. G., and Zhang, J. Z. H. (2008). Protein polarization is critical to stabilizing AF-2 and helix-2′ domains in ligand binding to PPAR-γ. J. Am. Chem. Soc. 130, 17129–17133. doi: 10.1021/ja807374x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kovalevsky, A. Y., Louis, J. M., Aniana, A., Ghosh, A. K., and Weber, I. T. (2008). Structural evidence for effectiveness of darunavir and two related antiviral inhibitors against HIV-2 protease. J. Mol. Biol. 384, 178–192. doi: 10.1016/j.jmb.2008.09.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamoureux, G., MacKerell, A. D. M. Jr., and Roux, B.T. (2003). A simple polarizable model of water based on classical Drude oscillators. J. Chem. Phys. 119, 5185–5197. doi: 10.1063/1.1598191

CrossRef Full Text | Google Scholar

Liu, X., Peng, L., Zhou, Y., Zhang, Y., and Zhang, J. Z. H. (2018). Computational alanine scanning with interaction entropy for protein–ligand binding free energies. J. Chem. Theory Comput. 14, 1772–1780. doi: 10.1021/acs.jctc.7b01295

PubMed Abstract | CrossRef Full Text | Google Scholar

Louis, J. M., Weber, I. T., Tözsér, J., Clore, G. M., and Gronenborn, A. M. (2000). HIV-1 protease: maturation, enzyme specificity, and drug resistance. Adv. Pharmacol. 49, 111–146. doi: 10.1016/S1054-3589(00)49025-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Yang, C. Y., and Wang, S. (2006). Binding free energy contributions of interfacial waters in HIV-1 protease/inhibitor complexes. J. Am. Chem. Soc. 128, 11830–11839. doi: 10.1021/ja058042g

PubMed Abstract | CrossRef Full Text | Google Scholar

Meher, B. R., and Wang, Y. (2012). Interaction of I50V mutant and I50L/A71V double mutant HIV-protease with inhibitor TMC114 (darunavir): molecular dynamics simulation and binding free energy studies. J. Phys. Chem. B 116, 1884–1900. doi: 10.1021/jp2074804

PubMed Abstract | CrossRef Full Text | Google Scholar

Menéndez-Arias, L., and Tözsér, J. (2008). HIV-1 protease inhibitors: effects on HIV-2 replication and resistance. Trends Pharmacol. Sci. 29, 42–49. doi: 10.1016/j.tips.2007.10.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, M., Jaskólski, M., Rao, J. M., Leis, J., and Wlodawer, A. (1989). Crystal structure of a retroviral protease proves relationship to aspartic protease family. Nature 337, 576–579. doi: 10.1038/337576a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Navia, M. A., Fitzgerald, P. M., McKeever, B. M., Leu, C.-T., Heimbach, J. C., Herber, W. K., et al. (1989). Three-dimensional structure of aspartyl protease from human immunodeficiency virus HIV-1. Nature 337, 615–620. doi: 10.1038/337615a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, D. T., and Case, D. A. (1985). On finding stationary states on large-molecule potential energy surfaces. J. Phys. Chem. 89, 4020–4026. doi: 10.1021/j100265a018

CrossRef Full Text | Google Scholar

Peterson, K., Jallow, S., Rowland-Jones, S. L., and de Silva, T. I. (2011). Antiretroviral therapy for HIV-2 infection: recommendations for management in low-resource settings. AIDS Res. Treat. 2011, 463704–463704. doi: 10.1155/2011/463704

PubMed Abstract | CrossRef Full Text | Google Scholar

Ponder, J. W., Wu, C., Ren, P., Pande, V. S., Chodera, J. D., Schnieders, M. J., et al. (2010). Current Status of the AMOEBA polarizable force field. J. Phys. Chem. B 114, 2549–2564. doi: 10.1021/jp910674d

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, L. Q., Yan, Y. N., Sun, Z. X., Song, J. N., and Zhang, J. Z. H. (2018). Interaction entropy for computational alanine scanning in protein–protein binding. Wiley Interdiscip. Rev. 8:e1342. doi: 10.1002/wcms.1342

CrossRef Full Text | Google Scholar

Rao, S. N., Singh, U. C., Bash, P. A., and Kollman, P. A. (1987). Free energy perturbation calculations on binding and catalysis after mutating Asn 155 in subtilisin. Nature 328, 551–554. doi: 10.1038/328551a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, P., and Ponder, J. (2003). Polarizable atomic multipole water model for molecular mechanics simulation. J. Phys. Chem. B 107, 5933–5947. doi: 10.1021/jp027815+

CrossRef Full Text | Google Scholar

Ryckaert, J. P., Ciccotti, G., and Berendsen, H. J. (1977). Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341. doi: 10.1016/0021-9991(77)90098-5

CrossRef Full Text | Google Scholar

Sanner, M. F., Olson, A. J., and Spehner, J. C. (1996). Reduced surface: an efficient way to compute molecular surfaces. Biopolymers 38, 305–320. doi: 10.1002/(SICI)1097-0282(199603)38:3<305::AID-BIP4<3.0.CO;2-Y

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J. N., Qiu, L. Q., and Zhang, J. Z. H. (2018). An efficient method for computing excess free energy of liquid. Sci. China Chem. 61, 135–140. doi: 10.1007/s11426-017-9106-3

CrossRef Full Text | Google Scholar

Soriano, V., Gomes, P., Heneine, W., Holguín, A., Doruana, M., Antunes, R., et al. (2000). Human immunodeficiency virus type 2 (HIV-2) in Portugal: clinical spectrum, circulating subtypes, virus isolation, and plasma viral load. J. Med. Virol. 61, 111–116. doi: 10.1002/(SICI)1096-9071(200005)61:1<111::AID-JMV18>3.0.CO;2-W

PubMed Abstract | CrossRef Full Text | Google Scholar

Srinivasan, J., Cheatham, T. E., Cieplak, P., Kollman, P. A., and Case, D. A. (1998). Continuum solvent studies of the stability of DNA, RNA, and phosphoramidate-DNA helices. J. Am. Chem. Soc. 120, 9401–9409. doi: 10.1021/ja981844+

CrossRef Full Text | Google Scholar

Straatsma, T., and Berendsen, H. (1988). Free energy of ionic hydration: analysis of a thermodynamic integration technique to evaluate free energy differences by molecular dynamics simulations. J. Chem. Phys. 89, 5876–5886. doi: 10.1063/1.455539

CrossRef Full Text | Google Scholar

Sun, H., Li, Y., Tian, S., Xu, L., and Hou, T. (2014). Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 16, 16719–16729. doi: 10.1039/C4CP01388C

PubMed Abstract | CrossRef Full Text | Google Scholar

Tannor, D. J., Marten, B., Murphy, R., Friesner, R. A., Sitkoff, D., Nicholls, A., et al. (1994). Accurate first principles calculation of molecular charge distributions and solvation energies from ab initio quantum mechanics and continuum dielectric theory. J. Am. Chem. Soc. 116, 11875–11882. doi: 10.1021/ja00105a030

CrossRef Full Text | Google Scholar

Tie, Y., Wang, Y. F., Boross, P. I., Chiu, T. Y., Ghosh, A. K., Tozser, J., et al. (2012). Critical differences in HIV-1 and HIV-2 protease specificity for clinical inhibitors. Protein Sci. 21, 339–350. doi: 10.1002/pro.2019

PubMed Abstract | CrossRef Full Text | Google Scholar

Tözsér, J., Bláha, I., Copeland, T. D., Wondrak, E. M., and Oroszlan, S. (1991). Comparison of the HIV-1 and HIV-2 proteinases using oligopeptide substrates representing cleavage sites in Gag and Gag-Pol polyproteins. FEBS Lett. 281, 77–80. doi: 10.1016/0014-5793(91)80362-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Triki, D., Billot, T., Visseaux, B., Descamps, D., Flatters, D., Camproux, A. C., et al. (2018). Exploration of the effect of sequence variations located inside the binding pocket of HIV-1 and HIV-2 proteases. Sci. Rep. 8:5789. doi: 10.1038/s41598-018-24124-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Vukovic, S., Brennan, P. E., and Huggins, D. J. (2016). Exploring the role of water in molecular recognition: predicting protein ligandability using a combinatorial search of surface hydration sites. J. Phys. Condens. Matter 28:344007. doi: 10.1088/0953-8984/28/34/344007

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Nguyen, P. H., Pham, K., Huynh, D., Le, T. B., Wang, H., et al. (2016). Calculating protein–ligand binding affinities with MMPBSA: method and error analysis. J. Comput. Chem. 37, 2436–2446. doi: 10.1002/jcc.24467

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Morin, P., Wang, W., and Kollman, P. A. (2001). Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA. J. Am. Chem. Soc. 123, 5221–5230. doi: 10.1021/ja003834q

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., and Kollman, P. A. (2000). Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model1. J. Mol. Biol. 303, 567–582. doi: 10.1006/jmbi.2000.4057

CrossRef Full Text | Google Scholar

Wang, W., Lim, W. A., Jakalian, A., Wang, J., Wang, J., Luo, R., et al. (2001). An analysis of the interactions between the Sem- 5 SH3 domain and its ligands using molecular dynamics, free energy calculations, and sequence analysis. J. Am. Chem. Soc. 123, 3986–3994. doi: 10.1021/ja003164o

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Liu, J. F., Zhang, L. J., He, X., and Zhang, J. Z. H. (2017). Computational search for aflatoxin binding proteins. Chem. Phys. Lett. 685, 1–8. doi: 10.1016/j.cplett.2017.07.024

CrossRef Full Text | Google Scholar

Warshel, A., Kato, M., and Pisliakov, A. V. (2007). Polarizable force fields: history, test cases, and prospects. J. Chem. Theory Comput. 3, 2034–2045. doi: 10.1021/ct700127w

PubMed Abstract | CrossRef Full Text | Google Scholar

Wlodawer, A., and Vondrasek, J. (1998). Inhibitors of HIV-1 protease: a major success of structure-assisted drug design. Annu. Rev. Bioph. Biom. 27, 249–284. doi: 10.1146/annurev.biophys.27.1.249

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, B., Shen, H., Zhu, X., and Li, G. (2011). Fast and Accurate computation schemes for evaluating vibrational entropy of proteins. J. Comput. Chem. 32, 3188–3193. doi: 10.1002/jcc.21900

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Sun, H., Li, Y., Wang, J., and Hou, T. (2013). Assessing the performance of MM/PBSA and MM/GBSA methods. 3. The impact of force fields and ligand charge models. J. Phys. Chem. B 117, 8408–8421. doi: 10.1021/jp404160y

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Y., Yang, M., Ji, C. G., and Zhang, J. Z. H. (2017). Interaction entropy for computational alanine scanning. J. Chem. Inf. Model. 57, 1112–1122. doi: 10.1021/acs.jcim.6b00734

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, D. W., and Zhang, J. Z. H. (2003). Molecular fractionation with conjugate caps for full quantum mechanical calculation of protein–molecule interaction energy. J. Chem. Phys. 119, 3599–3605. doi: 10.1063/1.1591727

CrossRef Full Text | Google Scholar