Introduction

Viral diseases like MERS-CoV, SARS-CoV, influenza, SARS-CoV-2, and HCV pose global threats1. Viral polymerases, crucial for genome replication and transcription, are classified into RNA-dependent RNA, RNA-dependent DNA, and DNA-dependent RNA polymerases2,3. HCV and SARS-CoV-2 have a positive-sense RNA-dependent RNA polymerase, RdRp. These enzymes have the most important roles in the replication cycle of viruses4,5. More importantly, HCV polymerase shares unique characteristic finger, palm, and thumb domains common to most polymerases, and its catalytic activity is associated with three key residues: Asp319, Asp318, and Asp220 (Fig.1A, B). It replicates HCV viral RNA from the positive-strand RNA4,6. The SARS-CoV-2 polymerase is composed of a core protein with the usual cofactors nsp7 and nsp8. It is an enzyme that plays an essential role in the replication and transcription of the viral genome. This polymerase contains two major domains, termed the nucleotidyl-transferase -NiRAN domain and the N-terminal domain. An interface to the C-terminal domain connects these two domains. The C-terminal domain is composed of three subdomains: finger, thumb, and palm. Furthermore, its catalytic functions are provided by four aspartate residues: Asp618 and Asp623 in motif A and Asp760 and Asp761 in motif C (Fig.1C, D).5,7

Fig. 1
figure 1

Representation of the structures of A) SARS-CoV-2 RdRp in complex with its non-structural protein cofactors (nsp7 and nsp8). The domain diagram of SARS-CoV-2 RdRp in complex with these two cofactors shows a large N-terminal, which contains two domains; NiRAN (1–250) and Interface (251–397), that connects to the polymerase (RdRp) domain (398–932). This polymerase also comprises three main subdomains: the fingers (397–581 and 627–687), palm (582–627 and 688–815), and thumb (816–932). The RdRp domain contains seven motifs: motif A (612–626), motif B (678–710), motif C (753–767), motif D (771–796), motif E (810–820), motif G (499–511), and motif F (544–560), B) Ribbon diagram of SARS-CoV-2 RdRp in complex with nsp7 and nsp8; NiRAN (colored as purple), Interface (colored as green), Fingers (colored in orange), Palm (colored in blue), and Thumb (colored in green). The active site contains four catalytic aspartate residues (Asp618, Asp623, Asp760, Asp761) and two metallic cofactors, Mg2+(A) and Mg2+(B), essential for enzymatic activities. C) the 3D illustration of HCV NS5B, containing the fingers (1–210, and 224–248) in red, palm (210–224, and 248–358) in yellow, thumb (358–530) in blue, and a linker (530–591) in gray. D) The active site contains two magnesium cofactors and three aspartate residues (Asp220, Asp318, Asp319). (The Figure was generated by Discovery Studio V 4.5, and images were prepared by Adobe Photoshop V 2022).

Full size image

Natural products and analogs are key in cancer, cardiovascular, and infectious disease therapies, offering fewer side effects, lower costs, and diverse activities like antioxidant anti-inflammation and antiviral effects8,9. Fungi are rich sources of bioactive compounds. Endophytic fungi, living symbiotically in plant tissues without causing disease, produce valuable natural products throughout part or all of their life cycle10. The order Sordariales, which belongs to the phylum Ascomycota, encompasses an extensive and diverse fungal group that produces secondary metabolites like terpenes, alkaloids, and polyketides. These might have promising biological properties, including being immunosuppressive11, cytotoxicity12, antifungal13, antioxidant14, antimicrobial15, and anti-inflammation16. Though the structures seemed very complicated, modern technologies enhance and enable their application in pharmaceutical and biotechnological approaches17.

Computational techniques, including molecular docking and molecular dynamics (MD) simulations, becomes handy in drug discovery and development, especially natural compounds with fewer side effects, turning them into vital inhibitors for targeting viral polymerases. These virus-encoded enzymes are important for genome replication and transcription; thus, they remain ideal therapeutic targets2,3,18. Previously, natural compounds were screened against the SARS-CoV-2 polymerase, e.g. in 2020, Ghibi et al. identified 36 compounds with stronger binding affinity than remdesivir, including alkaloids and flavonoids like cryptospirolepine and usararotenoid A19. In 2021, Ebrahimi and coworkers reported a few fungal metabolites, including methoxycytochalasin and pyricidin-A, as SARS-CoV-2 RdRp inhibitors20. More recently, sterostrein B and C have been identified by Zarei et al. (2023) as potent SARS-CoV-2 and HCV polymerase inhibitors, recommending further in vitro and in vivo studies to investigate their antiviral efficacy21. This study focuses on the antiviral potential of endophytic fungi from the ascomycete order Sordariales, including secondary metabolites17, against aspartyl polymerases of HCV and SARS-CoV-2 through an extensive virtual screening pipeline. The virtual screening workflow begins with blind docking, wherein the simulation box covers the total targets’ structures to find the ligands that bind to the active site. Then, in the targeted docking, Ribavirin was chosen as a multi-target standard drug22,23,24, along with those natural compounds with the stronger binding affinities were introduced to 300 ns MD simulations to shed more light on the effects of their interactions on the structures and dynamics of the viral targets. The outcome of this work might provide significant progress in identifying novel multi-target compounds against the aspartyl-based polymerases; HCV NS5B and SARS-CoV-2 RdRp.

Materials and methods

Molecular structure preparation

The PDB codes for the target enzymes, SARS-CoV-2 RdRp (7ED5) and the HCV NS5B RdRp (2XI3), were downloaded from the Protein Data Bank. ChemDraw software version 20.1.1 was used to draw the two-dimensional structures of fungal secondary metabolites. These structures were then transferred in the Avogadro software V1.2.0 to conduct energy minimization and optimization of the conformation via the steepest descent approach. Then, the saved 3D structures in Gaussian format were introduced to Gaussian for specifying geometry optimization at the B3LYP/6-31G(d) level of theory. The final geometry of the ligands were converted and saved in PDB format using GaussView for compatibility with molecular visualization and docking studies.

Docking studies

Molecular docking, as a sub-strategy in computer-assisted drug design, has opened up routes for studying the binding of ligands to targets25,26,27,28,29,30,31. These methods offer some advantages, the major ones being the simplification of the complexity of the ligand-target binding problem by reducing degrees of freedom and proper estimation of binding energies. Subsequently, the docking protocol was pursued using an AutoDock Vina and divided into two steps: first, a blind docking over the whole protein structure. Ligands showing binding energies less than −7 kcal/mol and binding into the active site were taken to the next step of targeted docking. The optimized chemical structures of all ligand in PDB format were introduced into Autodock Tools 1.5.6 for the incorporation of polar hydrogens, Gasteiger charges, and determination of torsion of freedom. Finally, they have been saved in PDBQT format. As for the receptors, their PDB formats were subjected to Autodock Tools 1.5.6 for merging non-polar hydrogens, addition of Kollman charges and polar hydrogens, finally saving the files in PDBQT format. In the blind docking process, the grid boxes were set large enough to cover all of the protein atoms. However, in the targeted molecular docking, the grid box dimensions were defined with the x, y, and z coordination centers of 128, 134, 127 for SARS-CoV-2 RdRp; and 8, 6.567, 13 for HCV NS5B, respectively. Moreover, the grid box sizes were determined at 40 × 40 × 40. Actually, the grid box sizes were confined to the active sites of the targets in the targeted docking. Finally, molecular docking studies were performed using Autodock Vina. Regarding the docking run process, the viral proteins were set to be fixed in the box, whereas the studied compounds were allowed to be flexible. Finally, the exhaustiveness was set to 24 to regulate the number of calculation iterations, with other parameters maintained at default settings. Ligands Emestrin K (Lig-1), Omega-hydroxyemodin (Lig-2), and Zopfiellamide A (Lig-3) with the lowest binding energies are selected for MD simulations.

MD simulations

The present study has utilized MD simulations to investigate the atomic-scale dynamics in ligand-receptor binding32. Three protein–ligand systems were selected for further MD simulations to investigate the interactions of these ligands on the structures of the studied polymerases using the Gromacs software package (version 2021). Topology data for the receptors were derived directly from the OPLS-AA force field33, while ligands with the best binding energies were taken from the LigParGen server34. For saturation, SPC/E water was used, and counterions such that the overall charge of the systems was neutralized were added to the simulation boxes35. Energy minimization was performed by the steepest descent method until a force threshold less than 10 kJ/mol−1·nm−1 was achieved36. Periodic boundary conditions were applied in all directions (X, Y, Z). Prior to running MD simulations, we performed a two-step equilibration process. First, under the NVT ensemble, systems were gradually heated to 310 K using a V-rescale thermostat over 100 ps, with temperature adjustments at 0.1 ps intervals. This was followed by 1 ns under the NPT ensemble, maintaining pressure at 1 bar through the Parrinello–Rahman algorithm with a pressure coupling constant of 0.2 ps. Van der Waals and electrostatic interactions were evaluated within a cut-off of 1 nm. Bond constraints applied on heavy atoms forced the application of the LINCS algorithm37. The MD simulations of 300 ns each were performed by the leapfrog algorithm38. Two-dimensional molecular interactions, along with graphical representations were generated using LigPlot+39 and Discovery Studio V4.5.0. Finally, the toxicity prediction of fungal metabolites was performed using the ProTox-II web server39, while pharmacokinetic properties were analyzed using the SwissADME server. Moreover, in this study, the binding free energy calculations for receptor-ligand complexes were performed using the MM/PBSA (Molecular Mechanics/Poisson-Boltzmann Surface Area) technique40.

Results and discussion

Molecular docking analysis provides valuable information about ligand binding and their respective spatial arrangements on targets such as proteins. From the total of 174 docked secondary metabolites, 76 low-toxic compounds with LD50 over 5,000 mg/kg in humans or animals showed potential candidates for the blind docking approach. This resulted in the selection of 10 compounds based on the binding energies below −7 kcal/mol across all viral polymerase targets for the second targeted docking stage. The results are summarized in Table 1. Upon the results of the targeted docking stage, three compounds named Emestrin K (Lig-1), Omega-hydroxyemodin (Lig-2), and Zopfiellamide A (Lig-3) showed the lowest binding energies of −9.11, −8.9, and −9.75 kcal/mol against HCV NS5B, and −9.05, −8.53, and −9.83 kcal/mol against SARS-CoV-2 RdRp. These three ligands with superior binding energies and good clustering ranking across the two viral receptors were then chosen for MD simulations. Table 2 summarizes the detailed biological properties of these compounds. The two-dimensional and three-dimensional interactions of three selected secondary metabolites with three viral protein targets are shown in Figs. 1 and 2.

Table 1 Targeted Docking Results for the viral receptors.
Full size table
Table 2 Secondary Metabolites with superior binding affinities to the two studied viral receptors in Molecular Docking Studies.
Full size table
Fig. 2
figure 2figure 2figure 2figure 2

3D representation of interactions of SARS-CoV-2 RdRp (A-D), HCV NS5B (E–H) and Ribavirin, Lig-1, Lig-2, and Lig-3.(The Figure was generated by Discovery Studio V 4.5, and images were prepared by Adobe Photoshop V 2022).

Full size image

Docking studies on SARS-CoV-2 RdRp

About the molecular docking results, Ribavirin was used as a control drug, which showed a binding affinity of −6.5 kcal/mol when it bound to the active site of SARS-CoV-2 RdRp. It has formed nine hydrogen bonds with residues from finger subdomain Ser501, Ala502 (motif G), Val560 (motif F), Thr565, Gly683, Asp684, and Thr687 (motif B) (Figs. 2A and 3). In comparison to the control drug, interestingly, all of the three ligands showed superior binding energies, among which Lig-3 is bound to the SARS-CoV-2 RdRp through the lowest binding energy of −9.8 kcal/mol, which was followed by Lig-1 and Lig-2 with binding affinities of −9.1 and −8.5 kcal/mol, respectively.

Fig. 3
figure 3

2D representations of interactions between viral polymerases (I: SARS-CoV-2 RdRp, II: HCV NS5B) and Lig-1, Lig-2, and Lig-3. (The Figure was generated by LigPlot + V 2.2.8, and images were prepared by Adobe Photoshop V 2022).

Full size image

Lig-1 formed two strong hydrogen bonds with residues Thr556 in the catalytic motif F and one strong H-bond with the catalytic residue Asp623 located in the catalytic motif A. This compound also interacts through hydrogen bonding with the Mg2+ cofactor (Figures. 2B, and 3). Strong hydrogen bonds prevent structural adjustments necessary for the correct enzymatic reaction; thus, the enzymatic function may be disturbed. Moreover, this secondary metabolite made hydrophobic contacts with key residues located at the catalytic motifs F (Arg553, Lys 545, Arg555), A (Arg624), and C (Ser759, and Asp761). These interactions are predicted to disturb the binding of incoming nucleotides, inhibit the metal-dependent catalytic mechanism, and decrease the structural stability of the enzyme. Apparently, such hydrophobic interactions might cause important disturbances in the catalytic performance of the enzyme45. Moreover, Lig-2 interacted with the viral target through two H-bonds with motif F (Thr556) and two H-bonds with the catalytic motif A (Cys622, Arg624) and B (Asn691). Hydrogen bonding helps stabilize the inhibitors in the active site of the enzyme and inhibit catalytic activity. Further, in the catalytic motifs A, B, and F, hydrophobic are also observed with Thr687, Ser682, Asp623, Arg553, Arg555, and Val557. Amino acids in these motifs regulate RNA synthesis, and their interaction with natural inhibitors may reduce viral replication. Thus, natural inhibitors interact only with threonine and other amino acids, which are essential in decreasing enzymatic activity and may disturb viral replication (Figures. 2C, and 3). Lig-3 was positioned within the F and G motifs of the finger subdomain, a region conserved among all RdRps and critical for the polymerization mechanism. It establishes two hydrogen bonds with Val560 (motif F) and one with Gly683 (motif B), stabilizing its positioning at the enzyme active site so that proper placing of RNA or rNTPs is impeded. Lig-3 also formed hydrophobic contacts with Asp684, Arg569, Thr565, Lys500, Ala558, Ser501, Glu573, Tyr689, Ala685, Leu576, and Lys577 of finger and palm subdomains. These interactions may change the binding and translocation of RNA within the replication channel. Consequently, such interruptions can influence viral replication by impairing proper nucleotide translocation and preventing RNA chain elongation (Figures. 2D, and 3)45,46. Armed with the finding of molecular docking, it can be concluded that the chosen ligands have shown strong polar interactions with catalytic motifs A, F, and G, which these conserved regions have large role in the template binding, NTP binding, and polymerization47,48. Thus, these hydrogen bonds may cause significant interferences in the entrance and suitable accommodation of NTP and template binding, ending up with important inefficiency in the enzyme’s catalytic function, resulting in an important disruption in the replication of the virus.

Molecular docking studies of HCV NS5B

The fungal secondary metabolites were evaluated for their inhibitory potential against HCV NS5B under the same computational protocol (Fig. 2, 3). Ribavirin bound to HCV NS5B with a binding energy of − 6.1 kcal/mol, forming five hydrogen bonds with Asn291, the catalytic Asp318, Gly449, and Ser556 (Figures. 2E, and 3). Additionally, the drug interacts with residues in the thumb domain (Phe193, Arg200, Thr287, and Ser288) and the palm domain (Cys366, Met414, and Tyr448).

It is worth mentioning that all three secondary metabolites demonstrated stronger binding energies than Ribavirin, with Lig-3 showing the strongest binding affinity (−9.6 kcal/mol), followed by Lig-1 (−9.0 kcal/mol) and Lig-2 (−8.9 kcal/mol). The three catalytic aspartate residues (Asp220, Asp318, and Asp319) located at the motif C of the palm subdomain are critical for the catalytic activity of this enzyme. It is essential to note that Lig-1 forms strong hydrogen bonds with the catalytic residues Asp318 and Cys316 (catalytic motif C), along with H-bonds with Arg386, and Ser556 located at the thumb subdomain. The interactions can stabilize the conformation of the compound in the active site of the HCV NS5B, and restricts the movement of these regions that are needed for the viral RNA synthesis. Furthermore, hydrophobic interactions with residues in the thumb subdomain (Gly410, Asn411, Met414, Tyr448, Glu446, Ile447, Arg394) and palm subdomain (catalytic Asp319) may further disturb the proper orientation of RNA towards the active site, leading to significant disruption in the polymerase function and obstructing viral replication (Figures. 2F, and 3). Surprisingly, Lig-2 formed seven strong hydrogen bonds with catalytic motifs A (Thr221, Asp225), motif C (Asp318), NTP tunnel (Tyr4, Arg158, and Leu159), and metallic cofactor Mg2+. More importantly, the H-bonding with these residues can stabilize the structure of the ligand at the catalytic pocket of the viral enzyme and also might disrupt the exact binding of the RNA template and incoming NTP to the catalytic pocket, disrupting catalytic functions of the polymerase (Figures. 2G, and 3). It is interesting to note that Lig-3 also made four hydrogen bonds; one with catalytic residues Asp225 (motif A), Asn291 (motif B), and two with Asp318 (motif C). Moreover, this secondary metabolite also formed hydrophobic contacts with Ser407, Gly410, Asn411, Met414, Glu446, Tyr448, and Ser556 in the thumb region and with Phe193 in the finger region (Figures. 2H, and 3). According to the molecular docking results, it can be asserted that the three selected ligands basically interacted with the catalytic motifs A and C. Since these two conserved regions have potential roles in the substrate binding and catalytic activities, thus, these polar interaction can strongly disrupt the exact binding of the substrate to the viral pocket, ending up with the inhibiting the viral replication.

MD simulations

MD simulations of SARS-CoV-2 RdRp

Molecular dynamics simulation is one important approach to investigating protein behaviors, and it provides considerable promise for accelerating drug discovery and development49. To investigate the dynamics of each ligand-receptor complex and their effects on viral building blocks, 300 ns MD simulations were added for each ligand-enzyme pair. The RMSD (root-mean-square deviations) of all protein atoms were monitored initially to verify the system had reached an equilibrated state (Fig. 4). Regarding the free SARS-CoV-2 RdRp, the RMSD exhibited a rapid increase beyond 0.22 nm during the initial 10 ns. Following this initial rise, only minor fluctuations were observed until the halfway point of the simulation (at 150 ns) (Fig. 4). After this, the RMSD plot displayed very slight upward fluctuations with the RMSD magnitude of 0.21 nm on average and plateaued at this value for the rest of simulation time. This suggests that the protein did not experience harsh oscillations or changes, and quickly attained its equilibrium, confirming that a 300 ns simulation provides adequate time for effective analysis of this system. The RMSD values for Ligands 1–3 and Ribavirin (the control drug) showed a rapid rise during the first 20 ns, reaching over 0.27 nm for Ribavirin, approximately 0.28 nm for Lig-2 and Lig-3, and around 0.30 nm for Lig-1. Notably, Lig-1 exhibited higher RMSD values compared to Ribavirin, while Lig-2 and Lig-3 displayed RMSD magnitudes similar or comparable to the standard drug. In total, the average RMSD magnitudes for free RdRp were 0.112 nm and 0.147 nm, 0.169 nm, 0.132 nm, and 0.139 nm, for the frameworks containing Ribavirin, Lig-1, Lig-2, and Lig-3 correspondingly. These results align well with the binding energies of the ligands, further highlighting their distinct effects on RdRp stability.

Fig. 4
figure 4

Changes in the RMSD, Rg, and SASA and time evaluation plots of H-bond of viral polymerases (I: SARS-CoV-2 RdRp, II: HCV NS5B) in free state and in complex with Ribavirin, Lig-1, Lig-2, and Lig-3. (The plots were generated by Excel 2013, and images were prepared by Adobe Photoshop V 2022).

Full size image

Radius of gyration (Rg) is a measure of the size of a protein, and its variation calculates the compressing/expanding of the protein. The compactness of an enzyme can potentially alter the substrate binding to the active site49. As shown in Fig. 4, the Rg magnitudes for free RdRp stood at the same level somehow during simulation with an average Rg value of 3.21 nm. Surprisingly, Rg magnitudes decreased for all of the ligand-bound complexes, e.g., magnitudes of Rg for the complexes comprising Ribavirin, Lig-1, Lig-2, and Lig-3 to be 3.20 nm, 3.198 nm, 3.207 nm, and 3.204 nm, respectively. The highest compression and compactness in the viral polymerase structure was observed after binding to Lig-1 and Lig-3 compared to that caused by the standard drug. Interestingly, this compression in the SARS-CoV-2 RdRp structure may destabilize the protein structure, leading to a remarkable negative impact on its enzymatic functions and preventing the substrate from binding. This hypothesis is supported by the solvent accessible surface area (SASA) analysis shown in Fig. 4. Overall, SASA magnitudes decreased for the complexes. However, it can be noticed that the magnitudes of Rg and SASA were observed to decrease more radically than those of Ribavirin, which could be attributed to the compactness of the protein after binding to the two ligands. Such compactness may reduce solvent exposure by disturbing internal water channels that could destabilize the overall structure, affecting enzymatic activity.

Another important contribution concerns H-bond interactions, which play a significant role in ligand stability at a protein’s active site50. Thus, the time evolution plots of the H-bonds involved in the molecular interactions between the ligands and the studied polymerases were further analyzed (Fig. 4). Over the simulation period, the number of hydrogen bonds oscillated to some point. On average, all the ligands maintained two hydrogen bonds, except for Lig-1, which established only one hydrogen bond within the catalytic pocket of the viral enzyme.

The root mean square fluctuation (RMSF) is analyzed to determine how each amino acid of protein fluctuates. These fluctuations of protein are changed when the proteins are bound with the ligands51. Ligands 1–3 and the control drug Ribavirin modified the RMSF of some regions in the SARS-CoV-2 RdRp enzyme. For the control drug, significant RMSF fluctuations were observed in residues 11–30, 34, and 227–230 within the N-terminal domain, residues 259–263 in the interface domain, residues 582–585 in the palm subdomain, and residues 821–826 and 829 in the thumb subdomain. However, the binding of Ribavirin led to reduced RMSF fluctuations in several regions of the protein. In the N-terminal domain, these decreases were noted at residues 41, 62, 107–109, 125–139, 163–182, and 203. In the interface domain, reductions occurred at residues 252, 260, 304–312, 337, 363–365, and 386–399. The finger subdomain displayed decreases at residues 404–413, 431, 444–461, 494–512 (motif G), 514–516, 543, and 554–565 (motif F). In the palm subdomain, fluctuations diminished at residues 607–611 (motif A), 643, 673–685 (motif B), 720–725, 767–771, 780–782, 795 (motif D), and 807–812 (motif D). Lastly, in the thumb subdomain, reduced fluctuations were detected at residues 845–848 and 903–905 (Fig. 9A). In comparison to the control drug, the binding of our ligands led to important changes in the RMSF magnitudes of residues responsible for substrate binding and catalytic functions.

For example, Lig-1 increased the RMSF values of residues 27, 41, 52, 62, 152, 169, 199, 228–230 (N-terminal domain), and 617–626 (catalytic motif A, palm subdomain). However, binding of Lig-1 led the RMSF values of the viral polymerase to experience significant fall near residues 62, 81–92, 135–159, 166–188, 213 (N-terminal domain), 262, 336–338, 363–365 (interface domain), 390–403, 408–424, 433–457, 493–496, 499–516 (motif G, fingers subdomain), 544–560 (NTP entry tunnel, motif F, fingers subdomain), 575–583 (fingers subdomain), 590, 595, 600–606, 645, 668–673, 679–682 (fingers subdomain), 751–766 (catalytic motif C, palm subdomain), 797, 804–810 (palm subdomain), 842, and 845–848, and 903 (thumb subdomain) (Fig. 5). As for Lig-2-SARS CoV-2 RdRp complex (Fig. 5), the RMSF magnitudes experienced some rise near residues 7, 16–24, 40–52, 76–92, 137, 155, 169–171, 191–203, 219–223, 230 located at the N-terminal domain, and around residues 257–265, 303–305, 321–323, 365, 384 (interface domain), 430–442, 470–485 (fingers subdomain), 593–597, 635–648, 643–648, 686–691 (fingers subdomain), 695–708 (motif B, palm subdomain), 726–735, 754–765 (catalytic motif C, palm subdomain), 790–794, 803–810 (palm subdomain), 854, and 876 (thumb subdomain). However, the RMSF magnitudes dropped around residues 60–62, 107, 178–181 (N-terminal domain), 337, 384–387, 403, 446–454, 499–514 (motif G, fingers subdomain), 547–560 (NTP entry tunnel, motif F, fingers subdomain), 561–565, 584, 611–624 (catalytic motif A, palm subdomain), 670 (fingers subdomain), 768–773(palm subdomain), and 841 (thumb subdomain). Upon binding Lig-3 to the viral polymerase, the RMSF plot saw remarkable decreases around residues 4–6, 14, 25, 107, 158, 174, 191, 260 (N-terminal domain), 337–339, 353, 360–366 (interface domain), 409–425, 451, 494, (fingers subdomain), 498–510 (motif G, fingers subdomain), 547–560 (NTP entry tunnel, motif F, fingers subdomain), 565–570, 577–580, 581–586 (fingers subdomain), 597–608, 610–621, 623, (catalytic motif A, palm subdomain), 643–644, 679–682, 759, 762, 807–810 (palm subdomain), 826, and 905 (thumb subdomain). However, the RMSF plot increased at some spots around residues 41, 51, 62–64, 92, 136, 228 (N-terminal domain), 261–272, 302, 389 (interface domain), 432 (fingers subdomain), 713–718, 731–740, 750, 759 (palm subdomain), 849–853, and 875 (thumb subdomain) (Fig. 4). In general, it should be noted that ligand binding caused significant changes in the fluctuations degrees of the residues located at the catalytic motifs A and F, which agree with the molecular docking results. These conserved regions are responsible for NTP binding and template binding, therefore these phenomena can cause important disruption on the catalytic function of SARS CoV-2 RdRp and hinder the replication of the virus.

Fig. 5
figure 5figure 5figure 5figure 5

Changes in the root-mean-square fluctuation (RMSF) of SARS-CoV-2 RdRp (A-D), and HCV NS5B (E–H) in their free state and in complex with Ribavirin, Lig-1, Lig-2, and Lig-3. (The plots were generated by Excel 2013, and images were prepared by Adobe Photoshop V 2022).

Full size image

In summary, it is plausible that in comparison to Ribavirin, the binding of the three natural ligands to the viral polymerase caused notable changes in the RMSF values of RdRp amino acids, particularly those located at the catalytic motifs A, G, C, F, and the NTP entry channel. These domains have important functions in binding RNA into the catalytic pocket and play a supporting role in the polymerase’s catalytic action. In contrast, this interaction may occlude the NTP entry channel and interfere with the RNA binding at the active site. This phenomenon thus blocks the enzyme’s catalytic activities and inhibits the formation of new infectious viral particles. To examine the stability of the selected hits (Ligands 1–3) within the viral polymerase catalytic pocket, the time evolution plots of the average distance between the ligands and the receptors’ active site were tracked (Fig. S1, Supporting Information). All of the studied ligands, except for Lig-3, were positioned at the mean distances between 0.31 and 0.39 nm from the catalytic site of the SARS-CoV-2 RdRp, with few minimal drifts during the simulation period. However, Lig-3 is situated at the average distance of 0.5 nm, showing a small sharp drift between 70 and 80 ns. Armed with these findings, it can be postulated that all of the three natural ligands are spatially will-situated at the catalytic pocket of the viral polymerase, forming stable complexes over the simulation period.

Principal component analysis (PCA) was applied to gain a two-dimensional model of the movements of SARS-CoV-2 RdRp, both in its unbound state and when bound to various ligands (Fig. 6). The unbound polymerase exhibited a wide range of positional clusters, extending from −6 to 6 on the PC1 plane (range = 12) and from −7 to 7 on the PC2 plane (range = 14); while ligand binding appeared to alter the movement range and pattern of the polymerase. For example, compared to Ribavirin, the RdRp complexes with Lig-1 and Lig-3 showed more restricted clustering, with Lig-1 ranging from −5 to 5 (PC1, range = 10), and −5 to 4 (PC2, range = 9), and Lig-3 ranging from −5 to 5 on both PC2 planes (range = 10 on two PCs). As for the Lig-2-RdRp complex, the positional clusters varied between −6 and 4.5 (PC1) and from −6 to 4 (PC2). These findings indicate that ligand binding may decrease the spatial ranges of SARS-CoV-2 RdRp, which aligns with the Rg results, suggesting that these ligands could destabilize the structure and dynamics of RdRp.

Fig. 6
figure 6

Changes in the PCA of viral polymerases (I: SARS-CoV-2 RdRp, II: HCV NS5B) in the free state and complex with Ribavirin, Lig-1, Lig-2, and Lig-3. (The plots were generated by Excel 2013, and images were prepared by Adobe Photoshop V 2022).

Full size image

MD simulations of HCV NS5B

As for free HCV NS5B, the RMSD magnitudes soared to approximately 0.22 nm after 25 ns, followed by a decline to 0.19 nm, with some variations until 80 ns (Fig. 4). Afterward, the values stabilized around an average of 0.09 nm, showing that equilibrium was achieved during the simulation period. The mean RMSD magnitudes for free viral polymerase were 0.093 nm, whereas these magnitudes for the complexes with the standard drug (Ribavirin), Lig-2, and Lig-3 soared and exhibited the RMSD average of 0.114 nm, 0.113 nm, and 0.116 nm, respectively. In contrast to SARS-CoV-2, Lig-3 and Lig-2 were revealed to have the largest RMSD values in HCV NS5B; these numbers were higher or even comparable than those of the Ribavirin. These findings are in line with the findings of molecular docking regarding with the well with the binding affinities of the ligands, confirming their different effects on HCV NS5B stability.

The Rg diagram of the free HCV NS5B and in complex with Ribavirin and Ligands 1–3 were also presented in Fig. 4. The Rg values for the free HCV NS5B moderately fluctuated and showed a small decreasing trend over the simulation time, with the Rg mean magnitude of 2.426 nm. However, the Rg magnitudes for complex composing the control drug fell, revealing the Rg mean values of 2.42 nm. Intriguingly, the Rg diagram for the other frameworks showed some significant drops, particularly in the framework comprising Lig-3 and Lig-2, with the mean Rg values of 2.411 nm and 2.417 nm, correspondingly. The Rg mean values in a complex containing Lig-1 (2.419 nm) were also comparable to the Ribavirin. These drops in the Rg numbers show the possible compression in the structure of HCV NS5B, especially in frameworks containing Lig-2 and Lig-3, destabilizing the polymerase structure and disturbing its enzymatic duties. The SASA values also saw the same decreasing trend and agreed with the Rg results (Fig. 4). The reduction in the protein’s SASA value indicates that the viral enzyme’s structure experiences conformational changes, becoming more compact after binding to these ligands. As a result, the internal parts of the protein are less exposed to the solvent. These compactness and conformational changes have a notable effect on internal water entry points, limiting water movement into the enzyme’s internal regions, which ultimately disrupts its polymerase activity. Moreover, the time evaluation plots of hydrogen bonds involved in the interaction between the studied ligands and viral polymerase were also assessed (Fig. 4). The standard drug (Ribavirin) formed two hydrogen bonds on average during the simulation time, whereas Lig-2 and Lig-3 are stabilized in the catalytic pocket of HCV NS5B through four and three H-bonds on average. Lig-1 also showed one H-bond on average during the simulation period.

The RMSF plots indicate that, consistent with the other enzyme systems examined, ligand binding led to notable alterations in the RMSF values of HCV NS5B residues (Fig. 5). For Ribavirin, increased fluctuations were observed in multiple regions of the finger, thumb, and palm subdomains. In the finger subdomain, notable increases were observed in residues 21–48, 59–79, 80–99, 104–126, 129–138, 150–152, 162–179, and 183–185. Similarly, the thumb subdomain exhibited escalated fluctuations in residues 194–205, 207–215, 221–230 (catalytic motif A), 244–251, 265–272, 283–293, 288–301, 301–310, 317–319 (catalytic motif C), 326–337, 352, and 360–362. The palm subdomain showed increases in residues 376–378, 385–398, 403–405, 431–449, 453–467, 468–483, and 512–525. Conversely, decreased RMSF values were detected for residues 8, 15–17, and 142 in the finger subdomain, along with residues 496 and 533–535 in the palm subdomain.

The RMSF plot for the complex containing Lig-1 saw some increases near residues 20, 26, 44–50, 60, 65, 84–86, 100–103, 109–117, 137, 142, 156–160, 166, 186–188 of the fingers subdomain, 206–209, 216–219 (palm subdomain), 220–225 (catalytic motif A), 247, 269–272 (palm subdomain), 280–291 (catalytic motif B), 304–310, 316–319 (catalytic motif C), 326–330, 341–355 (palm subdomain), 360–370, 387–393, 434–442, and 513–517 (palm subdomain). However, few decreases in the RMSF magnitudes were observed for residues located at the fingers subdomain (8, 19, 94, 131, 150, 243) and thumb subdomain (377, 402–405, 449, 470, 483–490, 494, and 533). Moreover, Binding of Lig-2 also caused the RMSF values to rise around residues 19–24, 29–36, 66–73, 77, 86, 111–113, 141–145 (fingers subdomain), 157–159 (NTP tunnel) 166, 195–201, 203–208, 211–215 (fingers subdomain), 219–226 (catalytic motif A), 229, 242, 248 (fingers subdomain), 271, 288–292 (catalytic motif B), 303–309, 316–319 (catalytic motif C), 327–343, 349–355 (palm subdomain), 359–364, 370–373, 388–395, 417–420, 449, 469, 479, and 515 (thumb subdomain). The RMSF magnitudes fell at residues 1, 8, 18, 41, 57, 95, 148–150, 152–154, 188 (fingers subdomain), 376, 404, 495, and 533 (thumb subdomain). After binding to Lig-3, the RMSF plots also experienced some soaring at residues 8, 26–35, 69–76, 83–95, 97–117, 124–134, 141–147, 163–183, 188–196, 202–207 (fingers subdomain), 219–225 (catalytic motif A), 255–263 (palm subdomain), 289–292 (catalytic motif B), 316–319 (catalytic motif C), 330–332, 354–358 (palm subdomain), 360–369, 388–396, 404, 436–462, 466–488, 510–519, and 532–535 (thumb subdomain), whereas the RMSF plots decrease at residues 2, 15, 43, 57, 149, 230 (fingers subdomain), 377, and 498 (thumb subdomain). In conclusion, ligand binding changed the RMSF fluctuations of the residues located at the catalytic motifs A and C, which these results align with the molecular docking findings and suggest that these ligands interfere with substrate recognition, impede the correct binding of the RNA template and positioning of NTPs, and ultimately disrupt RNA synthesis.

The stability of the studied ligands within the catalytic site of HCV NS5B was evaluated by tracking the time evolution of the average distance between the ligands and the enzyme’s catalytic pocket (Fig. S1, Supporting Information). Throughout the simulation, all ligands consistently maintained an average distance from the catalytic site ranging from 0.26 nm to 0.37 nm, with only a few minor-to-moderate fluctuations observed at 40 and 70–80 ns for Lig-1 and at 30 ns and 60 ns for Lig-2.

The 2D model of HCV NS5B movements, both in its free form and when complexed with the ligands, was obtained through PCA (Fig. 6). The movement behavior of the HCV NS5B enzyme, both in its unbound state and when bound to ligands, was analyzed using PCA (Fig. 6). The free HCV NS5B demonstrated positional clusters extending between −6 and 4 (PC1; range = 10) and from −4.5 to 8 on the PC2 plane (range = 12.5); whereas the binding of the studied ligands changed the mobility range and pattern of the enzyme. For example, after binding to Ribavirin, the mobility range of the protein changed from −5 to 3.5 (PC1; range = 8.5) and −5 to + 4.5 on the PC2 plane (range = 9.5). Compared to Ribavirin, the most severe deformation in the protein mobility range and pattern can be detected in complexes containing Lig-3, followed by Lig-2. For example, Lig-3 binding caused some changes in the mobility range between −5 and 4 (PC1; range = 9) and from −4 to 4 on PC2 (range = 8), while binding of Lig-2 altered the movement range from −4.5 to 5 and between −4 and 4 on PC1 and PC2 planes. Lig-1 also could change the protein movement range from −4.5 to 4 (PC1; range = 8.5) and from −5.5 to 5 (PC2; range = 10.5). The results indicated that ligand binding could reduce the spatial range of HCV NS5B movements. This observation aligns with the Rg analysis, implying that these ligands might disrupt the enzyme’s structural stability and dynamic behavior.

Finally, MM/PBSA approach was utilized to calculate the binding energies for the dynamic interactions between the Ribavirin and our studied ligands (Ligs 1–4) and the viral receptors (Table 3). In terms of the complexes with SARS-CoV-2 RdRp, Lig-1 and Lig-3 make the strongest complexes with the studied polymerase since they have shown the higher binding energies of − 143.90, − 129.64 kj/mol than that of the Ribavirin (− 103.06 kJ/mol). Moreover, it should be noted that the majority of the binding energies lies in van der Waals forces, and also the largest shares of van der Waals and electrostatic interactions stand for Lig-1 and Lig-3. In comparison to SARS-CoV-2 RdRp, the strongest binding energy in complex with HCV NS5B attributes to the Ribavirin (− 187.01 kj/mol), which is followed by Lig-3 (− 166.48 kj/mol) and Lig-2 (− 159.64 kj/mol). Furthermore, the main driving forces are van der Waals interactions in all of the systems. The biggest share of electrostatic interactions is related to the Ribavirin and Lig-3, while the largest contribution of van der Waals interactions are seen in the systems containing the Lig-2, Lig-3 and Ribavirin.

Table 3 The mean interaction energies of the drug Ribavirin and Ligs 1–3 and to the viral targets (i.e., SARS-CoV-2 RdRp and HCV NS5B) in a dynamic state.
Full size table

The pharmacokinetic properties of the three fungal secondary metabolites were analyzed using the SwissADME ProTox-II online webservers (Tables 4 and 5). Bioavailability is essential because drugs with low bioavailability fail to reach the bloodstream and cannot perform their intended functions effectively. It is important to note all of the three chosen natural compounds met the Lipinski rule of five, except for Lig-1 showing some deviations, demonstrating strong drug-likeliness characteristics and a bioavailability score of 0.55. Sufficient intestinal absorption also indicates high oral bioavailability. Interestingly, except for Lig-1, none of the three fungal compounds revealed any drug interactions with the other pharmaceuticals and would not play as a substrate for liver enzymes (Cytochrome P-450). None of these ligands can’t pass the blood–brain barrier and central nervous system. However, as far as the viral polymerases are not in the nervous system, it is not necessary for these inhibitors to cross the blood–brain barrier. Importantly by far, they have shown very low toxicities and revealed a high magnitude of LD50 (5000–10,000 mg/kg). To sum up, it can be concluded that Lig-3 and Lig-2 are the most suitable drug candidates among all studied fungal compounds, thereby warranting further investigations as multi-target agents against aspartyl polymerases; HCV NS5B and SARS-CoV-2 RdRp.

Table 4 Drug likeliness of the three selected fungal compounds (Swissadme).
Full size table
Table 5 Pharmacological features and toxicity prediction for the three selected fungal compounds.
Full size table

Conclusion

This study investigated the antiviral potential of 174 secondary metabolites from the Sordariales order against aspartyl polymerases; HCV NS5B and SARS-CoV-2 RdRp. A two-step virtual screening was conducted, beginning with blind docking to identify ligands that will bind at the active site, followed by targeted docking using Ribavirin and natural compounds with superior binding affinities. Blind docking revealed that 76 out of the studied ligands were bound to the catalytic pocket. Among these, 10 ligands showed binding energies below −7 kcal/mol, according to the targeted docking. Regarding the molecular docking results, Ligands 1–3 indicated superior binding affinities to both viral polymerases compared to Ribavirin as the control drug, and Lig-3 showed the lowest binding energies among the studied ligands. These ligands form hydrogen bonds and hydrophobic interactions with key residues located at the catalytic motifs A, C, and F in SARS-CoV-2 RdRp and with catalytic motifs A, B, and C in HCV NS5B, showing their high potential in disrupting the viral replication by inhibiting polymerases’ catalytic functions.

Finally, these three fungal metabolites were subjected to 300 MD simulations to further analyze their effects on the structure and dynamics of the viral polymerases. As for SARS CoV-2 RdRp, the chosen three ligands interacted with catalytic motifs A and F. Interestingly, RMSF analysis revealed that ligand binding caused significant changes in the fluctuations degrees of the residues located at the catalytic motifs A and F, which agree with the molecular docking results. These conserved regions are responsible for NTP binding and template binding, therefore these phenomena can cause important disruption in the catalytic function of SARS CoV-2 RdRp and hinder the replication of the virus. According to the molecular docking results for HCV NS5B, it can be asserted that the three selected ligands basically interacted with the catalytic motifs A and C. In this line, ligand binding changed the RMSF fluctuations of the residues located at these two catalytic motifs, which these results align with the molecular docking findings and suggest that these ligands interfere with substrate recognition, impede the correct binding of the RNA template and positioning of NTPs, and ultimately disrupt RNA synthesis.

The RMSD analysis also taken from the MD simulations of SARS-CoV-2 RdRp indicated that the rearrangements in the protein were very fast, and compared to Ribavirin, Lig-1 and Lig-3 may cause more significant deviations. Further, Rg and SASA analyses have indicated structural compression for all the complexes, particularly, in the presence of Lig-1 and Lig-3. PCA has further presented restricted movements in the ligand-bound complex, which further corroborates the hypothesis that these ligands destabilize the polymerase structure. As for HCV NS5B, RMSD analyses for free HCV NS5B indicated the equilibration of the protein within the simulation period, while complexes with Ribavirin, Lig-2, and Lig-3 showed higher RMSD values. With the binding of Lig-2 and Lig-3, there is a decrease in the values of Rg and SASA, indicating that there is significant structural compression and, thereby, destabilization of the HCV NS5B structure. Compared to Ribavirin, Lig-2 and Lig-3 have formed more hydrogen bonds, enhancing their stability in the catalytic pocket. From the RMSF and PCA analysis, the increased fluctuations with reduced protein mobility upon the addition of the ligands highlight the interference with RNA synthesis, disturbing the structural integrity and dynamics of the enzyme. MM/PBSA analysis revealed that for Lig-1 and Lig-3 demonstrated the strongest binding energies in complex with SARS-CoV-2 RdRp compared to Ribavirin. Most binding energies stem from van der Waals forces, with Lig-1 and Lig-3 exhibiting the highest contributions. In contrast, Ribavirin showed the strongest binding with HCV NS5B, followed by Lig-3 and Lig-2. Overall, van der Waals interactions were the main driving forces across all systems, with Ribavirin and Lig-3 contributing the most electrostatic interactions. Overall, compared to the Ribavirin, Lig-1 and Lig-3 demonstrated greater structural impact against SARS-CoV-2 RdRp. At the same time, Lig-3 and somehow Lig-2 showed stronger structural compression in the structure of HCV NS5B, suggesting they may block enzymatic function and viral replication more effectively than Ribavirin. Armed with the results of ADME/T properties, Lig-3 can be nominated as a multi-target antiviral agent aspartyl-polymerases of HCV and SARS-CoV-2, being on the call for future in vitro and in vivo studies.