Novel Hit Compounds as Putative Antifungals: The Case of Aspergillus fumigatus

The prevalence of invasive fungal infections has been dramatically increased as the size of the immunocompromised population worldwide has grown. Aspergillus fumigatus is characterized as one of the most widespread and ubiquitous fungal pathogens. Among antifungal drugs, azoles have been the most widely used category for the treatment of fungal infections. However, increasingly, azole-resistant strains constitute a major problem to be faced. Towards this direction, our study focused on the identification of compounds bearing novel structural motifs which may evolve as a new class of antifungals. To fulfil this scope, a combination of in silico techniques and in vitro assays were implemented. Specifically, a ligand-based pharmacophore model was created and served as a 3D search query to screen the ZINC chemical database. Additionally, molecular docking and molecular dynamics simulations were used to improve the reliability and accuracy of virtual screening results. In total, eight compounds, bearing completely different chemical scaffolds from the commercially available azoles, were proposed and their antifungal activity was evaluated using in vitro assays. Results indicated that all tested compounds exhibit antifungal activity, especially compounds 1, 2, and 4, which presented the most promising minimum inhibitory concentration (MIC) and minimum fungicidal concentration (MFC) values and, therefore, could be subjected to further hit to lead optimization.


Introduction
In recent decades, fungal infections have been one of the most common and serious health problems worldwide. Their importance is highlighted by the fact that they are responsible for over one million human deaths per year [1]. Additionally, the prevalence of invasive fungal infections has dramatically increased as the size of the immunocompromised population grows. Today more than 5 million fungal species have been identified and, among them, approximately 300 species induce adverse impacts on human life [2].

Ligand-Based Pharmacophore Model Generation
Due to the absence of a crystal structure of Aspergillus fumigatus CYP51A enzyme, a ligand-based pharmacophore model was generated based on the common features of commercially available active compounds against the tested fungus. Particularly, the training set of 22 active compounds was comprised of five commercially available antifungal drugs from the drug class of azoles (clotrimazole, ketoconazole, oxiconazole, fluconazole, and voriconazole) and 17 active compounds against the tested fungus (1.0 nmol•mL −1 < MIC < 305 nmol•mL −1 ), bearing characteristic imidazole and triazole rings ( Figure 2, Figure S1, and Table S1 see Supplementary Materials). According to this set, the common chemical features of the compounds were identified, and the pharmacophore model was created. The test set included 20 commercially available active compounds ( Figure 3, Figure S2, and Table S2 see Supplementary Materials) against Aspergillus fumigatus with structural and functional diversity. More specifically, these compounds belong to two major classes. The first group includes azole drugs (miconazole, econazole, bifonazole) and the second one contains compounds (19 nmol•mL −1 < MIC < 350 nmol•mL −1 ) with substituted heterocyclic rings (imidazole, triazole, and tetrazole).

Ligand-Based Pharmacophore Model Generation
Due to the absence of a crystal structure of Aspergillus fumigatus CYP51A enzyme, a ligand-based pharmacophore model was generated based on the common features of commercially available active compounds against the tested fungus. Particularly, the training set of 22 active compounds was comprised of five commercially available antifungal drugs from the drug class of azoles (clotrimazole, ketoconazole, oxiconazole, fluconazole, and voriconazole) and 17 active compounds against the tested fungus (1.0 nmol·mL −1 < MIC < 305 nmol·mL −1 ), bearing characteristic imidazole and triazole rings ( Figure 2, Figure S1, and Table S1 see Supplementary Materials). According to this set, the common chemical features of the compounds were identified, and the pharmacophore model was created. The test set included 20 commercially available active compounds ( Figure 3, Figure S2, and Table S2 see Supplementary Materials) against Aspergillus fumigatus with structural and functional diversity. More specifically, these compounds belong to two major classes. The first group includes azole drugs (miconazole, econazole, bifonazole) and the second one contains compounds (19 nmol·mL −1 < MIC < 350 nmol·mL −1 ) with substituted heterocyclic rings (imidazole, triazole, and tetrazole).     Then, the generation and assessment of a set of 10 pharmacophore hypotheses were carried out. Their examination indicated that the best hypothesis included six common pharmacophore features and all the test set compounds were fitted on them. Particularly, two hydrogen bond acceptors (HBA), three hydrophobic regions (H), one aromatic ring (AR), and 29 exclusion volumes constituted the pharmacophore model ( Figure 4A). Among the compounds of the training set, voriconazole fitted optimally (Pharmacophore-Fit Score = 62.23) to the features of the initial model ( Figure 4B). Due to the large number of pharmacophore features (6), the present model was considered as improper for practical applications, like virtual screening [30]. Consequently, it was subjected to further optimization based on pharmacophore feature modifications. For this scope, two hydrophobic regions (H) were discarded, while the rest features were preserved in order to keep the major features of azoles, which could play a key role in the binding at the active site of the enzyme. The optimization of the model was completed by changing the size and number (from 29 to 15) of the exclusion volumes to increase the degrees of freedom of the ligands during screening. In Figure 4C,D, the pharmacophore features of the optimized model and the pharmacophore fit of fluconazole (Pharmacophore-Fit Score = 45.21) to the optimum model are illustrated, respectively. In Table 1 the calculated coordinates of the pharmacophore features and the tolerance radii (Å) of the optimum model was presented. Then, the generation and assessment of a set of 10 pharmacophore hypotheses were carried out. Their examination indicated that the best hypothesis included six common pharmacophore features and all the test set compounds were fitted on them. Particularly, two hydrogen bond acceptors (HBA), three hydrophobic regions (H), one aromatic ring (AR), and 29 exclusion volumes constituted the pharmacophore model ( Figure 4A). Among the compounds of the training set, voriconazole fitted optimally (Pharmacophore-Fit Score = 62.23) to the features of the initial model ( Figure 4B). Due to the large number of pharmacophore features (6), the present model was considered as improper for practical applications, like virtual screening [30]. Consequently, it was subjected to further optimization based on pharmacophore feature modifications. For this scope, two hydrophobic regions (H) were discarded, while the rest features were preserved in order to keep the major features of azoles, which could play a key role in the binding at the active site of the enzyme. The optimization of the model was completed by changing the size and number (from 29 to 15) of the exclusion volumes to increase the degrees of freedom of the ligands during screening. In Figure 4C,D, the pharmacophore features of the optimized model and the pharmacophore fit of fluconazole (Pharmacophore-Fit Score = 45.21) to the optimum model are illustrated, respectively. In Table 1 the calculated coordinates of the pharmacophore features and the tolerance radii (Å) of the optimum model was presented.

Pharmacophore Model Validation
In general, the reliability of a pharmacophore model depends on the number of biologically active compounds which are retrieved from a structurally diverse compound database and on the number of inactive compounds which are discarded from the same database. For the above scope, the validation database includes three different sets of compounds, including actives, inactives, and decoys. The creation of the decoys set was necessary since the number of available inactive compounds in the literature was inadequate for the study. The enrichment of active compounds is quantified by the hit-rate results from the validation database [31].
The optimized pharmacophore model was assessed by the receiver operating characteristic (ROC) curve analysis in order to identify its ability to correctly classify a list of compounds as actives or inactives. The quantification of ROC-curve analysis was determined by the area under the curve (AUC) of the examined ROC, as well as by sensitivity (Se), false positive rate (1-Sp), and enrichment factor (EF) values. The ROC-curve is illustrated in Figure 5 and the quantitative key parameters are presented in Table 2. The validation of the model was completed by calculating statistical significance variables ( Table 3). The ROC-curve ( Figure 5) analysis clearly indicated that the derived model has a good selection score and is certainly better than random selection (AUC= 0.69>0.5). It presents a steep slope at the beginning of screening, indicating a high enrichment of actives among the top-ranked hit list compounds. This fact is amplified by sensitivity (Se) and false positive rate (1-Sp) values (Table 2), as well as the enrichment factor value (EF= 2.4). Additionally, the reliability of the model was confirmed by the fact that it succeeded to retrieve almost 50% of the active compounds and the percentage of retrieved inactives and decoys was by far lower than the actives (Table 3). According to the above results, the generated model was used as a reliable "filter" in virtual screening.

Pharmacophore Model Validation
In general, the reliability of a pharmacophore model depends on the number of biologically active compounds which are retrieved from a structurally diverse compound database and on the number of inactive compounds which are discarded from the same database. For the above scope, the validation database includes three different sets of compounds, including actives, inactives, and decoys. The creation of the decoys set was necessary since the number of available inactive compounds in the literature was inadequate for the study. The enrichment of active compounds is quantified by the hit-rate results from the validation database [31].
The optimized pharmacophore model was assessed by the receiver operating characteristic (ROC) curve analysis in order to identify its ability to correctly classify a list of compounds as actives or inactives. The quantification of ROC-curve analysis was determined by the area under the curve (AUC) of the examined ROC, as well as by sensitivity (Se), false positive rate (1-Sp), and enrichment factor (EF) values. The ROC-curve is illustrated in Figure 5 and the quantitative key parameters are presented in Table 2. The validation of the model was completed by calculating statistical significance variables ( Table 3). The ROC-curve ( Figure 5) analysis clearly indicated that the derived model has a good selection score and is certainly better than random selection (AUC= 0.69>0.5). It presents a steep slope at the beginning of screening, indicating a high enrichment of actives among the top-ranked hit list compounds. This fact is amplified by sensitivity (Se) and false positive rate (1-Sp) values (Table  2), as well as the enrichment factor value (EF= 2.4). Additionally, the reliability of the model was confirmed by the fact that it succeeded to retrieve almost 50% of the active compounds and the percentage of retrieved inactives and decoys was by far lower than the actives (Table 3). According to the above results, the generated model was used as a reliable ʹʹfilterʹʹ in virtual screening.

Virtual Screening and Docking Studies
The library "Drugs Now", a subset of the ZINC database (http://zinc.docking.org/) [29], containing~11 million compounds which satisfy the drug-likeness criteria [32], was selected for the pharmacophore-based virtual screening [33]. Hits with the top-ranked Pharmacophore-Fit score were subjected to further filtering based on the drug-likeness values of the commercially available antifungal class of azoles (Table S3, see Supplementary Materials). Molecular docking studies were performed on compounds (~1 million) which passed the filtering criteria and showed the highest fit to the pharmacophore model features (Pharmacophore-Fit score range from 41 to 46). Especially, high-throughput docking algorithms were utilized to study the binding mode of the selected compounds to the homology model of Aspergillus fumigatus CYP51A enzyme [24]. The lack of the heme ring and any inhibitor from the active site of the homology model was treated by utilizing the respective molecules from the human CYP51A crystal structure (PDB code: 3JUS). This crystal structure includes the heme ring containing a ferric iron atom (Fe 3+ ) and the inhibitor R-econazole [34]. The created model was subjected to energy minimization and molecular dynamics simulations to obtain a stable and low energy model.
The final selection of the most promising compounds was based not only on docking scores but also on Pharmacophore-Fit score values and on the visual inspection focusing on the presence of interactions with crucial amino acids. The qualitative criteria of visual inspection were: (a) The potential metal coordination between the heme ferric iron and the candidate inhibitor [20,28], (distance < 3 Å); and (b) The formation of π-π interaction with Tyr107 and Tyr121 aminoacids [35]. Finally, 8 compounds bearing different chemical scaffolds in comparison to the azole class were selected and acquired for further antifungal activity evaluation ( Figure 6). The docking criteria (Glide Score and Emodel), the Pharmacophore-Fit score and the distance from the heme ferric iron of each of the selected compounds are presented in Table 4. The docked poses and the crucial interactions of the above compounds are illustrated in Figures 7 and 8. In Figure S3    Additionally, the docking analysis suggested that the selected compounds may coordinate the heme metal ion as in the case of azoles binding (Figures 8-10). Specifically, compounds 1 ( Figure 7) and 7 ( Figure 8) develop π-π stacking with the aromatic rings of Tyr107 and Tyr121, as in the case of R-econazole ( Figure 9) and both coordinate ferric ion with the sulfonyl group oxygen atom. In the case of compound 2, the oxygen atom of the nitro group forms a salt bridge with the heme moiety and metal coordination with the iron atom. Moreover, the benzofuran and nitrobenzene rings of the examined compound develop π-π interactions with the side chains of Tyr121 and His296, respectively ( Figure 7). The sulfonyl group of compound 3 is responsible, for the formation of metal coordination, as in the case of compound 1 and 7. Moreover, the halobenzene ring and the pyrimidine ring interact via π-π stacking with the aromatic rings of Tyr121 and His296, respectively. Additionally, its carbonyl group forms a hydrogen bond with the side chain of Thr111 (Figure 7). The oxygen atom of the sulfonyl group of compound 4 develops metal coordination with the heme iron atom, as in the case of compound 1, 3, and 7. Additionally, a salt bridge is created between the Natom and the heme ring. The methoxy benzene moiety of this compound interacts via a π-π stacking with Tyr107 and its carbonyl group forms two hydrogen bonds with the backbone of Val120 and with  Additionally, the docking analysis suggested that the selected compounds may coordinate the heme metal ion as in the case of azoles binding (Figures 8-10). Specifically, compounds 1 ( Figure 7) and 7 (Figure 8) develop π-π stacking with the aromatic rings of Tyr107 and Tyr121, as in the case of R-econazole ( Figure 9) and both coordinate ferric ion with the sulfonyl group oxygen atom. In the case of compound 2, the oxygen atom of the nitro group forms a salt bridge with the heme moiety and metal coordination with the iron atom. Moreover, the benzofuran and nitrobenzene rings of the examined compound develop π-π interactions with the side chains of Tyr121 and His296, respectively (Figure 7). The sulfonyl group of compound 3 is responsible, for the formation of metal coordination, as in the case of compound 1 and 7. Moreover, the halobenzene ring and the pyrimidine ring interact via π-π stacking with the aromatic rings of Tyr121 and His296, respectively. Additionally, its carbonyl group forms a hydrogen bond with the side chain of Thr111 (Figure 7). The oxygen atom of the sulfonyl group of compound 4 develops metal coordination with the heme iron atom, as in the case of compound 1, 3, and 7. Additionally, a salt bridge is created between the N-atom and the heme ring. The methoxy benzene moiety of this compound interacts via a π-π stacking with Tyr107 and its carbonyl group forms two hydrogen bonds with the backbone of Val120 and with the side chain of His285 (Figure 7). The furan and the thiazole ring of compound 5 form π-π interactions with the aromatic rings of Tyr121 and Phe214. The approach towards iron is induced via the nitrogen atom of the imidazole ring of this compound (Figure 8). In the case of compound 6, the carbonyl group coordinates with the heme iron atom. Furthermore, the o-xylene group forms π-π interactions with Tyr121 and Phe214, as in the case of compound 5. Additionally, the formation of a hydrogen bond between the nitrogen of triazole ring and the side chain of Ser363 stabilizes the binding (Figure 8). For compound 8, the oxygen atom of the furan ring is shown to interact with the heme iron atom. Its halobenzene ring makes π-π stacking with the aromatic rings of Tyr107 and Tyr121. Additionally, the methylbenzene group forms a π-π interaction with Tyr121. Moreover, a hydrogen bond is formed with the side chain of His285 and a π-cation interaction is presented with Lys132 ( Figure 8). The interactions of the selected compounds in the active site of Aspergillus fumigatus CYP51A modeled enzyme, are presented in Table 5.

Antifungal Activity Evaluation
In vitro assays were performed to determine the antifungal activity of the selected compounds, against different strains of Aspergillus fumigatus. The minimal inhibitory concentrations (MIC) and minimal fungicidal concentrations (MFC) were identified for the selected compounds compared to reference azole drugs, econazole, and ketoconazole (Table 6). Experiments were performed in duplicate and repeated three times. Values are expressed as means ± SD. In each column, different letters mean significant differences between samples (p < 0.05).
The in vitro results indicated that all tested compounds displayed antifungal activity against the two strains of Aspergillus fumigatus. Especially, in the case of the clinical isolate the MIC and MFC values ranging from 0.12 to 0.27 µmol·mL −1 and from 0.40 to 0.54 µmol·mL −1 , respectively. Additionally, the observed MIC (0.05-0.26 µmol·mL −1 ) and MFC (0.10-0.54 µmol·mL −1 ) values of the tested compounds against the ATCC204305 strain possessed higher antifungal activity compared to clinical isolate strain. It is obvious that compounds 1 and 2 presented the highest antifungal activity against the ATCC204305 strain, whereas compound 4 achieved the highest inhibitory effect in the case of clinical isolate Aspergillus fumigatus. It should be mentioned that compound 7 exhibited the lowest antifungal activity against both examined strains. Additionally, the observed MIC (0.05-0.26 μmol•mL −1 ) and MFC (0.10-0.54 μmol•mL −1 ) values of the tested compounds against the ATCC204305 strain possessed higher antifungal activity compared to clinical isolate strain. It is obvious that compounds 1 and 2 presented the highest antifungal activity against the ATCC204305 strain, whereas compound 4 achieved the highest inhibitory effect in the case of clinical isolate Aspergillus fumigatus. It should be mentioned that compound 7 exhibited the lowest antifungal activity against both examined strains.   Additionally, the observed MIC (0.05-0.26 μmol•mL −1 ) and MFC (0.10-0.54 μmol•mL −1 ) values of the tested compounds against the ATCC204305 strain possessed higher antifungal activity compared to clinical isolate strain. It is obvious that compounds 1 and 2 presented the highest antifungal activity against the ATCC204305 strain, whereas compound 4 achieved the highest inhibitory effect in the case of clinical isolate Aspergillus fumigatus. It should be mentioned that compound 7 exhibited the lowest antifungal activity against both examined strains.    A general observation is that the MIC and MFC values of the hit compounds were remarkably different from those of econazole compared to the predictions of the molecular docking. Since the in vitro testing of compounds was performed directly to fungal strains (and not to CYP51 protein determining binding affinities) there may be several reasons explaining the different activities. These may include fungal cell wall permeability, cell membrane permeability, solubility of tested compounds, etc. Indeed, Table S4 presents the predicted values of physicochemical parameters such as lipophilicity and polar surface area. From these results we may notice that econazole is more lipophilic with a significantly lower polar surface area (27.05) compared to tested compounds (109.95 for compound 4 to 136.18 for compound 2) which may constitute a putative explanation.

Molecular Dynamics Simulations
In order to evaluate the stability of compounds 1, 2, and 4, compared to R-econazole, the docking poses of the four molecules were subjected to unconstrained molecular dynamics simulations for 1 μs each. In all cases, the RMSD of the structural model backbone was reasonable ( Figure 10A), as well as the ligand RMSD ( Figure 10B).
To further evaluate the interactions which contribute to complex stability, the averaged distances were measured from the molecular dynamics simulations between the heme iron and the center of mass of the two sulfonyl oxygens of compound 1 (cyan), the center of mass of the two nitro group oxygens of compound 2 (yellow),the center of mass of the two sulfonyl oxygens of compound 4 (magenta), and the imidazole nitrogen of R-econazole (orange).
Results indicate ( Figure 10C) that the imidazole nitrogen of R-econazole, the sulfonyl group oxygens of compounds 1 and 4 and nitro group oxygens of compound 2 were found in constant proximity to Fe 3+ . Although these interactions with heme iron remain stable for all tested compounds as shown in the distance profiles along their trajectories ( Figure 10C), we postulate that the forces maintaining such a stable binding profile are hydrophobic. This could be explained by the fact that the RMSD of these three compounds was stable along the trajectories as shown in Figure 10B.
Moreover, the three compounds under investigation share common binding characteristics. Reconazole forms steady hydrophobic contacts between the two-chlorine bearing aromatic rings with mainly Thr289, Leu110, and Val120 and the imidazole ring is stabilized by Ala293 ( Figure 10D). Compound 1, as a larger in size compound, besides hydrophobic contacts between its aromatic groups and mainly Thr289, Leu110, and Val120, bears an extra benzo-dioxole group, which penetrates deeper to interact with the two prolines Pro394 and Pro446, as well as Ile364 and Ser362 ( Figure 10E). Compound 2 forms similar hydrophobic interactions with Val120, Tyr121, as well as Tyr107. The benzyl group attached to the Fe interacting nitro group is also stabilized by Ala239 ( Figure 10F). Compound 4 has a similar interaction profile, forming aromatic-hydrophobic interactions with Thr289, Val120, and Ile364, and hydrophobic interactions with Met286 ( Figure 10G). A general observation is that the MIC and MFC values of the hit compounds were remarkably different from those of econazole compared to the predictions of the molecular docking. Since the in vitro testing of compounds was performed directly to fungal strains (and not to CYP51 protein determining binding affinities) there may be several reasons explaining the different activities. These may include fungal cell wall permeability, cell membrane permeability, solubility of tested compounds, etc. Indeed, Table S4 presents the predicted values of physicochemical parameters such as lipophilicity and polar surface area. From these results we may notice that econazole is more lipophilic with a significantly lower polar surface area (27.05) compared to tested compounds (109.95 for compound 4 to 136.18 for compound 2) which may constitute a putative explanation.

Molecular Dynamics Simulations
In order to evaluate the stability of compounds 1, 2, and 4, compared to R-econazole, the docking poses of the four molecules were subjected to unconstrained molecular dynamics simulations for 1 µs each. In all cases, the RMSD of the structural model backbone was reasonable ( Figure 10A), as well as the ligand RMSD ( Figure 10B).
To further evaluate the interactions which contribute to complex stability, the averaged distances were measured from the molecular dynamics simulations between the heme iron and the center of mass of the two sulfonyl oxygens of compound 1 (cyan), the center of mass of the two nitro group oxygens of compound 2 (yellow),the center of mass of the two sulfonyl oxygens of compound 4 (magenta), and the imidazole nitrogen of R-econazole (orange).
Results indicate ( Figure 10C) that the imidazole nitrogen of R-econazole, the sulfonyl group oxygens of compounds 1 and 4 and nitro group oxygens of compound 2 were found in constant proximity to Fe 3+ . Although these interactions with heme iron remain stable for all tested compounds as shown in the distance profiles along their trajectories ( Figure 10C), we postulate that the forces maintaining such a stable binding profile are hydrophobic. This could be explained by the fact that the RMSD of these three compounds was stable along the trajectories as shown in Figure 10B.
Moreover, the three compounds under investigation share common binding characteristics. R-econazole forms steady hydrophobic contacts between the two-chlorine bearing aromatic rings with mainly Thr289, Leu110, and Val120 and the imidazole ring is stabilized by Ala293 ( Figure 10D). Compound 1, as a larger in size compound, besides hydrophobic contacts between its aromatic groups and mainly Thr289, Leu110, and Val120, bears an extra benzo-dioxole group, which penetrates deeper to interact with the two prolines Pro394 and Pro446, as well as Ile364 and Ser362 ( Figure 10E). Compound 2 forms similar hydrophobic interactions with Val120, Tyr121, as well as Tyr107. The benzyl group attached to the Fe interacting nitro group is also stabilized by Ala239 ( Figure 10F). Compound 4 has a similar interaction profile, forming aromatic-hydrophobic interactions with Thr289, Val120, and Ile364, and hydrophobic interactions with Met286 ( Figure 10G).

Pharmacophore Model Generation
For the creation and the validation of the pharmacophore model, LigandScout 4.0 Advanced software was used (InteLigand, GmbH, Vienna, Austria [36,37]. Due to the absence of a crystal structure of fungi CYP51A, a ligand-based pharmacophore approach was applied. The compounds with known biological activity against Aspergillusfumigatus were extracted from the ChEMBL database [38] and were utilized for the creation of the training and the test set. The training set included five representative drugs from the class of azoles and 17 active compounds (1 nmol•mL −1 ≤ MIC ≤ 310 nmol•mL −1 ) with chemical diversity (Figure 2). According to this set, the common chemical features of the compounds were aligned, and the model was created. The test set contained 3 commercially available antifungal drugs (miconazole, econazole, and bifonazole) and 19 active compounds (19 nmol•mL −1 ≤ MIC ≤ 350 nmol•mL −1 ) (Figure 3). All compounds were converted from

Pharmacophore Model Generation
For the creation and the validation of the pharmacophore model, LigandScout 4.0 Advanced software was used (InteLigand, GmbH, Vienna, Austria [36,37]. Due to the absence of a crystal structure of fungi CYP51A, a ligand-based pharmacophore approach was applied. The compounds with known biological activity against Aspergillusfumigatus were extracted from the ChEMBL database [38] and were utilized for the creation of the training and the test set. The training set included five representative drugs from the class of azoles and 17 active compounds (1 nmol·mL −1 ≤ MIC ≤ 310 nmol·mL −1 ) with chemical diversity (Figure 2). According to this set, the common chemical features of the compounds were aligned, and the model was created. The test set contained 3 commercially available antifungal drugs (miconazole, econazole, and bifonazole) and 19 active compounds (19 nmol·mL −1 ≤ MIC ≤ 350 nmol·mL −1 ) (Figure 3). All compounds were converted from 2D to 3D and then were minimized in Maestro 2013 (Schrodinger, LLC, New York, NY, USA) software. Training set compounds were prepared in pH 8.0 ± 0.5 using the LigPrep module of MAESTRO software (Schrodinger, LLC). In the case of the test set, all compounds were prepared to retain their chiralities. The chemical structures and activities of training and test sets are presented in Tables S1 and S2 (see Supplementary Materials), respectively. The OMEGA algorithm was applied for the generation of conformations of the ligand set [39] and the maximum number of conformations per ligand was set to 50. Then, the ligand set was clustered according to the pharmacophore alignment score and 10 pharmacophore hypotheses were created. The performed scoring function was the "Pharmacophore-Fit and Atom Overlap", the selected pharmacophore type was the "merged" pharmacophore features, and the number of omitted pharmacophore features was 6 ( Figure 4).

Evaluation of Validation Results
For the validation process, three libraries were created by using the ChEMBL [38] database. The first one included 916 active compounds, the second, 117 inactive compounds and, due to the small number of inactive compounds available in literature, adecoy library was also created, consisting of 3095 compounds structurally similar to the active ones, but experimentally not tested for biological activity. The creation of this library was based on the application of an in-house protocol.
Sensitivity (Se) determines the percentage of truly active compounds selected during virtual screening and it is defined as: Se = Selected Actives/All Actives (1) Specificity (Sp) describes the percentage of truly inactive compounds rejected during virtual screening and it is defined as: Sp = Discarded (Inactives + Decoys)/All (Inactives + Decoys) Enrichment factor measures how much better the method is than a random selection by considering the fraction of found actives compounds in the top of the ranked hit list, compared to the ration of compounds in the entire database: The ROC-curve is a graphical plot of Se versus 1-Sp and the area under the curve (AUC) is measured by:

Pharmacophore-Based Virtual Screening
Pharmacophore-based virtual screening was applied to the ZINC database (http://zinc.docking. org/) subset "Drugs Now", which includes approximately 11 million compounds. The idbgen-tool of LigandScout [36,37] was used to convert this pool of compounds to the appropriate screening format. Twenty-five conformers per compound were created by the OMEGA algorithm [39].
The higher pharmacophore-fit score values (43)(44)(45)(46) were utilized as the main filter of the hit list. The next filter was based on the values of the drug-likeness criteria of the azole drug class (molecular weight, lipophilicity, number of hydrogen-bond donors, number of hydrogen-bond acceptors, polar surface area and number of rotatable bonds) (Table S3, see Supplementary Materials). All drugs were prepared in the optimum pH = 8.0 ± 0.5 by the LigPrep [41] version 2.7 (Schrodinger, LLC) module and their physicochemical properties were calculated by the QikProp version 3.9 module of MAESTRO (Schrodinger, LLC). One out of nine million compounds which were screened fit to the pharmacophore model features and satisfy the drug-likeness criteria of azoles.

Molecular Docking Studies
The derived compounds were subjected to further filtering through molecular docking studies at the binding site of the Aspergillus fumigatus CYP51A constructed model [24]. This model was based on the human CYP51A model ortholog, providing a high-quality model. The heme ring and the inhibitor econazole (R) were transferred from the human crystal structure (PDB: 3JUS) and were aligned to the constructed model. The hits were docked with the Standard Precision mode of Glide, version 6.0 (Schrodinger, LLC) [42,43]. The resulting poses (5000 compounds) were visually inspected and the most promising ones were further docked with the Induced-Fit protocol 2013-2, Glide, version 5.9, Prime, version 3.2 (Schrodinger, LLC) [44].

Antifungal ActivityEvaluation
For the antifungal bioassays, two strains of Aspergillus fumigatus were used: Aspergillus fumigatus (ATCC 204305), and Aspergillus fumigatus (human clinical isolate). The organisms are deposited at the Mycological Laboratory, Department of Plant Physiology, Institute for Biological Research "Siniša Stankovic," Belgrade, Serbia.
The micromycetes were maintained on malt agar and the cultures were stored at 4 • C and sub-cultured once a month. The antifungal assay was carried out by a modified microdilution technique [45,46] in order to determine the minimum inhibitory (MIC) and minimum fungicidal concentrations (MFC) of the examined compounds. Briefly, the fungal spores were washed from the surface of agar plates with sterile 0.85% saline containing 0.1% Tween 80 (v/v). The spore suspension was adjusted with sterile saline to a concentration of approximately 1.0×10 5 in a final volume of 100 µL per well. The examined compounds were dissolved in 5% of DMSO, serially diluted in broth malt medium after which fungal inoculum was added. The microplates were incubated for 72 h at 28 • C. The lowest concentrations without visible growth (in a binocular microscope) were defined as MICs. The fungicidal concentrations (MFCs) were determined by serial sub cultivation of 2 µL of the wells content into microtiter plates containing 100 µL of broth per well and further incubation for 72 h at 28 • C. The lowest concentration with no visible growth was defined as MFC indicating 99.5% killing of the original inoculum. The commercial antifungals econazole and ketoconazole were used as positive controls.

Molecular Dynamics
All MD simulations were performed using GROMACS 2018 v3 (University of Groningen& Royal Institute of Technology, Uppsala University, Groningen, Netherlands) [47]. The structural model, described by Fraczeket al. [24] in complex with heme (covalently) and compounds R-econazole, compound 1, compound 2, and compound 4, in each case as a result of the docking previously described, was inserted in a pre-equilibrated box containing water and a 0.15M concentration of Na and Cl ions. The latest AMBER99SB-ILDN [48] force field was used for all the dynamics simulations along with the TIP3P water model. Force field parameters for the ligands were generated using the general Amber force field (GAFF) and HF/6-31G*-derived RESP atomic charges [49]. AMBER-compatible heme parameters were derived from elsewhere [50]. Each system consisted of the protein, the ligand, 15,000water molecules, and~160 ions in a 10 × 10 × 10nm simulation box. The model systems were energy minimized and subsequently subjected to a 20 ns MD equilibration, with positional restraints on protein Cα atoms. These restraints were released and 1000 ns MD trajectories were produced at a constant temperature of 300K using separate v-rescale thermostats for the protein, the ligand, and solvent molecules. A time step of 2 fs was used, and all bonds were constrained using the LINCS algorithm. Lennard-Jones interactions were computed using a cutoff of 10 Å, and the electrostatic interactions were treated using PME with the same real-space cutoff as described before [51].

Conclusions
In the present study, a combination of in silico techniques and methodologies were implemented to identify novel scaffolds with antifungal activity against Aspergillus fumigatus targeting the fungal cytochrome P450 dependent lanosterol 14-α demethylase (CYP51A) enzyme. For this scope, a carefully designed and validated ligand-based pharmacophore model was used to screen the ZINC database. Compounds that passed the pharmacophore filtering were further examined using molecular docking studies. In total, eight compounds were selected, six of which have different azole structural motifs. The MIC and MFC values from the antifungal activity tests showed that all compounds present antifungal activity in comparison to commercially available antifungals and could serve as starting materials for further modifications. Especially, compounds 1, 2, and 4 exhibited the highest activity and were further studied to explore their mechanism of action on the CYP51A enzyme by means of molecular dynamics simulations. All three compounds were found to share common binding characteristics. Importantly, the imidazole nitrogen of R-econazole, the sulfonyl group oxygens of compounds 1 and 4, and nitro group oxygens of compound 2, were found in constant proximity to Fe 3+ . The role of hydrophobic interactions with amino acids Thr289, Leu110, Val120, Tyr107, and Tyr121to the stability of the complexes was also observed.
In conclusion, compounds 1, 2, and 4 are the most promising scaffolds that differentiate from classical azoles and could serve as starting points for further hit-to-lead optimization.