Revealing the drug resistance mechanism of saquinavir due to G48V and V82F mutations in subtype CRF01_AE HIV-1 protease: molecular dynamics simulation and binding free energy calculations

Abstract Human immunodeficiency virus-1 (HIV-1) protease is one of the important targets in AIDS therapy. The majority of HIV infections are caused due to non-B subtypes in developing countries. The co-occurrence of mutations along with naturally occurring polymorphisms in HIV-1 protease cause resistance to the FDA approved drugs, thereby posing a major challenge in the treatment of antiretroviral therapy. In this work, the resistance mechanism against SQV due to active site mutations G48V and V82F in CRF01_AE (AE) protease was explored. The binding free energy calculations showed that the direct substitution of valine at position 48 introduces a bulkier side chain, directly impairing the interaction with SQV in the binding pocket. Also, the intramolecular hydrogen bonding network of the neighboring residues is altered, indirectly affecting the binding of SQV. Interestingly, the substitution of phenylalanine at position 82 induces conformational changes in the 80’s loop and the flap region, thereby favoring the binding of SQV. The V82F mutant structure also maintains similar intramolecular hydrogen bond interactions as observed in AE-WT. Communicated by Ramaswamy H. Sarma


Introduction
HIV-1 protease enzymes belong to the family of aspartyl proteases and play a crucial role in the life cycle of HIV by cleaving the Gag and Gag-pol polyproteins into the production of structural (MA, CA, NC, and P6) and functional virions (protease, reverse transcriptase, and integrase) that help in viral assembly (Sundquist & Kr€ ausslich, 2012). HIV-1 protease has broad substrate specificity and cleaves the peptide bonds between proline residue and aromatic residues (phenylalanine/tyrosine) in polypeptides leading to viral maturation (Brik & Wong, 2003;Potempa et al., 2018). The maturation of the virus particle is essential for infectivity rendering the protease enzyme a major therapeutic target for the treatment of HIV infection (Wlodawer, 1998).
The diversity in HIV-1 protease has given rise to different subtypes (A-D, F-H, and J-K) and recombinant strains (CRF01_AE and CRF01_AG). These subtypes carry naturally occurring polymorphisms, some of which are associated with drug resistance in subtype B protease (A. F. Santos & Soares, 2010). Although subtype B is more prevalent in the western world, the circulation of non-B subtypes in North America, South America, and Western Europe seems to be increasing presumably (De Paschale et al., 2011). Limited data are confirming the effectiveness of existing HIV-1 protease drugs against the non-B subtypes and CRFs (Bessong, 2008;Clemente et al., 2006;Velazquez-Campoy et al., 2001).
The high genetic variability in HIV-1 subtypes may influence the resistance pathways to the currently available antiretroviral therapy. In HIV-1 protease, it has been shown that when drug pressure-induced mutations are concurrent with these polymorphisms, the binding affinity of the PIs is dramatically decreased which is one of the major challenges in AIDS therapy (Huang et al., 2014). Currently, subtype CRF01_AE (AE) has formed a global pandemic and is expanding rapidly, infecting a large proportion of people, throughout the world (Giallonardo et al., 2021;Sun et al., 2020). The polymorphisms of AE protease (I13V, E35D, M36I, S37N, R41K, H69K, and L89M) differ by $10% of amino acid sequence from the first identified and characterized clade subtype B protease (Ode et al., 2007).
The most common mutations associated with antiretroviral drug resistance in subtype B are G48V and V82F mutations, and they are considered as unique mutations characteristically selected by saquinavir (SQV). In subtype B HIV-1 protease, G48V is a unique mutation characteristically selected against SQV therapy. Maschera et al. (1996) reported that mutations G48V, L90M, and G48V/L90M arise in HIV-1 protease in the presence of SQV and indicated that the binding affinity of SQV was much reduced for mutant proteases compared to NFV, IDV, and APV. Various studies showed that prolonged exposure of SQV in patients resulted in G48V, G48V/L10I, V82A, and L90M mutations (Boucher, 1996;Deeks et al., 1997;Schapiro et al., 1996).
In subtype B protease, the effectiveness of SQV is hindered because SQV shifts way from the protease to accommodate the side chain of valine at position 48, thereby resulting in larger gaps between P3 quinoline ring of SQV and flap region of the enzyme, further leading to disruption of hydrogen bonding interactions (Hong et al., 1997(Hong et al., , 2000Saen-oon et al., 2007;Wittayanarakul et al., 2005). Aruksakunwong et al. (2006) emphasized the importance of protonation of catalytic residue in the G48V mutant complex. Stoica et al. (2008) performed MD simulations and binding free energy calculations and showed that G48V mutant result in unfavorable balance of electrostatic and polar desolvation energies.
Atomic resolution crystal structures of V82A mutant complexed with SQV show that the 80's loop has undergone conformational changes to accommodate the large decahydroisoquinoline ring of SQV at P1subsite (Tie et al., 2007). Mosebi et al. (2008) reported that in subtype C protease, the binding affinity of SQV with V82A mutant was due to loss in entropy. Ahmed et al. (2013) performed MD simulations and binding free energy calculations to understand the impact of active site mutations V82A and V82F/I84V and showed that the both the single and double mutation had unfavorable impact on the binding of SQV. Tzoupis et al. elucidated the SQV resistance mechanism through MM-PBSA and thermodynamic integration calculations. Their results show that G48V confers major resistance due to lack of hydrogen bond interactions with active site and flap residues (Tzoupis et al., 2013).
Similar to subtype B protease, in AE protease also, the mutations G48V and V82F are more commonly present. Clemente et al. showed that the V82F mutation in AE protease exhibit increased binding to SQV relative to AE WT and also subtype B protease (Clemente et al., 2006). Studies by Hsu et al. and Nukoolkarn et al. show that the mutation of Glycine to Valine at position 48 (G48V) has been observed in patients undergoing SQV therapy in AE protease (Hsu et al., 2005;Nukoolkarn et al., 2004). These experimental and clinical studies prompted to investigate the mechanism behind higher affinity of SQV in the V82F mutant and resistance mechanism of G48V mutant in AE protease through molecular dynamics and binding free energy calculations. The positions of SQV resistance mutations G48 and V82, the catalytic triad (Asp25-Gly27) is shown in Figure 1(a), and the chemical structure of SQV is shown in Figure 1(b). The findings from this study would hint in designing more effective inhibitors against subtype AE protease.

Materials and methods
A detailed flowchart in Figure 2 shows the methodology followed for molecular dynamics simulation and binding free energy calculations of AE-WT and mutants complexed with SQV.

Homology modeling
The modeling and preparation of the protease systems was followed using the same protocol established in our previous research work (Vasavi et al., 2017(Vasavi et al., , 2019. The sequences of AE-protease exhibited > 90% identity with the template structure of HIV-1 subtype B protease complexed with SQV (PDB id: 1HXB). In the present study, the AE-WT structure was modeled by incorporating the polymorphisms (I13V, E35D, M36I, S37N, R41K, H69K, and L89M) within the framework of subtype B protease, as there are no experimentally resolved protein structures for AE protease complexed with SQV. Further, the mutations (G48V and V82F) were introduced in AE protease to create the respective mutant structures.
Each initial structure of AE-WT and mutants in complex with SQV were modeled using Modeller 9.17 (Sali & Blundell, 1993) with atomic coordinates of template structure along with its crystallographic water molecules. Discrete Optimized Protein Energy (DOPE) score was used to evaluate the generated models. Further, the stereochemical quality of the modeled structures was assessed using SAVES server (Laskowski et al., 1993), QMEAN (Benkert et al., 2009), and ANOLEA (Melo et al., 1997). The Ramachandran plots were generated using Discovery Studio 3.5 Visualizer (San Diego: Accelrys Software Inc., 2012). The generated models with lowest DOPE score validated by SAVES server were considered for energy minimization and molecular dynamic simulations.

System preparation and molecular dynamics simulation
Energy minimization (G€ otz etal., 2012) and MD simulations were performed using AMBER16 package (Case et al., 2016). For the protein residues, the AMBER force field ff14SB parameters (Maier et al., 2015) were used to describe the system. The force field parameters for the inhibitor SQV were assigned by the updated General AMBER Force Field (GAFF2) (Case et al., 2005;Wang et al., 2004). As the atomic partial charges play an important role in electrostatic and other non-covalent interactions, the following procedure was used to calculate the atomic partial charges for SQV. The LEap module of AMBER16 (Duan et al., 2003) was used to add the missing hydrogen atoms of SQV, extracted from crystal structure (PDB id: 1HXB), and the resulted structure was then optimized using (HF)/6-31G Ã basis set of Gaussian 09 package (Frisch et al., 2010). Further, these electrostatic potentials were used to obtain RESP charges using Antechamber module. Each complex system was then solvated by placing it in a TIP3P water cubic box with at least 10 Å of solvent on each face of the AE protease. Chloride counter ions were added to neutralize the positively charged system. A three-step energy minimization was performed to allow the system to reach an energetically favorable conformation. In each of the three steps, energy minimization was executed by the steepest descent method for the first 2500 cycles and the conjugate gradient method for the subsequent 2500 cycles. In the first energy minimization step, the restraints were placed on all heavy atoms of the AE-proteases with a harmonic force constant of 10 kcal mol À 1 Å -2 . In the second step, restraints were applied only to the backbone nitrogen, oxygen, and carbon atoms, and 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. The temperature of each energy minimized system was gradually  raised from 0 K to 300 K in NVT ensemble over a period of 100ps. Following heating, 100 ps of constant-pressure equilibration was carried out at 1 atm. During heating and constant-pressure equilibration process, the solute heavy atoms were restrained with a harmonic force constant of 10 kcal mol À1 Å -2 . The force constant on the restrained atoms was then reduced in the following order 8, 4, and 2 kcal mol À1 Å -2 , and the equilibration had been repeated for 50 ps at each value. Then, unrestrained constant-pressure equilibration was extended for 500 ps. Following the equilibration, a total length of 100 ns MD simulations were performed in the NPT ensemble at a temperature of 300 K by using Langevin dynamics temperature scaling with collision frequency 2 ps À1 , and the coordinates of the trajectories were saved at every 10 ps for each simulation. The simulation time step was set to 2 fs. All the bond lengths involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977). The particle mesh Ewald (PME) method was used to handle long-range electrostatic interactions and a cutoff of 12.0 Å was used for MD simulations (Castro et al., 2021). The trajectories obtained from the last 25 ns were used to analyze the structure in detail. The plots were constructed using the Xmgrace program (Turner, 2005), and VMD program (Humphrey et al., 1996) was used for analyzing the trajectories. Further, to understand the energetic contribution of each residue in binding of SQV, the per-residue energy decomposition was computed using the MM-GBSA decomposition method. Also, the hydrogen bond was analyzed using PTRAJ module of AMBER16 program. The methodology followed for this work is similar to our previous studies on understanding the drug resistance mechanism for NFV and IDV using MD simulations and Binding Free energy calculations (Vasavi et al., 2017(Vasavi et al., , 2019).

Determination of protonation state
As protonation state of the Asp dyad varies depending upon different PIs and they remarkably affect the binding mechanism of protease-inhibitor complex, a proper determination of the protonation state of catalytic aspartates Asp25/Asp25ˊis important. In this work, the following protonation states of catalytic aspartates have been considered: The aspartate in the A chain protonated and B chain unprotonated. The aspartate B chain is protonated and chain A is unprotonated. The protonation states of AE-WT, G48V, and V82F structures complexed with SQV were determined by calculating the binding energies using MM/GBSA method (Rastelli et al., 2010). Further, for the favorable protonation state, the simulation time length is increased to understand the resistance mechanism in detail.

Binding free energy calculation
The binding free energy of AE protease and SQV was calculated using MM-GBSA method, and the binding free energy of a system is expressed in terms of enthalpic and entropic contributions (Equation (1)): The enthalpy (DH) of binding is given by where DE ele and DE vdw in Equation (2) define the interaction energy between the protease and SQV computed with the molecular mechanics method, and DG sol denotes the change in solvation free energy upon ligand binding. Further, the solvation energy in Equation (2) is defined by the combination of polar (DG pol ) and non-polar (DG np ) contributions (Equation (3)).
The polar solvational energy (DG polar ) in Equation (3) was approximated using generalized Born (GB) model of AMBER16 program, and the non-polar contributions was determined using the linear combinations of pairwise overlaps (Weiser et al., 1999) (LCPO) method as described in Equation (4).
For the surface tension g and offset b, the default values of 0.00542 kcal mol À1 Å -2 and 0.92 kcal mol À1 were considered, respectively. The dielectric constant for the solute and solvent was set to 1.0 and 80.0, respectively (Cruz et al., 2019;C. B. R. Santos et al., 2021).
The conformational entropy upon ligand binding was estimated by performing the normal-mode analysis using nmode module of AMBER16. The contributions of entropy (TDS) to binding free energy result from changes in the translational, rotational and vibrational degrees of freedom.

Per-residue energy decomposition
Per-residue energy decomposition (Ara ujo et al., 2020;Leão et al., 2020) was computed using the MM-GBSA decomposition method to examine the energetic contribution of each residue in binding of SQV. The SQV-residue interaction is approximated by three terms: van der Waals (DE vdw ) and electrostatic interactions (DE ele ) between SQV and each residue in gas phase and polar contribution to solvation energy (DG PB ) was calculated using the GB model, and the parameters were developed by Onufreiv and Bashford (Onufriev et al., 2004).

Hydrogen bond criterion
The hydrogen bonds were analyzed using the PTRAJ module of AMBER16 program. The formation of hydrogen bonds was defined in terms distance and angle cutoff as follows: (a) the distance between donor D and acceptor A atoms was 3.5 Å, and (b) the angle between D-hydrogen-A was !120 .

Homology modeling of AE-WT and mutant structures complexed with SQV
The sequences of AE-protease exhibited > 90% identity with the template structure of HIV-1 subtype B protease complexed with SQV (PDB id: 1HXB) (Krohn et al., 1991). The AE-WT wild type was constructed by substituting the following amino acids I13V, E35D, M36I, S37N, R41K, H69K, and L89M in subtype B protease sequence. The mutations G48V Aand V82F were introduced in AE-WT to create respective mutant structures. The amino acid sequence alignment of subtype B, AE-WT, and mutant structures G48V and V82F, shown in Figure 3, was generated using ESPript (Robert & Gouet, 2014). The polymorphisms and mutations are highlighted in yellow color. The Ramachandran plot for the modeled structures indicated that $98% of the residues were present in the favored region for AE-WT, G48V, and V82F protease models. 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 Figure S1. The QMEAN score ranges between 0 and 1 with higher value to reflect a better quality of the input model. The QMEAN score for the AE-WT, G48V mutant is 0.81 and for V82F mutant it is 0.83. Another tool ANOLEA was also used to examine the local model quality, and the results were in good agreement with the results by QMEAN, indicating that majority of the residues in the computational model were located in the favorable regions.

Comparing the stability of trajectories from RMSD: AE-WT versus mutant
To explore the effect of mutations G48V and V82F on the conformational stability of the AE-WT, the RMSD values for the Ca backbone coordinates relative to the initial structures were calculated in Figure 4. The RMSD indicate that the conformations of AE-WT, G48V, and V82F complexes were in the good equilibrium. All the three trajectories go almost parallel to each other till 38 ns; after 38 ns till 48 ns, trajectories show the maximum fluctuation of 1.8 Å in AE-WT compared to mutant structures. From 71 ns until 80 ns, the fluctuation of mutant structures seems to be less compared to AE-WT. According to the RMSD average values of the three complexes, G48V has a mean (1.18 Å) than the AE-WT and V82F structures (1.16 Å), with respective standard deviations of 0.15, 0.14, and 0.11 Å. The result signifies that all the three complexes deviate to a quite similar extent from their starting structures with values around 1.0 to 1.7 Å during simulation, ensuring stable trajectories. This suggests that the stabilities of the three complexes were reliable and could provide a suitable basis for the subsequent analyses.

Comparing the atomic fluctuations: AE-WT versus mutant
To understand the average fluctuations taken over the MD simulation, detailed analyses of RMSF of backbone atoms versus residue number for AE-WT and mutant complexes were calculated. Although Figure 5 shows that all the structures show similar RMSF distributions, the AE-WT and G48V complexes have slightly larger RMSF values than the V82F complex. The average residual atomic fluctuations for the V82F (0.78 Å) complexes seem to be slightly lesser than AE-WT (0.79 Å) and G48V (0.80 Å) complexes. The residues Ile15 and Gly16 in the G48V complex exhibit large fluctuations than the other two structures. In AE-WT and G48V complexes, the residues in the flap region seem to be more flexible compared to the V82F complex. This might affect the change in significant interactions between the flaps and as a result leads to less effective entrapment of SQV into the binding cavity. The residues 79-82 in the 80's loop show less flexibility in the V82F complex compared to AE-WT and G48V structures. The difference in RMSF values of AE-WT and mutant complexes is shown in Figure 6(a). The residues Gly16 (fulcrum), Asp35/35ˊ-Asn37/37ˊ(flap elbow) of G48V complex, and the residues Cys67 (cantilever), Ile15-Gln18 (fulcurum) of V82F complex show absolute difference larger than 0.5 Å and are considered to be highly fluctuating residues. The sites that show greatest difference in RMSF have been highlighted in Figure 6(b).

Distance between flap tip to catalytic asp (Ile50 2 Asp25)
The extent of flap opening was measured by the distance between the flap tip (Ile50C a /Ile50ˊC a ) and the catalytic aspartate (Asp25C a /Asp25ˊC a ). The Ile50-Asp25 and Ile50-Asp25ˊdistances for AE-WT and mutant complexes (G48V and V82F) were measured from the MD trajectories. The histogram distributions are shown in Figure 7. The mean distances of flap tip-catalytic Asp in chain A were 15.42, 16.03, and 14.25 with SD values of 0.56, 0.61, and 0.34, respectively, for AE-WT, G48V, and V82F, respectively (Figure 7(a)). This indicates that the distance between flap tips and the catalytic aspartate is smaller in the V82F complex by $ 1.0 Å than AE-WT and in G48V complex the distance is stretched by $ 0.6 Å. For chain B, as seen in Figure 7(b), the mean distances and SD values AE-WT (G48V) are 13.66 (15.25) Å and 0.74 Å (0.41), respectively, whereas for the V82F complex, the mean and SD were 14.35 and 0.30, respectively. The mean of distribution for G48V and V82F complexes was greater by $1.60 Å and $ 0.7 Å, respectively, than AE-WT. In general, the flap tip-catalytic aspartate distance is different for the AE-WT and G48V complexes in both chains A and B, whereas in the V82F complex, the distance is similar in both the monomers. It was demonstrated that in chain A of V82F mutant complex, the distance between flap tips and the active site is smaller, whereas in G48V, the distance is farther than AE-WT. From this analysis, it seems that for the G48V complex, the average distance from flap tip to the active site is longer in both the chains compared to AE-WT, making the active site volume marginally larger. The residue Gly48 lies close to the flap tip. Therefore, mutations in this position are likely to affect the overall structure of the protein. Studies have reported that substitution of Val, a hydrophobic residue at position 48, causes a specific structural change at the tip of the flap (Ile50/Ile50ˊ) (Saen-oon et al., 2007;Tie et al., 2007;Tzoupis et al., 2013). We observed that the residue Ile50 of both the chains have a large conformational change in G48V structure which would have probably induced changes in Ile50C a /Ile50ˊC a -Asp25C a /Asp25ˊC a distances thereby deforming the active site cavity. On examining the V82F mutant structure, large conformational changes were observed for Ile50 and Asp25 residues. This is because the 80's loops (residues 79-84) lie in close contact with Ile50á nd Asp25. The changes in the conformation of catalytic and flap residue would probably make the distance smaller in volume, thereby enhancing the binding affinity of SQV.

Total binding free energy: AE-WT versus mutants
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 (G48V, V82F) complexes using molecular mechanics GB surface area method in AMBER. A summary of the contributions of the binding interactions between SQV and all the complexes is given in Figure 8 and Table 1. As shown in Table 1, the calculated binding free energies of AE-WT, G48V, and V82F complexes are -32.1, -29.2, and -41.0 kcal/mol, respectively. This shows that the V82F complex shows an enhanced binding affinity and G48V mutant demonstrate drug resistance to SQV. According to the binding free energy components, the major contributions to SQV binding were from van der Waals interaction (DE vdw ), closely followed by electrostatic interactions (DE ele ). Non-polar solvation energies (DG np ), resulting from the burial of solvent accessible surface area upon SQV binding provides a slightly favorable contribution to the binding energy, indicating that the cavity region in all the three systems are well packed. Conversely, the polar solvational energies (DG pol ) and entropy components (ÀTDS) produced unfavorable contributions to the binding energy. The loss in binding affinity of SQV in the G48V mutant is due to a decrease in van der Waals interaction and the electrostatic interaction. One of the major reason could be substitution of Val at position 48/ 48ˊin the flap region of the enzyme introduces a bulkier side chain into the binding pocket directly impairing the interaction with SQV. On the other hand, owing to the close contact of 80's loop with flaps and active site cavity, probably, the V82F mutant induces conformational change thereby resulting in favorable interactions with SQV.
To understand the drug resistance mechanism, the binding energies of AE-WT and mutant complexes were decomposed into inhibitor-residue pairs to create an inhibitorresidue spectrum as shown in Figure 9. The per-residue contributions were calculated using the MM-GBSA method. This per-residue decomposition method is especially useful to explain the drug resistance mechanism at the atomistic level and also helps to locate the individual residue contribution to protein-inhibitor interactions. The energy contributions of both backbone and side-chain atoms for each active site residues are shown in Tables 2 to 4. From Figure 9, it is clear that the interaction spectra were a little different for all the three complexes. Active site residues that obtained a total energy value below -1kcal/mol were named in each complexes, as they play an important role in binding of SQV.  Overall, the major contributions come from a few groups around Asp25/Asp25ˊ, Ala28/Ala28ˊ, Ile50/Ile50ˊ, Pro81ˊ, Val/ Phe82ˊ, and Ile84/Ile84ˊ. The energy decomposition spectra showed that the binding mode of SQV to chain A is slightly different from that of chain B, this is due to the asymmetrical nature of SQV. The contributions of most residues to binding affinities were mainly due to van der Waals and electrostatic energy components.

Hydrogen bonds between SQV and protease: AE-WT versus mutants
Hydrogen bond analysis was performed on the trajectories from MD simulation to complement the energetic analysis. Hydrogen bond interactions between the AE-WT and mutant complexes are given in Tables 5 to 7. The binding of SQV to AE-WT is mainly mediated by hydrogen bonds from Asp25, Asp29, Asp30, and Gly48. Similar hydrogen bond interactions were observed between SQV and subtype B protease (F. Liu et al., 2008;Tie et al., 2007;Tzoupis et al., 2013). The catalytic aspartates Asp25/Asp25ˊin the active site region plays a key role in binding of the inhibitors. In AE-WT, the side chains of Asp25 directly interacted with central hydroxyl group of SQV. Table 5 shows that the interaction between Asp25 and SQV is mainly driven by electrostatic energy, whereas for Asp25ˊ, the interaction with SQV is favored only by van der Waals energy. The results in Table 5 show that the O4 hydroxyl group in hydroxyethylamine moiety of SQV formed hydrogen bond with carboxyl oxygen of Asp25 with occupancy of 99%. In G48V and V82F mutant complexes, both the aspartates Asp25/Asp25ˊformed hydrogen bond interaction with SQV. As seen from Tables 6 and 7, Asp25 of both the mutants formed hydrogen bond with almost a similar occupancy as observed in AE-WT and Asp25ˊformed hydrogen bond interactions with O4 atom of SQV with an occupancy of 28%. In V82F complex, Asp25ˊformed hydrogen bond with a higher occupancy of 99.8%. On comparing with subtype B wild-type protease complexed with SQV, both the aspartates of AE-V82F mutant structure show strong hydrogen bond interactions with a higher frequency. The results from per residue decomposition in Tables 2 to 4 show that the contribution of electrostatic interaction from Asp25 in binding of SQV is decreased by 3.09 kcal/mol and 8.35 kcal/ mol in G48V complex and V82F complex, respectively, compared to AE-WT. Interestingly, the V82F mutation induces conformational change of Asp25ˊside chains and positional deviation of SQV. These changes make the Asp25ˊto align in close proximity with the O4 atom, thereby favoring strong electrostatic interaction by -7.6 kcal/mol and -8.3 kcal/mol compared to AE-WT and G48V complexes, respectively.
The aspartates at position 29 and 30 in the active site also play a key role in binding of SQV in AE-WT and G48V   complexes. These interactions have been observed even in subtype B wild-type complexes. From Tables 5 and 6, it could be seen that the residue Asp29 formed hydrogen bond interaction with SQV with a higher occupancy of 62.5% compared to AE-WT, whereas the occupancy of hydrogen bond interaction between Asp30 and SQV seems to be slightly higher in AE-WT compared to G48V complexes. The per-residue decomposition analysis in Tables 2 and 3 shows  All energies are reported in kcal/mol. that the van der Waals energy between Asp29 and SQV directs the binding in both AE-WT and G48V complexes, and the contributions toward van der Waals energy were mainly derived from the side-chain atoms. The electrostatic contributions were observed to be unfavorable toward binding affinity. Moreover in AE-WT, the increase in van der Waals energy contribution by -0.23 kcal/mol is accompanied by an unfavorable polar solvational energy. In G48V complex, the conformational changes of the residue Asp29 and slight deviation of O1 atom of SQV favors van der Waals contribution. On analyzing the Asp30 residue, it was observed that the binding affinity is mainly driven by van der Waals contribution and polar solvational energy in AE-WT. In G48V complex, though there was favorable electrostatic contribution of -1.17 kcal/mol, it was counterbalanced due to unfavorable contribution of polar solvational energy thereby leading to decrease in total energy contribution by 0.13 kcal/mol. The G48V mutation causes a complete shift of residue Asp30 and large positional deviation of SQV, such that the side-chain atoms of Asp30 align near the P3 subsite of SQV, thereby resulting in strong electrostatic interactions. Moreover, in V82F complex, though both the aspartates Asp29 and Asp30 are conformationally similar as observed in AE-WT, the positional deviation of SQV makes the P3 subsite move far away from these aspartates, thereby resulting in unfavorable binding affinity.
The residue Gly48 located in the flap region is crucial for effective inhibition of SQV. Additionally, it has been reported in subtype B complexes that Gly48 also plays a role in shaping the active site and also stabilizes the enzyme-substrate complex. In AE-WT, we observed a similar interaction with occupancy of 15%. Several studies have provided insights into SQV resistance due to G48V mutation in subtype B protease (Aruksakunwong et al., 2006;F. Liu et al., 2008F. Liu et al., , 2008  All energies are reported in kcal/mol. All energies are reported in kcal/mol. at position 48, introduces a bulkier side chain into the binding pocket, thereby leading to resistance of SQV. In AE-protease, the substitution of Gly at position 48 to a more hydrophobic residue also results in loss of interaction with SQV. Though there is a loss of hydrogen bond interactions, the substitution of valine results in more favorable van der Waals interaction compared to electrostatic interaction. It is observed from Table 6 that the increase in van der Waals interaction by -1.79 kcal/mol in the G48V complex is mainly due to the positional deviation of SQV, which aligns P2 and P3 subsites close to the substituted valine. In the V82F complex, as seen in Table 7, Gly48 forms strong hydrogen bond interaction with P2 subsite of SQV with an occupancy of 82% and also with N2 atom with an occupancy of 32% compared to AE-WT. Studies by Tzoupis et al. (2013) show that this strong hydrogen bond interaction is observed only in subtype B wild-type and not in any other mutant structures. Though the interactions by aspartates (Asp25ˊ/Asp29) eventually stabilized the binding of SQV, the absence of interactions between the flap residues and SQV shows that the drug is not strongly anchored in the active site, which partially explains the decrease in binding affinity for SQV in subtype B protease. Interestingly, in subtype AE protease, the V82F mutation favors the formation of hydrogen bond interactions between the flap residue Gly48 and P2 subsite of SQV. The per-residue decomposition analysis also shows that the total energy contribution of Gly48 is increased by -2.06 kcal/mol compared to AE-WT. The structural analysis shows that residue Gly48 is completely shifted away in V82F complex and large positional deviation was observed in the P3 subsite. These structural changes make Gly48 in the V82F mutant to have strong electrostatic interaction and the electrostatic contribution is increased by -3.26 kcal/mol compared to AE-WT. Another important residue that plays a key role in the binding of SQV in V82F is Ile50/Ile50ˊ(I. T. Weber & Harrison, 1999). Interestingly, in the V82F complex, the water molecule (WAT 265) bridges by tetrahedrally coordinating the nitrogen atoms of Ile50 and Ile50ˊwith SQV. A similar interaction has been observed in the crystal structure of subtype B wild-type protease complexed with SQV (Tie et al., 2007). The V82F mutation induces a shift in the conformation of residue Ile50/Ile50ˊand, along with the positional deviation of SQV, favors the watermediated interaction. The van der Waals contribution from Ile50 in the binding of SQV seems to be more favorable than electrostatic contribution in AE-WT and G48V complexes. The side chains of residue Ile50ˊis aligned near O2 atom of SQV, such that it favors strong van der Waals and electrostatic contribution compared to mutant structures.

Contributions from individual residues to binding free energy: AE-WT versus mutants
Apart from the hydrogen bond and water-mediated interactions, the residues Gly27, Ala28/Ala28ˊ, Ile47/Ile47ˊ, Gly49/   Phe82ˊO 2.77 67.2 a The hydrogen bonds are determined by acceptor … donor atom distance of 3.5 Å. and acceptor … H-donor angle of ! 120 . b Occupancy in fraction to evaluate the stability and strength of the hydrogen bonds.
Gly49ˊ, Val82/Val82ˊ, Ile84/84ˊ, Arg8ˊ, Leu23ˊ, and Pro81f avor van der Waals interactions with SQV in AE-WT. These residues have been observed to contribute to the binding of SQV even in subtype B protease (Tie et al., 2007). In addition to hydrogen bond interactions, the van der Waals interactions also play a significant role in stabilizing the inhibitor and HIV-1 protease interactions.
The van der Waals energy between Gly27 and SQV directs the binding in AE-WT and G48V complexes, and the contributions toward van der Waals energy were mainly derived from the backbone atoms. The electrostatic contributions were observed to be unfavorable toward binding affinity. In V82F complex, the binding affinity of SQV with Gly27 was increased by -0.83 kcal/mol and -0.30 kcal/mol compared to AE-WT and G48V. Structural analyses of the V82F mutant showed that Gly27 lies in close proximity to quinoline ring at P3 subsite of SQV compared to AE-WT and G48V mutant structures, thus favoring the contributions from electrostatic and van der Waals energy.
The binding affinity of Ala28 with SQV in AE-WT and mutant structures is favored by both van der Waals and electrostatic contribution. In G48V complex, though the mutation induces only slight conformational changes in Ala28, a large positional deviation of SQV makes Ala28 to lie closer to central hydroxyl group of hydroxyethylamine moiety and O1 atom, thus favoring strong van der Waals contribution compared to AE-WT and V82F structures by -0.31 kcal/mol and -0.25 kcal/mol, respectively. For residue Ala28ˊ, the binding affinity with SQV with AE-WT and mutant structures is mainly due to van der Waals and polar solvational energy. The electrostatic contribution seems to be unfavorable. Both the mutant structures show positional deviations at P2ˊsubsite, but in V82F complex, the positional deviation of P2ˊof SQV favors van der Waals interactions with the back bone atoms of Ala28ˊthan AE-WT and G48V complexes by -0.68 kcal/mol and 0.39 kcal/mol, respectively.
The residue Ile47 is oriented toward P2 subsite of SQV favoring van der Waals interaction compared to electrostatic interactions in AE-WT and G48V complexes. Though conformational changes of Ile47 and positional deviation of SQV was observed in both the mutant structures, Ile47 is shifted far from the P2 subsite, thereby leading unfavorable electrostatic contribution in V82F complex. Also, the backbone van der Waals contribution is reduced in the V82F complex, thereby leading to a decrease in total binding affinity of Ile47 with SQV by $0.5 kcal/mol compared to AE-WT and G48V complexes.
The amino acid Gly49 showed less contribution to the binding affinity in AE-WT and G48V mutant complexes compared to V82F. It was interesting to note that the orientation of Gly49 toward SQV is different in each AE-WT and mutant structures. The two main reasons for such differences were the conformational changes of Gly49 compared to AE-WT and the large positional deviation of SQV in the mutant structures. In AE-WT, we see that the residue Gly49 is oriented toward quinoline ring at P3 subsite of SQV. Although there is a slight conformational change of Gly49 in both the structures, they seem to move away compared to the AE-WT complex and due to the positional deviation of SQV, Gly49 lies near P1 subsite and P2 subsite of SQV in G48V and V82F complexes, respectively. These changes led to an increase in distance between Gly49 and SQV, thereby leading to a loss in van der Waals and electrostatic contribution in AE-WT and G48V complexes. On the other hand, Gly49ˊoriented toward decahydroisoquinoline ring at P1ˊsubsite also undergoes conformational changes in mutant structures. In the G48V mutant structure, the decahydroisoquinoline ring has deviated, such that the distance of Gly49ˊis increased, thereby leading to a loss in van der Waals and electrostatic contribution by 0.77 kcal/mol and 0.95 kcal/ mol, respectively. Comparing AE-WT and V82F complexes, in the P1 subsite of SQV, slight conformational changes were observed, but the conformational changes of Gly49í n the V82F complex make it to move farther from P1 subsite thereby leading to loss of both van der Waals and electrostatic contribution by $ 0.5 kcal/mol.
Substitution of phenylalanine at position 82 results in an increase in van der Waals interaction with SQV by -0.58 kcal/mol and -0.47 kcal/mol, respectively, compared to AE-WT and G48V complexes. A possible reason for this could be loss in van der Waals interaction due to conformational changes and an increase in distance between Val82 and decahydroisoquinoline ring at P1ˊof SQV. The electrostatic contribution and polar solvational energy seem to be unfavorable for Val/Phe82 in both the chains A and B of AE-WT and mutant structures. The van der Waals energy between Val/Phe82ˊand SQV directs the binding in AE-WT, and mutant complexes and the contributions toward van der Waals energy were mainly derived from the side-chain atoms with an increase in total van der Waals contribution by -0.33 kcal/mol and -0.21 kcal/mol compared to AE-WT and G48V complexes, respectively. This is because Phe82ˊlies nearby between quinoline ring at P3 and P1 subsites, whereas Val82ˊin AE-WT and G48V complexes are aligned close to the P1 subsite. However, the favorable van der Waals contribution in the V82F complex is counterbalanced by more unfavorable polar solvational energy.
The residue Ile84 showed less van der Waals contribution in AE-WT and G48V complexes. The van der Waals contribution was decreased by 0.27 kcal/mol in AE-WT and by 0.18 kcal/mol in the G48V complex. There were small favorable electrostatic contributions in AE-WT, whereas in mutant complexes, the electrostatic contribution was unfavorable. On analyzing the V82F structure, conformational shift of alkyl groups in Ile84 and deviation of P1ˊsubsiteprobably led Ile84 in the V82F complex to align near the decahydroisoquinoline ring of SQV, favoring van der Waals contribution. The binding affinity of Ile84ˊin all the complexes is mainly driven by van der Waals contribution from backbone atoms and a slightly favorable electrostatic contribution, while in AE-WT, the electrostatic contribution was unfavorable. The conformational changes of Ile84ˊwere observed in all the mutant structures, but the positional deviation of SQV in the V82F complex makes the P1 subsite to align closer toward Ile84ˊ, thereby favoring van der Waals interaction by -0.77 kcal/mol and electrostatic interactions compared to AE-WT and G48V complexes.
The guanidinium moiety of Arg8ˊis oriented toward the quinoline ring at P3 subsite, thereby favoring strong van der Waals interaction in both AE-WT and G48V complexes. In the V82F complex, the conformational shift of Arg8ˊand positional deviation of SQV leads to an increase in distance between the guanidinium moiety of Arg8ˊand P3 subsite, thereby leading to a loss in the side chains van der Waals contribution by $1.00 kcal/mol. Also, the contribution of electrostatic interaction was unfavorable in the V82F complex. The isobutyl side chain of Leu23ˊis shifted away in mutant structures. This makes both the mutant structures to show less contribution of side-chain van der Waals interaction, whereas in AE-WT, the P1 subsite is aligned closer to Leu23ˊ, resulting in an increased binding affinity by -0.39 kcal/mol and -0.19 kcal/mol than G48V and V82F complexes, respectively. The P1 subsite of SQV is deviated farther leading to a tremendous decrease in binding affinity in the G48V complex.
The residue Pro81ˊlies close to the P1 subsite for all the complexes. The pyrrolidine group of Pro81ˊin the V82F complex was shifted closer to P1 subsite, thus resulting in favorable van der Waals interactions from both backbone and side-chain atoms, thus enhancing the total binding affinity of SQV by -0.41 kcal/mol and -0.24 kcal/mol in AE-WT and G48V complexes, respectively. The electrostatic interactions were favorable in all the three complexes. Although AE-WT shows more favorable electrostatic contributions by $ -0.20 kcal/mol than mutants, a large decrease in polar solvational energy and less contribution of van der Waals interaction makes Pro81ˊto show a decreased binding affinity than V82F complex. The binding of key residues with SQV in AE-WT and mutant structure are shown in Figure 10.

Impact of mutations on the binding affinity of SQV
For each subtype of HIV-1 protease, the molecular pathway to drug resistance differs and is yet to be completely characterized for non-B subtypes. Many studies have been carried out to understand the SQV resistance mechanism of mutations in subtype B (G48V and V82A) and subtype C proteases (V82A), very little is known about the resistance mechanism in AE-protease (Ahmed et al., 2013;Aruksakunwong et al., 2006;F. Liu et al., 2008;Saen-oon et al., 2007;Tie et al., 2007;Tzoupis et al., 2013). However, the effect of G48V and V82F mutations in the binding of SQV in AE-protease is still unknown. In this study, we performed MD simulations on AE-protease to understand how the G48V and V82F mutations confer resistance against SQV. The simulations indicate that each mutation caused a change in binding free energy compared to AE-WT. The G48V complex led to a decrease in binding free energy, whereas the binding free energy is increased for V82F. As observed by Hao et al. (2010), G48V, a high resistance mutation in subtype B resulted in a decrease in enthalpy (A-type mechanism) leading to resistance of SQV, whereas for V82A, a low-level resistance mutation the decrease in enthalpy was compensated by an increase in entropy (E-type mechanism). The V82A mutation in subtype C protease resulted in a decrease of both enthalpic and entropic contributions (C-type mechanism), leading to resistance of SQV (Ahmed et al., 2013). The resistance mechanism of G48V in AE protease is almost similar to subtype B protease, where the entropy is gained and the decrease in free energy of binding is due to less enthalpy contribution. The loss in enthalpy contribution is mainly from the residues Gly49/Gly49ˊ, Ile50/Ile50ˊ, Leu23ˊ, and Asp25ˊ. Contrary to the resistance mechanism of V82F mutation in subtype B and C proteases, the V82F mutation in AE-protease shows favorable binding affinity of SQV and the increase is mainly due to gain in enthalpy contribution. This gain in enthalpy is mainly due to residues Gly27, Gly49, Ile50, Asp25ˊ, and Ile84ˊ.

Influence of mutations on intramolecular hydrogen bond interactions
To clarify the effect of mutations, the intramolecular hydrogen bond interactions were analyzed in AE-WT and mutant protease structures. It has been reported in subtype B protease that a salt bridge interaction is present between Glu35-Arg57. It has been reported in various crystal structures that  E35D mutation removes the salt bridge. This kind of variation has been observed in different subtypes of HIV-1 protease. Moreover, substitutions at E35/R57 modulate the interaction between flap and hinge regions that can consequently affect the rigidity and conformational flexibility of the protein backbone (Huang et al., 2014;Z. Liu et al., 2016). A recent study has also shown that mutations in the hinge region at position 35 eliminates the ion pair interaction between Glu35-Arg57 (Kneller et al., 2020). In AE-WT, we observed that the hinge loop, where residue Asp35 is located, aligns in closer proximity to Arg57 at a distance favored to form a salt bridge. This Arg57/Arg57ˊ-Asp35/Asp35ˊsalt bridge was observed even in the V82F complex. But in the G48V complex, the Asp35/Asp35ˊ-Arg57/Arg57ˊsalt bridge was disrupted. This is because the conformation of Asp35 in the hinge loop is observed to move away from Arg57 in both the monomers compared to AE-WT and V82F complexes. This makes the Asp35/Asp35ˊ-Arg57/Arg57ˊdistance too large, unfavorable to form a salt bridge. Instead in the G48V complex, Asp35 formed hydrogen bond interaction with Leu33 and Asp35ˊformed interactions with Gly78ˊ, as seen in Table 6. Loss of Asp35/Asp35ˊ-Arg57/Arg57salt-bridge interaction may result in a more flexible mode of flap motion in the G48V complex, thereby affecting the binding of SQV.
It was interesting to observe that Thr80-Val82, highly conserved hydrogen bonds, important for stabilizing the P1 loop was retained in the V82F mutant structure. Usually, structural changes in the P1 loop have been observed when binding of inhibitors and drug resistance mutations at Val82 and Ile84 (Cameron et al., 1994;Foulkes et al., 2006). However, the substitution of phenylalanine at position 82 favored the formation of hydrogen bonds Thr80-Phe82 in both chain A and chain B, with an occupancy of 41.5% and 67.2%, respectively. In G48V complex, the Thr80-Val82 occupancy is reduced in both the monomers and in AE-WT the occupancy of Thr80-Val82 is reduced to 15%.

Analysis of intraflap and interflap interactions: AE-WT versus mutants
Hydrogen bonding analysis between SQV and the flap region reveals the significance of the flaps in the process of ligand binding. Hence, it is extremely important to analyze the hydrogen bonding interactions in the flap region. Analysis of hydrogen bonds in the flap region (residues 47/47ˊ-54/54ˊ) in AE-WT and mutant structures is given in Tables 8 to 10. It has been reported that hydrogen bonds formed at the tip of the flaps Gly49/Gly49ˊ-Gly52/Gly52ˊhelp in inhibitor association (Fern andez, 2005). A similar kind of interaction has been observed in each monomer for all the complexes. However, the substitution of Valine at position 48 reduces the occupancy of Gly49/Gly49ˊ-Gly52/Gly52ˊhydrogen bond interaction compared to AE-WT and V82F complexes. The residue Ile47 from both the monomers forms hydrogen bond interactions with the adjacent residue Ile54 in all the complexes. In the V82F complex, Ile47ˊ-Ile54ˊforms strong hydrogen bond interaction with a higher occupancy of 84.1%. It is evident from Tables 9 and 10 that both interflap and intraflap hydrogen bonding interactions are stabilized by both G48V and V82F complexes; however, there are no interflap hydrogen bond interactions formed in AE-WT, and it is stabilized only by intraflap hydrogen bond interactions. The interflap hydrogen bond interactions help in keeping the flaps in the closed conformation (Badaya & Sasidhar, 2020). In subtype B proteases, it has been observed that G48V mutation leads to structural changes in flaps thereby abolishing the interflap interactions between Ile50-Ile50ˊ (Tzoupis et al., 2013). Contrary to subtype B, in AE protease, we observed that the interflap interactions between Ile50-Ile50ˊis present only in G48V complexes. In the V82F complex, the interflap interactions were formed between Gly51 and Ile50ˊ. However, the interaction of Gly51 with Ile50ˊis of importance because they indirectly help in holding SQV in the binding pocket through a water molecule in the V82F complex.

Remarks on binding affinity
The MD simulation results show that there is high structural fluctuation involving the residues Gly16 (fulcrum), Asp35/35-Asn37/37ˊ(flap elbow) in G48V complex, and high fluctuations were observed for the residues in the flap region (44-55)/(49ˊ-55ˊ) for both AE-WT and G48V complexes. The residue Gly48 lies close to the flap tip and a mutation in this position induces conformational changes of Ile50/Ile50. These conformational changes could probably stretch the flap tip-active site distances in the G48V complex, probably making the active site volume a little wider. As a result, the increased flexibility in the flaps and active site region of the mutants leads to a decrease in binding affinity of the enzyme for SQV. On the other hand, in the V82F mutant structure, the residues Ile50ˊand Asp25 lie in close contact with the 80's loop. Mutation in this 80's loop region induces conformational changes for the residues Asp25/Asp25ˊand Ile50/Ile50ˊ. These conformational changes probably induce changes in Ile50C a /Ile50ˊC a -Asp25C a /Asp25ˊC a distances, smaller in volume, thereby enhancing the binding affinity of SQV. The energy decomposition analysis suggests that the decrease in binding affinity of the G48V mutant was mainly attributed to the loss in enthalpy component. In the V82F mutant, an increase in binding affinity was due to favorable enthalpy components. The results are in agreement with the experimental studies reported by Clemente et al. (2006) that SQV shows effective inhibition against the V82F mutant.

Conclusion
The molecular basis for SQV resistance due to active site mutations G48V and V82F shows that the G48V mutation reduces the binding affinity of SQV compared to AE-WT. The V82F mutation shows increased binding affinity due to favorable enthalpy contribution. An extensive comparison of dynamics and energetics for all the systems indicated the following reasons for SQV resistance in the G48V complex. As the residue Gly48 is directly responsible for entrapment of SQV into the active site, substitution of bulkier side chain impairs the interaction with SQV. Moreover, the G48V mutation also influences the intramolecular hydrogen bonding network of the neighboring residues, thereby leading to reduced stability and altered interactions, indirectly affecting the binding of SQV. On the other hand, in V82F complex, the substitution of phenylalanine at position 82 induces conformational changes in the 80's loop, flaps, and active site cavity, thereby resulting in an enhanced binding affinity with SQV. The intramolecular hydrogen bond analysis shows that in the V82F mutant salt bridge and P1 loop is stabilized similar to AE-WT. The findings from the present study indicated that the preservation of (i) hydrogen bonds between Asp25/ Asp25ˊ, Gly48, and SQV and (ii) water-mediated interaction between Ile50/Ile50ˊand SQV is crucial for effective inhibition.

Disclosure statement
No potential conflict of interest was reported by the authors.