Understanding The Mechanism of Mutant HIV-1 Protease Resistance Against Darunavir Using Molecular Dynamic Study


 The human immunodeficiency virus type 1 protease (HIV-1 PR) is an important enzyme in life cycle of the HIV virus. It cleaves inactive pre-proteins of the virus and changes them into active proteins. Darunavir suppresses the wild type HIV-1 PR (WT-Pr) activity, but can’t inhibit the mutant resistant forms (MUT-Pr). Increasing knowledge about the resistance mechanism can be helpful for designing of more effective inhibitors. In this study, the mechanism of resistance of Ile47val and Ile54Met MUT-Pr strain against Darunavir was investigated. For this purpose, complexes of Darunavir with WT-Pr (WT-Pr-D) and MUT-Pr (MUT-Pr-D) were simulated for 200 ns and structure, dynamic and energetic properties of both simulations were investigated based on essential dynamics (principal component analysis (PCA)), root mean square fluctuation (RMSF), radial distribution function (RDF), molecular mechanics/Poisson Boltzman surface area (MM/PBSA) energies and etc. Our data revealed that mutations increased the flap tips flexibility and increased the active-site space, probably due to the reduction in hydrophobic forces. So, the protease’s conformation changed from closed state to semi-open state. Formation of semi open structure along with a reduction in van der Waals interactions and hydrogen bonds with Darunavir facilitates disjunction of Darunavir from the active-site in MUT-Pr-D.


Introduction
The human immunode ciency virus type 1 protease (HIV-1 PR) is an enzyme that cleaves the HIV polyproteins (the Gag and Gag-pol), helping the virus to reach maturity at the last stage of its life cycle 1 . So, the enzyme inhibitors may play an important role in the battle with the acquired immunode ciency syndrome (AIDS) disease 2 . However, the strains with a mutant protease demonstrate a great resistance against inhibitors 3 . Therefore, understanding the molecular mechanisms underlying the resistance of HIV-1 PR antiviral drugs is crucial for the design of highly potent inhibitors against the drug-resistant strains.
HIV-1 PR is a member of aspartyl protease family, which contains the aspartate residue in its active-site 4 .
Use of drugs against the HIV-1 PR leads to increased diversity of its mutant strains 38 . Recently ~ 50 mutations have been discovered at 30 different sites in the HIV-1 PR gene 39 . The mutations fall into two types: (1) mutations that occur in the active-site, which directly reduce the drug-protease interactions; (2) mutations that occur distant from the active-site, which indirectly reduce the HIV-1 PR a nity to drug by affecting the conformational dynamics of the enzyme 36,40 . However, some mutations may exert both effects.
This study was devoted to compare the HIV-1 PR structure, dynamics and interactions with Darunavir in mutant and wild-type strains. Our data revealed that mutations could increase the exibility of the aptips, make them separated and rotated with respect to each other and change active-site space, which probably help to conversion of protease's conformation from closed state to semi-open and curled states that consequently facilitate disjunction of Darunavir from the active-site, leading to the protease resistance against the Darunavir. The mechanistic knowledge from this study may provide clues for designing new inhibitors against the HIV-1 protease.  44,45 . Darunavir charges were obtained by using Restrained Electrostatic Potential (RESP) method 46 . General Amber force eld (GAFF) 47 parameters and the RESP charges were determined for Darunavir using the Antechamber module in the AMBER20 package 48 . All missing hydrogen atoms were added using the LEaP module. The WT-Pr-D and MUT-Pr-D systems were solvated with the TIP3P 49 water model in the periodic boxes of size 7.71, 7.71, 5.45 (x, y, z: nm) and 8.14, 8.14, 5.75 (nm), respectively. Distance between water molecules and box surface was determined as 10 Å. Then, chloride ions were added in order to neutralize the positive charges of systems. For WT complex system, 9697 TIP3P water molecules and 7 Cl ions, and for the mutant complex system, 11533 water molecules and 3 Cl ions were added to the simulation boxes. Then, topology les for the WT-Pr and MUT-Pr systems were built by GROMACS version 2019 and the ff99SB force eld 50 .

Molecular dynamics simulation
At rst step, the WT-Pr-D and MUT-Pr-D in water box were minimized by using of steepest descent minimization algorithm for 150000 steps. Constant temperature and pressure conditions were applied for both simulations 51,52 . The systems were equilibrated at the NVT ensemble for 150000 steps by using leap-frog integrator by time step as 2 fs and temperature coupling was set by V-rescale algorithm. At last, temperature of systems was equilibrated at 310 K and followed by NPT ensemble equilibrating for another 150000 steps. The Pressure coupling and its type were determined to be Berendsen and isotropic. The pressure of systems was equilibrated at 1 bar. For all energy minimization, NVT ensemble, NPT ensemble and MD production steps, neighbor list search was performed by using grid with 1.2 nm as cut-off. Lennard-Jones interactions cut-off were determined 1.2 nm. Long-range electrostatic interactions were calculated by particle mesh ewald (PME) method 53 . Short range electrostatic interactions cut-off was determined as 1. Comparison of the exibility of the ap-tips In order to evaluate the residual exibility, the root mean square uctuations (RMSF) of C α atoms were calculated for chain-A (Fig. 1a) and chain-B (Fig. 1b) of WT-Pr-D and MUT-Pr-D. The average RMSF and standard deviation (SD) values for ap-tip and active-site binding region of the MUT-Pr-D and WT-Pr-D were 0.97 ± 0.41 Å and 0.83 ± 0.21Å, respectively. As seen, mutations increased the exibility of residues that are around Darunavir, but the residues Asp25-Glsdy26-Thr27(chain-A/B) from catalytic site remained rigid, which is in accordance with the previous experimental and theoretical studies [59][60][61][62] . Figure 1c and 1d show the RMSF difference between the MUT-Pr-D and WT-Pr-D for chain A and B, respectively. Residues with an RMSF difference more than 0.50 Å were considered to be highly uctuating, which represent signi cant mutation-induced conformational changes. As seen, in MUT-Pr-D, there are signi cant decrease in the exibility of residues Lys14-Arg20 (fulcrum of chain-A), Ile36-Leu38 ( ap-elbow of chain-A) and Val63-Ile64 (cantilever of chain-A) (Fig. 1c). These changes can help forming of semi-open conformation in mutant strain as mentioned in Sect. 3-3. In contrast, the exibility of the ap-tip region in chain-B (residues Gly51-Phe53) of MUT-Pr-D was remarkably increased compared with that of WT-Pr-D. It has been shown that increasing the ap-tip exibility facilitate opening of the active-site gate, which consequently facilitate releasing the inhibitor from active-site 10,12,63 . The highly affected residues as a result of mutations in both WT-Pr and MUT-Pr complexes were illustrated by superposed 20 trajectory structures, with 10 ns intervals ( Fig. 1e and 1f). According to these data, it can be suggested that the ap-elbow, fulcrum and cantilever regions in chain-A and the ap-tip region in chain-B experience the most conformational changes as a result of mutations.

PCA analysis
In order to extract essential dynamic motions, PCA analyses were performed 54 . The rst principal components (pc1) was calculated by projection of both trajectories on their rst eigenvector (ev1). The PCA of trajectories during 200 ns from 20000 frames is illustrated in Fig. 2, in which the cones length represents amplitude of movements and their direction represent the direction of movements.
As seen in Fig. 2, the ap-tips, cantilever (chain-A) and fulcrum (chain-A) regions of WT-Pr-D and the aptips, ap-elbow (chain-B) and interface regions of MUT-Pr-D have larger essential dynamics than other parts that is consistent with the RMSF data ( Fig. 1a and 1b). Both RMSF and PCA data showed that cantilever (chain-A) and fulcrum (chain-A) in WT-Pr-D have a higher dynamic with respect to MUT-Pr-D. As seen in Fig. 2, the movements of these regions in WT-Pr-D are directed toward the protease active-site. However, in MUT-Pr-D these movements are negligible and directed outward of the active-site (Fig. 2).
Also, in WT-Pr-D, ap-elbow (chain-A) that was more exible than that of MUT-Pr-D showed upward movement directed to closed conformation. However, in MUT-Pr-D it showed negligible movement (Fig. 2). These data showed that in wild type complex, regions surrounding active-site made its closed conformation more stable. Also, ap-tips residues in WT-Pr-D move closer to each other and to active-site (Fig. 2

Structural analysis
In order to investigate the structural differences between WT-Pr-D and MUT-Pr-D complexes, the radial distribution function (RDF) and ap-tip distance (between C α of residues Ile50 (chain-A)-Ile50(chain-B) ( Fig. 3) were examined in both complexes. The RDF and vertical distance between C α of residues Asp25(chain-A)-Ile50(chain-A) and between C α of residues Asp25(chain-B)-Ile50(chain-B) (Fig. 3) were also investigated as a greet metric to estimate the extent of ap opening 12,69 .

Flap tips to active-site distance measurements
The intra-chain distances between Ile50 ( ap-tip) and Asp25 (active-site) were measured during the 200 ns MD simulations (see Fig. 3). The mean ± SD distances between Asp25 and Ile50 in chain-A were 14.22 ± 0.55 and 13.94 ± 0.43 Å for MUT-Pr-D and WT-Pr-D, respectively. The most probable distances (RDF) were 14.0-14.5 Å and 13.5-14.0 Å for MUT-Pr-D and WT-Pr-D, respectively (Fig. 4a). So, in chain-A of the MUT-Pr-D, the average distance between the ap-tip and active-site was increased by ~ 0.28 Å as a result of mutation. This increase was even greater in chain-B. The mean ± SD distance values between Asp25 and Ile50 in chain-B were 16.07 ± 0.95 and 15.31 ± 0.72 Å for MUT-Pr-D and WT-Pr-D, respectively. The most probable distances were 15.8-16.5 Å and 14.5-15.5 Å for MUT-Pr-D and WT-Pr-D, respectively (Fig. 4b). Thus, the average ap-tip to active-site distance in chain-B of the MUT-Pr-D was increased averagely by ~ 0.76 Å. As seen, the upward movements of the ap-tips were increased as a result of mutations that can help to conversion of closed conformation to semi-open conformation and releasing inhibitor from its active-site binding region 12,70 . Flap-tip to ap-tip distance measurements The Average ± SD distance between ap-tips of chain-A and chain-B were 6.20 ± 0.71 and 7.62 ± 0.73 Å in MUT-Pr-D and WT-Pr-D, respectively. The most probable distances were 5.8-6.5 Å and 7.5-8 Å for MUT-Pr-D and WT-Pr-D, respectively (Fig. 4c). Based on previous reports, the ap-tips separation more than10 Å was considered as open conformation, from 6 to 10 Å as semi-open conformation and the separation of less than 6 Å as the closed structure 8,10,12 . As seen, the ap-tip to ap-tip average distance was about 1.42 Å smaller in MUT-Pr-D. Because Ile50 in chains-A and B were located in the cap of ap-tips, when the ap-tips got far away from each other along the longitudinal axis, these residues become closer to each other (Fig. 3). However, in our nano second duration simulation, ap-tips conversion from closed conformation to semi-open have not been shown directly because it acquires micro second to millisecond time to occur as NMR experiments have revealed [22][23][24].

Hydrophobic interactions stabilize the closed conformation of ap-tips
In ligand-bonded WT-Pr, hydrophobic interactions between Ile50(chain-A) and its adjacent residues Ile47(chain-B)/Ile54(chain-B) as well as hydrophobic interactions between Ile50(chain-B) and its adjacent residues Ile47(chain-A)/Ile54(chain-A) (Fig. 5) play an important role in retention of the active-site in closed conformation, which consequently traps Darunavir in the active-site [71][72][73] . Previous experimental analysis on the extensive statistical community that suffer from AIDS disease have revealed that Ile47val and Ile54Met mutations are most prevalent in resistant strains to Darunavir 74 . So, this type of Mut-Pr was selected and the hydrophobic interactions during 200 ns MD simulations were compared in MUT-Pr-D and WT-Pr-D. For this purpose, the Radial Distribution Function and distance between Cα atoms of Ile50(chain-B) and Ile47(chain-A)/ Ile54(chain-A) (Fig. S3 and Fig. 5) and distance between Ile50(chain-A) with Ile47(chain-B)/ Ile54(chain-B) in WT-Pr and MUT-Pr-D, in which Ile47 was substituted by Val and Ile54 was substituted by Met, were calculated ( Fig. S4 and Fig. 5) The average distance ± SD between Ile50(B) and Val47(A) in MUT-Pr-D and between Ile50(chain-B) and Ile47(chain-A) in WT-Pr-D were calculated 9.25 ± 0.88 Å and 8.82 ± 0.51 Å, respectively. The most probable distances were calculated as 8.1-10 Å and 8.3-9.2 Å for MUT-Pr-D and WT-Pr-D, respectively ( Fig. S3-a). As seen, in MUT-Pr-D, the average distance was increased ~ 0.43 Å compared to that of WT-Pr-D. The Average Distance ± SD between Ile50(chain-B) and Met54(chain-A) in MUT-Pr-D and between Ile50(chain-B) and Ile54(chain-A) in WT-Pr-D were calculated to be 7.02 ± 0.55 and 7.48 ± 0.63 Å, respectively. The most probable distances were calculated as 6.5-7.2 Å and 7.3-7.8 Å for MUT-Pr-D and WT-Pr-D, respectively. So, in MUT-Pr-D, the average distance decreases about 0.46 Å relative to WT-Pr-D (Fig. S3-b).
The Average distance ± SD between Ile50(chain-A) and Val47(chain-B) in MUT-Pr-D and between Ile50(chain-A) and Ile47(chain-B) in WT-Pr-D were calculated to be 9.29 ± 0.91 and 8.10 ± 0.53 Å, respectively. The most probable distances were calculated as 8.3-8.7,9.9-10.3 Å and 7.7-8.2 Å for MUT-Pr-D and WT-Pr-D, respectively ( Fig. S4-a), showing that average distance increased about 1.19 Å with respect to WT-Pr-D due to mutations. The Average distance ± SD between Ile50(chain-A) and Met54(chain-B) in MUT-Pr-D and between Ile50(chain-A) and Ile54(chain-B) in WT-Pr-D were calculated to be 8.40 ± 0.75 Å and 7.41 ± 0.35 Å, respectively. The most probable distances were calculated as 8-8.5 Å and 7.2-7.6 Å for MUT-Pr-D and WT-Pr-D, respectively. So, this distance increased about 1 Å respect to WT-Pr-D due to mutations (Fig. S4-b). however, distance between Ile50(chain-B) and Met54(chain-A) in MUT-Pr-D was decreased but this can't increase hydrophobic interaction between ap-tips because MET isn't a hydrophobic residue (Fig. 5 and Fig. S3-b). Indeed, except lle50 (chain-B) and Met54(chain-A), distances between other mentioned three-pair hydrophobic residues increased in MUT-Pr-D compared to that in WT-Pr-D. So, the ap-tips separation relative to each other occurs due to the reduction of hydrophobic interactions that consequently forms semi-open conformation (Fig. 5) 72 . As a result, binding site gateway could open more easy in MUT-Pr-D than WT-Pr-D [10,12].

Comparing the interactions in WT-Pr-D and MUT-Pr-D systems
In order to compare the binding strength between Darunavir and WT-Pr or MUT-Pr, the binding energies were calculated by MM-PBSA approach 75 (see data in Table 1). According to the data presented in Table 1, the part of binding energy related to electrostatic interactions (ΔE electrostatic ) were almost the same in both complexes. the Polar solvation energy (∆G Polar−solvation ) and non-polar solvation energy (∆G SASA ) showed little changes, so that in the WT-Pr-D, polar interactions were a little more unstable and non-polar interactions were a little more stable than that in MUT-Pr-D. Thus, these interactions had same effect on total binding energy in both complexes too. However, van der Waals interactions between the protease and Darunavir have been decreased in MUT-Pr-D. So, the protease-Darunavir binding energy in MUT-Pr-D was smaller than that of WT-Pr-D by about 20 kJ/mol (Table 1). This indicates that the mutations have reduced the binding strength between the protease and Darunavir in the MUT-Pr-D system by decreasing van der waals interactions. Table 1 The binding energy components between the protease and Darunavir in WT-Pr-D and MUT-Pr-D systems.

Radius of gyration
In order to assess the compactness of the protease active-site cavity, the radius of gyration (g(r)) of backbone atoms was calculated for residues 23-32, 46-54 and 78-87 in both chains (Fig. S5). As seen, the average ± SD g(r) values during rst 110 ns simulation time for MUT-Pr-D and WT-Pr-D were 10.50 ± 0.09 Å and 10.62 ± 0.10 Å respectively. The average g(r) in WT-Pr-D was ~ 0.12 Å larger than that of MUT-Pr-D. So, during the rst 110 ns, the active-site of MUT-Pr-D was more compact than of WT-Pr-D. However, after 110 ns, the average g(r) and standard deviation for MUT-Pr-D and WT-Pr-D were calculated as 10.54 ± 0.07 Å and 10.52 ± 0.06 Å, respectively. Thus, the compactness of the MUT-Pr-D active-site decreased a little, while in WT-Pr-D it became more compact.
So, at last 90 ns, respect to the beginning of simulation, the active-site cavity in the WT-Pr-D was decreased about 0.12 Å but it was increased about 0.04 Å in MUT-Pr-D. As a result, compactness of WT-Pr-D become increased but compactness of MUT-Pr-D become a little decreased these results are in agreement with other structural and dynamical analysis ( Fig. S2-a).

Conclusion
The principal components and uctuation analyses showed that in MUT-Pr-D, the ap-tips get away from each other, curl upwardly and become more exible in chain-B. Also, movement correlation between Flaptip(chain-B) and Flap-Elbow(chain-B) in MUT-Pr-D was seen. All these evidences con rm that as a result of mutation, the formation of semi-open conformation become more probable and the closed conformation become unstable. Based on the distance and RDF analysis, it was also shown that the average distance of ap-tips from each other and from the active-site increased in MUT-Pr-D, leading to conversion from closed to a semi-open conformation more probably. Increase of ap-tip(chain-A) to aptip (chain-B) distance occurred as a result of decreasing hydrophobic interactions in MUT-Pr-D. PCA and g(r) analysis showed that the volume of active-site in MUT-Pr-D becomes a little larger. This is because in WT-Pr-D, the fulcrum, cantilever and ap-tips move toward the active-site core, while in MUT-Pr-D, these movements were slightly outward from the core of active-site. Thus, the MUT-Pr-D conformation become near to semi-open conformation. At last, hydrogen bonds analysis indicated that the hydrogen bonds between Darunavir and the active-site and catalytic residues of protease were weaker in MUT-Pr-D.
Binding energy calculations also approved that binding strength between Darunavir and MUT-Pr decrease as a result of mutations. Brie y, it can be concluded that resistance of MUT-Pr to Darunavir occurred due to changing HIV-1 PR structure and dynamic as a result of mutation.