Abstract

Today, fungal infection has become more common disease especially in some cases, such as AIDS, cancer, and organ transplant which the immune system is suppressed. On the other hand, due to the increasing resistance to current antifungal drugs, more and more options for design of novel more efficient compounds with higher resistance are needed. In this study, a series of a fluconazole analogues were subjected to quantitative structure-activity relationship analysis to find the structure requirements for modeling adequate candidate. The best multiple linear regression equation was achieved from GA-PLS and MLR modeling. Subsequently, in silico screening study was applied to found new potent lead compounds based on the resulted model. The ability of the best designed compounds for antifungal activity was investigated by using molecular dynamic (MD) and molecular docking simulation. The results showed that compound F13 can efficiently bind to lanestrol 14-α demethylase target similar to other antifungal azoles. The molecular docking studies revealed an interesting binding profile with very high receptor affinity to the CYP51 active site. The triazole moiety of ligand F13 pointed to HEM group in lanestrol 14-α demethylase site and coordinated to Fe of HEM through its N4 atom. Also, there was a convenient relevance between QSAR and docking results. With the compound F13 which demonstrated the most promising minimum inhibitory concentration (MIC) values, it can be concluded that F13 is appropriate candidate for the development as antifungal agent.

1. Introduction

Recently, fungal infections have become more prevalent and play a considerable role into human morbidity and mortality. Currently, there are five categories of medicines for the treatment of fungal diseases worldwide, whether orally or injectable such as polyenes, pyrimidine analogues, allylamines, azoles, and the echinocandins [13]. Fluconazole was the first triazole antifungal agent which inhibits ergosterol biosynthesis and has been developed by Pfizer [4]. However, with the relatively limited spectrum of activity against Candida species and as some species of fungi are resistant against fluconazole [57], the demand for development of new fluconazole analogues is increased [1, 8].

Computational techniques could have benefits such as time-saving for early design of compounds and reducing costs in the new drug production process [9]. The quantitative structure activity relationships (QSAR) studies is a way to achieve more appropriate structures and modify the previous structure by using in silico modeling [1012]. This method could be used to predict the activity of nonsynthesized compounds. Due to the lack of information from molecular docking simulation about biomolecules and ligand complexes, the molecular dynamics (MD) simulations can compensate for this deficiency. Recently, one of the most distinguished tools for achieving data about chemical features of the protein, nucleic acid, and biomolecule-ligand complexes is MD simulations.

In this study, we utilized QSAR and molecular modelling techniques including MD simulation and molecular docking to achieve the meaningful relationship between fluconazole-like structures and their antifungal activities. Also, to access the exact molecular connection models, molecular modelling techniques were performed on all the data set components and designed compounds. Then, in silico study were applied to introduce novel compounds based on our QSAR and docking results as effective candidates toward fungal infection. Finally, molecular docking simulation was run to survey our designed compounds.

2. Material and Methods

2.1. Quantitative Structure Activity Relationship (QSAR)
2.1.1. Data Set

In order to provide a source and prior preparation, 97 fluconazole derivatives were collected with similarity in chemical structure and acceptable MIC [1322]. To classify the data set into calibration and test sets, the kenardston algorithm was used. The softwares including Hyperchem v8.0 (Hypercube Inc., Gainesville, FL, USA) (descriptor calculation), Gaussian98 software (Gaussian 09 Revision-A.02-SMP) (descriptor calculation and optimized the structure), Dragon (version 5.5) (descriptor calculation), MGL tools 1.5.6 (docking), MODELLER 9.17 (docking), Ligplot v.2.2, Chimera 1.9 viewer (visualize the interaction), GROMACS 2019.6 package (simulation), Avogadro software (version 1.2.0) (operate the 3D structure of ligands), Swiss-PDB viewer 4.10, and Chimera 1.10.2 software (visualize the interaction) were used. The structural features and antifungal activity details of these derivatives are given in Table 1.

2.1.2. Molecular Descriptor

After drawing chemical structure of compounds by using Hyperchem v8.0 (Hypercube Inc., Gainesville, FL, USA), the optimization process was carried out by AM1 and MM+ algorithm. In the next step, molecular descriptors (molecular volume (V), hydration energy (HE), hydrophobicity (LogP), molecular surface area (SA), and molecular polarizability (MP) were computed using Hyperchem v8.0. The highest occupied molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO) energies, the foremost positive, the most negative net atomic charges, the average absolute atomic charge, and also the molecular dipole moment were calculated by Gaussian98 software (Gaussian 09 Revision-A.02-SMP). Different topological, geometrical, empirical, and constitutional descriptors; 2D autocorrelations; atom-centered fragments; and functional groups were also calculated using the Dragon (version 5.5) for every molecule.

Three different regression methods including multiple linear regression with stepwise variable selection (MLR), MLR with factor analysis as the data preprocessing step for variable selection (FA-MLR), and genetic algorithm–partial least squares (GA-PLS) were applied to obtain QSAR equation. Statistical parameters such as standard error of the regression (SE), coefficient of determination (2), Fisher test (F) at specified degrees of freedom, and leave-one-out cross-validation correlation coefficient (2) were exploit for the validity of regression equation [23]. The predictive value of a QSAR model has been evaluated on an external set of data (2p).

2.2. Molecular Docking Protocols

An in-house batch script (DOCK-FACE) was applied to running docking procedure [24]. As an initial step, each ligand structure was prepared and optimized. Then, the Gasteiger charges were implemented, nonpolar hydrogen of compounds were merged, and rotatable bonds were assigned to obtain PDBQT format by MGL tools 1.5.6 [25]. The three-dimensional crystal structures of lanestrol 14-α demethylase (PDB ID: 1EA1) was obtained from the Protein Data Bank (PDB data base; http://www.rcsb.org) [26]. All water molecules and cocrystal ligand were removed, and the missing hydrogens were added. The missing atoms of PDBs were then corrected with MODELLER 9.17 [27]. A grid box of 70 × 70 × 70 points in , , and directions was built and centred on the ligand in the complex with a spacing of 0.375 Å for 1EA1. The cocrystal ligand (fluconazole) was excluded and treated the same as examined ligands for internal validation. All the docking procedures were done on validated procedure with a root mean square deviation (RMSD) value below of 2 Å. All interactions were visualized and evaluated on the basis of docking results using Ligplot v.2.2 and Chimera 1.9 viewer.

2.3. Molecular Dynamics Simulations Protocols

The molecular dynamic simulation was done through the GROMACS 2019.6 package. The crystal structure of cytochrome P450 enzyme as a (PDB ID: 1EA1) was achieved from http://www.rcsb.com/ The 3D structure of ligand was operated from Avogadro software (version 1.2.0). The Gaussian 98 software (09 W v7.0) was applied to the minimization of the 3D structure by the DFT, B3LYP. All interactions were performed using Swiss-PDB viewer 4.10 and Chimera 1.10.2 software. The amber tools package was used to obtain the parameters of ligand 32 and F13, fluconazole, and also heme parameters. The ACPYPE script was applied to convert the amber Gromacs format. The bond length, bond constant force, angle, and its constant force were defined. Firstly, the cys-heme-ligand complex was separated from the origin one; then, it was minimized by the B3LYP and basic system 6-31G. Then, the minimized-complex was utilized to calculate bond constant force of Fe-S and Fe-N through Seminario as well as derivation of constant force from hess2ff. For the calculation of partial charges, heme-cys was utilized by the restrained electrostatic potential via applying coulomb potential to calculate the electrostatics interactions. Although at the amber force fields, the HF/6-31G level is often employed for parametrization. The input structures were prepared by ff99SB + ILDN as the force field. The correct situation of Histidine hydrogen was defined for all histidine inside of protein and disulfide boned. In the next level, the surface charges of the structure were neutralized by adding Na and Cl. For defining the box grid, the complex was posed in an octahedron box; this stage was operated via a steep descent algorithm during 1000 steps in order to remove the Van der Waals interactions and construct hydrogen bonds between water and complex. Then, the temperature of the system was increased gradually from 0 to 310 k over 200 ps (volume is constant). Finally, it was equilibrated in the constant pressure. MD simulation was carried in the mentioned temperature over 100 ns.

3. Result and Discussion

3.1. QSAR Analysis

In this study, three chemometrics methods such as (stepwise MLR, FA-MLR, and GA-PLS) were applied to achieve a relationship between the biological activity and molecular descriptors. As shown in Table 2, based on the MLR results, G (CL..CL), n piperazine, nR06, nROR, and nX have a positive effect on the antifungal activity of the compounds. This means that compounds with more number of rings, especially hexagonal rings such as piperazine along with halogen substitution increases antifungal activity. In the other words, azole derivatives with ether group and 1,2,4-triazole ring could improve their antifungal activity, while the negative sign of SIC2, MAXDP, Mor05u, and HATS7m descriptor denotes the reverse effect of these parameters on the pMIC values of the compounds. The calculated statistical parameters are shown in Table 2. The FA-MLR model (Eq. 2) demonstrates the positive effect of nR06 and n piperazine and the negative effect of nArNH2 for the compounds tested toward antifungal activity. As it is observed, about 69% of variances in the original data matrix could be explained by the selected three parameters. The third equation was GA-PLS model. In the PLS model, the latent variables were used to make modeling, so the collinearity between descriptor was omitted. The GA-PLS model was selected as the best equation of our compounds. The best GA-PLS model with the best fitness contained eight descriptors, seven of which were obtained through the final MLR modelling (Figure 1). As it can be seen, a set of functional, constitutional, and geometrical descriptors are the most effective descriptors in the final GA-PLS model.

VIP was employed to understand the relative importance of the variable in the PLS model based on Erikson et al. VIP >1.0 and VIP <0.8 mean profoundly and limited influential, respectively. Also, 0.8< VIP <1.0 means moderately influential. The results of the VIP analysis on the final PLS model are illustrated in Figure 2. VIP indicates the SIC2 and n1,2,4-triazole descriptors as the most influential and significant descriptors in the QSAR equation obtained in the GA-PLS analysis, and the nX, MAXDP descriptors stands second.

In overall, we concluded that structure with pentagonal (1,2,4-triazole) and hexagonal (phenyl) rings represented better antifungal activity. Also, the presence of piperazine or piperidine rings and ether group on the structure can positively impact on the antifungal activity. Besides, the existence electronegative substitution such as chlorine and fluorine, especially in positions 2 and 4 of the phenyl ring, will lead to enhance the biological activity. Electron withdrawing groups such as NO2 can be also placed at the phenyl ring. On the other hand, we avoid to put the amine substitutions on the phenyl ring, as well as alcohol (10) as substitutes on the chain containing amine substitution, too.

3.1.1. Robustness and Applicability Domain of the Models

The leverage value was calculated to estimate the applicability domain which this range could be predict appropriately and reliably for those compounds that are placed in that location. The warning leverage (h) is another criterion for interpretation of the results. The warning leverage is defined with , where is the number of training compounds and is the number of model parameters. The leverage greater than warning leverage means that the model may not be reliable. As seen in Table 1(s), all of the compounds have leverage () values lower than the threshold leverages (). So, all of the compounds were assessment identically and within the suggested models applicability domain.

3.1.2. Y-Randomization Test

Y-randomization test is a method which could evaluate the model and make sure that they did not achieve accidentally. Subsequently, PMIC distributes randomly as a dependent factor for several time; if statistical parameters were obtained less than the origin model, it is considered an accidental model. The results of Y-randomization are presented for each model separately (Table 2(s)). As it is known, it was repeated ten times, and each time, 2LOOCV and 2 were lower than the values presented in the QSAR models.

3.2. Design of New Antifungal Compounds

Computational study was used to discover novel compounds with improved efficacy as antifungal agents based on the QSAR model. As shown in Table 3, new compounds were designed, and their predicted activities for antifungal activity based on GA-PLS equations and their docking binding energies were achieved. Among these molecules (Table 3), compounds F7 and F11-15 showed the best activity. These compounds can be considered as good candidates for becoming antifungal agents.

3.3. Molecular Docking Studies

As it was mentioned, fluconazole-like compounds illustrated their antifungal effects through inhibition of cytochrome P450 14α-demethylase (CYP51) which is a crucial enzyme for sterol biosynthesis in the fungal cell membrane. Molecular docking simulation was performed on all studied (ligand 1-97) and designed compounds (F1-F15) to predict the pattern of the interaction between ligands and receptors through their active sites. This approach can describe the behavior of ligand in the binding site of the target proteins. In the validity evaluation step of docking process, redocking of fluconazole was done; RMSD of docking for fluconazole in comparison with its coordination in the crystal structure was 0.45. The docking binding energies of all studied compounds and designed compounds are shown in Tables 1 and 3, respectively. As shown in these tables, compounds 32 and F13 with high antifungal activities represented stronger energies in binding to lanestrol 14-α demethylase site compared to the other compounds. The binding mode of the fluconazole as standard drug is displayed in Figure 3. A hydrogen bond interaction exists between the C3 atom of the 1,2,4-triazole ring with Ala 256 and Ser 252. The triazole of ligands 32, F13, and fluconazole pointed to HEM group in CYP51 active site and coordinated to Fe of HEM through its N4 atom, and there is also existed some hydrophobic interaction with Leu 321, Tyr 76, and Phe 83.

As it was depicted in Figure 4, in the cytochrome P450 14α-demethylase (CYP51) binding mode of 32, the fluorine of 2,4-di-fluorine phenyl rings is involved in the halogen bond interaction with Gln 109 residues. The piperazine and 2,4-di-fluorine phenyl rings are involved in the carbon-hydrogen bond interactions with Cys 151, Leu 152, and Gly 112 residues, respectively. The side chains consisting of amino acid residues Ala 104, Ala 103, and Met 225 create van der Waals contacts with ligand 32.

As it was shown in Figure 5, the most important residues in binding of F13 to 1EA1 receptor are three hydrogen bonds including a hydrogen bond interaction between 4-fluore phenoxy moiety with Ile 322, a hydrogen bond acceptor interaction between 2,4,6-triflourine phenyl ring with Arg96 and a hydrogen bond interaction between 4-fluore phenoxy moiety and Met 433. The 1,2,4-triazole and 2,4,6-triflourine ring are involved in the pi-cation and pi-alkyl interaction with HEM group and, on the other hand, the F13 compound coordinate to Fe of HEM through its N4 atom. These interactions stabilize the compound in the proper orientation within the active pocket of cytochrome P450 14α-demethylase (CYP51) enzyme. There also exists a halogen bond interaction with Met 433, Leu 321, Pro 320, Phe 255, Ser 252, and Met 99. Meanwhile, the groups of Met 253, Phe 83, ILE 323, Val 435, Val 434, His 259, Phe 78, and Thr 260 involved to create a pocket around the F13 by van der Waals forces, which lead to stabilization of F13 compound to introduce as potential antifungal agents.

3.4. MD Simulation

The simulation results were obtained from Gromacs software as root mean square deviations (RMSD), root mean square oscillations (RMSF), the radius of gyration and secondary structure, and potential energy. Calculating the RMSD between the obtained structures and a reference structure is the most important method for evaluating the stability of molecular dynamics simulations over time. To analyze the main motions in the system, the RMSD was obtained for alpha carbon atoms. The structure flexibility was also measured based on the RMSF of alpha carbon atoms. The radius of gyration, which is one of the determining parameters in the density of the structure, was obtained during the simulation period [28, 29].

3.4.1. RMSD

According to Figure 6, a comparison was carried out between fluconazole, 32, and F13. RMSD for all ligands was relatively equal over the simulation. Whereas the crystal curve as the control reached stability around 0.16 after 5000 ps, the 32 and F13 curves levelled out after 3000 ps and 2000 ps around 0.21 and 0.19, respectively. This finding can present the superiority of F13 in comparison with 32 and, as well, the proximity to the control curve.

According to Figure 7, the compound 32 structure is greater than fluconazole and F13. So it could be a reasonable justification for differences between their RMSD even though they stand at an acceptable quantity (around 0.2 nm).

Moreover, based on Figure 7, the 32 has more rotatable bonds than the crystal as well as F13 which has less rotation compared to the crystal. This matter increases the flexibility of this ligand at the binding site, resulting in protein instability. So, fluconazole and 32 have the most and the least stability, respectively. The MD simulated interaction of fluconazole and ligands 32 and F13 are shown in Figure 8. The interaction with HEM was the main bonding forces between our studied ligands and lanestrol 14-α demethylase enzyme. The binding interaction pattern during 100 ns MD simulation was concordance with the docking results.

3.4.2. RMSF

RMSF was utilized in our study to explore those goals (main movement and flexibility ratio), and the results re shown in Figure 9.

As shown in Figure 9, in most areas, the degree of protein flexibility in all simulations is the same, and only in a few limited regions has the vibration changed. For compound 32, in five regions, the protein flexibility was increased while in only one region did the flexibility shrinks. The areas with risen flexibility include amino acids 15-22 (green), 65-73 (purple), 85-93 (red), 225-232 (blue), and 325-330 (light purple). As it can be seen according to Figure 10, four of these areas are located around the ligand junction, which can be attributed to the increase in flexibility in these areas. Also, as mentioned, only in the area of 285-290, the amount of flexibility has decreased, which is shown with an arrow in Figure 10 of this area. According to this figure, the area is much further away from the ligand junction.

RMSF for F13 was similar to the crystal more or less, whereas there can be a range of changes in the following regions, 65-70, 159-161, 230-240, 271-272, and 210-220, that the flexibility increased and decreased in order.

3.4.3. Radius of Gyration

As shown in Figure 11, the radius of gyration in the presence of all ligands does not display a notable difference. The only variation between the two graphs is that in the presence of 32 the amount of radius of gyration fluctuation is slightly risen, while in the presence of F13 and fluconazole, the amount of radius of gyration to the end of the simulation is quite constant. In terms of the amount of radius of gyration according to this diagram in all proteins, the amount of radius is equal to about 2.27 nm.

3.4.4. Potential Energy

Based on Figure 12, both system under study and references system reached the acceptable stability during the simulation. According to this diagram, it can be concluded that in the presence of F13, the amount of potential energy of the system is higher, which indicates that it is more unstable than the other two simulations (the point to be noted is the difference in the number of water molecules in all three simulations. In fluconazole, 32, and F13 simulations, the number of water molecules is 16226, 16205, and 16048, respectively, which affects the potential energy of the system).

3.4.5. Secondary Structure

As shown in Figure 13, alterations in the secondary structures of the proteins are observed in 3 regions in the presence of the 32 ligand. These areas are amino acids 90-110, 220-230, and 285-290. In the 90-110 region, the protein structure turns from Bend (green) to turn (yellow). In the region 220-230, the protein is converted from the second alpha-helix structure (blue) to the Bend structure. In the 285-290 region, the protein structure is slightly shifted from the alpha-helix state to the Bend state. Furthermore, changes in the secondary structures of the protein in the presence of F13 are similar to the crystal ligand. Table 4 and Figure 13 also show the content of protein secondary structures in all three simulations.

4. Conclusion

Herein, 2D-QSAR, molecular docking, and MD simulation studies were done on a series of fluconazole-like compounds as antifungal agents. The QSAR results indicated that the presence of 1,2,4-triazole rings as well as hexagonal rings which contain the electron withdrawing groups such as chlorine and fluorine at positions 2 and 4 can increase the inhibitory power of the compounds. The best regression equation was obtained from GA-PLS method, which predicts 84% of variances. Subsequently, an in silico screening study was conducted, and the newly designed molecules (compounds F7, F11-15) were predicted to possess higher activity than the data set compounds. The molecular docking results represented that the triazole ring of ligand F13, 32, and fluconazole perched to HEM group in lanestrol 14-α demethylase site and coordinated to Fe of HEM through its N4 atom. According to docking result, there was an appropriate interaction between F13 as the most active compound and enzyme 1EA1. Hydrogen bond interaction with Ile 322, Arg 96, and Met 433 were shown critical for a good inhibitory activity. The groups of Met 253, Phe 83, ILE 323, Val 435, Val 434, His 259, Phe 78, and Thr 260 were involved to create a pocket around the F13 ligand by van der Waals forces, which lead to stabilization of F13 compound to introduce as potential antifungal agents. Two ligands and fluconazole matched up perfectly during the 100 ns of simulation according to RMSD, radius of gyration, potential energy, and secondary structure. The binding interaction pattern during 100 ns MD simulation was concordance with docking results. F13 can be considered a desire candidate for future studies.

Data Availability

There is no data to support the findings of this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Farnaz Salehi and Leila Emami contributed equally to this work.

Acknowledgments

Financial assistance from the Shiraz University of Medical Sciences by way of grant numbers 24523 and 24244 are gratefully acknowledged.

Supplementary Materials

Table 1(s). Leverage (h) of the external test set molecules for different models. Table 2(s). 2 and R2LOOCV after Y-randomization for each Model. Figure 1(s). The detail step of docking study. Figure 2(s). The detail step of QSAR. (Supplementary Materials)