Virtual screening-driven drug discovery of SARS-CoV2 enzyme inhibitors targeting viral attachment, replication, post-translational modification and host immunity evasion infection mechanisms.

The novel coronavirus SARS-CoV2, the causative agent of the pandemic disease COVID-19, emerged in December 2019 forcing lockdown of communities in many countries. The absence of specific drugs and vaccines, the rapid transmission of the virus, and the increasing number of deaths worldwide necessitated the discovery of new substances for anti-COVID-19 drug development. With the aid of bioinformatics and computational modelling, ninety seven antiviral secondary metabolites from fungi were docked onto five SARS-CoV2 enzymes involved in viral attachment, replication, post-translational modification, and host immunity evasion infection mechanisms followed by molecular dynamics simulation and in silico ADMET prediction (absorption, distribution, metabolism, excretion and toxicity) of the hit compounds. Thus, three fumiquinazoline alkaloids scedapin C (15), quinadoline B (19) and norquinadoline A (20), the polyketide isochaetochromin D1 (8), and the terpenoid 11a-dehydroxyisoterreulactone A (11) exhibited high binding affinities on the target proteins, papain-like protease (PLpro), chymotrypsin-like protease (3CLpro), RNA-directed RNA polymerase (RdRp), non-structural protein 15 (nsp15), and the spike binding domain to GRP78. Molecular dynamics simulation was performed to optimize the interaction and investigate the stability of the top-scoring ligands in complex with the five target proteins. All tested complexes were found to have dynamic stability. Of the five top-scoring metabolites, quinadoline B (19) was predicted to confer favorable ADMET values, high gastrointestinal absorptive probability and poor blood-brain barrier crossing capacities.Communicated by Ramaswamy H. Sarma.

Introduction

The family Coronaviridae comprised of a large group of viruses is associated with a number of human respiratory infections (Fehr & Perlman, 2015). In 2019, the novel coronavirus SARS-CoV2 (previously known as 2019-nCov) emerged and is currently causing an alarming global pandemic. The virus is highly infectious resulting in over two million confirmed cases worldwide catapulting human death tolls to more than one hundred twenty thousand mortality cases (Dong et al., 2020). The present situation is aggravated by the unavailability of viral therapeutics and biointerventions such as specific drugs and vaccines (Liu et al., 2020). Given the rippled effect of this pandemic, an urgent call to discover effective therapeutic agents to combat the SARS-CoV2 virus is needed.

Modern drug discovery propelled by advances in bioinformatics and computational modelling has enabled virtual screening of biologically active compounds for hit identification and lead optimization (Aanouz et al., 2020; Abdelli et al., 2020; Abraham Peele et al., 2020; Adeoye et al., 2020; Al-Khafaji et al., 2020; Basit et al., 2020; Bhardwaj et al., 2020; Boopathi et al., 2020; Das et al., 2020; Elasnaoui & Chawki, 2020; Elfiky, 2020a, 2020b, 2020c; Elfiky & Azzam, 2020; Elmezayen et al., 2020; Enayatkhani et al., 2020; Enmozhi et al., 2020; Gupta et al., 2020; Gyebi et al., 2020; Hasan et al., 2020; Hendaus, 2020; Islam et al., 2020; Joshi et al., 2020; Khan et al., 2020a, 2020b, 2020c; Kitchen et al., 2004; Kumar et al., 2020a, 2020b; Lobo-Galo et al., 2020; Mahanta et al., 2020; Mittal et al., 2020; Muralidharan et al., 2020; Nejadi Babadaei et al., 2020a, 2020b;; Paniri et al., 2020; Pant et al., 2020; Sarma et al., 2020; Sinha et al., 2020; Sk et al., 2020; Umesh et al., 2020; Wahedi et al., 2020). Thus, interests to uncover SARS-CoV2 pathogenesis that may facilitate discovery of potential enzyme inhibitors is rapidly growing. As a response to the COVID-19 pandemic, several studies utilized molecular docking technologies for virtual screening and repurposing pre-existing drugs and natural products (Arya et al., 2020; Chen et al., 2020; Elfiky, 2020a). However, efforts to discover multi-targeting, selective, druggable and less toxic compounds are also equally important and thus, warranting further investigations.

As a contribution to this global scientific undertaking, our study focused on the exploitation of fungal secondary metabolites with profound antiviral activity on a range of known pathogenic viruses such as the human immunodeficiency virus (HIV), influenza virus, herpes simplex virus (HSV), and hepatitis C virus (HCV) – as potential drug prototypes against the SARS-CoV2 virus. In this paper, we disclose the in silico aided discovery of two fumiquinazoline alkaloids scedapin C (15) and quinadoline B (19), and polyketide isochaetochromin D1 (8) exhibiting strong in silico binding activity against five proteins involved in SARS-CoV2 infection such as viral attachment (spike binding domain to GRP78), replication (RdRp), post-translational modification (PLpro and 3CLpro), and evasion of host immunity (non-structural protein 15, nsp15) mechanisms. In addition, a computational ADMET was also determined to define the potential of these multi-targeting fungal natural products for lead optimization and drug development.

Materials and methods

Structural properties of selected target enzymes

Papain-like protease (PLpro)

SARS-CoV2 PLpro is a multi-domain protein and is the largest replicase unit. It is responsible for the cleavages located at the N-terminus of the replicase polyprotein. Structurally, it shares 83% sequence identity with PLpro from SARS-CoV (Liu et al., 2020; Wu et al., 2020).

Chymotrypsin-like protease (3CLpro)

The 3CLpro enzyme functions by cleaving the polyprotein at eleven distinct sites to generate other non-structural proteins vital in viral replication. Thus, drugs that target this protease play a critical role in inhibiting the replication of the virus. The 3CLpro monomer has three domains, wherein the active site, housing the catalytic dyad (Cis145 and His41), is located between domains I and II. The third domain is connected to the rest of the protein structure by a long loop and this gap contains a larger pocket relative to the active site (Arya et al., 2020; Wu et al., 2020).

Rna-directed RNA polymerase (RdRp)

The RdRp of SARS HCoV is the closest strain to SARS-CoV2. Thus, the SARS-CoV2 RdRp model, consisting of 801 amino acids, was designed using the SARS HCoV RdRp as a homolog. The active site of RdRp is highly conserved and contains aspartate residues that protrude from the beta-turn joining the β15 and β16. With a vital role in viral genome replication, the active site residues of RdRp and most of the 5-Å-region surrounding it are surface accessible and can bind to free nucleotides, including ATP, GTP, CTP, and UTP (Elfiky, 2020a). Similar to the viral proteases, RdRp is regarded as a molecular target for the development of new drugs against SARS-CoV2.

Non-structural protein 15 (nsp15)

Nsp15 is an endoribonuclease that cleaves 3′ of uridylates through a ribonuclease A (RNase A)-like mechanism (Ortiz-Alcantara, 2010). It is not required for viral replication to occur but it inhibits detection of the host immune system. It does this by preventing simultaneous activation of host cell dsRNA sensors (Kindler et al., 2017). The structure of SARS-CoV2 nsp15 is quite similar to other coronavirus nsp15 monomers, which consists of three distinct regions: the N-terminal domain, a subsequent middle domain, and a catalytic C-terminal nidoviral RNA uridylate-specific endoribonuclease domain. Generally, for coronaviruses, the C-terminal domain houses the active site of the enzyme, specifically the proposed catalytic triad His235, His250, and Lys290 (Kim et al., 2020).

S protein (spikes)

The spike protein is a homo-trimer protruding from the outer surface of the virion and serves as the major driving force for host cell recognition and entry. It is among the four structural proteins encoded by the viral RNA in addition to nsp 1-16. Although the primary entry of the virus is mediated by ACE2 receptors forming complexes with spike proteins, previous study demonstrated via molecular docking that there might be four regions in the spike that may act as binding sites for the glucose-regulated protein 78 (GRP78) (Ibrahim et al., 2020). Being a master chaperone protein, GRP78 is found in the lumen of the endoplasmic reticulum (ER) where it regulates important cellular processes, such as protein folding, cell death, and cell differentiation (Quinones et al., 2008). During cellular stress, however, GRP78 is overexpressed which increases the chance for it to be translocated from the ER to the cell membrane. Once in the cell membrane, GRP78 is susceptible to virus recognition by its substrate-binding domain (SBD) facilitating viral invasion (Ibrahim et al., 2020).

Target enzyme preparation

Five target proteins crucial for the infectivity of SARS-CoV2 were selected as molecular targets (Figure 1). The homology model of the catalytic domains of the target proteins was generated in the SWISS-MODEL workspace (swissmodel.expasy.org/workspace) using as templates the high-resolution crystal structure of 3CLpro (PDB ID: 6LU7), PLpro (PDB ID: 6W9C), RdRp (PDB ID: 6M71), nsp15 (PDB ID: 6VWW), and the spike’s protein binding domain to GRP78 (PDB ID: 6VXX). The three-dimensional structures of the enzymes were added in PDB formats to the UCSF Chimera platform (Pettersen et al., 2004).

Figure 1.

Structures of SARS-CoV2 target proteins: (A) papain-like protease, PLpro, (B) chymotrypsin-like protease, 3CLpro, (C) RNA-directed RNA polymerase, RdRp, (D) non-structural protein 15, nsp15, and (E) spike binding domain to GRP78.

Ligand selection and preparation

A library of ninety seven secondary metabolites 1–97 from marine and terrestrial fungi (Figure 2; Table S1) previously reported to inhibit human immunodeficiency virus (HIV), herpes simplex virus (HSV), hepatitis B virus (HBV), hepatitis C virus (HCV), influenza virus, tobacco mosaic virus, and enteroviruses were screened against the aforementioned target proteins (Bashyal et al., 2014; Gao et al., 2013; He et al., 2013; Huang et al., 2017a, 2017b; Liu et al., 2010; Nong et al., 2014; Pang et al., 2017; Peng et al., 2013; Sandargo et al., 2019; Singh et al., 1998, 2003; Tan et al., 2019; Zhao et al., 2017; Zhu et al., 2014). The ligands (fungal secondary metabolites) were added to the UCSF Chimera platform, rendered from SMILES notation or SYBYL mol2 file. Avogadro (version 1.1.1) an open-source molecular builder, was used to prepare the ligands (Hanwell et al., 2012).

Figure 2.

Structures of fungal secondary metabolites 1–25 with in silico antagonism to SARS-CoV2.

Molecular docking

Molecular docking experiments were performed on the UCSF Chimera platform (Pettersen et al., 2004). The three-dimensional structures of the protein were added to the docking platform as PDB formats. Each protein crystal structure was processed by removing existing co-crystallized ligands and water molecules. Meanwhile, the ligands were added to the docking platform, rendered from SMILES notation or added as SYBYL mol2 file. Minimization and dock-prepping of ligand and protein structures were done by adding the missing hydrogen atoms and appropriate charges to the structures employing the Gasteiger charge method computed using Amber’s Antechamber module (Wang et al., 2006). The docking procedure was done using flexible ligand into flexible active site protocol, where the ligand was allowed to be flexible and torsion within a grid box encompassing the ligand-binding cavity of each enzyme, predicted using COACH algorithms (Yang et al., 2013). With all docking parameters maintained at default values, molecular docking simulation was performed following the Broyden-Fletcher-Goldfarb-Shanoo (BFGS) algorithm of AutoDock Vina (version 1.1.2) (Trott & Olson, 2010). After each run, Autodock Vina provides a set of docking poses for each ligand with calculated binding affinities in which the docking pose with best affinity was chosen to represent the set and subjected to post-dock analysis. Visualization and analysis of the enzyme–ligand complex conformation were carried out using UCSF Chimera and Biovia Discovery Studios (version 4.1).

Molecular dynamics simulation

To understand the dynamics behavior of all the top scoring complexes, molecular dynamics simulation was employed using Amber 18 (Case et al., 2005). Drug topologies were generated using antechamber (Wang et al., 2001). With TIP3P water box and Na+ ions, all the systems were solvated and neutralized. Gentle energy minimization in two steps was performed to relax all the systems followed by heating. PME (Particle Mesh Ewald) algorithm for short and long range interactions was also used. A total of 50 ns simulation for each system was performed (Toukmaji et al., 2000). CPPTRAJ and PTRAJ (Roe & Cheatham III, 2013) were used to evaluate the trajectories for stability (RMSD) and residual flexibility (RMSF).

Binding free energy calculation

Trajectories from molecular dynamics simulations were subjected to binding free energy calculations using the MMPBSA.py script (Sun et al., 2014). This method is widely implemented by different studies to calculate the binding energy of a ligand-protein, protein-protein and protein-nucleic acid complexes (Khan et al., 2018, 2019; Wang et al., 2019). The equation below was used to calculate the total binding energy (ΔG.

To further understand each energy term such as electrostatics energy, Van der Waal forces, polar and non-polar interactions which may contribute to the total energy (G), the equation below was used.

Drug-likeness, ADME, and toxicity prediction

The computational prediction of the absorption, distribution, metabolism, and excretion (ADME) properties of hit compounds were done using SwissADME software. Pharmacokinetic profiles were evaluated according to Lipinski’s ‘rule-of-five’ which analyzes the biochemical features of a drug that may influence its absorption and permeation across cell membranes. Lipinski’s criteria states that for a compound to exhibit drug likeness, at least three of the following criteria must be fulfilled: molecular weight < 500 Daltons (Da), calculated lipophilicity (Log P) < 5, number of hydrogen-bond acceptors < 10, and number of hydrogen-bond donors < 5 (Macabeo et al., 2020). Furthermore, OSIRIS property explorer program was employed for in silico toxicity prediction which takes into account the potential mutagenicity, tumorigenicity, irritant effects, and reproductive toxicity of the hit compounds (Phukhamsakda et al., 2019). Using the same software, the solubility (Log S) of hits which is a determining factor for absorption and bioavailability was also predicted. As previously elucidated, values of Log S equivalent or greater than -4 isindicative of favorable solubility (Escobedo-González et al., 2017).

Results

Owing to the potential of cysteine proteases (3CLpro and PLpro), RdRp, nsp15, and the spike proteins as drug targets for SARS-CoV2 (Figure 3), ninety-seven fungal natural products 1–97 previously reported to exhibit antiviral properties were virtually screened for in silico inhibitory activity (Figure 2; Table S1). Using the UCSF Chimera platform, a molecular docking approach was performed to simulate and visualize protein-ligand interactions of the different fungal metabolites against these five SARS-CoV2 proteins. Of the ninety-seven metabolites screened, twenty-five compounds 1–25 (Bashyal et al., 2014; Gao et al., 2013; He et al., 2013; Huang et al., 2017a, 2017b; Liu et al., 2010; Nong et al., 2014; Pang et al., 2017; Peng et al., 2013; Sandargo et al., 2019; Singh et al., 1998, 2003; Tan et al., 2019; Zhao et al., 2017; Zhu et al., 2014) (Figure 2) showed high binding affinities to at least one of the five target proteins (Table 1).

Figure 3.

Pathways inhibited by fungal natural products cytochalasin Z8 (4), isochaetochromin D1 (8), 11a-dehydroxyisoterreulactone A (11), scedapin C (15), scequinadoline A (16), quinadoline B (19), norquinadoline A (20) and rubrolide S (23) predicted in silico during SARS-CoV2 host-cell infection.

Table 1.

Binding energy values of fungal metabolites 1–25 against five SARS-CoV2 proteins.

Compound	Viral Target	Binding Energy (kcal/mol)
		PLpro	3CLpro	RdRp	nsp15	S Protein
Altertoxin V (1)	HIV	−9.4	−7.2	−8.5	−8.1	−7.1
Altertoxin II (2)	HIV	−9.2	−7.5	−8.3	−7.6	−7.6
Penicillixanthone A (3)	HIV	−9.5	−8.2	−8.1	−6.9	−8.1
Cytochalasin Z8 (4)	HIV	−10.3	−7.9	−9.6	−8.5	−8.1
Stachybotrosin D (5)	HIV	−9.5	−6.5	−7.8	−7.8	−7.1
Chloropupukeanolide A (6)	HIV	−10.0	−7.5	−7.7	−6.8	−8.1
Phomasetin (7)	HIV-1 Integrase	−9.6	−7.2	−6.9	−6.9	−8.3
Isochaetochromin D1 (8)	HIV-1 Integrase	−9.9	−7.9	−8.4	−9.1	−9.5
Aspergilol H (R) (9)	Herpes simplex	−9.6	−7.9	−7.6	−8.2	−9.1
Aspergilol H (S) (10)	Herpes simplex	−9.2	−8.3	−7.7	−8.1	−9.2
11a-Dehydroxyiso- terreulactone A (11)	Herpes simplex	−8.9	−8.9	−9.1	−8.4	−10.2
Arisugacin A (12)	Herpes simplex	−10.0	−8.0	−8.5	−8.5	−9.2
Aspernolide A (13)	Herpes simplex	−9.2	−7.1	−7.5	−7.5	−9.3
Rhodatin (14)	Hepatitis C	−9.0	−7.1	−7.9	−7.9	−9.9
Scedapin C (15)	Hepatitis	−10.9	−8.6	−9.2	−8.9	−9.4
Scequinadoline A (16)	Hepatitis C	−9.1	−8.7	−9.1	−7.9	−8.5
14S-Oxoglyantrypine (17)	Influenza	−9.6	−8.4	−8.3	−7.4	−7.6
Deoxynortryptoquivaline (18)	Influenza	−9.6	−7.6	−7.4	−7.1	−8.0
Quinadoline B (19)	Influenza	−10.6	−8.3	−9.8	−9.1	−10.5
Norquinadoline A (20)	Influenza	−10.9	−8.1	−8.5	−8.6	−8.3
Asperterrestide A (S) (21)	Influenza	−8.9	−8.3	−7.6	−7.3	−8.7
Asperterrestide A (R) (22)	Influenza	−8.9	−8.4	−7.6	−7.4	−9.7
Rubrolide S (23)	Influenza	−8.7	−7.7	−7.4	−7.8	−10.3
Isoaspulvinone (24)	Influenza	−7.7	−7.4	−7.1	−8.6	−9.7
Pulvic acid (25)	Influenza	−7.9	−6.6	−7.5	−8.5	−9.3

Papain-like protease (PLpro) and chymotrypsin-like protease (3CLpro)

Molecular docking (PLpro)

Through in-silico-based screening, thirteen fungal secondary metabolites displaying high affinity towards PLpro with binding energies (BE’s) ranging from -9.5 to -10.9 kcal/mol (Table 1 and S2) were noted. Among these compounds, the fumiquinazoline marine alkaloids scedapin C (15) and norquinadoline A (20) demonstrated highest affinity to PLpro, both with BE’s of -10.9 kcal/mol. Another quinazoline alkaloid quinadoline B (19) also showed high affinity with BE’s of -10.6 kcal/mol (Figure 4). Alkaloid 15 was isolated from the marine-derived fungus Scedosporium apiospermum F41 − 1 and displayed significant antiviral activity against hepatitis C (Huang et al., 2017a). During the post-dock analysis, this fumiquinazoline alkaloid was noted to be bound to the putative binding site of PLpro through H-bonding with two ketones of the quinazolinedione core against Arg712. Additionally, the indole moiety is hinged to the enzyme’s pocket via pi-cation intermolecular bonding with Lys711, pi-pi stacking with His342, and pi-alkyl interaction with Ala579. On the other hand, alkaloids 19 and 20 were previously isolated from the mangrove-derived fungus Cladosporium sp. PJX-41 and showed activity against influenza A (H1N1) (Peng et al., 2013). The top-scored pose of 20 is tightly nestled to the PLpro’s binding site, stabilized mostly by van der Waals forces. Two H-bonds further strengthened its interaction with the enzyme, particularly between the carbonyl oxygens of the pyrazinoquinozolinedione core with Lys711 and Arg712. Similarly, 19 interacted with PLpro’s binding site via H-bonds with Lys711 and Arg355. The presence of a fused phenyl ring in the latter interacted with amino acid residues Leu557, Ala579, and Ile580 through pi-alkyl interactions. For this study, we test out the docking mechanism of lopinavir, an antiretroviral protease inhibitor used as a trial drug for COVID-19. Interestingly, all thirteen reported metabolites exhibited a much better binding mechanism against PLpro. Lopinavir, with a binding energy of -6.9 kcal/mol, likewise formed H-bonds with Lys711 and Arg558. Comparison of the docking behavior of the top three inhibitors for SARS-CoV2 PLpro revealed that these alkaloids – scedapin C (15), quinadoline B (19) and norquinadoline A (20) – all contain a quinazolinone core primarily involved in interacting with PLpro’s putative binding site. Between compounds 19 and 20, the latter gave a higher binding energy probably due to the presence of two ketone groups in the quinazolinone core responsible for two intermolecular hydrogen bonds and the lack of spiro junction which is present in alkaloids 15 and 20, allowing a more relaxed conformational flexibility for enhanced non-covalent interactions with SARS-CoV2 PLpro catalytic residues.

Figure 4.

Docked poses of scedapin C (15), quinadoline B (19), and norquinadoline A (20) against SARS-CoV2 PLpro shown as (A) ribbon representation and (B) molecular surface representation.

Molecular docking (3CLpro)

To assess the binding activity of the various fungal metabolites to 3CLpro, each was docked directly to the active site of the enzyme. Fourteen compounds were found to display relatively high affinity with BE’s ranging from -7.9 to -8.9 kcal/mol (Table 1 and S3). All fourteen top-ranked compounds exhibited a higher BE compared to the reference drug lopinavir (BE = −7.6 kcal/mol). Of which, the most notable inhibitor is the anti-HSV metabolite from the marine-derived fungus Aspergillus terreus (Nong et al., 2014), 11a-dehydroxyisoterreulactone A (11) with a BE of -8.9 kcal/mol (Figure 5). The results of the docking analysis indicated that 11 is bound to the catalytic residue His41 through pi-pi stacking with the pyranone ring fused in its polycyclic core. The enzyme-ligand complex is further stabilized by H-bonding and pi-alkyl interactions with its methoxyphenyl substituent. An unbiased (blind) docking was also performed for this compound in order to assess the preferential attachment of the ligand to the active site as opposed to other binding pockets. Interestingly, 11 (unbiased BE = −8.9 kcal/mol) adapts to the same docking mechanism within the binding pocket of the catalytic dyad. Alkaloids scedapin C (15) and scequinadoline A (16) also showed significant binding affinities towards 3CLpro with BE’s of -8.6 and -8.7 kcal/mol, respectively (Figure 5). Both compounds were isolated from the marine fungus Scedosporium apoispermum and displayed significant antiviral activity against hepatitis C (Huang et al., 2017a). In the case of compound 16, it is bound to the catalytic residue, His41, via two H-bonds with a carbonyl oxygen of the quinazolinone core and a hydroxyl group. Two other H-bonds were observed for Asn142 and Gly143 as well as pi-sulfur bond between Met49 and the indole moiety. The quinazoline core of 15, on the other hand, interacted with 3CLpro through pi-sulfur bonding with the catalytic amino acid Cys145 and pi-pi stacking with Met165. Its indole moiety is linked to His41 and Met49 by means of pi-pi stacking and pi-alkyl interaction, respectively. The presence of an indole moiety and a pyramidoquinazolinone core in both scedapin C (15) and scequinadoline A (16) is generally responsible for these two to exhibit notable in silico inhibitory activity against 3CLpro through strong interaction with the enzyme’s catalytic dyad. Although both alkaloids have comparable binding energies, 16 gave a slightly stronger attachment (higher BE) because of four H-bonds (mostly by the ketone groups) with the amino acid residues as opposed to 15 with mostly pi interactions.

Figure 5.

Docked poses of 11-dehydroxyisoterreulactone (11), scedapin C (15), and scequinadoline A (16) against SARS-CoV2 3CLpro shown as (A) ribbon representation and (B) molecular surface representation.

Rna-directed RNA polymerase (RdRp)

Molecular docking

Twelve secondary metabolites from fungi were noted to exhibit high affinity for RdRp with BE’s ranging from -8.1 to -9.8 kcal/mol (Table 1 and S4). Interestingly, the same group of secondary metabolites that conferred strong binding affinities to cysteine proteases was also potent against RdRp. Quinadoline B (19) possesed the highest binding affinity for RdRp with a BE of -9.8 kcal/mol (Figure 6). Compound 19’s attachment to RdRp’s binding site is due to the participation of two particular moieties. In its quinazolinone core, H-bonding can be observed between Gln73 as well as pi-cation and pi-alkyl interactions with Arg569 and Ala686, respectively. The indoline moiety of 19 is affixed to Tyr689, Ala580 and Ala688 residues through pi-pi stacking and pi-alkyl interactions. Other notable secondary metabolites are cytochalasin Z8 (4) and scedapin C (15) with RdRp BE’s of -9.6 and -9.2 kcal/mol, respectively (Figure 6). For compound 4, there were two conventional hydrogen bonds demonstrated by the hydroxy substituent of the macrocylic core (with Ser682) and the pyrrolidinone carbonyl oxygen (with Arg569). Several pi-pi stacking and pi-alkyl interactions were also observed against Tyr689, Ala580, Lys577, Leu576, Lys500, and Ile495. Compound 15 showed two conventional hydrogen bonds between a hydroxy substituent against Lys593 and Cys813. Also, in its quinazolinone core, the following interactions were observed: pi-alkyl with Ile589 and Leu758 and pi-sulfur with Cys813. Interestingly, these fungal metabolites showed higher BE’s to RdRp than the reference prodrug favipiravir, docked in its active form favipiravir-ribofuranosyl-5′-triphosphate (favipiravir-RTP) with a BE of -7.6 kcal/mol. The docked poses of compounds 4, 15 and 19 are shown in Figure 6. Since RdRp is an essential protein encoded in all RNA-containing viruses, such as SARS-CoV2, responsible for catalyzing the synthesis of the RNA complementary strand, we expect that compounds with antiviral activities against RNA-based viruses would be good candidates for potential RdRp inhibition. The results of this study, expectedly, revealed that the top three RdRp inhibitors: quinadoline B (19), cytochalasin Z8 (4), and Scedapin C (15) were previously known to exhibit activities against RNA-based viruses, specifically HIV, influenza, and hepatitis C, respectively.

Figure 6.

Docked poses of cytochalasin Z8 (4), scedapin C (15), and quinadoline B (19) against SARS-CoV2 RdRp shown as (A) ribbon representation and (B) molecular surface representation.

Non-structural protein 15 (nsp15)

Based on the results of our screening, we identified twelve fungal metabolites with high affinity to nsp15 with BE’s ranging from -8.1 to -9.1 kcal/mol (Table 1 and S5). Two secondary metabolites namely isochaetochromin D1 (8) and quinadoline B (19) exhibited remarkable binding mechanisms to nsp15, both having BE of -9.1 kcal/mol (Figure 7). Compound 8 from Fusarium sp. was previously reported to interfere with HIV-integrase activity (Wu et al., 2020). The best docking pose of 8 is attached to the binding domain via H-bonding with Val292 and catalytic residue His250. Several pi interactions were observed due to the presence of a polyaromatic core in 8. Furthermore, van der Waals (VDW) forces enhanced the stabilization of the ligand-enzyme complex, particularly against His235 and Lys290 along with other residues in the catalytic zone. Compound 19, on the other hand, showed notable interaction with the catalytic triad through pi-cation intermolecular bonding with Lys290 and VDW affinity with His235 and His250. Also, 19 exhibited H-bonding with Val292 and pi-pi/pi-alkyl interactions with Tyr343, Lys345, and Leu346. Another metabolite with interesting binding affinity towards nsp15 is the quinazoline alkaloid scedapin C (15), with BE = −8.9 kcal/mol (Figure 7). Similar to the 8 and 19, 15 hinged itself to the catalytic triad via pi-pi stacking of its indole moiety against His235 and VDW affinity with His250 and Lys290. Other interactions of 15 include H-bonding with Thr341 and pi-pi stacking with Tyr343, both directed at the quinazolinone core. These results were compared to the binding behavior of benzopurpurin B, a potent inhibitor against SARS-CoV nsp15 endoribonuclease activity in vitro (Ortiz-Alcantara, 2010). The said standard displayed a relatively good binding affinity with a BE of -9.4 kcal/mol. Although the BE of the standard is slightly higher compared to the three aforementioned fungal metabolites, all shared a similar binding pattern, hinging on nsp15’s catalytic triad of the C-terminal domain through H-bonding, pi interactions, and VDW stabilization.

Figure 7.

Docked poses of isochaetochromin D1 (8), scedapin C (15), and quinadoline B (19) against SARS-CoV2 nsp15 shown as (A) ribbon representation and (B) molecular surface representation.

S protein (spikes)

Using fungal natural products 1–97, the binding between Region IV (C480- C488) of the spike protein model and GRP78 was challenged. Via molecular-docking-based screening, thirteen secondary compounds from fungi were identified to exhibit favorable binding to the spike protein with BE’s ranging from -9.1 to -10.5 kcal/mol (Table 1 and S6). Of these, quinadoline B (19), rubrolide S (23), and 11a-dehydroxyisoterreulactone A (11) were the most promising with BE’s of -10.5, −10.3, and -10.2 kcal/mol, respectively (Figure 8). It is worth noting that quinadoline B (19) as a multi-targeting compound antagonizes not only the spike proteins, but also the cysteine proteases, RdRp and nsp15. Meanwhile, 11a-dehydroxyisoterreulactone A (11), a territrem derivative previously isolated from the marine-derived fungus Aspergillus terreus showed to inhibit HSV-1 in vitro (Nong et al., 2014). The major ligand-enzyme interactions of compound 19 at the supposed GRP78-binding region of the trimeric spike protein include the following: pi-sulfur bonding with Cys454, pi-alkyl with Ala444, pi-pi stacking with Phe430, and pi-anion interaction with Asp441. At the same putative binding site, compound 11 exhibited a conventional hydrogen bond between the pyranone moiety and Ala449. The phenol substructure of 11 displayed pi-pi/pi-anion interplay with Phe430, Asp441, and Asn343. Another potential inhibitor of SARS-CoV2 spikes is rubrolide S (23) isolated from the marine-derived fungus Aspergillus terreus OUCMDZ-1925 with known inhibitory property against Influenza A (H1N1) (Zhu et al., 2014). Compound 23’s attachment to Region IV of the spike protein is due to H-bonding between the carbonyl oxygen of its furanone moiety against Glu445 and Tyr447. Moreover, the phenol substituent is stabilized by pi-anion attraction with Asp441 while the rest of the structure is stabilized by VDW forces. Compounds 11, 19 and 23 showed a more promising binding affinity than the reference drug umifenovir, an inhibitor of viral spikes, with a BE of only -7.1 kcal/mol. Thus, the following results generated via molecular docking open another probable use of fungal secondary metabolites. To the best of our knowledge, this is the first molecular docking experiment that targeted the binding domain of the spike protein to GRP78.

Figure 8.

Docked poses of 11-dehydroxyisoterreulactone (11), quinadoline B (19), and rubrolide S (23) against SARS-CoV2 spike protein shown as (A) ribbon representation and (B) molecular surface representation.

Molecular dynamics simulation of top-ranked ligands in complex with viral target enzymes

Molecular dynamics (MD) simulation is the most widely implemented method to study protein-ligand binding at atomic level. In this study, we used MD simulation to understand the dynamics behaviour of our top scoring compounds. Post-simulation analysis for stability, as well as the root mean square deviation (RMSD), of all the systems were calculated. As shown, the systems in the bound state attained dynamic stability (Figure 9). In the case of 3CLPro and spike complexes, acceptable fluctuation was observed between 0 and 10 ns. While for the rest of the simulation time, the two systems were observed to be stable. The average RMSD for 3CLPro was observed to be 2 Å, while for the spike protein the average RMSD was 2.5 Å. In the case of nsp15, an acceptable fluctuation was observed. However, the overall RMSD remained uniform between 2 to 3 Å. Systems RdRp and PLpro (norquinadoline A [20] and scedapin C [15]) were found to have incremental increase in the RMSD value but overall, the graphs are uniform and therefore we speculate that the increment increase in the RMSD value is due to the large systems. For larger systems, relatively more time is required to straighten the graph. Hence, these results support the stability of the ligand bound systems.

Figure 9.

Root mean square deviation (RMSD, Å) of top-scoring protein-ligand complexes as determined during 50-ns molecular dynamics simulation.

In addition, we also calculated the residual flexibility of each system. Differential residual fluctuations were observed due to secondary structure differences and effect of the bound ligand (Figure 10). The average (root mean square fluctuation) RMSF for all systems except in the case of the spike protein was observed to be in between 0.5 to 1.5 Å. For the spike protein, fluctuation was observed to be high with an average of 2.6 Å.

Figure 10.

Root mean square fluctuation (RMSF) of selected top-scoring protein-ligand complexes. The x-axis shows the total residues while RMSF value in Å is given on y-axis.

To further understand the binding of the top scoring compounds, the binding free energy was calculated using MMPBSA approach. Along with the total binding energy (Gtotal), other energy terms such as vdW and electrostatic energies were also calculated. All the ligands possessed strong inhibitory effect against the specific target (Table 2). The total binding energy for 3CLPro is -11.68 kcal/mol, for the spike protein -10.47 kcal/mol, for nsp15 -7.03 kcal/mol, for RdRp -8.47 kcal/mol, and for PLpro -9.17 (norquinadoline A, 20) and -8.63 (scedapin C, 15) kcal/mol. These results strongly require the in vitro and in vivo experimental testing of these compounds against their respective target viral proteins.

Table 2.

Binding energy (kcal/mol) results of the top scoring compounds 8, 11, 15, 19 and 20.

Target	vdW	Electrostatic	Gtotal
11a-Dehydroxyisoterreulactone A (11) - 3CLpro	−54.39	−16.40	−11.68
Quinadoline B (19) – spike protein	−51.32	−36.48	−10.47
Isochaetochromin D1 (8) – nsp15	−49.99	−12.23	−7.03
Quinadoline B (19) – RdRp	−36.54	−2.36	−8.47
Norquinadoline A (20) – PLpro	−32.48	−13.32	−9.17
Scedapin C (15) – PLpro	−44.28	−11.46	−8.63

In silico ADMET and drug-likeness

The twenty-five high dock-scoring compounds were also submitted for in silico ADME (absorption, distribution, metabolism, excretion) screening using Swiss ADME (Daina et al., 2017) to predict their overall pharmacokinetic behavior. The druggability of each compound was assessed based on Lipinski’s Rule of Five, which is primarily focused on the prediction of the following descriptors: molecular weight, lipophilicity (reported as octanol-water partition coefficient), and the number of hydrogen bond donors and acceptors. The corresponding results are summarized in Table 3. Likewise, the BOILED-Egg (brain or intestinal estimated permeation predictive model), which is an intuitive graphical plot of the functions of lipophilicity and apparent polarity as described by the parameters WLOGP (atomistic octanol-water partition coefficient) and TPSA (topological polar surface area), respectively, was used to predict passive intestinal absorption and brain permeation of the compounds (Figure 11). Compounds located in the yellow region (yolk) have a high probability of BBB penetration while those in the white region have the propensity for passive absorption through the GI tract. Of the five metabolites, 11a-dehydroxyisoterreulactone A (11), quinadoline B (19), and norquinadoline A (20) were predicted to have good gastrointestinal (GI) absorption. They also exhibited low blood-brain barrier (BBB) penetration as reflected by their high topological surface area (TPSA). Despite having good binding affinity to all viral proteins, scedapin C (15) and isochaetochromin D1 (8) were predicted to be poorly absorbed in the GI tract hence explaining their location outside the white portion of the BOILED-egg model. Since the compounds are considered for oral delivery, they were also analyzed using the Lipinski’s rule of 5 – a rule that qualitatively assesses the drug-likeness of the candidate compounds (Daina et al., 2017). Of these, fungal compounds 11, 19, and 20 exhibited good bioavailability and drug-likeness by fulfilling Lipinski’s criteria. Apart from fulfilling the criteria, it is important to highlight that these fungal metabolites have good solubility which reflects good absorption, movement in the bloodstream and better elimination by the urinary tract (Escobedo-González et al., 2017). Quinadoline B (19), having an S log value of -3.87, has the best solubility among the multitargeting metabolites.

Table 3.

Lipinski’s Rule of Five for ADME analysis of compounds 1–25.

Compound	MW (<500)	#H acceptors (<10)	#H-bond donors (<5)	Lipophilicity MLogP (<5)	Lipinski #violations	Drug Likeness
1	318.32	4	1	1.58	0	yes
2	350.32	6	3	0.47	0	yes
3	638.57	14	6	−1.36	3	no
4	465.58	5	3	2.83	0	yes
5	428.56	5	2	2.89	0	yes
6	643.08	11	2	2.09	2	no
7	413.55	4	2	2.39	0	yes
8	544.51	10	6	0.48	2	no
9	598.6	10	6	1.26	2	no
10	598.6	10	6	1.26	2	no
11	468.54	7	1	2.78	0	yes
12	496.55	8	2	1.57	0	yes
13	424.44	7	2	2.09	0	yes
14	398.45	6	3	1.69	0	yes
15	564.61	9	1	0.97	2	no
16	487.55	6	3	1.63	0	yes
17	358.35	4	2	1.84	0	yes
18	516.55	8	1	2.61	1	yes
19	439.47	5	1	2.00	0	yes
20	471.51	6	3	1.36	0	yes
21	480.56	5	4	0.86	0	yes
22	480.56	5	4	0.86	0	yes
23	348.39	4	1	3.29	0	yes
24	296.27	5	3	1.32	0	yes
25	308.28	5	2	1.76	0	yes

Figure 11.

Prediction of gastrointestinal (GI) tract and brain permeation of the top five multitargeting hit fungal metabolites – isochaetochromin D1 (8), 11a-dehydroxyisoterreulactone A (11), scedapin C (15), quinadoline B (19), and norquinadoline A (20) – by brain or intestinal estimated permeation predictive model (BOILED-Egg) method.

A prediction of the toxicities, such as mutagenicity, tumorigenicity, irritant effect and reproductive toxicity, of the top-ranked compounds was also performed using OSIRIS Property Explorer. The toxicity risks and significant solubility of the compounds are presented in Table 4. Of the five multitargeting compounds, 11a-dehydroxyisoterreulactone A (11), scedapin C (15), and norquinadoline A (20) did not demonstrate any form of toxicity. Isochaetochromin D1 (8), however, was predicted to exhibit tumorigenicity and reproductive toxicity. Considering its good oral biovailability and druggability, quinadoline B (19), despite being a high-risk irritant, still remains a promising drug candidate for SARS-CoV2 because it did not confer mutagenic, tumorigenic, and reproductive toxicities.

Table 4.

Predicted toxicity parameters and solubility of compounds 1–25.

Compound	Toxicity				Solubility (logS)
	Mutagenicity	Tumorigenicity	Irritant Effect	Reproductive Toxicity
1	None	None	None	None	−2.56
2	Medium Risk	None	High Risk	None	−5.12
3	None	None	None	None	−5.35
4	None	None	None	None	−4.21
5	None	None	High Risk	None	−4.8
6	High Risk	High Risk	High Risk	High Risk	−5.17
7	None	None	None	None	−3.38
8	None	High Risk	None	Medium Risk	−7.92
9	High Risk	None	High Risk	None	−8.32
10	High Risk	None	High Risk	None	−8.32
11	None	None	None	None	−4.71
12	None	None	None	None	−4.19
13	None	None	None	High Risk	−3.69
14	None	None	None	None	−3.95
15	None	None	None	None	−4.37
16	None	None	None	None	−4.05
17	None	None	None	None	−3.31
18	None	None	None	None	−4.48
19	None	None	High Risk	None	−3.87
20	None	None	None	None	−4.0
21	None	None	None	None	−3.4
22	None	None	None	None	−3.4
23	None	None	None	High Risk	−4.51
24	None	None	None	High Risk	−3.02
25	None	Medium Risk	None	None	−2.82

Discussion

Secondary metabolites from fungi have been successfully developed into several drug medications that include, but not limited to, antineoplastics, antibiotics, and antivirals (Zaiyou et al., 2017). Due to the immense range of habitats that fungi inhabit and their need to compete against a diverse array of organisms, fungal organisms have developed numerous ways to survive (Hyde et al., 2019) and offers great potential as tomorrows bioresource for the discovery of biologically active natural products (Kuephadungphan et al., 2019; Macabeo et al., 2020; Mapook et al., 2020; Notarte et al., 2017, 2018, 2019; Phukhamsakda et al., 2018).

In drug discovery and development programs, understanding the pathogenesis of the disease is a matter of importance. For COVID-19, the pathogenesis of SARS-CoV2 begins when the viral spike proteins interact with host cell receptors (Rothan & Byrareddy, 2020). The spike protein is a homo-trimer protruding from the outer surface of the virion and serves as the major driving force for host cell recognition and entry. Although the primary entry of the virus is mediated by ACE2 receptors forming complexes with these proteins, previous study demonstrated that the chaperone protein glucose-regulated protein 78 or GRP78 localized in the lumen of the endoplasmic reticulum, can be translocated to the cell membrane during cellular stress (Ibrahim et al., 2019). Upon integration in the cell membrane, GRP78 mediates the secondary entry of SARS-CoV2 virus. Thus, inhibiting the interaction between the SARS-CoV2 spike protein and the host cell receptor GRP78 offers the possibility of decreasing the rate of viral infection (Ibrahim et al., 2020). In our study, we demonstrated that quinadoline B (19) showed the highest binding affinity against the GRP78-binding domain of spike protein with a BE of -10.5 kcal/mol. Compound 19 is an alkaloid with known activity against influenza A (H1N1) (Peng et al., 2013).

Knowledge of the subsequent steps in the viral life cycle (Figure 3) prompted us to target other viral proteins as well. Following attachment via spike protein, endocytosis and uncoating of the virus ensues inside the host cell (Rothan & Byrareddy, 2020). The uncoated viral RNA is consequently used as template to directly translate polyprotein 1a/1ab (pp1a/pp1ab) whose function is to encode nsp 1-16. This includes the cysteine proteases, PLpro and 3CLpro, which are associated with nsp3 and nsp5, respectively. These proteases cleave and activate other nsps such as nsp12 and 15. Nsp12, also known as RdRp, is important in viral replication while nsp15 evades host dsRNA sensors (Chen et al., 2020). Our study revealed that quinadoline B (19) and scepadin C (15) exhibited highest binding affinities to PLpro, 3CLpro, RdRp and nsp15 as high as -10.6 kcal/mol. Their in vitro activities were also tested in previous studies; 19 having activity against influenza A (H1N1) (Peng et al., 2013) while 15 against hepatitis C (Huang et al., 2017b). Influenza A and hepatitis C are RNA viruses but are not related to SARS-CoV2 family, Coronaviridae. This shows that 19 and 15 are active against a range of viruses with the potential to be developed as broad-spectrum antiviral agents – providing additional protection from emerging and re-emerging viral diseases such as SARS-CoV2 (Andersen et al., 2020).

Out of the twenty five top-ranked fungal metabolites, the best scoring ligands for each protein target, as determined by molecular docking studies, were subjected to molecular dynamics simulation. In particular, MD simulation was done on the following protein-ligand complexes: norquinadoline A (20)-PLpro, scedapin C (15)-PLpro, 11a-dehydroxyisoterreulactone A (11)-3CLpro, quinadoline B (19)-RdRp, isochaetochromin D1 (8)-nsp15, and quinadoline B (19)-spike protein. The RMSD of these top-ranked complexes obtained during a 50-ns MD simulation revealed that these ligand-bound systems attained dynamic stability. Furthermore, the computed total binding free energy of these complexes using MMPBSA approach provides credence on the strong binding capacity of the ligands to their respective target proteins.

Additionally, the pharmacokinetics of the biologically active fungal natural products were predicted through the elucidation of their ADMET (absorption, distribution, metabolism, excretion and toxicity) profiles in silico. Notably, the secondary metabolites 11a-dehydroxyisoterreulactone A (11), quinadoline B (19) and norquinadoline A (20) conferred high gastrointestinal (GI) absorption, poor blood-brain barrier penetrability and high drug-likeness as per Lipinski’s rule of five. As per toxicity prediction, 11 and 20 did not demonstrate any form of toxicity. Despite being a high-risk irritant, 19 still remains a promising drug candidate for SARS-CoV2 because it did not confer mutagenic, tumorigenic, and reproductive toxicities.

Conclusion

The rapid transmissibility and the lack of specific drugs and vaccines against the COVID-19 disease inspired the virtual screening of fungal secondary metabolites. Although RdRp and the cysteine proteases, 3CLpro and PLpro, have long been targeted in silico, our study focused on utilizing fungal secondary metabolites as potential inhibitors for these enzymes. Our study reports for the first time inhibitors that could target the spike-binding domain to GRP78 which mediates the secondary entry of the virus. Additionally, it is also one of the few studies that have identified inhibitors against nsp15 – a non-structural protein essential for evading host immunity. Of the fungal natural products that were virtually screened, two fumiquinazoline alkaloids quinadoline B (19), scedapin C (15), and the polyketide isochaetochromin D1 (8) were identified to have strong in silico inhibitory activity against the five target proteins in SARS-CoV2. Finally, the in silico ADMET prediction showed quinadoline B (19) to confer high gastrointestinal (GI) absorption, poor blood-brain barrier penetrability and high drug-likeness as per Lipinski’s rule of five. Molecular dynamics simulation and subsequent total free energy calculation revealed that the complexes formed between the top-ranked ligands and their respective protein targets are dynamically stable with high total free energy of binding. Therefore, these compounds can be utilized as templates or prototypes for further development of multitargeting ligands against SARS-CoV2. The present work may be used as a framework for future discovery and development of novel multitargeting compounds against the COVID-19 disease.