Discovery of inhibitors against SARS-CoV-2 associated fungal coinfections via virtual screening, ADMET evaluation, PASS, molecular docking, dynamics and pharmacophore studies

Abstract The black fungus is a SARS-CoV-2-associated fungal co-infection. Mechanisms by which antifungal drugs act against Mucormycosis are by inhibition of enzymes implicated in ergosterol biosynthesis or cell wall formation or the drug interaction directly with ergosterol. Since it has been shown that the reduction of ergosterol biosynthesis prevents this infection, we carried out a virtual screening with molecular docking analysis and tested a dataset of 880 molecules for its effects on an enzyme involved in ergosterol biosynthesis (CYP51) and fungal infection (lipase), as well as ADMET and toxicity properties of the studied molecules. We calculated the stability of the complex’s inhibitor-enzyme throughout the time using molecular dynamics simulations. From this, initial screening indicated that 600 out of a total of 880 compounds have a good affinity towards sterol 14α-demethylase. The two compounds ZINC95486139 and ZINC33833639 were determined as safe and have good drug-likeness properties, and have binding energies of −75.07 and −83.89 kcal/mol, respectively by MMGBSA and create strong and stable interactions with the residues of the active site of the enzyme. To the best of our knowledge, ZINC95486139 and ZINC33833639 have not yet been described in the literature as anti-sterol 14α-demethylase inhibitors or evaluated for their anti-fungal activity.

Mucormycosis is a rare fatal disease that impacts humans (Elfiky, 2019;Kazak et al., 2013), correlated with high morbidity and mortality ). It appears mainly in patients with malignancies, uncontrolled diabetes, undergoing organ transplantation, it also appears in immunocompromised patients due to the immunosuppressive drugs use and patients affected with COVID-19 or recovering from it (Elfiky, 2019;Kazak et al., 2013;Rahman, Islam, & Bhuiyan, 2021). Mucormycosis has mortality rates range from 40% to 100% depending on the associated diseases, sites of infection, and underlying conditions Elfiky, 2019;Kazak et al., 2013;Roden et al., 2005), reaching 100% in haematological malignancies (Elfiky, 2019;Lee et al., 1999). Rhinocerebral mucormycosis cases are the most cases of mucormycosis, however, pulmonary, gastrointestinal, cutaneous, and disseminated infections are probably as well observed (Kazak et al., 2013). The mortality rate increases with the delay in the initiation of the therapy . Besides the surgical techniques, antifungal drugs are also considered one of the main available strategies to treat mucormycosis . Currently, the antifungal drugs with anti-Mucorales activity are restricted to triazole posaconazole, amphotericin B, and isavuconazole (Chakrabarti & Singh, 2020;Cornely et al., 2014;Marty et al., 2016;Schwarz, Cornely, & Dannaoui, 2019). Amphotericin B deoxycholate is the most chosen drug for decades (Chakrabarti et al., 2009;Pagano et al., 2004;Roden et al., 2005;Skiada et al., 2011), Although, it is effective against Mucorales, its usages are limited due to its toxicity Ullmann et al., 2006;Walsh et al., 1999). Terbinafine and azoles (in particular fluconazole) are other antifungals drugs, that block the synthesis of ergosterol, by inhibiting the fungal squalene epoxidase and lanosterol 14a-demethylase, respectively (Ryder, 1992;Schwarz et al., 2019;Van Daele et al., 2019). Despite the uses of azoles as available treatments, it has been observed a Mucorales inherent resistance against short-tailed azoles like fluconazole, this resistance is related to an evolutionary substitution of amino acid in the target enzyme (14a-demethylase) (Chakrabarti & Singh, 2020;Schwarz et al., 2019). Other mutations have been reported in FKS1 and FKS2 genes that coded to the enzyme 1.3b-glucan synthetase, which lead to resistance to echinocandins (da N obrega Alves et al., 2020;del Valle, 2015). Lanosterol 14a-demethylase (microsomal P450 (P45014DM)) is a key enzyme in the process of ergosterol production (Desai et al., 2016;Pasqualotto & Falci, 2013;Van Daele et al., 2019), the important sterol of the cell membrane of the fungus (Van Daele et al., 2019). It is confirmed that blocking this enzyme is an effective strategy for the elimination of viral infections.
Secondary metabolites are organic compounds produced by living organisms such as bacteria, fungi, plants, and animals (James, 2017). To date, about 500,000 secondary metabolites have been described (Bills & Gloer, 2016). In plants, they can be found in the leaves, flowers, stems, roots, or bark (Anulika, Ignatius, Raymond, Osasere, & Hilda, 2016). Terpenes, phenolic compounds, and nitrogen-containing chemicals like alkaloids are the three primary types of secondary metabolites found in plants (Anulika et al., 2016;Guerriero et al., 2018). Secondary metabolites play a crucial role in powerful and effective defense systems, plant secondary metabolites have long been used in traditional medicine due to their significant biological activity (Bourgaud, Gravot, Milesi, & Gontier, 2001), and their potent defense mechanisms (Bennett & Wallsgrove, 1994). There are many published papers, patents, and patent applications on natural products used for patents in antifungal agent discovery (Aldholmi, Marchand, Ourliac-Garnier, Le Pape, & Ganesan, 2019;Heard, Wu, & Winter, 2021;Negri et al., 2014;Santo, 2008). This study aims to discover new natural molecules that have anti-fungal effects to develop drugs against fungal infection.

Construction of molecules dataset and ligands' preparation
The list and structures of approved drugs were obtained from AfroDb a database of natural products from African medicinal plants: https://zinc12. docking.org/catalogs/afronp) (Ntie-Kang et al., 2013). The version of the database downloaded on 27 July 2021 contained 880 purchasable compounds. Gasteiger charges were added to each molecule, and structural optimization was carried out using the MMFF94 force field and Obabel. The structures were then converted to pdbqt format for use in AutoDock Vina docking simulations.

Enzyme preparation
We have chosen two enzymes that are involved in fungal pathology function and were admitted as pharmacological targets for the development of new anti-fungal drugs. Those enzymes are sterol 14ademethylase and Lipase. Virtual screening was implemented by exploiting the published crystal structure of the enzyme's sterol 14a-demethylase and Lipase from Candida albicans. The PDB files of the targets (PDB IDs: 5FSA and 1LPP) were obtained from the Protein Data Bank (PDB) (Berman, Henrick, & Nakamura, 2003). The targets were carefully prepared using AutoDock tools (ADT) (version 1.5.4) (Morris et al., 2009) as follows: Water and other nonspecific molecules were removed. The proteins were protonated by adding polar hydrogens. The co-crystallized ligand was used as a centre of the docking grid box of 40 Ã 40 Ã 40 Å size, the centre of the rectangular box has set 195.506 Ã -5.037 Ã 38.928 (Serseg, Benarous, & Yousfi, 2020). All heteroatoms, ligands, co-crystallized solvent, and water molecules were removed. Partial charges and polar hydrogens were added to the structure.

Virtual screening
AutoDock Vina was used for the initial virtual screening. The African natural products from AfroDb were docked into the active site cavity of sterol 14ademethylase (5FSA) and Lipase (1LPP). According to the docking score, the output results were ranked. The obtained docking results were analyzed to include a comparison with fluconazole. All compounds with predicated DG less than -9Kcal/mol were selected to evaluate their ADMET properties.

ADMET and drug-likeness evaluation
The vast majority of drug candidates fail due to their toxicity. ADME/Tox research plays an important role in drug candidate success (Li, 2001). In order to select the ideal candidates to identify drug-like compounds ADME/Tox properties have been evaluated. The drug-likeness and toxicity prediction of the studied compounds have been carried out using web servers: admetSAR 2.0 (http://lmmd.ecust.edu. cn/admetsar2/) (Yang et al., 2019) and SwissADME server (http://www.swissadme.ch/index.php) (Daina, Michielin, & Zoete, 2017). The canonical SMILES of the compounds used for ADMET properties prediction were obtained from PubChem Database (Kim et al., 2021). The compounds that follow Lipinski's 'rule of five' and demonstrate good predicted ADME profiles were selected to conduct molecular docking experiments.

Molecular docking simulation
Medicinal chemistry research frequently uses the molecular docking methods (Ram ırez & Caballero, 2018), it is used for predicting the affinity and binding mode of small molecules to targets and calculate the binding energy of complexes (Moraga-Nicol as et al., 2018;Taghipour, Zakariazadeh, Sharifi, Ezzati Nazhad Dolatabadi, & Barzegar, 2018). A further docking experiment was done, using GOLD (Genetic Optimization for Ligand Docking) program (Alvarez & Shoichet, 2005;Jones, Willett, Glen, Leach, & Taylor, 1997) to predict the selectivity of the selected molecules. GOLD program uses a genetic algorithm for docking and performs an automated docking with partial flexibility of protein and ligand (John, Thangapandian, Sakkiah, & Lee, 2011). The water molecules were removed, the co-crystallized inhibitor posaconazole was removed from the binding site and the selected molecules were docked into the active site of sterol 14a-demethylase. The early termination option was set to 3 with a maximum of 10 conformations for each molecule, The GoldScore fitness function was used to evaluate the obtained conformations. The standard set parameters were used in all the calculations. The top conformation of each molecule was saved for analysis of the prediction of the binding mode (Karatas et al., 2017). Finally, the generated files were visualized and analyzed using Discovery Studio visualizer, v 4.0.

Molecular dynamics simulation
We have conducted a molecular dynamics simulation to analyze the stability of ligand-protein interactions in the complex and check its dynamic features (Haq, Abro, Raza, Liedl, & Azam, 2017). The ligand-protein complexes resulting from molecular docking were used in this simulation. Ligand molecules parameters were generated using CGenFF web server (Vanommeslaeghe et al., 2009). GROMACS 2021 (Abraham et al., 2015) was used for the molecular dynamics simulation of ligandprotein complexes CYP51-ZINC95486139 and CYP51-ZINC33833639. Then the system was neutralized by adding Clions and submerged into a rectangular box containing three-point transferable intermolecular potential (TIP3P) water molecules, the energy was minimized, and the system was heated gradually to 300 K at 1 atm pressure (NPT conditions) for 100 picoseconds. Finally, The MD simulation was run for 110 ns, and the recording was made every 10 ps (Lemkul, 2019;Linani, Benarous, Bou-Salah, Yousfi, & Goumri-Said, 2022;Lindahl, Hess, & Spoel, 2022). The resulting trajectories were analyzed, then and plotted using XMGRACE (Turner, 2005) to evaluate the root mean square fluctuation (RMSF), root mean square deviation (RMSD), intermolecular hydrogen bond interactions, and radius of gyration (Rg) for the elucidation of structural, conformational and compactness of the ligand protein.

Biological activity prediction (PASS)
In order to check the possible biological activities for ZINC95486139 and ZINC33833639 the best hits, we have used the PASS web server (http://way2drug. com/passonline) (Goel, Singh, Lagunin, & Poroikov, 2011). PASS stands for Prediction of Activity Spectra for Substances . This server predicts the different pharmacological effects of different compounds including multilevel neighbourhoods of atoms (MNA) descriptors . This method is based on structural activity relationship (SAR) analysis of a training set containing more than 205,000 compounds with more than 3750 types of biological activities (Garg et al., 2021;Goel et al., 2011;Jamkhande, Wattamwar, Pekamwar, & Chandak, 2014;Khurana, Ishar, Gajbhiye, & Goel, 2011). It also predicts the toxicity of these molecules. The results are expressed as Pa (probable activity) and Pi (probable inactivity) (Garg et al., 2021;Goel et al., 2011;Jamkhande et al., 2014;Khurana et al., 2011), If the probability Pa value is more than Pi of a compound, then it is believed to be potential for a special pharmacological activity (Garg et al., 2021;Goel et al., 2011;Jamkhande et al., 2014;Khurana et al., 2011). The canonical SMILES of the two investigated compounds are extracted from ZINC database (Kim et al., 2021, Kawsar et al., 2022.

Pharmacophore study
Energy-optimized pharmacophore (e-pharmacophore) model development can consider a target active site's inherent flexibility. According to reports, pharmacophores with high energy efficiency allow us to screen millions of compounds in a short period of time.
The e-Pharmacophore approaches combine ligand properties with structure-based methodologies to investigate the intrinsic flexibility of residues of binding pockets and ligand binding (Salam, Nuti, & Sherman, 2009). An E-pharmacophore model was constructed using the protein-ligand complexes from the Phase module. The Phase module defaulted to the negatively ionizable region (N), aromatic ring (R), hydrogen bond donor (D), hydrogen bond acceptor (A), hydrophobic group (H), and positively ionizable (P) pharmacophore attributes.

Binding free energies calculation (MM-GBSA)
The molecular mechanics-generalized born surface area (MM-GBSA) technique was used to quantify the binding affinity between the protein and small ligands. The theoretical binding free energies of the prospective inhibitors were calculated using Schr€ odinger's prime module. The MM-GBSA (Molecular Mechanical-Generalized Born Surface Area) study was performed on the hits received as a result of docking (Cappel et al., 2016).

Virtual screening
A subset of 880 small molecules initially selected from AfroDb database was used in a virtual docking screening targeting the binding site of sterol 14ademethylase and lipase.
The virtual screening was carried out using AutoDock Vina. The results with DG >-8.0 Kcal/mol were filtered out, keeping 600 compounds (Figure 1; Supplementary Table 1). ZINC95486182 is predicted to have the best binding energy (-13,4 kcal/mol) among the natural products classed in the Africa database (880 molecules), followed by ZINC95486333, ZINC13411589 respectively. While the score for ZINC95485962, ZINC95486139, and ZINC33833639 was predicted to be À11.5, À11.4, and À10.8 Kcal/mol, respectively. The only poses with DG <-8.0 Kcal/mol were considered to have a good affinity towards the enzyme in comparison with the approved drug fluconazole (-7.3 Kcal/mol). The selected molecules were evaluated using ADMET and drug-likeness servers.

ADMET and drug likeness evaluation
ADMET and drug-likeness evaluation predict the properties of a drug based on its chemical structure. This method allows the selection of drugs with good pharmacological properties such as toxicological potential, drug-drug interaction potential, metabolic stability, and intestinal permeability. This method reduces unnecessary time and costs in drug development.
The oral bioavailability of a drug is the most property evaluated, as it is the most common route of medication rout. Therefore, it is important to check the absorption effectively through the intestinal epithelium of the developed drugs, the used model to evaluate this property is Caco-2, a derived cell from human colon adenocarcinoma. In addition, the Inhibitor-Substrate property, which is called Drugdrug interaction, the metabolic stability of these drugs could be affected if there is an interaction between two drugs. Drug toxicity is another crucially important drug property, which generally examines hepatotoxicity of the developed drug (Li, 2001).
Lipinski, Lombardo, Dominy, & Feeney, 1997 proposed the 'Rule of Five,' which is a filter that predicts if a molecule is orally absorbed well or not. This rule provided that, a molecule is could not be orally active if it violates two or more of the following criteria: (i) Molecular weight (MW) should be less or equal to 500 Da. (ii) Octanol/water partition coefficient (A Log P) should be less than or equal to 5. (iii) The amount of hydrogen bond donors (HBDs) should be less than or equal to 5. (iv) The number of hydrogen bond acceptors (HBAs) should be less than or equal to 10 (Lipinski et al., 1997). As shown in Table S2, only a few molecules (59 of 600 molecules) passed all drug-likeness criteria successfully and are predicted to be orally bioavailable, they respect Lipinski rules, have high human intestinal absorption, they are not toxic nor carcinogens or Table 1. Docking results of natural products from AfroDb database in the active pocket of the target enzyme (sterol 14ademethylase) using GOLD program.

Molecular docking
Molecular docking study is a computational simulation that predicts the binding affinity of a small molecule to a protein target, using scoring functions (Jiao et al., 2021). In this work, the inhibitory effect of the selected molecules which passed through the ADMET and toxicity evaluation successfully with sterol 14a-demethylase was investigated by docking experiments and compared with anti-fungal approved drugs. The results are presented in Table 1 and Figure 2. We have chosen, only two molecules to discuss their binding mode towards 5FSA: ZINC95486139 and ZINC33833639. Molecule ZINC95486139 formed 18 hydrophobic interactions with Tyr118, Leu121, Phe228, Phe233, Leu376, Met508, and Hem group of type Alkyl, Pi-Alkyl, and Pi-Sigma. The alkyl groups and rings in this molecule are responsible for the hydrophobic interactions with the amino acids of the active pocket, where we observe that the benzene groups participate in the formation of 6 hydrophobic interactions with amino acids Try118, Leu121, Leu376, and Met508. We do not observe any hydrogen bonds between ZINC95486139 and the amino acids of the sterol 14a-demethylase pocket ( Figure  3a,b). While ZINC33833639 forms 5 hydrogen bonds with residues: His377 (2.97 Å), His377 (2.19 Å), Ser378 (1.96 Å), Ser378 (1.69 Å), and Met508 (1.62 Å). These hydrogen bonds are supported by 17 hydrophobic interactions of type Alkyl, Pi-Alkyl, and Pi-Lone pair, with Tyr64, Leu87, Leu88, Pro 230, Phe233, Leu376, and Phe380. ZINC33833639 is tied in the rim of the binding pocket and does not form any interactions with the Hem group, even though, it could block the rim of the binding pocket impeding the substrate from reaching the catalytic site. From the docking results, both ZINC95486139 and ZINC33833639 show good affinities À11.4 and 10.8 Kcal/mol, respectively, towards sterol 14a- demethylase, and form a high number of interactions with the amino acids of the active pocket, in comparison with the reference molecule fluconazole, the current approved anti-fungal drug, where the affinity of fluconazole was À7.4 Kcal/mol and forms only seven hydrophobic interactions and one hydrogen bond with Tyr64 (2.84 Å). According to these observations, we suggest that ZINC95486139 and ZINC33833639 are good inhibitors for sterol 14ademethylase, and could be promising molecules to develop new drugs against black fungus.

Molecular dynamics simulation
Molecular dynamics simulation was carried out to evaluate the stability of the most potent compounds in the binding pocket of sterol 14a-demethylase, and the effect of these compounds on the enzyme structure. Accordingly, molecular dynamics simulation (110 ns) was run for the following complexes 5FSA-ZINC95486139 and 5FSA-ZINC33833639 ( Figure  4) using Gromacs2021. The trajectory files were analyzed for Root-mean-square deviation (RMSD), Root-mean-square fluctuation (RMSF), and radius of gyration (Rg). The graphs related to RMSD, RMSF, and Rg properties are commonly used in the evaluation of the complex protein-ligand structural stability (Falsafi-Zadeh, Karimi, & Galehdari, 2012).
3.4.1. RMSD ZINC95486139 and ZINC33833639 exhibited a stable binding to sterol 14a-demethylase, as confirmed by recording minor changes in the positional rootmean-square deviation (RMSD) of both the ligands and residue backbone of the protein (Figure 4a). RMSD of both ZINC95486139 and ZINC33833639 was lower than 0.8 Å, which means that they are stable throughout the 100 ns simulation. Whereas RMSD of the backbone increased gradually in the first 40 ns, then remained stable at 3 Å till the end of the simulation in both cases (Figure 4b). Thus, the two compounds were firmly fixed to the active site of sterol 14a-demethylase, with the hydrophobic pocket formed by Tyr64,Leu87,Leu88,Tyr118,Leu121,Phe228,Pro230,Phe233,Leu376,Phe380,and Met508, as confirmed by molecular docking study. The backbone RMSD plot suggests that the (sterol 14a-demethylase þ ZINC95486139) and (sterol 14a-   demethylase þ ZINC33833639) show the same affinity and behaviour. Thus, these natural molecules do not disturb the structural stability of the studied enzyme and remained tightly bound to sterol 14ademethylase active site.

RMSF
RMSF analysis of sterol 14a-demethylase (5FSA) backbone shows a notable fluctuation at residues 329 to 341 that gives a peak at 9 Å, while the rest of the residues recorded low fluctuations less than 2 Å (Figure 4c), the high values of RMSF of the residues 329 to 341 correspond to residues that belong to a loop against the active site. These residues do not participate in the ligand binding, so we think that they do not have any effect on the stability of the inhibitor on the binding pocket of the studied enzyme. It has been confirmed that conformational changes can often be observed in the loops (Lee, Shin, Park, Shin, & Seok, 2012;Teague, 2003).

The radius of gyration (Rg)
The radius of gyration (Rg) shows the impact of ligand on the structure of the protein, by measuring the compactness of a protein (Lobanov, Bogatyreva, & Galzitskaya, 2008). Sudden fluctuations of the radius of gyration values indicate the instability of protein folding, whereas a constant radius of gyration value shows that the ligands do not affect the protein folding (Khan et al., 2021;Mosquera-Yuqui, Lopez-Guerra, & Moncayo-Palacio, 2022).   The Rg graphs show small changes between 22.75 and 23.5 Å during the 110 ns simulation (Figure 4d) for both complexes. There were no abrupt fluctuations in the Rg plot. These observations indicate that both ZINC95486139 and ZINC33833639 preserve the folding behaviour of sterol 14a-demethylase and this enzyme remained in its native structure during the simulation.

Hydrogen bonds
The intermolecular hydrogen bonds have been checked to confirm the binding affinity of molecules towards sterol 14a-demethylase. The more hydrogen bonds in the complex, the more affinity and strong binding (Men endez, Accordino, Gerbino, & Appignanesi, 2016).
Both ZINC95486139 and ZINC33833639 form a stable hydrogen bond with the sterol 14a-demethylase, and the number of hydrogen bonds in the complexes changes during 110 ns simulation. In some points (e.g. 24-33 ns and 73-75 ns for ZINC33833639, and 18-19ns for ZINC95486139), the ligands form more than one hydrogen bond (Figure 4e). The stability of ZINC95486139 in the binding pocket is resulting from the 18 hydrophobic interactions shown in the molecular docking analysis. However, the stability of the studied molecules over time can be explained by the hydrogen bonds formed in the complexes alongside the high number of hydrophobic interactions formed between ZINC95486139/ZINC33833639 and sterol 14ademethylase. In general, concerning side effects, non-covalent inhibitors show multiple advantages compared with irreversible inhibitors (Pillaiyar, Manickam, Namasivayam, Hayashi, & Jung, 2016) (Figure 5).

Prediction of biological activity (PASS)
After obtaining the results from the server, we noted that ZINC95486139 is predicted to have 462 biological activities (PBA), among them five PBA having a value of Pa up to 0.9. The sterol 14a-demethylase inhibition is not predicted for both compounds, which allows us to suggest adding this activity to the online PASS server. For the second compound ZINC33833639, we have found a total of 435 PBA with five PBA having a value of Pa more than 0.9.  The five PBA with values of Pa more than 0.9 of the two compounds are presented in Table 2. The other PBA are in the supplementary data.

Pharmacophore study
Schr€ odinger's PHASE module was used to produce the E-pharmacophore hypothesis of the chosen target in order to develop an energetically optimized pharmacophore model. The pharmacophore hypothesis (ADHH) was generated for the 5FSA-ZINC33833639 complex and (HHH) for the 5FSA-ZINC95486139 complex. One hydrogen bond donor, one hydrogen bond acceptor, and two hydrophobic groups make up the former pharmacophore theory, whereas the latter has three hydrophobic groups (Figures 6, 7, 8, & 9).

MM-GBSA calculation
We have calculated the binding free energies using the MM-GBSA technique. MM-GBSA analysis was used to predict binding energies and find the most promising inhibitors for compounds with high glide scores. The compounds ZINC95486139 and ZINC33833639 showed binding energies of À75.07 kcal/mol and À83.89 kcal/mol, respectively. These energies indicate the high binding affinity of compounds for the receptor.

Conclusion
Mucormycosis is a fatal infectious disease that could affect patients with immunodeficiency caused by diabetes, cancer, COVID-19, or after organ transplantation procedures. The presented work suggests that the natural products could be successful candidates to inhibit fungal growth by inhibiting sterol 14ademethylase (an enzyme involved in ergosterol biosynthesis), where the major studied compounds (600 out of 880 compounds), were presented a free energy between À13.4 kcal/mol and -9Kcal/mol using AutoDock Vina. The presented molecular docking study shows that ZINC95486139 and ZINC33833639 can bind to sterol 14a-demethylase tightly by hydrogen bonds and different types of hydrophobic interactions. The binding energy obtained by MMGBSA and molecular dynamics simulations of the complexes proves the high affinity of ZINC95486139 and ZINC33833639 towards sterol 14a-demethylase, the high stability of the complexes through the time. The pharmacophore features were also predicted using protein-ligand complexes of both compounds. These findings open the door to testing the natural compounds from the AfroDb database as potent inhibitors against mucormycosis and to developing new safe drugs.

Disclosure statement
No potential conflict of interest was reported by the authors.