N88S and L90M mutations in CRF01_AE HIV-1 protease: Molecular dynamics simulations and binding free energy calculations

Accepted Manuscript Title: Drug Resistance Mechanism of L10F, L10F/N88S and L90M mutations in CRF01 AE HIV-1 protease: Molecular dynamics simulations and binding free energy calculations Authors: C.S. Vasavi, Ramasamy Tamizhselvi, Punnagai Munusami PII: DOI: Reference:

S1093-3263(17)30101-8 http://dx.doi.org/doi:10.1016/j.jmgm.2017.06.007 JMG 6938

To appear in:

Journal of Molecular Graphics and Modelling

Received date: Revised date: Accepted date:

13-2-2017 2-6-2017 4-6-2017

Please cite this article as: C.S.Vasavi, Ramasamy Tamizhselvi, Punnagai Munusami, Drug Resistance Mechanism of L10F, L10F/N88S and L90M mutations in CRF01 AE HIV-1 protease: Molecular dynamics simulations and binding free energy calculations, Journal of Molecular Graphics and Modellinghttp://dx.doi.org/10.1016/j.jmgm.2017.06.007 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Drug Resistance Mechanism of L10F, L10F/N88S and L90M mutations in CRF01_AE HIV-1 protease: Molecular dynamics simulations and binding free energy calculations

C.S. Vasavia, Ramasamy Tamizhselvia , Punnagai Munusamib*,

a

School of Biosciences and Technology, VIT University, Vellore 632014, Tamilnadu, India

b

Center for Computational Natural Sciences and Bioinformatics, International Institute of Information Technology, Hyderabad 500 032, India

*Corresponding author. E-mail address:[email protected]

Graphical abstract

Highlights  

 

The molecular basis of nelfinavir resistance due to strongly linked L10F/N88S mutation and L90M mutation in AE-protease was analyzed. The studies reveal that due to spatial deviation of nelfinavir, a hydrogen bond that is crucial for effective inhibition is lost in L10F and L10F/N88S complexes. When L10F co-occurs with N88S, an intramolecular hydrogen bond network is altered thus impairing the interaction of nelfinavir. In L90M mutant, the hydrogen bond interactions between nelfinavir and the active site are maintained similar to subtype AE wildtype and it does not affect the binding affinity.

Abstract HIV-1 protease plays a crucial role in viral replication and maturation, which makes it one of the most attractive targets for anti-retroviral therapy. The majority of HIV infections in developing countries are due to non-B subtype. Subtype AE is spreading rapidly and infecting huge population worldwide. The mutations in the active site of subtype AE directly impair the interactions with the inhibitor. The non-active site mutations influence the binding of the inhibitor indirectly and their resistance mechanism is not well understood. It is important to design new effective inhibitors that combat drug resistance in subtype AEprotease. In this work, we examined the effect of non active site mutations L10F, L10F/N88S and L90M with nelfinavir using molecular dynamics simulation and binding free energy calculations. The simulations suggested that the L10F and L10F/N88S mutants decrease the binding affinity of nelfinavir, whereas the L90M mutant increases the binding affinity. The formation of hydrogen bonds between nelfinavir and Asp30 is crucial for effective binding. The benzamide moiety of nelfinavir shows large positional deviation in L10F and L10F/N88S complexes and the L10F/N88S mutation changes the hydrogen bond between the side chain atoms of 30th residue and the 88th residue. Consequently the hydrogen bond interaction between Asp30 and nelfinavir are destroyed leading to drug resistance. Our present study shed light on the resistance mechanism of the strongly linked mutation L10F/N88S observed experimentally in AE subtype.

Keywords: HIV-1 protease; Nelfinavir; Subtype AE; Non-active site mutations; Resistance mechanism

Introduction The spread of Human Immuno defeciency virus (HIV) has proven to be a major clinical challenge and is responsible for Acquired immuno defeciency syndrome (AIDS). It has been more than three decades after its recognition, AIDS still has no cure, but treatment with combination of antiretroviral drugs has improved the health of people with longer life expectancy [1, 2]. HIV is subdivided into two major types HIV-type 1 (HIV-1) and HIV-type

2 (HIV-2) [3]. HIV-1 is responsible for majority of the infection and is spread worldwide. HIV-2 is largely restricted to countries like West Africa. HIV-1 is categorized into different groups, subtypes, sub-subtypes, circulating recombinant forms (CRFs) and unique recombinant forms. Each subtype (A–D, F–H, and J–K), exhibits a unique set of naturally occurring polymorphism [4]. The life cycle of HIV presents numerous potential targets for anti retroviral therapy. The promising targets include distinct stages starting from viral entry, reverse transcription, integration and maturation [5]. HIV-1 protease is responsible for cleaving gag and gag-pol polyproteins into structural proteins which are necessary for viral maturation. Since inhibition of this enzyme prevents the maturation of viral proteins and the replication of the virus, HIV-1 protease has been considered as an important target in the treatment of HIV infection and AIDS [6, 7]. HIV-1 protease is a homodimeric aspartic protease. Each subunit is made up of 99 amino acids. The active site with the catalytic triads (Asp25/25΄-Gly27/27΄) is interfaced between subunits and is covered by two flexible β-hairpin flaps (residues 43-58/43΄-58΄). These flaps are important for binding of the substrate and product release [8]. The movement of these flaps is essential for activity and specificity of the protease [9]. Protease inhibitors (PIs) have emerged as the most effective drugs in the treatment of HIV. PIs are competitive active site inhibitors that prevent HIV-1 protease from cleaving the polyproteins into mature infectious virions. Currently there are nine FDA approved PIs available for the treatment of HIV: saquinavir (SQV), indinavir(IDV), ritonavir(RTV), nelfinavir(NFV), amprenavir(APV), lopinavir(LPV), atazanavir(ATV), tipranavir(TPV) and darunavir(DRV) [10-12].

The

currently available protease inhibitors (PIs) were primarily developed based on subtype B isolates, although they account for approximately 11% of HIV-1 infection world-wide. Though the therapeutic intervention with these PIs has greatly improved the survival of patients infected with HIV, it has led to rapid emergence of drug resistance [13]. The infection of circulatory recombinant form 01 A/E (AE) is expanding rapidly and infecting large proportion of people throughout the world, more predominantly in south asian countries [14]. The consensus AE sequence has polymorphisms at positions I13V, E35D, M36I, S37N, R41K, H69K and L89M and they differ by ~10% from subtype B protease. Some of these naturally occurring polymorphisms are associated with drug resistance in subtype B protease [15]. It has been shown that when the drug resistance mutations occur along with natural polymorphisms, the binding affinity of the PIs are dramatically decreased [16]. NFV has been widely used for several years for the treatment of HIV. The molecular

mechanism of NFV resistance caused by active-site and non-active site mutations has been well studied in subtype B protease. There are two possible pathways for acquisition of mutations related to NFV resistance. Resistance to NFV most often arises with the selection of active site mutation D30N in subtype B protease [17, 18]. The other possible pathway leading to NFV resistance involves the non-active site mutation L90M [19]. The N88S is a rarely selected during NFV and IDV therapy in subtype B protease. Instead the N88D mutation is often observed in patients who fail NFV therapy. In addition, a combination of N88D and D30N mutations were seen in subtype B and C patients treated with NFV [20]. Molecular dynamics (MD) simulation was performed by Ode et al. to understand the resistance mechanism of NFV due to D30N, N88D, D30N/N88D and L90M in subtype B protease. They reported that N88D and D30N/N88D mutants modulated the hydrogen bond network and thus leading to resistance to NFV in subtype B [21]. Computational simulations of L90M with NFV by same group showed that L90M mutation caused large conformational change at the active site, flap and 80's loop region of the protease, thus leading to displacement of the NFV from its active site [21, 22]. NFV resistance pathway of AEprotease differs from that of subtype B. Ariyoshi et al. reported the pattern of drug resistance mutations that differ between subtype B and AE. They observed that the mutations L10F, K20I, L33I and N88S were seen more frequently in patients infected with AE-protease. They reported that in AE-protease, the mutations L10F, L33I and N88S were significantly associated with NFV resistance, whereas the K20I mutation was not associated with the use of any particular PIs. Interestingly, they observed that the N88S mutation was strongly linked with L10F mutation in AE patients, whereas in subtype B patients the N88D mutation was strongly associated with D30N mutation. They also observed that none of the subtype B patients had both L10F and N88S mutations. Conversely, none of the AE patients had D30N, A71V, or N88D mutations, which were often found in PI-resistant subtype B patients [2324]. Bandaranayake et al. performed biochemical and biophysical methods to determine the effect of polymorphisms in AE-protease by measuring enzyme activity and inhibitor binding [25]. They reported that inherent weaker affinity due to natural polymorphism in AE and the single non-active site mutation at position 88 lead to NFV resistance. The mechanism of resistance against NFV due to N88S mutation in subtype AE-protease was investigated by Ode et al. through MD simulation. They reported that the N88S mutation led to conformational changes at D30 reducing the interaction between D30 and NFV [26]. The L90M mutation which is specific to NFV resistance is also reported in subtype AE and other subtypes C, F and G [27,28]. Yam et al, performed genotypic resistance test to determine

mutation patterns in subtypes AE and B patients receiving antiretroviral therapy. They observed that L90M mutation was commonly occuring in both subtype B and AE patients who had been exposed to treatment with NFV [29]. From experimental and theoretical studies it is apparent that the combined effects of naturally occurring polymorphisms and drug resistant mutations affect the viability of currently available PIs. Although studies have been reported in detail about the drug resistance mechanism of L90M mutation in subtype B, N88D/D30N in subtype B and C protease [30-32], the resistance mechanism of strongly linked mutation L10F/N88S with NFV is still unknown. In the present work, we carried out MD simulations and binding free energy calculations to investigate the resistance mechanism against NFV due to L10F, L10F/N88S and L90M mutations (Fig.1). To quantitatively describe the influence of these non-active site mutations, the binding free energies were calculated using the MM-PBSA methods. The free energies were further decomposed on a per-residue basis. Our computational simulations indicate the co-occurrence of L10F with N88S mutation enhances the resistance mechanism against NFV. It was interesting to note that L90M shows increase in binding free energy compared to AE-WT. The large spatial deviation of benzamide moiety of NFV and the loss in hydrogen bond interactions between Asp30 and NFV lead to drug resistance in L10F and L10F/N88S complexes. The findings from this study would hint in designing more effective PIs against subtype AE-protease.

Materials and Methods Homology modeling Homology modeling can act as a powerful tool to model the three-dimensional structure of the proteins if the target sequence and the template structures show high similarity. The three-dimensional structure of a given protein sequence is predicted based on its sequence similarity to one or more experimentally resolved protein structures [33-34]. In several studies, researchers have modeled the three-dimensional structures of HIV-1 protease subtypes to understand the drug resistance mechanism conferred by mutations [21, 22, 26, 30, 35-37]. Currently, there are no experimentally determined structures available for AEprotease complexed with NFV. The sequences of AE-protease exhibited > 90% identity with the template structure of HIV-1 subtype B protease complexed with NFV (PDB id: 1OHR) [38]. In the present study the wild type AE-WT was constructed by substituting the following amino acids I13V, E35D, M36I, S37N, R41K, H69K, and L89M in subtype B protease

sequence. The mutations L10F, L10F/N88S and L90M were introduced in AE-WT to create respective mutant structures. The sequence alignment of subtype B, AE-WT and mutant structures L10F, L10F/N88S, L90M generated using ESPript (http://espript.ibcp.fr) [39] are shown in Fig.2.

Each initial structure for the AE-protease in complex with NFV was modeled from the atomic coordinates of the template structure along with its crystallographic water molecules. The structures were modeled using Modeller 9.17 [40]. The generated models were further evaluated energetically based on its Discrete Optimized Protein Energy (DOPE) score and the stereo-chemical quality of the modeled structures was validated using RAMPAGE server [41]. The Ramachandran plots were generated using Discovery Studio 3.5 Visualizer [42]. The models with the lowest energy DOPE score, validated by RAMPAGE were considered for energy minimization and molecular dynamic simulations. System setup and Molecular Dynamic Simulations Energy minimizations and MD simulations were carried out using the Sander module of AMBER 12 [43, 44]. The force field parameters for the protease residues were adopted from the AMBER database (ff03). The force field parameters for NFV were obtained from General AMBER Force Field (GAFF) [45]. The atomic partial charges of NFV were obtained using following procedure. The atomic coordinates of NFV was extracted from the crystal structure (PDB id: 1OHR) and the missing hydrogen atoms were added by using the LEap module of AMBER12 [46]. The resulted molecular geometry was then optimized using (HF)/6-31G* basis set of Gaussian 09 package [47] to obtain the electrostatic potential. The Antechamber module [48] of the AMBER 12 package was used to obtain RESP charges from the electrostatic potential [49]. The explicit solvent model was chosen for all MD simulations. Each structure was solvated in a cubic box filled with TIP3P water molecules to allow at least 10 Å of solvent on each face of the AE-protease. The positively charged system was neutralized by adding chloride ions to form a charge-neutral system. Energy minimizations were performed to allow the system to reach an energetically favorable conformation. First, the restraints were placed on all heavy atoms of the AE-proteases with a harmonic force constant of 10 kcal mol−1 Å−2. Then the restraints were applied only to the backbone nitrogen, oxygen, and carbon atoms. The strength of the restraint was maintained as 10 kcal mol−1 Å−2. Finally, the restraints were turned off, and all atoms were allowed to relax.

In each step, energy minimization was executed by the steepest descent method for the first 2500 steps and the conjugate gradient method for the subsequent 2500 steps. After energy minimization, the temperature of each system was gradually raised from 0 K to 300 K under constant volume (NVT), over a period of 100 ps. Following heating, 100 ps of constantpressure equilibration was carried out at 1 atm. A harmonic force constant of 10 kcal mol−1 Å−2 was applied to restrain the solute heavy atoms during heating and constant-pressure equilibration processes. The force constant on the restrained atoms was then reduced in the following order 8, 4, and 2 kcal mol−1 Å−2, the equilibration had been repeated for 50 ps at each value. Then unrestrained constant-pressure equilibration was extended for 500 ps. Following the equilibration, MD simulations were performed in the NPT ensemble at a temperature of 300K by using Langevin dynamics temperature scaling with collision frequency 2 ps-1. The coordinates of the trajectories were saved at every 10 ps for each simulation. Adequate cutoff sizes are necessary to obtain reliable MD results; hence a cutoff of 12.0 Å was used for MD simulations. The particle mesh Ewald (PME) method was used to handle long range electrostatic interactions. The SHAKE algorithm [50] was employed to constrain all bond lengths involving hydrogen atoms. The simulation time step was set to 2 fs. The MD simulations were carried out using periodic boundary conditions (PBC). The trajectories obtained from the last 5 ns were used to analyze the structure in detail. The plots were constructed using the Xmgrace program and VMD program was used for analyzing the trajectories. Protonation State Determination The protonation states of Asp25/Asp25΄ play a key role in catalytic mechanism and inhibitor recognition. Protonation state of these two aspartic groups (Asp25/Asp25΄) vary according to different PIs and they remarkably affect the mechanism of binding of the protease-inhibitor complex. Therefore, a proper determination of the appropriate protonation state of Asp25/Asp25΄ is important. Several studies have reported that the mono-protonated states favor the ligand binding [51, 52]. Two different protonation states of Asp residues have been considered. (i) The aspartate in the A chain protonated and B chain unprotonated. (ii) The aspartate B chain is protonated and chain A is unprotonated. In order to determine the protonation states of L10F, L10F/N88S and L90M complexed with NFV, the binding energies of mono-protonated states of each complex were computed and compared from the trajectories obtained during 4 ns MD simulation. The binding energies were calculated using MM/PBSA method [53].

Binding Free Energy Calculation The binding free energy between the AE-protease and NFV was calculated using the MM/PBSA methodology according to the equation: Gbind = Gcomplex – (Greceptor + Gligand)

(1)

Where Gcomplex, Greceptor and Gligand in Eq.1 denotes the free energies of complex, receptor and ligand averaged over the entire MD trajectories. In MM/PBSA, the binding free energy of a system is expressed as Gbind =H –TS

(2)

H is enthalpy and –TS in Eq.2 is the contribution of entropy to binding. The enthalpy component (H) in Eq.3 is further decomposed into sum of electrostatic (Eele), van der waals interaction (Evdw) and solvation free energy (Gsol). H = Eele + Evdw + Gsol

(3)

Gsol =Gpol +Gnp

(4)

Gnp = SASA + β

(5)

The electrostatic and van der Waals components between the protein and ligand are calculated using molecular mechanics. The solvation free energy contribution includes the polar (Gpol) and non-polar (Gnp) contributions. The polar portion (Gpol) in Eq.4 was calculated using PBSA module of AMBER12 program. The dielectric constant for the solute and solvent was set to 1.0 and 80.0, respectively. The nonpolar contribution (Gnp) in Eq.4 ) to the solvation free energy was estimated from linear relation to solvent accessible surface area (SASA) using the LCPO method as described in Eq.5, in which the surface tension  and offset β were set with default values of 0.00542 kcal mol−1 Å−2 and 0.92 kcal mol−1 respectively. Per-residue decomposition To examine the energetic contribution of individual residues, the interaction between NFV and each AE-protease residue were computed using the MM-GBSA decomposition method in AMBER12. The binding interaction of each NFV-residue pair consists of three terms: Electrostatic (Eele), van der Waals interaction (Evdw) contribution in gas phase and polar solvation contribution (GPB) and are computed as follows: GNFV_residue = Evdw + Eele + GPB (6)

The polar solvation contribution was calculated using the Generalized Born model and the parameters were developed by Onufreiv and Bashford [54]. All energy components of Eq.6 were calculated from the last 5 ns of MD simulation. Hydrogen Bond Criterion The PTRAJ module of AMBER was used for analyzing the hydrogen bonds. The formation of hydrogen bonds was defined in terms of distance and angle cutoff as follows: (a) the distance between donor D and acceptor A was ≤ 3.5Å (b) the angle between D-hydrogenA was ≥120˚. Results and discussion Homology modeling of AE-WT and mutant structures The three-dimensional structures of AE-protease were modeled as a step to understand the resistance mechanism against NFV.

Ramachandran plot represents the

distribution of phi and psi angles for each residue in the protein. The results of the Ramachandran plot for the modeled structures indicated that ~98% of the residues were present in the favored region for AE-WT, L10F/N88S, L90M protease models and ~97% of residues for L10F. For all the modeled structures ~2% of the residues were located in allowed regions and no residues were restrained to the disallowed region. The Ramachandran plots for the modeled structures are shown in Fig.S1. Ideally a good quality model is preferred to have more than 90% of the residues in the favored regions of the Ramachandran plot [55]. Hence the quality of the modeled structures is valid and could be further used to understand the resistance mechanism against NFV through MD simulation and binding free energy calculations. Comparing the Stability and flexibility of complexes and structures from MD trajectories: AE-WT vs. Mutant To explore the effect of mutations on the conformational stability of the AE-WT, the root mean square displacement (RMSD) values for the Cα coordinates relative to the initial structures were calculated and shown in Fig.3. The RMSD of each mutant structure with respect to AE-WT is given in Fig.S2. The RMSD plots indicate that the conformations of AE-WT, L10F, L10F/N88S and L90M complexes were in good equilibrium. The trajectories of the L10F mutant complex fluctuate more than the other three complexes with values around 1.70 to 2.00 Å from 2 ns till the end of simulation period. The L10F/N88S and AE-

WT complexes go parallel to each other till 7.5 ns, after 7.5 ns a slight deviation was observed in L10F/N88S complex reaching a maximum of 1.7 Å. In L90M complex, the structure started deviating to a greater extent after 9 ns, with values fluctuating around 1.6 to 2.00 Å. According to the RMSD average values of the four complexes, L10F has a higher mean (1.54 Å) than the L90M (1.45 Å), L10F/N88S (1.35 Å) and AE-WT (1.20 Å) with respective standard deviations of 0.22, 0.17, 0.15 and 0.13 Å. The result signifies that all the four complexes deviate to a quite similar extent from their starting structures with values around 1.0 to 1.7 Å during the course of simulation, ensuring stable trajectories. This suggests that the stabilities of the four complexes were reliable and could provide a suitable basis for the subsequent analyses.

The detailed analyses of root-mean-square fluctuation (RMSF) of backbone atoms versus residue number for all complexes are illustrated in Fig.4. The RMSF plots of each mutant structure with respect to AE-WT are given in Fig.S3. As shown in Fig.4, although the RMSF shape of the L10F, L10F/N88S and L90M mutants were similar to AE-WT, the mutated complexes have slightly larger RMSF values than that of AE-WT complex. The catalytic Asp25 stays relatively rigid in both the monomers except L10F/N88S complex. The active site residues 24΄-29΄ of L10F/N88S complex showed higher flexibility than other complexes. Though in all the three mutant complexes fluctuations were observed in the flaps, the L90M complex shows high flexibility in the flap region. The amino acids Thr91-Gly94 showed additional mobility in L10F complex compared to AE-WT. The L10F/N88S mutant resulted in higher fluctuations in the region 61΄-65΄ and 85΄-88΄. Differences in RMSF values of the AE-WT and mutant complexes are shown in Fig.5. The residues Pro1 (N-terminal), 14-16/16΄ (fulcurum), Gly40-Trp42 (Flap elbow), Asp25΄ (catalytic triad), Leu38΄(Flap elbow), Lys69΄, Lys70΄(cantilever), Asn88΄ (alpha-helix) and Phe99΄ (C-terminal) show absolute difference of larger than 0.4Å and were considered as the highly fluctuating residues.

Flap tip to catalytic Asp (Ile50−Asp25) distances The extent of flap opening was measured by the distance between the flap tip (Ile50Cα/ Ile50΄ Cα) and the catalytic aspartate (Asp25Cα/Asp25΄Cα). The measure of this distance was believed to be more reasonable than I50-I50΄ distance because both flap tip curling and flap symmetry can affect the flap tip-tip distances. In Ile50-Asp25 and Ile50΄-Asp25΄ distances for

mutant complexes (L10F, L10F/N88S and L90M) were measured from the MD trajectories the histogram distributions are shown in Fig.6. It was found that Ile50΄-Asp25΄ distribution (chain B) was much overlapped compared with to Ile50-Asp25 distribution (chain A). The mean distances of flap tip-catalytic aspartate in chain A were 15.08, 14.89, 15.77 and 13.76 with SD values of 0.41, 0.44, 0.63 and 0.42 for AE-WT, L10F, L10F/N88S and L90M complexes respectively. This indicates that in L10F/N88S complex the distance between catalytic aspartate and flap tips was stretched. The mean o distribution for L90M was smaller by ~ 1.3Å, whereas the L10F complex is smaller by ~ 0.2Å from AE-WT complex. For chain B, the mean distances and SD values of AE-WT (L10F) are 14.63 (14.52) Å and 0.70Å (0.61), respectively; whereas for the L10F/N88S complex the mean and SD are 13.94 and 0.51 respectively. The respective mean distance and SD of L90M complex was 13.72 and 0.51. The mean of distribution of L10F/N88S was smaller by ~0.7Å than AE-WT complex and ~0.6Å from L10F complex. In general for all complexes (AE-WT, L10F, L10F/N88S and L90M) the average flap tip-catalytic aspartate distance in chain A is longer than that of chain B. In L10F/N88S complex, the distances from flap regions to the catalytic aspartates are farther than the AE-WT. However the flap tip-catalytic aspartate distance of L10F is slightly less, whereas in L90M complexes the distance is smaller in both the monomers by ~1.0Å compared to AE-WT. This change in distances possibly affects or enhances the binding of NFV. To further clarify this, the binding energy of the complexes were calculated.

Total binding free energies In order to gain insight into individual contributions to total binding free energy, the electrostatic interaction, van der Waals interaction, solvation energy and entropy were calculated for the AE-WT and mutant (L10F, L10F/N88S and L90M) complexes using molecular mechanics generalized Born surface area method in AMBER12 package. MMPBSA method is a versatile tool used in the calculation of binding free energy of a given protein-inhibitor complex. A lot of studies have reported that MM-GBSA approach provides acceptable correlation with experimental values even when carried for lower nanosecond simulation [56]. A summary of the contributions of the binding interactions between NFV and all the complexes are given in Fig.7 and Table 1. The calculated binding energy of AEWT was consistent with the binding free energy reported by Ode et al., [26]. As shown in Table 1, the calculated binding free energies of the AE-WT, L10F, L10F/N88S and L90M

complexes were −39.8, −35.5, −34.4 and −43.2 kcal/mol respectively, suggesting weaker binding affinity of NFV to mutants L10F and L10F/N88S, whereas the single mutant L90M exhibited increase in binding affinity to NFV compared to AE-WT. The most prominent contributions were from van der Waals (Evdw) interactions, closely followed by the electrostatic (Eele) contributions to the total binding free energy. Non-polar solvation energies (Gnp), resulting from the burial of solvent accessible surface area upon NFV binding also contributed to the binding energy a bit favourably. The polar solvation energies and entropy components produced unfavourable components to the binding free energy. The L10F and L10F/N88S complexes showed favourable van der Waals interaction, but there is decrease in electrostatic contributions compared to AE-WT. The L10F showed decrease in electrostatic contribution. When N88S mutation occurs in combination with L10F (L10F/N88S), it resulted in dramatic decrease in electrostatic contribution relative to AE-WT. In L10F, the loss in binding energy is solely attributed to the change in entropy compared to enthalpy components. Furthermore it was observed that the entropy for L10F vs AE-WT was unfavourable. In L10F/N88S the loss in binding energy was mainly due to the impact of enthalpy components relative to the AE-WT. The L90M mutant resulted in binding enhancement in terms of electrostatic energy. The entropic component of L90M complex showed similar contribution as AE-WT. Contributions from individual residues to the binding free energy To understand the drug resistance mechanism, the binding energies of AE-WT and mutant complexes were decomposed into inhibitor-residue pairs to create an inhibitor-residue spectrum as shown in Fig.8. The per-residue contributions were calculated using MM-GBSA method. The energy contributions from backbone and side chain atoms of individual residues are given in Tables 2-5. This per-residue decomposition method is especially useful to explain the drug resistance mechanism at atomistic level and also helpful to locate the individual residue contribution to protein-inhibitor interactions. From Fig.8, it was clear that the interaction spectra were little different for all the four complexes. Overall the major contributions come from a few groups around Asp25/Asp25΄, Ala28/Ala28΄, Ile50/Ile50΄ and Ile84/Ile84΄. The energy decomposition spectra showed that the binding mode of NFV to chain A is slightly different from that of chain B, this is due to asymmetrical nature of NFV. The contributions of most residues to binding affinities were mainly due to van der Waals and electrostatic energy components. Most of the residues showed unfavourable polar solvational energy contribution towards binding affinity. The amino acids Asp25/Asp25΄,

Gly27, Ala28/Ala28΄, Asp30, Ile47, Gly49/Gly49΄, Ile50/Ile50΄, Pro81΄, Val82΄and Ile84/Ile84΄ showed binding energy of more than −1.0 kcal/mol in AE-WT.

To complement the energetic analysis, hydrogen bond analysis was performed on the trajectories obtained from MD simulation. Hydrogen bonds are ubiquitous in nature and play an essential role in many biological and chemical processes that includes protein folding, protein-ligand interaction and enzyme catalysis [57, 58]. The drug resistance mechanisms of HIV-1 protease with several PIs are mainly due to disruption of hydrogen bond interactions. Hydrogen bond interactions between the AE-WT and mutant complexes are given in Table 6-9 and Table S1-S4. b

Occupancy percentage to evaluate the stability and strength of the hydrogen bonds.

The catalytic aspartates (Asp25/Asp25΄) play a crucial role in binding of the inhibitors to HIV-1 protease. The interaction between Asp25΄ and NFV in AE-WT was driven mainly by electrostatic energy, but it showed unfavourable van der Waals energy contribution towards binding affinity. As shown in Table 6, the catalytic Asp25΄ and central hydroxyl group (P1 subsite) of NFV formed hydrogen bond in AE-WT with occupancy of 99.8%. Similarly the L10F and L10F/N88S complexes were also dominated by the electrostatic contributions from catalytic Asp25/Asp25΄, but the contribution of electrostatic energy from Asp25΄ in L10F and L10F/N88S complexes were reduced by 5.69 and 9.43 kcal/mol respectively compared to AE-WT. A favourable contribution of polar solvation energy from Asp25΄ (−0.19 kcal/mol) was observed in L10F/N88S mutant, whereas other complexes showed unfavourable contribution of polar solvation energy. The results in Tables 7 and 8 showed that the occupancy of Asp25 and NFV hydrogen bond in L10F was 99.8%, whereas in L10F/N88S mutant, the Asp25(OD1/OD2) occupancies were reduced to ~60%. Another important residue in the active site is Asp30 which was involved in hydrogen bond formation with various PIs and they play a vital role in binding of NFV. The residue based energy decomposition of AE-WT clearly highlighted that Asp30 as a major contributor towards binding energy in terms of electrostatic interactions. The electrostatic energy in AEWT, originates from the hydrogen bond interactions between m-phenol group (P2 subsite) of NFV and carboxyl oxygen of Asp30 with occupancy of 95%. Interestingly, no hydrogen bonds were observed between Asp30 and NFV for the mutants L10F and L10F/N88S which further led to decrease in electrostatic energy contribution by ~6.0kcal/mol. But there was

increase in van der Waals energy contribution by −0.49 and −0.71 kcal/mol respectively compared to the AE-WT. The total energy contribution of residues Asp25, Asp25΄, and Asp30 were increased in L90M approximately by −0.30, −0.70 and −0.50 kcal/mol respectively. Similar to AE-WT, the electrostatic interaction were the key driving force for binding of residues Asp25 and Asp30 in L90M. The central hydroxyl group (P1 subsite) of NFV forms hydrogen bond with both the carboxyl oxygens OD1 and OD2 of Asp25 with occupancies of ~30% and 99.8% respectively. However the occupancy of hydrogen bond between Asp30 and m-phenol group in NFV was reduced to 61% in L90M. Interestingly, the OD2 atom of Asp30 and O4 atom of NFV also showed weak hydrogen bond formation in L90M (Fig.9). The van der Waals energy between Gly27 and NFV directs the binding in AE-WT and mutant complexes and the contributions are mainly from the side-chain atoms of Gly27. In L10F, the P1 subsite lies distant to Gly27 thus showing less contribution to van der Waals energy, whereas in other complexes they tend to lie in close proximity.

Like hydrophobic residue Gly27, the binding between NFV and Ala28/Ala28΄, was governed by van der Waals interactions. Both backbone and side-chain of Ala28 showed almost equal contribution towards van der Waals interaction in AE-WT, whereas in mutant complexes only the backbone atoms were involved. In Ala28΄ the van der Waals interactions were mainly derived from the backbone atoms for all the complexes and the electrostatic contributions were observed to be unfavourable towards binding affinity. In AE-WT and mutant complexes, Ala28 was oriented towards P2 subsite of NFV favouring van der Waals interaction compared to electrostatic interactions. Due to substitution of phenylalanine in the place of leucine in 10th position, a slight shift in the alkyl group of Ala28 was observed. This led to the increase in distance between Ala28 and P2 subsite of NFV thus reducing the side chain van der Waals contribution by 0.79kcal/mol. Structural analyses of the L10F/N88S and L90M showed that the P2 subsite of NFV is oriented in a similar manner as in AE-WT favouring van der Waals interaction with Ala28. The binding affinity of Ala28΄ with NFV was increased by −0.21 and −0.25 kcal/mol respectively for mutants L10F and L90M. This may be due to the close proximity of Ala28΄ and P2 subsite of NFV, thus resulting in increased van der Waals interaction than AE-WT. The L10F/N88S mutant almost showed similar binding affinity as AE-WT. For residue, Ile47, though the main driving force is van der Waals energy, there was also favourable

contribution from electrostatic interactions. The mutants L10F and L90M showed more contribution to binding affinity compared to AE-WT complex with differences of ~−0.6 and −0.3 kcal/mol respectively. This may be due to conformational change in side chain of Ile47 making Ile47 align in close proximity to P2 subsite. In L10F/N88S complex the P2 subsite is oriented in a similar way as AE-WT showing same van der Waals contribution towards binding affinity. The residue Gly49 showed less contribution towards binding affinity in all mutant complexes compared to AE-WT. A possible reason for this could be partial loss of electrostatic and van der Waals interaction due to increased distance between Gly49 and P1 subsite. The residue Gly49΄ of L10F/N88S showed increase in contribution to binding affinity by −0.27kcal/mol, whereas other mutants showed similar trend as AE-WT. The Ile50/Ile50΄ at the flap tip plays a key role in providing stabilization to the protease dimer. In L10F, the conformational changes of Ile50 makes the backbone atoms align little closer to P2΄ in NFV, resulting in favourable electrostatic contribution of −1.08 kcal/mol, whereas in AE-WT there was unfavourable contribution from the backbone atoms. In L10F/N88S and L90M complexes, the Ile50 showed less van der Waals and electrostatic contribution, thus resulting in decrease in binding affinity. The residue Ile50΄ considerably decreased the binding affinity of L10F and L90M by 0.17 and 0.34 kcal/mol respectively. This is because Ile50΄ is shifted away from P1΄ of NFV, which would probably lead to decrease in binding affinity. The Ile50΄ of L10F/N88S showed similar binding affinity to NFV as in AE-WT. For the residue Ile84/Ile84΄, the van der Waals energy is the key contribution for binding energy in all the complexes. The P1΄ subsite of L90M is moved away from Ile84, resulting in decreased contribution towards binding affinity by 0.52 kcal/mol, whereas other complexes showed binding affinity similar to AE-WT in the range of ~−1.5 kcal/mol. The AE-WT and L90M showed favourable polar solvational energy contributions towards binding affinity. The conformational changes of Ile84΄ in L10F, L10F/N88S and L90M complexes led to decrease in binding affinity by 0.30, 0.16 and 0.34 kcal/mol respectively. The amino acids Pro81΄ and Val82΄ were mainly driven by van der Waals contribution towards binding affinity. The distance between Pro81΄ and P1 subsite of NFV was increased, leading to decrease in binding affinity in L10F, L10F/N88S and L90M complexes by 0.30, 0.38 and 0.13 kcal/mol respectively. In L90M, the binding affinity of Val82΄ was decreased by 0.13kcal/mol, whereas other mutant complexes maintained similar

binding affinity with NFV as observed in AE-WT. Further the binding of Val82΄ with NFV in all complexes were favoured by contributions from polar solvational energies.

Effects of mutations on the binding affinity: Molecular pathways to drug resistance vary for each subtype and are yet to be fully characterized in non-B subtypes. Though lot of studies have been carried out to understand the drug resistance mechanism in HIV-1 protease, very little is known about the resistance mechanism in AE-protease. Also it is difficult to speculate the resistance mechanism caused due to non-active site mutations. In this study, we performed MD simulations on AE-protease to understand how the non-active site mutations L10F, L10F/N88S and L90M confer resistance against NFV. The simulations indicate that each mutation caused change in binding free energy compared to AE-WT. The L10F and L10F/N88S complexes lead to decrease in binding free energies, whereas the binding free energy increased for L90M. In L10F even though the enthalpy is gained, there was decrease in free energy of binding due to unfavourable entropy contribution. In L10F, the contributions from residues Ala28΄, Asp30 and Ile47 favoured the increase in enthalpy. In L10F/N88S, the loss in electrostatic energy is due to less contribution from residues Ile47, Ile50 and Gly49΄. This loss further weakened the enthalpy contribution towards the binding free energy. The amino acids Asp25, Asp30 and Ile47 of L90M mutant showed favourable contributions to enthalpy, thus enhancing the binding of NFV to the active site.

Intramolecular hydrogen bond interactions: To clarify the effect of mutations, the intramolecular hydrogen bond interactions were analysed in AE-WT and mutant protease structures. Our hydrogen bond analysis indicated that the catalytic Asp25 of both the chains form hydrogen bond interactions in AE-WT, but this interaction was lost in mutant complexes. Apart from the catalytic triad present in aspartic proteases, the Gly86-Arg87-Asn88 triad is specific to retroviral protease and is present close to the active site loop and dimer interface but does not form any direct interaction with the inhibitor. It has been reported that a salt bridge between Asp29-Arg87 is crucial for stabilizing the dimer interface of the protease [9, 59]. This salt bridge was observed in both chains of all the complexes. In L10F/N88S, the Asp29-Arg87 salt bridge was observed only in chain A of the protease. It shows that the substitution of serine at position 88 thereby results in loss of hydrogen bond interaction between Asp29 and Arg87 in

chain B, thus weakening the dimer stability. In L90M, a salt bridge between Asp29 and Arg8΄ was observed providing additional stability. According to Ode et al. [26] the residues Thr31, Thr74 and Asn88 formed a hydrogen bond network in AE-WT. We observed hydrogen bond between Thr31- Asn88 and Thr31΄Asn88΄. In L10F, hydrogen bond between Thr31 and Asn88 was observed, but showed a reduced occupancy than AE-WT. When L10F co-occurs with N88S, the hydrogen bond of Asn88 with Thr31 is lost; instead the Ser88 formed hydrogen bond interaction with Asp30 in the active site. These correlated changes from distal site to active site induce resistance by altering the intrasubunit hydrogen bond network, thereby reducing the interaction between Asp30 and NFV. The amino acids Lys45, Asp30 and Asn88 are in the vicinity of the binding pocket and the side chains of these residues are spatially adjacent to each other. Any mutation at one of these residues will affect the other two sites [60]. We observed that Asp30 and Lys45 formed hydrogen bonds only in L90M and L10F, whereas other structures lack this hydrogen bond interaction. The Asp60 and Thr74 which connects the β chains of the residues 52-66 and 69-78 [61] formed hydrogen bond in all complexes, but in L10F, the hydrogen bond occupancy is strongly reduced.

Hydrogen Bonds between protease and NFV: The formation of hydrogen bonds between protein and ligand appears to be crucial to the specificity of ligand recognition by proteins. The protease inhibitors bind with protease via both hydrogen bond and hydrophobic interactions [62]. The lowest energy structures calculated from MD simulation were considered to understand the changes in geometries of the key residues due to mutation. The binding of key residues with NFV in AE-WT and mutant structures are shown in Fig.10. The hydrogen bond interactions are mainly observed between the active site residues of HIV-1 protease and PIs. The combined changes in hydrogen bond and hydrophobic interactions are associated with drug resistance in HIV-1 protease. Ode et al.,[26] observed that in AE wild-type protease and N88S mutant, the side chains of catalytic residues Asp25 and Asp25΄ interact with the central hydroxyl group of NFV. They reported that the side chain of Asp30 formed hydrogen bond interaction with NFV in AE wild-type protease. They investigated the effect of N88S mutation on resistance mechanism. They reported that the hydrogen bond interaction between the side chains of Asp30 and Ser88 altered the location of Asp30. Though the location of Asp30 is altered, the

N88S mutation still retained the hydrogen bond between Asp30 and NFV, but the interaction of Asp30 with NFV was reduced.

We also observed similar hydrogen bond interaction between the side chains of either of the Asp25/Asp25΄ in AE-WT, L10F and L90M complex. Interestingly, in L90M, Asp30 also formed hydrogen bonds with NFV similar to AE-WT. Though in subtype B, NFV resistance occurs via L90M pathway; our results show that in subtype AE, the L90M mutant as such does not affect the binding affinity. In L10F/N88S mutant, the occupancy between Asp25 and NFV was reduced compared to AE-WT. The drug resistance mechanism of L10F and L10F/N88S against NFV is due to the following reasons. In AE-WT, the side chain of Asp30 formed a strong hydrogen bond with the benzamide moiety of NFV. In L10F and L10F/N88S, the P2 subsite of NFV showed a larger positional deviation, thus the hydrogen bond formation between Asp30 and benzamide moiety of NFV was lost. When L10F cooccurs with N88S, the hydrogen bond between residues at position Thr31 and Asn88 was disturbed; instead the side chain of Ser88 formed hydrogen bond with Asp30 which further impaired the interaction of Asp30 with NFV which is in agreement with the results reported by Ode et al. [26] for N88S mutation. Thus the co-occurrence of L10F with N88S mutation, enhances the resistance of NFV against AE-protease

Conclusions: The molecular basis of NFV resistance due to non-active site mutations L10F, L10F/N88S and L90M in AE-protease was studied using MD simulation, binding free energy calculation and per-residue decomposition. The L10F and L10F/N88S mutations reduce the binding affinity of NFV compared to AE-WT. The L90M mutation show increased binding affinity due to favourable enthalpy and entropic contributions. Extensive comparison of dynamics and energetics for all systems indicated the following reasons for NFV resistance in L10F and L10F/N88S complexes: The entropy component is unfavourable towards binding free energy in L10F complex. The hydrogen bond between m-phenol group of NFV and Asp30 is lost due to large positional deviation of NFV in L10F and L10F/N88S complexes. At the same time, due to L10F/N88S mutation the hydrogen bond is formed between Asp30 and Ser88 which causes Asp30 to lose its interaction with P2 subsite of NFV. Further the loss in electrostatic energy contribution from residues Ile47, Ile50 and Gly49΄ also decreases the enthalpy contribution towards binding free energy. Our simulation results could explain the

reason for the strong link between L10F and N88S mutation observed by Ariyoshi et al., experimentally. The findings from the present study have revealed the resistance mechanism of non-active site mutations in AE-protease and indicated that the preservation of hydrogen bonds between NFV and Asp30 is crucial for effective inhibition. This knowledge can give a clue in designing more effective protease inhibitor for drug resistant variants of a targeting enzyme.

Declaration of interest The authors declare no conflict of interest.

Acknowledgment: PM thanks the Department of Science and Technology (DST)-Science and Engineering Research Board (SERB), New Delhi, for financial assistance through the Fast-Track (YSS/2015/000955) project and IIIT for support. RT and CSV are grateful to the management of VIT University to carry out this study.

References [1] Gayle HD, Hill GL (2001) Global Impact of Human Immunodeficiency Virus and AIDS. Clin Microbiol Rev 14:327–335 [2] Troncoso A (2016) Will there be a cure for HIV/AIDS? Making the dream a reality. Asian Pacific Journal of Tropical Biomedicine 6:441–442 [3] Sharp PM, Hahn BH (2011) Origins of HIV and the AIDS Pandemic. Cold Spring Harb Perspect Med 1:a006841 [4] Kear JL, Blackburn ME, Veloro AM, et al (2009) Subtype Polymorphisms Among HIV-1 Protease Variants Confer Altered Flap Conformations and Flexibility. J Am Chem Soc 131:14650–14651 [5] Reeves JD, Piefer AJ (2005) Emerging drug targets for antiretroviral therapy. Drugs 65:1747–1766 [6] Sundquist WI, Kräusslich H-G (2012) HIV-1 Assembly, Budding, and Maturation. Cold Spring Harb Perspect Med 2:a006924 [7] Brik A, Wong C-H (2003) HIV-1 protease: mechanism and drug discovery. Org Biomol Chem 1:5–14 [8] Yang H, Nkeze J, Zhao RY (2012) Effects of HIV-1 protease on cellular functions and their potential applications in antiretroviral therapy. Cell Biosci 2:32 [9] Ishima R, Ghirlando R, Tözsér J, et al (2001) Folded Monomer of HIV-1 Protease. J Biol Chem 276:49110–49116 [10] Hughes PJ, Cretton-Scott E, Teague A, Wensel TM (2011) Protease Inhibitors for Patients With HIV-1 Infection. P T 36:332–345 [11] Ali A, Bandaranayake RM, Cai Y, et al (2010) Molecular Basis for Drug Resistance in HIV-1 Protease. Viruses 2:2509–2535 [12] Konvalinka J, Kräusslich H-G, Müller B (2015) Retroviral proteases and their roles in virion maturation. Virology 479–480:403–417 [13] Kantor R, Katzenstein DA, Efron B, et al (2005) Impact of HIV-1 Subtype and Antiretroviral Therapy on Protease and Reverse Transcriptase Genotype: Results of a Global Collaboration. PLoS Med. 2 [14] Taylor BS, Sobieszczyk ME, McCutchan FE, Hammer SM (2008) The Challenge of HIV-1 Subtype Diversity. N Engl J Med 358:1590–1602 [15] Manosuthi W, Butler DM, Pérez-Santiago J, et al (2010) Protease polymorphisms in HIV-1 subtype CRF01_AE represent selection by antiretroviral therapy and host immune pressure. AIDS 24:411–416 [16] Huang X, Britto MD, Kear-Scott JL, et al (2014) The role of select subtype polymorphisms on HIV-1 protease conformational sampling and dynamics. J Biol

Chem 289:17203–17214 [17] Colonno R, Rose R, McLaren C, et al (2004) Identification of I50L as the Signature Atazanavir (ATV)-Resistance Mutation in Treatment-Naive HIV-1Infected Patients Receiving ATV-Containing Regimens. J Infect Dis 189:1802– 1810 [18] Ragland DA, Nalivaika EA, Nalam MNL, et al (2014) Drug Resistance Conferred by Mutations Outside the Active Site through Alterations in the Dynamic and Structural Ensemble of HIV-1 Protease. J Am Chem Soc 136:11956–11963 [19] Weber IT, Agniswamy J (2009) HIV-1 Protease: Structural Perspectives on Drug Resistance. Viruses 1:1110–1136 [20] Santos AF, Soares MA (2010) HIV Genetic Diversity and Drug Resistance. Viruses 2:503–531 [21] Ode H, Ota M, Neya S, et al (2005) Resistant mechanism against nelfinavir of human immunodeficiency virus type 1 proteases. J Phys Chem B 109:565–574 [22] Ode H, Neya S, Hata M, et al (2006) Computational Simulations of HIV-1 ProteasesMulti-drug Resistance Due to Nonactive Site Mutation L90M. J Am Chem Soc 128:7887–7895 [23] Ariyoshi K, Matsuda M, Miura H, et al (2003) Patterns of point mutations associated with antiretroviral drug treatment failure in CRF01_AE (subtype E) infection differ from subtype B infection. J Acquir Immune Defic Syndr 33:336– 342 [24] Nukoolkarn S, Pongthapisith V, Panyim S, Leelamanit W (2004) Sequence variability of the HIV type 1 protease gene in thai patients experienced with antiretroviral therapy. AIDS Res Hum Retroviruses 20:1368–1372 [25] Bandaranayake RM, Kolli M, King NM, et al (2010) The Effect of Clade-Specific Sequence Polymorphisms on HIV-1 Protease Activity and Inhibitor Resistance Pathways. J Virol 84:9995–10003 [26] Ode H, Matsuyama S, Hata M, et al (2007) Mechanism of drug resistance due to N88S in CRF01_AE HIV-1 protease, analyzed by molecular dynamics simulations. J Med Chem 50:1768–1777 [27] Naicker P, Sayed Y (2014) Non-B HIV-1 subtypes in sub-Saharan Africa: impact of subtype on protease inhibitor efficacy. Biol Chem 395:1151–1161 [28] Rhee S-Y, Taylor J, Wadhera G, et al (2006) Genotypic predictors of human immunodeficiency virus type 1 drug resistance. Proc Natl Acad Sci U S A 103:17355–17360 [29] Yam WC, Chen JHK, Wong KH, et al (2006) Clinical utility of genotyping resistance test on determining the mutation patterns in HIV-1 CRF01_AE and subtype B patients receiving antiretroviral therapy in Hong Kong. J Clin Virol 35:454–457

[30] Soares RO, Batista PR, Costa MGS, et al (2010) Understanding the HIV-1 protease nelfinavir resistance mutation D30N in subtypes B and C through molecular dynamics simulations. J Mol Graph Model 29:137–147 [31] Liu X, Dai Q, Li L, Xiu Z (2011) Resistant mechanism against nelfinavir of subtype C human immunodeficiency virus type 1 proteases. Journal of Molecular Structure 986:30–38 [32] Kožíšek M, Bray J, Řezáčová P, et al (2007) Molecular Analysis of the HIV-1 Resistance Development: Enzymatic Activities, Crystal Structures, and Thermodynamics of Nelfinavir-resistant HIV Protease Mutants. Journal of Molecular Biology 374:1005–1016 [33] Xiang Z (2006) Advances in Homology Protein Structure Modeling. Curr Protein Pept Sci 7:217–227 [34] Fiser A (2010) Template-based protein structure modeling. Methods Mol Biol 673:73–94 [35] Lockhat HA, Silva JRA, Alves CN, et al (2016) Binding Free Energy Calculations of Nine FDA-approved Protease Inhibitors Against HIV-1 Subtype C I36T↑T Containing 100 Amino Acids Per Monomer. Chem Biol Drug Des 87:487–498 [36] Ahmed SM, Kruger HG, Govender T, et al (2013) Comparison of the molecular dynamics and calculated binding free energies for nine FDA-approved HIV-1 PR drugs against subtype B and C-SA HIV PR. Chem Biol Drug Des 81:208–218 [37] Antunes DA, Rigo MM, Sinigaglia M, et al (2014) New Insights into the In Silico Prediction of HIV Protease Resistance to Nelfinavir. PLOS ONE 9:e87520 [38] Kaldor SW, Kalish VJ, Davies JF, et al (1997) Viracept (Nelfinavir Mesylate, AG1343):  A Potent, Orally Bioavailable Inhibitor of HIV-1 Protease. J Med Chem 40:3979–3985 [39] Robert X, Gouet P (2014) Deciphering key features in protein structures with the new ENDscript server. Nucleic Acids Res 42:W320–W324 [40] Sali A, Blundell TL (1993) Comparative protein modelling by satisfaction of spatial restraints. J Mol Biol 234:779–815 [41] Lovell SC, Davis IW, Arendall WB, et al (2003) Structure validation by Calpha geometry: phi,psi and Cbeta deviation. Proteins 50:437–450 [42] Accelrys Discovery Studio (2012) San Diego: Accelrys Software Inc. Accelrys Discovery Studio modeling environment, release 3.5 [43] Case DA, Cheatham TE, Darden T, et al (2005) The Amber biomolecular simulation programs. J Comput Chem 26:1668–1688 [44] Götz AW, Williamson MJ, Xu D, et al (2012) Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J Chem Theory Comput 8:1542–1555

[45] Wang J, Wolf RM, Caldwell JW, et al (2004) Development and testing of a general amber force field. J Comput Chem 25:1157–1174 [46] Duan Y, Wu C, Chowdhury S, et al (2003) A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J Comput Chem 24:1999–2012 [47] Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, et al(2009) Gaussian 09, Revision E.01. Gaussian, Inc.; Wallingford, CT, USA [48] Jakalian A, Jack DB, Bayly CI (2002) Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J Comput Chem 23:1623–1641 [49] Vanquelef E, Simon S, Marquant G, et al (2011) R.E.D. Server: a web service for deriving RESP and ESP charges and building force field libraries for new molecules and molecular fragments. Nucleic Acids Res 39:W511–W517 [50] Ryckaert JP, Ciccotti G, Berendsen HJC(1977) Numerical-integration of cartesian equations of motion of a system with constraints – molecular-dynamics of nalkanes. J Comput Phys 23:327–341 [51] Yang M, Jiang X, Jiang N (2014) Protonation state and free energy calculation of HIV-1 protease–inhibitor complex based on electrostatic polarisation effect. Molecular Physics 112:1659–1669 [52] Wittayanarakul K, Hannongbua S, Feig M (2008) Accurate prediction of protonation state as a prerequisite for reliable MM-PB(GB)SA binding free energy calculations of HIV-1 protease inhibitors. J Comput Chem 29:673–685 [53] Rastelli G, Del Rio A, Degliesposti G, Sgobba M (2010) Fast and accurate predictions of binding free energies using MM-PBSA and MM-GBSA. J Comput Chem 31:797–810 [54] Onufriev A, Bashford D, Case DA (2004) Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 55:383–394 [55] Laskowski RA, MacArthur MW, Moss DS, Thornton JM (1993) PROCHECK: a program to check the stereochemical quality of protein structures. J Appl Cryst, J Appl Crystallogr 26:283–291 [56] Srivastava HK, Sastry GN (2012) Molecular dynamics investigation on a series of HIV protease inhibitors: assessing the performance of MM-PBSA and MM-GBSA approaches. J Chem Inf Model 52:3088–3098 [57] Huang YM, Kang M, Chang CA (2014) Switches of hydrogen bonds during ligand-protein association processes determine binding kinetics. J Mol Recognit 27:537–548 [58] Chen D, Oezguen N, Urvil P, et al (2016) Regulation of protein-ligand binding affinity by hydrogen bond pairing. Sci Adv

[59] Ishima R, Gong Q, Tie Y, et al (2010) Highly conserved glycine 86 and arginine 87 residues contribute differently to the structure and activity of the mature HIV-1 protease. Proteins 78:1015–1025 [60] Mahalingam B, Louis JM, Hung J, et al (2001) Structural implications of drugresistant mutants of HIV-1 protease: high-resolution crystal structures of the mutant protease/substrate analogue complexes. Proteins 43:455–464 [61] York DM, Darden TA, Pedersen LG, Anderson MW (1993) Molecular dynamics simulation of HIV-1 protease in a crystalline environment and in solution. Biochemistry 32:1443–1453 [62] Yu Y, Wang J, Shao Q, et al (2015) Effects of drug-resistant mutations on the dynamic properties of HIV-1 protease and inhibition by Amprenavir and Darunavir. Scientific Reports 5:10517

Fig. 1. (a)The positions of NFV resistance mutations (L10, N88, and L90), the catalytic triad (Asp25- Gly27), and the residue Asp30 crucial for NFV binding in AE-protease are highlighted in red. All residues and substitutions occur in both the subunits (b) Chemical structure of NFV with only the sulphur, nitrogen and oxygen atoms labelled.

Fig. 2. Amino acid sequence alignment of AE-WT, L10F, L10F/N88S and L90M against the subtype B protease. The polymorphisms and mutations are highlighted in white color.

Fig. 3. Root mean square displacement of the Cα coordinates in AE-WT, L10F, L10F/N88S and L90M complexes as a function of time.

Fig. 4. Root mean square fluctuation (RMSF) of backbone atoms N, Cα and C for every residue in AE-WT, L10F, L10F/N88S and L90M complexes.

Fig. 5. The differences of root mean square fluctuation (RMSF) values of every residue in AE-WT, L10F, L10F/N88S and L90M complexes.

(a)

(b)

Fig. 6. Distributions of the (a) Cα−Ile50/Cβ−Asp25 and the (b) Cα−Ile΄50/Cβ−Asp΄25 distances for AE-WT, L10F, L10F/N88S and L90M complexes.

Fig. 7. Energy components (kcal/mol) for the binding of NFV to the AE-WT, L10F, L0F/N88S and L90M complexes. Evdw , van der waals energy; Eele , electrostatics energy in gas phase; Gpol, polar solvation energy; Gnp, nonpolar solvation energy; H, entropy contribution; TSconf, , entropy contribution, Gbind = H−TS.

(a)

(b)

(c) (d) Fig. 8. Decomposition of G on a per-residue basis for the (a) AE-WT, (b) L10F, (c) L10F/N88S and (d) L90M complexes.

Fig. 9. Distances between O4 atom of NFV and OD2 atom of Asp30 in AE-WT, L10F, L10F/N88S and L90M complexes.

Fig. 10. The key residues interacting with NFV in (a) AE-WT, (b) L10F, (c) L10F/N88S and (d) L90M. The key residues and NFV are shown in stick representation. The hydrogen bond interactions and the distances between NFV and proteases are shown with a dashed line.

Table 1. Decomposed contributions to the free energy of binding of NFV to AE-WT,L10F, L10F/N88S and L90M complexes. All energies are given in kcal/mol. Varian ts AEWT AEWT§ AEN88S§ L10F L10F/ N88S L90M

Evdw

Eele

Gpol

Gnp

Gsol

H

−65.7 (3.6) −70.8 (4.0) −67.8 (4.0) −68.2 (3.8) −69.8 (3.9) −63.6 (4.2)

−45.2( 4.1) −12.8 (1.7) −10.6 (1.4) −36.5 (4.7) −26.6 (4.3) −51.7 (4.6)

52.6 (3.8)

−8.8 (0.3)

−

−

−

−

44.8 (3.3) 41.9 (3.8) 53.7 (3.8)

−8.9 (0.2) −8.8 (0.3) −8.5 (0.3)

43.8 (3.8) 15.1 (1.4) 12.8 (1.9) 35.9 (3.3) 33.2 (3.7) 45.2 (3.9)

−67.2 (4.2) −68.6 (3.7) −65.6 (3.9) −68.8 (3.8) −63.3 (4.0) −70.1 (3.9)

−TSc onf

Gbind

Gbin d

27.4 (5.9)

−39.8

−

−

−

−

3.0

33.3 (5.0) 28.9 (3.5) 26.9 (6.8)

−35.5

4.3

−34.4

5.4

−43.2

−3.4

§

Results from Ref [26]., Gbind=Gmutant−GAE The values in parentheses are the standard deviation Table 2. Decomposition of binding energies for the AE-WT complex on a per-residue basis. All energies are reported in kcal/mol. AE−WT Residue Svdw Bvdw Tvdw Sele Bele Tele SPB BPB TPB TGBTOT −0.83 0.00 −0.35 0.00 −0.35 0.15 0.00 0.15 −1.03 Asp25 −0.83 −0.04 −1.31 −1.35 0.22 −0.58 −0.36 −0.40 0.85 0.45 −1.26 Gly27 −0.99 −1.17 −2.16 −0.02 −0.44 −0.45 0.08 0.43 0.51 −2.10 Ala28 0.39 −0.43 −0.04 −7.91 −0.15 −8.06 6.39 0.41 6.80 −1.30 Asp30 −0.80 −0.18 −0.98 −0.60 0.13 −0.47 0.51 −0.10 0.40 −1.05 Ile47 b −0.08 −0.87 −0.96 −1.88 0.55 −1.34 1.79 −0.38 1.41 −0.89 Gly49 −1.53 −0.61 −2.14 −1.14 0.37 −0.77 0.97 −0.03 0.95 −1.96 Ile50 −1.39 −0.13 −1.52 0.03 0.09 0.12 0.00 −0.12 −0.12 −1.52 Ile84 Asp25΄ 0.32 −0.12 0.19 −9.12 −0.45 −9.57 6.51 0.60 7.11 −2.27 Ala28΄ −0.62 −1.08 −1.71 −0.03 0.90 0.87 −0.04 −0.15 −0.19 −1.03 Gly49΄ b −0.15 −1.19 −1.34 −1.35 −0.02 −1.36 1.26 0.53 1.80 −0.90 Ile50΄ −1.61 −0.86 −2.47 −1.03 −0.67 −1.70 0.99 0.93 1.92 −2.25 Pro81΄ −0.77 −0.22 −0.99 −0.75 0.18 −0.57 0.70 −0.16 0.54 −1.02 Val82΄ b −0.52 −0.20 −0.72 −0.08 −0.01 −0.09 −0.09 0.04 −0.06 −0.87 Ile84΄ −1.58 −0.17 −1.75 −0.02 −0.03 −0.05 0.01 0.01 0.02 −1.78 a Energies are shown as contributions from van der Waals energy (vdw), electrostatics energy (elec), polar solvational energy (PB) of side-chain atoms (S), backbone atoms (B), and sum of them (T) of AE-WT-NFV complex. Residues of |G|≥1.0kcal/mol were listed. bResidueswith energies that did not meet |G|≥1.0kcal/mol.

Table 3. Decomposition of binding energies for the L10F complex on a per-residue basis.All energies are reported in kcal/mol. L10F Residue Svdw

Bvdw

Asp25

0.63

−0.12

Gly27

−0.06

Ala28

Tvdw

Sele

Bele

Tele

SPB

BPB

TPB

0.51

−8.17

−0.29

−8.46

6.98

0.35

7.33

−0.62

−0.82

−0.87

0.20

0.07

0.27

−0.27

0.17

−0.10

−0.70

−0.47

−1.26

−1.73

0.17

−0.51

− 0.34

−0.09

0.52

0.43

−1.64

Asp30

−0.28

−0.25

−0.53

1.05

−2.98

−1.93

−0.70

1.38

0.68

−1.78

Ile47

−1.15

−0.38

−1.53

−0.41

−0.18

−0.59

0.38

0.11

0.49

−1.63

Gly49

−0.07

−0.84

−0.91

−1.15

0.22

−0.94

1.14

0.19

1.34

−0.51

Ile50

−1.55

−0.69

−2.24

−0.82

−0.26

−1.08

0.72

0.57

1.29

−2.03

Ile84

−1.42

−0.17

−1.59

0.11

−0.10

0.01

−0.04

0.03

−0.01

−1.59

Asp25΄

−0.08

0.00

−0.08

−3.88

0.00

−3.88

1.39

0.00

1.39

−2.57

Ala28΄

−0.73

−1.19

−1.92

−0.08

0.66

0.58

0.01

0.06

0.07

−1.27

Gly49΄

−0.14

−1.42

−1.57

−1.59

0.26

−1.32

1.44

0.55

1.99

−0.90

Ile50΄

−1.40

−0.89

−2.28

−0.72

−0.56

−1.28

0.70

0.78

1.48

−2.08

Pro81΄

−0.57

−0.16

−0.73

−0.33

0.06

−0.26

0.32

−0.05

0.27

−0.72

Val82΄

−0.61

−0.16

−0.77

−0.01

0.05

0.04

−0.15

0.03

−0.12

−0.85

Ile84΄ −1.42 −0.11 −1.53 −0.01 −0.04 −0.05 0.03 0.07 0.10 Energies are shown as contributions from van der Waals energy (vdw), electrostatics energy (elec), polar solvational energy (PB) of side−chain atoms (S), backbone atoms (B), and sum of them (T) for L10F-NFVcomplex.

TGBTOT

−1.48

Table 4. Decomposition of binding energies for the L10F/N88S complex on a per-residue basis. All energies are reported in kcal/mol. L10F/ N88S Residue Svdw Bvdw Tvdw Sele Bele Asp25 −0.13 −0.16 −0.29 −8.59 −0.28 Gly27 −0.01 −1.39 −1.40 0.06 0.16 Ala28 −0.74 −1.41 −2.15 0.13 −0.46 Asp30 −0.25 −0.50 −0.75 0.42 −2.61 Ile47 −0.88 −0.11 −0.99 −0.35 0.04 Gly49 −0.10 −0.87 −0.97 −1.84 0.73 Ile50 −1.43 −0.50 −1.93 −0.54 −0.11 Ile84 −1.40 −0.16 −1.56 0.15 −0.08 Asp25΄ −0.82 0.00 −0.82 −0.14 0.00 Ala28΄ −0.70 −1.19 −1.89 −0.08 0.84 Gly49΄ −0.26 −1.35 −1.62 −1.32 0.37 Ile50΄ −1.56 −0.93 −2.49 −0.71 −0.81 Pro81΄ −0.49 −0.15 −0.65 −0.44 0.10 Val82΄ −0.63 −0.18 −0.81 −0.04 0.06 Ile84΄ −1.53 −0.15 −1.68 −0.05 − 0.01 Energies are shown as contributions from van der Waals solvational energy (PB) of side−chain atoms (S), backbone NFV complex.

Tele SPB BPB −8.87 7.08 0.42 0.22 −0.24 0.45 −0.33 −0.09 0.61 −2.19 −0.19 1.51 −0.31 0.34 −0.07 −1.11 1.69 −0.37 −0.64 0.57 0.51 0.07 −0.09 0.02 −0.14 −0.19 0.00 0.75 −0.02 0.16 −0.95 1.35 0.03 −1.52 0.77 0.97 −0.34 0.43 −0.07 0.02 −0.12 0.00 −0.07 0.04 0.09 energy (vdw), electrostatics atoms (B), and sum of them

TPB TGBTOT 7.50 −1.66 0.20 −0.97 0.52 −1.95 1.32 −1.62 0.26 −1.04 1.32 −0.76 1.08 −1.49 −0.07 −1.56 −0.19 −1.15 0.13 −1.00 1.39 −1.17 1.74 −2.27 0.35 −0.64 −0.12 −0.91 0.13 −1.62 energy (elec), polar (T) for L10F/N88S-

Table 5. Decomposition of binding energies for the L90M complex on a per-residue basis. All energies are reported in kcal/mol. L90M Residue Svdw Bvdw Tvdw Sele Bele Tele SPB BPB TPB Asp25 0.59 −0.11 0.48 −8.67 −0.24 −8.92 6.84 0.30 7.15 Gly27 −0.04 −1.19 −1.24 −0.01 −0.07 −0.09 −0.21 0.37 0.16 Ala28 −0.82 −1.40 −2.23 0.05 −0.17 −0.12 −0.01 0.33 0.31 Asp30 0.37 −0.56 −0.18 −7.70 −0.49 −8.20 5.87 0.70 6.58 Ile47 −0.88 −0.31 −1.20 −0.68 −0.05 −0.73 0.55 0.05 0.60 Gly49 −0.06 −0.66 −0.73 −1.08 0.28 −0.79 1.05 −0.11 0.94 Ile50 −1.41 −0.47 −1.88 −0.89 0.54 −0.35 0.81 −0.03 0.78 Ile84 −1.15 −0.11 −1.26 0.11 0.11 0.23 −0.08 −0.15 −0.24 Asp25΄ −0.04 0.00 −0.04 −4.64 0.00 −4.64 1.76 0.00 1.76 Ala28΄ −0.76 −1.28 −2.05 −0.09 0.74 0.65 0.01 0.11 0.12 Gly49΄ −0.19 −1.23 −1.43 −1.61 0.19 −1.42 1.47 0.42 1.89 Ile50΄ −1.39 −0.83 −2.22 −0.66 −0.70 −1.36 0.75 0.92 1.67 Pro81΄ −0.69 −0.17 −0.86 −0.68 0.15 −0.53 0.62 −0.12 0.50 Val82΄ −0.50 −0.16 −0.66 −0.07 0.02 −0.05 −0.05 0.02 −0.03 Ile84΄ −1.37 −0.11 −1.48 −0.05 −0.03 −0.09 0.09 0.04 0.13 Energies are shown as contributions from van der Waals energy (vdw), electrostatics energy (elec), polar solvational energy (PB) of side−chain atoms (S), backbone atoms (B), and sum of them (T) of L90M-NFV complex.

TGBTOT −1.28 −1.16 −2.03 −1.80 −1.33 −0.58 −1.45 −1.27 −2.92 −1.28 −0.96 −1.91 −0.89 −0.74 −1.44