Simulation stability
The root mean square deviation (RMSD) for backbone atoms of proteases were analyzed in wild-type and mutant complexes to find out if the structures reached their global minimum and stable state in the 200 ns simulation time. The RMSD of 200 trajectory structures with 1 ns interval were analyzed (Fig. S2). The Mean RMSD and SD for backbone atoms were calculated as 0.12 ± 0.01 and 0.11 ± 0.01 for the WT-Pr-D and MUT-Pr-D, respectively. As seen, both of the mutant and wild type proteases immediately reached their stable states at RMSD value of almost 0.1 nm and RMSD fluctuations were small during simulation. The Darunavir/MUT-Pr-D complex showed a big structural deviation in 110 ns and then reached a steady state at a mean RMSD of 0.127 ± 0.01 nm with small fluctuations. Darunavir/WT-Pr-D complex also reached a stable state at RMSD of 0.126 ± 0.009 nm. It also showed big structural deviation in earlier steps of simulation (see Fig. S2). These data indicated that both complexes reached their stable states during simulations.
Comparison of the flexibility of the flap-tips
In order to evaluate the residual flexibility, the root mean square fluctuations (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 flap-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 flexibility 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–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 fluctuating, which represent significant mutation-induced conformational changes. As seen, in MUT-Pr-D, there are significant decrease in the flexibility of residues Lys14-Arg20 (fulcrum of chain-A), Ile36-Leu38 (flap-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 flexibility of the flap-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 flap-tip flexibility 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 flap-elbow, fulcrum and cantilever regions in chain-A and the flap-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 first principal components (pc1) was calculated by projection of both trajectories on their first 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 flap-tips, cantilever (chain-A) and fulcrum (chain-A) regions of WT-Pr-D and the flap-tips, flap-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, flap-elbow (chain-A) that was more flexible 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, flap-tips residues in WT-Pr-D move closer to each other and to active-site (Fig. 2), but in MUT-Pr-D the flap-tip and flap-elbow of chain-B got far away from flap-tip of chain-A, downwardly on longitudinal axis (Fig. 2). The correlation between this flap-elbow and flap-tip downward movement lead to increasing of flap-tips opening and facilitate change to semi-open and open-like conformation 10,64−67. Indeed, flap-tips in MUT-Pr-D show the rotation that change their conformation from two parallel planes to curled planes described as curled conformation, facilitating formation of open-like conformation (Fig. 2) 6,25,68. These results indicate that the gate of active-site in MUT-Pr-D becomes susceptible to formation of open-like states compared to WT-Pr-D. In addition, the interface in MUT-Pr-D move upwardly to flap-tips and looks like it can push Darunavir outside of active-site and may help flap-tips to open (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 flap-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 flap opening 12,69.
Flap tips to active-site distance measurements
The intra-chain distances between Ile50 (flap-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 flap-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 flap-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 flap-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 flap-tip distance measurements
The Average ± SD distance between flap-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 flap-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 flap-tip to flap-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 flap-tips, when the flap-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, flap-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–24].
Hydrophobic interactions stabilize the closed conformation of flap-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–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 flap-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 flap-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 (ΔEelectrostatic) were almost the same in both complexes. the Polar solvation energy (∆GPolar−solvation) and non-polar solvation energy (∆GSASA) 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.
System
|
∆Evdw
(Kcal/mol)
|
∆Eelectrostatic (Kcal/mol)
|
∆GPolar−solvation (Kcal/mol)
|
∆GSASA
(Kcal/mol)
|
∆GMMPBSA (Kcal/mol)
|
WI-Pr-D
|
-280.71+/-26.0
|
-88.52+/-14.9
|
252.40+/-17.3
|
-26.258 +/- 0.723
|
-143.08 +/- 22.4
|
MUT-Pr-D
|
-260.91+/-16.3
|
-88.82+/-27.7
|
250.59+-30.1
|
-24.44+/- 0.714
|
123.583+/-18.2
|
All of the energies are in Kcal/mol. Total binding energies (∆GMMPBSA: Molecular mechanics/Poisson-Boltzmann Surface Area) obtained from sum of the ∆Evdw (van der Waals interaction energy), ∆Eelectrostatic (Electrostatic interaction energy), ∆GPolar−solvation (Polar solvation energy) and ∆GSASA (Non-Polar solvation energy) energies. |
Hydrogen bonds between Darunavir and HIV-1 proteases
Hydrogen bonds pattern between Darunavir and HIV-1 PR are important and can affect inhibition 16,76,77. In order to investigate the effect of mutation on these interactions, hydrogen bonds between Darunavir and HIV-1 PR from all 20000 frames in both trajectories were analyzed by Hbonanza 78 program. Figure 6 shows the hydrogen bonds data, shown by color gradient, in which white color dashed lines represent the hydrogen bonds appearing at least in 10 percent of whole simulation time and green color dashed lines represent those appearing in 100 percent of entire 200 ns simulation time. Only hydrogen bonds in which the distance between hydrogen donor and acceptor atoms were 3.5Å and angle of Donor-H-acceptor atoms were between 0–30° were considered. Figure 6 shows the comparison of hydrogen bond pattern in WT-Pr-D with MUT-Pr-D. As seen, the Asp25(chain-A/B) are very important residues in the catalytic site, locating near each other, forming strong hydrogen bond interactions with each other in whole simulation time in both WT-Pr-D and MUT-Pr-D. Both sidechain oxygen atoms of Asp25(chain-A/OD1, OD2) are involved in interaction with two Darunavir’s atoms (O18, N20) in WT-Pr-D by hydrogen bond interactions, while in MUT-Pr-D only one oxygen atom (OD2) of Asp25(chain-A) is involved in hydrogen bond interaction with only one Darunavir atom (O18). This is because of Darunavir rotation in the active-site of MUT-Pr-D (by 180°) with respect to WT-Pr-D. Also, in WT-Pr-D, hydrogen bond between N21 atom of Darunavir and Asp30 (chain-B) in backbone of protease was formed in whole trajectory. However, in MUT-Pr-D, hydrogen bond between N21 atom of Darunavir and backbone of Asp30(chain-A) was formed almost half of trajectory’s frames. These results indicate that in WT-Pr-D, strongest hydrogen bond interactions are involved in binding of Darunavir to the active-site and catalytic residues.
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 first 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 first 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).