Molecular Dynamics Investigations of Binding Mechanism for Triazoles Inhibitors to CYP51

The sterol 14α demethylase enzyme (CYP51) is an important target of fungal infections. However, the molecular mechanism between triazoles inhibitors and CYP51 remains obscure. In this study, we have investigated the binding mechanism and tunnel characteristic upon four triazoles inhibitors with CYP51 based on the molecular docking and molecular dynamics simulations. The results indicate the four inhibitors stabilize in the binding cavity of CYP51 in a similar binding mode. We discover a hydrophobic cavity (F58, Y64, Y118, L121, Y132, L376, S378, S506, S507, and M508) and the hydrophobic interaction is the main driving force for inhibitors binding to CYP51. The long-tailed inhibitors (posaconazole and itraconazole) have stronger binding affinities than short-tailed inhibitors (fluconazole and voriconazole) because long-tailed inhibitors can form more hydrophobic interactions with CYP51. The tunnel 2f is the predominant pathway for inhibitors ingress/egress protein, which is similar to the other works of CYP51. This study could provide the theoretical basis for the development of efficient azoles inhibitors and may lead a better insight into structure–function relationships of CYP51.

GRAPHICAL ABSTRACT

Introduction

Life-threatening infections caused by fungi have increased rapidly, especially for patients that have immunocompromised diseases, such as AIDS, cancer, and organ-transplantation (Bongomin et al., 2017; Lee et al., 2020). It is reported that fungal diseases kill 1.5 million people per year, whose number is almost equal to the death of tuberculosis, and nearly three times as that of malaria (Brown et al., 2012; Bongomin et al., 2017). Among all kinds of fungal pathogens, Candida albicans is the most general fungi, which leads to candidemia (Antinori et al., 2016; Lee et al., 2020). To deal with the serious effects of life-threatening fungal infections, the development of antifungal agents has become a widespread concern (Lepesheva et al., 2018). In fungal cell membrane, the sterol 14α demethylase enzyme (CYP51) is responsible for catalyzing the lanosterol 14α methylation to produce ergosterol which can regulate the integrity, fluidity and permeability of the cell membrane (Balding et al., 2008; Lass-Floerl, 2011; Ene et al., 2012; Hargrove et al., 2016). Thus, influencing the growth and replication of fungi by inhibiting CYP51 has become a strategy for the development of antifungal agents (Lepesheva et al., 2008; Choi and Roush, 2017).

As we all know, there are four clinical triazoles agents targeting CYP51 for the treatment of systemic fungal infections: fluconazole (Flu), itraconazole (Itc), voriconazole (Vor), and posaconazole (Pos) (Lepesheva et al., 2018; Lee et al., 2020). The structures of four inhibitors are shown in Figure 1, they can be divided into short-tailed agents (ST: Flu and Vor) and long-tailed agents (LT: Itc and Pos) (Keniya et al., 2018). As the drug resistance mutations in CYP51 of Candida albicans increased, the effectiveness of four existing inhibitors is limited (Warrilow et al., 2019; Nishimoto et al., 2020). Therefore, elaborating the molecular mechanism of the existing drugs is very positive for further design and development of new drugs.

FIGURE 1

Figure 1. The structures of (A) fluconazole, (B) voriconazole, (C) itraconazole, and (D) posaconazole.

The successful resolution of the crystal structure of Candida albicans CYP51 by Lepesheva group inspired us to explore the molecular mechanisms between inhibitors and CYP51 (Hargrove et al., 2017). Molecular dynamics (MD) simulations are widely used to study the molecular mechanisms of inhibitors and enzyme (Watson et al., 2017; Sun et al., 2018; Xiao et al., 2020). Thus, we employed MD simulations and molecular docking to explore the binding mechanism between inhibitors and CYP51 in the present study. The results might offer insights into the structure–function relationships of CYP51 and provide the molecular basis for the rational design of new azoles inhibitors.

Materials and Methods

Preparation of Molecular Systems

The three-dimensional structure of Candida albicans CYP51 enzyme was obtained from the Protein Data Bank (PDB code: 5FSA) (Hargrove et al., 2017). Moreover, the structures of the four inhibitors: Flu (Compound CID: 3365), Vor (Compound CID: 71616), Itc (Compound CID: 55283), and Pos (Compound CID: 468595) were obtained from the PubChem database. After removing the ligand from the complex, the apo protein (Supplementary Figure S1) was saved by the Discovery Studio 3.1 (Studio, 2011). The CDOCKER protocol of Discovery Studio 3.1 (Studio, 2011) was employed to build the four complex structures, including Flu-CYP51 complex, Vor-CYP51 complex, Itc-CYP51 complex, and Pos-CYP51 complex.

Molecular Docking

CDOCKER (Studio, 2011) is a grid-based molecular docking method by CHARMm force field. The geometry optimization of four ligands was performed by Gaussian 09 (Frisch et al., 2009) at the B3LYP/6-31G (d) level. The CHARMm force field was applied to ligands and receptors. The receptor was held rigid, whereas ligands were allowed to flex during the refinement. The input site sphere was defined as a radius of 12 and 18 Å for short-tailed inhibitors (Flu and Vor) and long-tailed inhibitors (Itc and Pos), respectively. The other docking information was set as the default value. A conformational search of the ligands was performed by a grid-based simulated annealing method. The ligands were firstly heated to 700 K (2000 steps) and then annealed to 300 K (5000 steps). The value of the grid extension was set as 8 Å. The top 20 poses of each complex were saved for comparison and analysis. Finally, combined with the information of binding site in the literature (Nair et al., 2016; Hargrove et al., 2017) and docking score, the bested conformation of each system was chosen as the initial structure for the subsequent MD simulations.

Molecular Dynamics Simulations

For each system, geometry optimization of four inhibitors was performed by Gaussian 09 (Frisch et al., 2009) with the ab initio calculation method at the B3LYP/6-31G (d) level (Lee et al., 1988; Andersson and Uvdal, 2005). For charge derivation, the restrained electrostatic potential (RESP) fitting procedure was used (Bayly et al., 1993). The force field parameters of inhibitors were supplied by the general AMBER force field (GAFF) in the Antechamber module of AMBER 16 package (Wang et al., 2004, 2006; Case et al., 2016). The force field parameters developed by Shahrokh et al. were assigned to heme (Shahrokh et al., 2012). Finally, all missing atoms and hydrogen atoms were added using the t-leap procedure in the AMBER 16 package (Case et al., 2016). MD simulations were performed by the AMBER 16 package using the ff14SB force field (Maier et al., 2015; Case et al., 2016). To ensure the overall neutrality of the systems, counterions were added. All systems were subjected to MD simulations in explicit solvent, and all systems were solvated with TIP3P water box with 10 Šbetween the solute boundary (Yoo and Xantheas, 2011). First, protein and inhibitor were fixed with a 500 kcal/mol/Ų, and minimized the energy of all water molecules and counterions for 10000 steps of steepest descent (SD) followed by 10000 steps of the conjugate gradient. Subsequently, to remove conflicting contacts, the entire system was repeated for 12000 steps of SD minimization and 8000 steps of CG minimization. Next, the system was gradually heated up to 310 K in the NVT ensemble, thereby applying harmonic restraints with a force constant of 10.0 kcal/mol/Ų on the solute atoms, and equilibration was performed three times with 3000 ps using a force constant of 5.0 kcal/mol/Ų (Loncharich et al., 1992). Finally, 200 ns MD simulations were performed using the NPT ensemble without restraints. We used the Particle mesh Ewald (Darden et al., 1993) technique with a non-bonded cutoff of 12.0 Što limit the direct space sum to treat the long-range electrostatic interactions. The SHAKE (Krautler et al., 2001) algorithm was used to constrain bonds involving hydrogen atoms. The time step of MD simulation was set to 2 fs, and sampling was performed every 10 ps into the MD file. Cluster analysis was performed by employing average linkage as the clustering algorithm (Shao et al., 2007). All figures in this contribution were generated by PyMOL (DeLano, 2002). The hydrogen bonds and hydrophobic interactions between the inhibitors and CYP51 were studied using LigPlot + v.2.2 (Laskowski and Swindells, 2011).

Free Energy Calculations

The Molecular Mechanics Generalized Born Surface Area (MM-GB/SA) method implemented in AMBER 16 package was performed to calculate the binding free energy (Sun et al., 2014; Case et al., 2016), as well as the free energy decomposition of the four inhibitors systems. For each system, 2000 snapshots were selected from the last 100 ns MD trajectories to calculate the relevant energies. The binding free energy (ΔG_bind) in MM-GB/SA between enzyme and ligand was summarized as the follows:

$$\mathrm{\Delta }{G}_{\mathit{bind}}=G_{\mathit{complex}}-(G_{\mathit{receptor}}+G_{\mathit{ligand}})(1)$$

$$\mathrm{\Delta }G=\mathrm{\Delta }{E}_{\mathit{MM}}+\mathrm{\Delta }{G}_{\mathit{sol}}-T\mathrm{\Delta }S(2)$$

$$\mathrm{\Delta }{E}_{\mathit{MM}}=\mathrm{\Delta }{E}_{\mathit{int}}+\mathrm{\Delta }{E}_{\mathit{ele}}+\mathrm{\Delta }{E}_{\mathit{vdW}}(3)$$

$$\mathrm{\Delta }{G}_{\mathit{sol}}=\mathrm{\Delta }{G}_{\mathit{GB}}+\mathrm{\Delta }{G}_{\mathit{SA}}(4)$$

In equation (2), ΔE_MM, ΔG_sol, and TS represent molecular mechanics components in the gas phase, the stabilization energy due to solvation, and a vibrational entropy term, respectively. ΔE_MM represents the summation of ΔE_int, ΔE_ele, and ΔE_vdW which are the internal, coulomb, and van der Waals interaction terms, respectively. ΔG_sol represents the solvation energy, which is divided into the electrostatic solvation free energy (ΔG_GB) (Hawkins et al., 1996) and the non-polar solvation free energy (ΔG_SA). ΔG_GB can be obtained by using the generalized Born method, and ΔG_SA is calculated as follows:

$$\mathrm{\Delta }{G}_{\mathit{SA}}=\mathrm{\gamma }\mathit{SASA}+\mathrm{\beta }(5)$$

Here, γ and β, two empirical constants, were set as 0.0072 and 0.00 kcal/mol/Ų, respectively, and SASA (Weiser et al., 1999) represents the solvent accessible surface area determined by a probe radius of 1.4 Å. To estimate the change in conformational entropy (TΔS) for all atoms, the normal-mode analysis was performed using the nmode module of AMBER 16 package (Case et al., 2016). 100 snapshots from the last 100 ns MD trajectories were used to calculate the entropic contribution.

Energy decompositions were performed to identify the important residues within the systems. Here, only per-residue decomposition was included, which was used to separate the energy contribution of each residue from the combination of enzyme with the inhibitor into three terms: van der Waals contribution (ΔE_vdW), electrostatic contribution (ΔE_ele), and solvation contribution (ΔG_GB + ΔG_SA).

Analysis of Access Tunnels

CAVER (Chovancova et al., 2012) is a famous software to explore routes leading from buried cavities (active sites) to enzyme surfaces. The starting point for the tunnel search was located in the position between heme and the inhibitor. The CAVER algorithm (Petrek et al., 2006) divides three-dimensional space into a grid and calculations are based on grid points. During calculations, the probe radius and the clustering threshold were set to 0.8 and 4.5 Å, respectively. A total of 200 frames of each system were extracted from the last 100 ns MD simulations trajectories. Other parameters were maintained at their default settings. Subsequently, tunnels were visualized by using PyMOL (DeLano, 2002).

Results and Discussion

Determination of the Optimal Binding Pose of Inhibitor by Docking Analysis

The CDOCKER protocol (Studio, 2011) is a CHARMm-based docking method, which was carried out to obtain an optimum initial model of the complex. To determine the reliability of this docking method, co-crystallized ligand (Pos) was firstly re-docked into defined cavity with the CDOCKER protocol. It has been reported that the distance between the nitrogen atoms on the triazole ring of ligands and the iron atom of heme (N6_Flu-Fe, N1_Vor-Fe, N4_Pos-Fe; N7_Itc-Fe) is less than 5 Å (Nair et al., 2016; Hargrove et al., 2017). The optimal conformation was the one with the best score among the structures that satisfy the above distance condition. The root-mean-square deviation (RMSD) value between the docking and initial conformation of Pos was 1.18 Å, which suggested that the CDOCKER protocol was suitable for docking in this work. The ligands Flu, Vor, and Itc were successively docked into the receptor, and the optimal binding pose was selected for the further MD analyses according to the above criteria. For the sake of clarity, Flu-CYP51 complex, Vor-CYP51 complex, Itc-CYP51 complex, and Pos-CYP51 complex was referred to as Flu system, Vor system, Itc system, and Pos system, respectively.”

The Structural Stability and Dynamics Properties of the Inhibitor-CYP51 Systems

In the 200 ns MD simulations of five systems, the root-mean-square deviation (RMSD) value of backbones atom of protein, binding cavity residues, and heavy atom of inhibitors were calculated to investigate the structural stability of CYP51. As shown in Figures 2A–C, each system gradually reached equilibrium, which remained quite stable during the last 100 ns. Thus, all subsequent analysis was performed on the last 100 ns of the simulation trajectories. Comprehensively considered these RMSD values of systems, the binding of the inhibitors reduced the perturbation of the protein to some extent. To further explore the effect of inhibitors binding on fluctuations of a certain residue, the root-mean-square-fluctuation (RMSF) of backbones Cα atoms in CYP51 was calculated (Figure 2D). The comparison of RMSF between the Apo system and four inhibitors systems showed that the fluctuations were mostly similar except for the structural elements of F-F′′ loop, F′′ helix, and F′′-G loop (Supplementary Figure S1). This region is also called F-G loop in P450 enzyme, which may affect the channel characteristics of the enzyme (Cojocaru et al., 2007). The F- F′′ loop, F′′ helix, and F′′-G loop showed large RMSF values in Apo system (Figure 2D). Compared with Apo system, the RMSF values of Flu, Vor, and Itc systems were reduced, while that of Pos system was increased. Further analysis showed that due to the 2-hydroxypentan of Pos is close to the F-G loop, the instability of 2-hydroxypentan (Supplementary Figure S2) caused the F-G loop to change greatly, thus increasing the RMSF value of the F-G loop. These results indicated that inhibitors binding may affect protein tunnel characteristics by influencing the conformations of F- F′′ loop, F′′ helix, and F′′-G loop. These results suggested that inhibitors’ binding increased the stability of CYP51. Local conformational changes of F-G loop caused by inhibitors binding may affect the protein tunnel characteristics.

FIGURE 2

Figure 2. (A) RMSDs of backbone atoms of protein, (B) backbone atoms of binding cavity residues, and (C) heavy atoms of the inhibitor as a function of time. Black, red, blue, pink, and green for Flu, Vor, Itc, Pos, and Apo system, respectively. (D) RMSF values of the backbone Cα atom of each residue. The region of F-F′′ loop, F′′ helix, and F′′-G loop is shown in cyan rectangle.

Analysis Inhibitor Binding Mode

To explore the binding mode of inhibitor in the binding cavity of CYP51, clustering analysis was used to extract the representative conformation. As shown in Figures 3A, 4A, the shared triazole ring of four inhibitors located above the heme was coordinated with the heme. The halogenated phenyl was pointed toward to the crack between the I helix and the B′-C loop. The rest side chain oriented toward the entrance of the binding cavity of CYP51. These results indicated the four inhibitors maintained a similar binding pattern. The 2D diagrams displayed the interaction between inhibitor and protein. As shown in Figures 3B, 4B, the hydrophobic interaction was the main driving force for inhibitors binding to CYP51. Based on molecular shapes and scaffolds, the four inhibitors were divided into two kinds of inhibitors: ST inhibitors and LT inhibitors (Keniya et al., 2018). ST inhibitors (Flu and Vor) formed hydrophobic interaction with the shared residues Y118, F126, Y132, F228, G303, G307, and T311 (Figure 3B). For LT inhibitors (Itc and Pos), they formed hydrophobic interaction with the common residues F58, Y64, Y118, L121, Y132, L376, S378, S506, S507, and M508 (Figure 4B). The important and obvious discrepancy of two kinds of inhibitors was the length of the side chain. Comparing with ST inhibitors, the extended side chains of LT inhibitors (Itc and Pos) provided additional points in contact with the azole target CYP51. Thus, LT inhibitors can form more hydrophobic interactions with CYP51 than ST inhibitors, which may demonstrate LT inhibitors have stronger binding affinities with CYP51.

FIGURE 3

Figure 3. (A) The representative structures of Flu system and Vor system. (B) The 2D diagrams of the detailed binding information of Flu system and Vor system. The protein is shown in white cartoon, the inhibitor and heme are displayed in orange sticks. The molecular interactions show hydrophobic interactions as semi-arcs with red eyelashes.

FIGURE 4

Figure 4. (A) The representative structures of Itc system and Pos system. (B) The 2D diagrams of the detailed binding information of Itc system and Pos system. The protein is shown in white cartoon, the inhibitor and heme are displayed in orange sticks. The molecular interactions show hydrogen bonds as green dashed lines and hydrophobic interactions as semi-arcs with red eyelashes.

The analysis of inhibitor binding mode suggested four inhibitors hold a similar binding pattern and the hydrophobic interactions were the dominant driving force for inhibitors’ binding to CYP51. For both types of inhibitors, LT inhibitors can form hydrophobic interactions with more residues due to the characteristics of their long side chains. This may indicate that LT types of inhibitors are more suitable for targeted CYP51.

Rational Ranking of Binding Ability by Binding Free Energy Calculations

To gain energic information about the four inhibitors systems, the binding free energy calculations were performed by the MM-GB/SA method, and the entropy contributions were also considered. As presented in Table 1, the total ΔG_bind values of Flu, Vor, Itc, and Pos systems were −8.50, −15.08, −52.52, and −48.53 kcal/mol, respectively. LT inhibitors (Pos and Itc) had stronger binding affinities when compared with ST inhibitors (Flu and Vor), which confirmed our previous speculation. Our results showed the rational ordering of binding affinities in different systems which were consistent with that of the experimental inhibitory effects (Hargrove et al., 2017). Further analysis of the binding free energy indicated that the contributions of van der Waals interactions (ΔE_vdW), electrostatic energy (ΔE_ele), and non-polar solvation energy (ΔG_non–polar) were favorable for the formation of the inhibitors’ complexes. As listed in Table 1, non-polar interaction (ΔG_non–polar) was mainly responsible for the formation of the Flu system (−47.97 kcal/mol), Vor system (−54.83 kcal/mol), Pos system (−101.34 kcal/mol), and Itc system (−100.09 kcal/mol). In comparison to the non-polar interaction, polar interaction of the four systems had an unfavorable contribution. The entropy change values of the four systems were less than zero, which were detrimental to the binding of the inhibitors to the CYP51. The results of free energy analysis elaborated that non-polar interaction was the key factor for the binding of inhibitors and CYP51.

TABLE 1

Table 1. Binding free energies (kcal/mol) and its components of four inhibitors systems.

The total binding free energy was decomposed into residues to identify key residues for inhibitors binding to CYP51. Essential residues with the binding free energy value below −1.0 kcal/mol were listed in Table 2. The number of residues meeting to the criterion were 4, 3, 14, and 12 in four systems, respectively, which also indicated that the LT inhibitors were tightly bound to CYP51. We found that Y118 and L376 had significant contributions in inhibitors binding of all four inhibitors systems, Y64, L87, L88, P230, F233, F380, and M508 made outstanding contributions during the LT inhibitors binding.

TABLE 2

Table 2. Decomposition of binding free energy (kcal/mol) on per residue basis for Flu system, Vor system, Pos system, and Itc system.

Tunnel Analysis

Illustrating the tunnel characteristic of CYP51 is beneficial to develop new inhibitors and understand the structure–function relationships of CYP51 (Yu et al., 2016; Fischer et al., 2018). In this work, 200 frames were extracted from the last 100 ns trajectories to classify and analyze the characteristics of access pathways in four inhibitors systems and Apo system. The tunnels were clustered by the average-link algorithm according to the pairwise distances of tunnels. On the basis of spatial and secondary structure, the nomenclature of these tunnels is defined systematically by Wade group (Cojocaru et al., 2007). The five highest ranked tunnels of five systems were all displayed in Figure 5, and the characteristics of these tunnels were summarized in Table 3. As shown in Figure 5, the locations of five tunnels were marked with different color spheres (Flu system: 2f, W, S, 1, and 2e; Vor system: 2f, 2a, W, S, and 2e; Itc system: 2f, S, W, 1, and 2a; Pos system: 2f, 2a, W, S, and 2ac; Apo system: 2f, S, W, 2a, and 2ac). Tunnels 2f, 2a, 2e, and 2ac are subclasses of tunnel 2. Tunnel 2f locates between the F-G loop, Pro-rich loop, and A helix, whereas tunnel 2a locates between the F-G loop and B′-C loop (Supplementary Figure S1). Tunnel 2e egresses through the B′-C loop, and tunnel 2ac egresses between the B′ helix and the G helix (Supplementary Figure S1). Tunnel W (water tunnel) egresses at the base of the B′-C loop near the C-terminus of the B helix, and tunnel S (solvent tunnel) runs between the E, F, and I helices and β5 sheet (Supplementary Figure S1). As listed in Table 3, tunnel 2f was the most frequently identified collective pathway and had the highest bottlenecks radius in five systems. Thus, tunnel 2f was regarded as the predominant tunnel for inhibitors ingress/egress from the active site to the surface of CYP51, which was similar as the other works of CYP51 (Monk et al., 2014; Yu et al., 2016; Gao et al., 2018). The occurrence of tunnel 2f was different slightly in five systems, which may be related to the size of the inhibitors. In the LT systems (Itc and Pos), the inhibitor is long and bulky, and its binding mode analysis showed that its long side chain extended to the entrance of the tunnel 2f, resulting in a fully opened tunnel 2f. In the ST systems (Flu and Vor), the inhibitor is in small size, and the inhibitor was submerged in the binding cavity of CYP51, causing decreased slightly in the opening frequency of tunnel 2f, which was in line with the RMSF analysis that inhibitors’ binding affected the tunnel characteristic. Further, we determined the essential residues lining the dominant tunnel 2f, and all residues located within the 3 Å distance from the tunnel surface will be regarded as tunnel-lining residues (Supplementary Table S1). The key residues (F58, Y64, Y118, L121, Y132, L376, S378, S506, S507, and M508) determined based on the binding mode analysis and per-residue binding free energy decomposition analysis also belong to the tunnel-lining residues. Most of the tunnel-lining residues were hydrophobic residues, which formed a stable hydrophobic cavity and provided hydrophobic interactions that play an indispensable role in inhibitor stabilization. Thus, when designing more efficient inhibitors, the interactions between inhibitors and these residues should be rationally increased and the new inhibitors should be hydrophobic ligands.

FIGURE 5

Figure 5. Access tunnel identified from the average structures of (A) Flu system, (B) Vor system, (C) Itc system, (D) Pos system, and (E) Apo system, respectively. The protein backbone is shown in cartoon. Tunnels are shown in sphere for blue (2f), green (2a), magenta (2e), orange (2ac), yellow (s), red (w), and cyan (1).

TABLE 3

Table 3. Characteristics of the five top ranked tunnels of Flu system (A), Vor system (B), Itc system (C), Pos system (D), and Apo system, respectively.

Conclusion

The sterol 14α-demethylase enzyme (CYP51) belongs to cytochrome P450 family essential in sterol biosynthesis, which is the target for fungal infections. In this work, molecular docking and molecular dynamics simulations were employed to investigate the binding mechanism and tunneling characteristics between four inhibitors and CYP51, so as to provide useful information for inhibitors design. The results show that four inhibitors bind CYP51 in similar binding mode and hydrophobic interactions are the main driving force for inhibitors binding to CYP51. Due to long-tailed inhibitors (posaconazole and itraconazole) can contact with more residues of CYP51 by hydrophobic interactions than short-tailed inhibitors (fluconazole, voriconazole), long-tailed inhibitors have stronger binding affinities. Tunnel analysis showed that tunnel 2f is the predominant pathway for inhibitors ingress/egress from the active site to the surface of CYP51. We discover a hydrophobic cavity and identify the key residues (F58, Y64, Y118, L121, Y132, L376, S378, S506, S507, and M508) which are responsible for anchoring the inhibitors binding to CYP51. Therefore, in order to enhance the binding affinity of inhibitors to CYP51, we should focus on strengthening hydrophobic interactions of inhibitors and these residues, while longer inhibitors are probably best suited to target CYP51. Taken together, the results obtained in this study will be valuable for designing potent azoles inhibitors and improve the understanding of the structure–function relationships of CYP51.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

Author Contributions

NS performed the research, analyzed the data, and wrote the manuscript. QZ and HZ designed the research and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 21273095) and the Fundamental Research Funds for the Central Universities, JLU.

Conflict of Interest

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

Supplementary Material

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

References

Andersson, M. P., and Uvdal, P. (2005). New scale factors for harmonic vibrational frequencies using the B3LYP density functional method with the triple-xi basis set 6-311+G(d,p). J. Phys. Chem. A 109, 2937–2941. doi: 10.1021/jp045733a

PubMed Abstract | CrossRef Full Text | Google Scholar

Antinori, S., Milazzo, L., Sollima, S., Galli, M., and Corbellino, M. (2016). Candidemia and invasive candidiasis in adults: a narrative review. Eur. J. Intern. Med. 34, 21–28. doi: 10.1016/j.ejim.2016.06.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Balding, P. R., Porro, C. S., McLean, K. J., Sutcliffe, M. J., Marechal, J. D., Munro, A. W., et al. (2008). How do azoles inhibit cytochrome P450 enzymes? A density functional study. J. Phys. Chem. A 112, 12911–12918. doi: 10.1021/jp802087w

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Bongomin, F., Gago, S., Oladele, R. O., and Denning, D. W. (2017). Global and multi-national prevalence of fungal diseases-estimate precision. J. Fungi 3:57. doi: 10.3390/jof3040057

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, G. D., Denning, D. W., Gow, N. A. R., Levitz, S. M., Netea, M. G., and White, T. C. (2012). Hidden killers: human fungal infections. Sci. Transl. Med. 4:165rv13. doi: 10.1126/scitranslmed.3004404

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D., Betz, R., Cerutti, D. S., Cheatham, T., Darden, T., Duke, R., et al. (2016). Amber 2016. San Francisco, CA: University of California.

Google Scholar

Choi, J. Y., and Roush, W. R. (2017). Structure based design of CYP51 inhibitors. Curr. Top. Med. Chem. 17, 30–39. doi: 10.2174/1568026616666160719164933

PubMed Abstract | CrossRef Full Text | Google Scholar

Chovancova, E., Pavelka, A., Benes, P., Strnad, O., Brezovsky, J., Kozlikova, B., et al. (2012). CAVER 3.0: a tool for the analysis of transport pathways in dynamic protein structures. PLoS Comput. Biol. 8:e1002708. doi: 10.1371/journal.pcbi.1002708

PubMed Abstract | CrossRef Full Text | Google Scholar

Cojocaru, V., Winn, P. J., and Wade, R. C. (2007). The ins and outs of cytochrome P450s. Biochim. Biophys. Acta Gen. Subj. 1770, 390–401. doi: 10.1016/j.bbagen.2006.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: an N ⋅ log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397

CrossRef Full Text | Google Scholar

DeLano, W. (2002). PyMOL Molecular Graphics System. San Carlos, CA: DeLano Scientific.

Google Scholar

Ene, I. V., Adya, A. K., Wehmeier, S., Brand, A. C., MacCallum, D. M., Gow, N. A. R., et al. (2012). Host carbon sources modulate cell wall architecture, drug resistance and virulence in a fungal pathogen. Cell. Microbiol. 14, 1319–1335. doi: 10.1111/j.1462-5822.2012.01813.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fischer, A., Don, C. G., and Smiesko, M. (2018). Molecular dynamics simulations reveal structural differences among allelic variants of membrane-anchored cytochrome P450. J. Chem. Inform. Model. 58, 1962–1975. doi: 10.1021/acs.jcim.8b00080

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M., Trucks, G., Schlegel, H., Scuseria, G., Robb, M., Cheeseman, J., et al. (2009). Gaussian 09. Wallingford, CT: Gaussian Inc.

Google Scholar

Gao, P., Cui, Y. L., and Wu, R. L. (2018). Molecular dynamic modeling of CYP51B in complex with azole inhibitors. J. Biomol. Struct. Dyn. 36, 1511–1519. doi: 10.1080/07391102.2017.1328315

PubMed Abstract | CrossRef Full Text | Google Scholar

Hargrove, T. Y., Friggeri, L., Wawrzak, Z., Qi, A., Hoekstra, W. J., Schotzinger, R. J., et al. (2017). Structural analyses of Candida albicans sterol 14 alpha-demethylase complexed with azole drugs address the molecular basis of azole-mediated inhibition of fungal sterol biosynthesis. J. Biol. Chem. 292, 6728–6743. doi: 10.1074/jbc.M117.778308

PubMed Abstract | CrossRef Full Text | Google Scholar

Hargrove, T. Y., Friggeri, L., Wawrzak, Z., Sivakumaran, S., Yazlovitskaya, E. M., Hiebert, S. W., et al. (2016). Human sterol 14 alpha-demethylase as a target for anticancer chemotherapy: towards structure-aided drug design. J. Lipid Res. 57, 1552–1563. doi: 10.1194/jlr.M069229

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawkins, G. D., Cramer, C. J., and Truhlar, D. G. (1996). Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem. 100, 19824–19839. doi: 10.1021/jp961710n

CrossRef Full Text | Google Scholar

Keniya, M. V., Sabherwal, M., Wilson, R. K., Woods, M. A., Sagatova, A. A., Tyndall, J. D. A., et al. (2018). Crystal structures of full-length lanosterol 14 alpha-demethylases of prominent fungal pathogens Candida albicans and Candida glabrata provide tools for antifungal discovery. Antimicrob. Agents Chemother. 62:e01134-18. doi: 10.1128/aac.01134-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Krautler, V., Van Gunsteren, W. F., and Hunenberger, P. H. (2001). A fast SHAKE: algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22, 501–508. doi: 10.1002/1096-987x(20010415)22:5<501::aid-jcc1021>3.0.co;2-v

CrossRef Full Text | Google Scholar

Laskowski, R. A., and Swindells, M. B. (2011). LigPlot+: multiple ligand-protein interaction diagrams for drug discovery. J. Chem. Inform. Model. 51, 2778–2786. doi: 10.1021/ci200227u

PubMed Abstract | CrossRef Full Text | Google Scholar

Lass-Floerl, C. (2011). Triazole antifungal agents in invasive fungal infections a comparative review. Drugs 71, 2405–2419. doi: 10.2165/11596540-000000000-00000

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C., Yang, W., and Parr, R. G. (1988). Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B Condensed Matter 37, 785–789. doi: 10.1103/PhysRevB.37.785

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y., Puumala, E., Robbins, N., and Cowen, L. E. (2020). Antifungal drug resistance: molecular mechanisms in Candida albicans and beyond. Chem. Rev. [Epub ahead of print]. doi: 10.1021/acs.chemrev.0c00199

PubMed Abstract | CrossRef Full Text | Google Scholar

Lepesheva, G. I., Friggeri, L., and Waterman, M. R. (2018). CYP51 as drug targets for fungi and protozoan parasites: past, present and future. Parasitology 145, 1820–1836. doi: 10.1017/s0031182018000562

PubMed Abstract | CrossRef Full Text | Google Scholar

Lepesheva, G. I., Hargrove, T. Y., Kleshchenko, Y., Nes, W. D., Villalta, F., and Waterman, M. R. (2008). CYP51: a major drug target in the cytochrome P450 superfamily. Lipids 43, 1117–1125. doi: 10.1007/s11745-008-3225-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Loncharich, R. J., Brooks, B. R., and Pastor, R. W. (1992). Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N’-methylamide. Biopolymers 32, 523–535. doi: 10.1002/bip.360320508

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Monk, B. C., Tomasiak, T. M., Keniya, M. V., Huschmann, F. U., Tyndall, J. D. A., O’Connell, J. D. I. I. I., et al. (2014). Architecture of a single membrane spanning cytochrome P450 suggests constraints that orient the catalytic domain relative to a bilayer. Proc. Natl. Acad. Sci. U.S.A. 111, 3865–3870. doi: 10.1073/pnas.1324245111

PubMed Abstract | CrossRef Full Text | Google Scholar

Nair, P. C., McKinnon, R. A., and Miners, J. O. (2016). Cytochrome P450 structure-function: insights from molecular dynamics simulations. Drug Metab. Rev. 48, 434–452. doi: 10.1080/03602532.2016.1178771

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishimoto, A. T., Sharma, C., and Rogers, P. D. (2020). Molecular and genetic basis of azole antifungal resistance in the opportunistic pathogenic fungus Candida albicans. J. Antimicrob. Chemother. 75, 257–270. doi: 10.1093/jac/dkz400

PubMed Abstract | CrossRef Full Text | Google Scholar

Petrek, M., Otyepka, M., Banas, P., Kosinova, P., Koca, J., and Damborsky, J. (2006). CAVER: a new tool to explore routes from protein clefts, pockets and cavities. BMC Bioinformatics 7:316. doi: 10.1186/1471-2105-7-316

PubMed Abstract | CrossRef Full Text | Google Scholar

Shahrokh, K., Orendt, A., Yost, G. S., and Cheatham, T. E. III (2012). Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle. J. Comput. Chem. 33, 119–133. doi: 10.1002/jcc.21922

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, J., Tanner, S. W., Thompson, N., and Cheatham, T. E. III (2007). Clustering molecular dynamics trajectories: 1. Characterizing the performance of different clustering algorithms. J. Chem. Theory Comput. 3, 2312–2334. doi: 10.1021/ct700119m

PubMed Abstract | CrossRef Full Text | Google Scholar

Studio, D. (2011). Version 3.1. San Diego, CA: Accelrys.

Google Scholar

Sun, D. R., Wang, Z. J., Zheng, Q. C., and Zhang, H. X. (2018). Exploring the inhibition mechanism on HIF-2 by inhibitor PT2399 and 0X3 using molecular dynamics simulations. J. Mol. Recogn. 31:e2730. doi: 10.1002/jmr.2730

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, H., Li, Y., Shen, M., Tian, S., Xu, L., Pan, P., et al. (2014). Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 16, 22035–22045. doi: 10.1039/c4cp03179b

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006). Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graph. Model. 25, 247–260. doi: 10.1016/j.jmgm.2005.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. M., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

Warrilow, A. G., Nishimoto, A. T., Parker, J. E., Price, C. L., Flowers, S. A., Kelly, D. E., et al. (2019). The evolution of azole resistance in Candida albicans sterol 14 alpha-demethylase (CYP51) through incremental amino acid substitutions. Antimicrob. Agents Chemother. 63:e02586-18. doi: 10.1128/aac.02586-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, G. M., Lucas, W. A. H., Gunzburg, M. J., and Wilce, J. A. (2017). Insight into the selectivity of the G7-18NATE inhibitor peptide for the Grb7-SH2 domain target. Front. Mol. Biosci. 4:64. doi: 10.3389/fmolb.2017.00064

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiser, J., Shenkin, P. S., and Still, W. C. (1999). Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 20, 217–230. doi: 10.1002/(sici)1096-987x(19990130)20:2<217::aid-jcc4>3.0.co;2-a

CrossRef Full Text | Google Scholar

Xiao, F., Song, X., Tian, P., Gan, M., Verkhivker, G. M., and Hu, G. (2020). Comparative dynamics and functional mechanisms of the CYP17A1 tunnels regulated by ligand binding. J. Chem. Inform. Model. 60, 3632–3647. doi: 10.1021/acs.jcim.0c00447

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoo, S., and Xantheas, S. S. (2011). Communication: the effect of dispersion corrections on the melting temperature of liquid water. J. Chem. Phys. 134:121105. doi: 10.1063/1.3573375

CrossRef Full Text | Google Scholar

Yu, X., Nandekar, P., Mustafa, G., Cojocaru, V., Lepesheva, G. I., and Wade, R. C. (2016). Ligand tunnels in T-brucei and human CYP51: insights for parasite-specific drug design. Biochim. Biophys. Acta Gen. Subj. 1860, 67–78. doi: 10.1016/j.bbagen.2015.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar