1.
Introduction
It has been more than three years since the global community has been negatively impacted in many aspects of life as a result of COVID–19, and the future of this disease appears bleak. COVID–19's eradication may be threatened if a new strain emerges that can overcome the current vaccines, much like the flu viruses can. COVID–19's many diseases have prompted requests for complementary therapies to address them. For a successful COVID–19 treatment, the drug or drug regimen must overcome not only inhibiting virus proliferation but also complications due to infection and, in some cases, its accompanying maladies, such as adverse events following immunization (AEFI), which has been observed in several cases of COVID–19 vaccination [1],[2] and also in numerous historical vaccines usage [3]–[6].
Stopping virus replication is only half the story; the other half is minimizing the development of cytokine storms and thrombocytopenic events. Furthermore, it is critical to restore platelet levels after a thrombocytopenic event, which individuals with severe COVID–19 are more likely to experience than those with milder forms of the disease, to prevent hypercoagulation and other thrombotic consequences. Furthermore, reports on spike protein toxicity [7]–[12] may complicate spike–protein–based vaccines and vaccine–induced immune thrombotic thrombocytopenia (VITT) [13]–[15] have slowly trickled in [14]–[18] meaning a battery of drugs/treatment may be needed to address all of these issues. Because each component of the medication would have its own set of safety and efficacy concerns, the overall issue presents a mountain of difficulties.
There is currently an increase in the number of clinical trials using repurposed drugs to supplement vaccine use. Recent advances in the repurposing of drugs originally approved for other uses as COVID–19 treatments have introduced drugs, such as ivermectin [19], chloroquine, hydroxychloroquine [20] and cyclosporin [21]. In addition, the highly popularized use of the repurposed drug ivermectin has been recently found in a large-scale trial in Brazil to have no benefit to COVID–19 patients [22]. Furthermore, to complement current and future pharmaceutics–based approaches, there has been an increasing call to mobilize the use of antioxidants, functional foods, natural bioactive compounds [23],[24] and herbal–based antivirals to combat COVID–19 [25]–[27]. Due to this, the use of other approaches, including natural–based therapies, to complement vaccine research is needed. Because of their high absorption and low toxicity, phytocompounds may be the most effective therapeutic candidate available at the moment, particularly compounds derived from edible plants or plants that have been traditionally used for medicine for hundreds of years. SARS-CoV-2 virus inhibition is one of the potential targets of phytocompounds. Since its inception, the area of computer–aided drug design (CADD) has been dominated by methods that allow compounds to bind to protein targets [28]–[31] and even to understand the molecular mechanisms of several biological processes in the cell [32]–[35]. Most of the screening works to find potential SARS–CoV–2 inhibitors nowadays are carried out through the in silico approach [36]–[38].
SARS–CoV–2 has a genome that is entirely made up of single–stranded RNA and produces numerous polypeptides. The genes encode for structural proteins such as N–protein that binds the viral RNA genome to form the ribonucleoprotein complex, RNA–dependent RNA polymerase (RdRp) that is important in viral proliferation, spike protein that mediates viral entry into the host, matrix, envelope proteins, and chymotrypsin–like (MPro) proteases, such as 3CLpro, which are involved in viral replication and can cleave a variety of non–structural proteins that are essential for viral replication [38]–[40]. The most often viral targets for antiviral purposes are RNA–dependent RNA polymerase or RdRp, MPro and spike protein [27],[41]–[45].
Before the onset of the COVID–19 pandemic, one of the top virus death cases in Malaysia is caused by the dengue virus, which is a positive–sense single–stranded RNA virus similar to SARS–CoV–2 in COVID–19. More than 80,000 cases are reported in 2019 with more than 100 deaths reported [46]. One of the top herbal remedies that emerge as a prime weapon for the dengue viral infection is Carica papaya. Research has demonstrated that papaya leaves extract (PLE) decreases dengue complication through another route by inhibiting viral production. A significant lowering of the intracellular viral load supports papaya leaves extract antiviral activity [47]. To date more than 1000 dengue patients in several countries have been included in clinical and pilot studies to test the therapeutic effect of PLE in dengue will an overwhelmingly positive response on the efficacy of PLE as far as reducing dengue complications is concerned [48]–[60].
The rationale for this work covers the interesting observation that dengue and COVID–19 share many similarities and that the use of bioactive compounds from PLE, a traditional herbal preparation for dengue that have been proven effective in several clinical studies, might be useful for COVID–19. In addition, compared to synthetic drugs, herbal preparations rarely cause adverse effects or mortality, except in certain cases where the herbs have been misidentified [61]. The SARS-CoV-2 and the dengue virus (DENV) both enter the body through separate entry points, but both cause a systemic infection and have some of the same clinical manifestations, including fever, headache, myalgia, and gastrointestinal symptoms. Similarities in risk factors for severe illness, endothelial dysfunction, cytokine storms, and multi-organ failure follow the early similarities in clinical presentation [62]. A proinflammatory immune response and a delayed and weakened type I IFN response are hallmarks of both illnesses. Low platelet counts, which are common in dengue illness, are shown to be linked with a more than fivefold heightened risk of severe COVID–19 infection or in vaccinated people [63]. In a subset analysis, platelet counts are shown to be substantially lower in critically sick patients, and this was related to higher mortality [64]–[66]. Platelets expressed the ubiquitous ACE2, which is a major viral entry receptor for SARS–CoV–2, as well as TMPRSS2, a serine protease that is needed in the priming process of the spike protein. Aside from an increase in the fibrin degradation product D–dimer, COVID–19 patients also exhibit prolonged prothrombin time and decrease antithrombin, which is also seen in dengue infection [67]. Both lead to an increased risk of thromboembolic disease and bleeding with the development of disseminated intravascular coagulation (DIC) is reported in as many as 71% of non-survivors of COVID–19 [66]. There have been calls for the use of anti–TNF therapy as a possible treatment for COVID–19 [68],[69]. The possible use of papaya leaves extract (PLE) in dengue virus inhibition and in alleviating cytokine storm has been demonstrated in dengue infection in mice models [70]. The role of PLE as an inhibitor to TNF–alpha [71] may play a bigger role than its ability to inhibit IL–6 as it is more upstream than IL–6 in the inflammation pathway cascade [72]. The fact that PLE can alleviate the severity of dengue in clinical studies through possibly antiviral, anti–TNF–alpha and also increasing platelet count may come in handy in dealing with co–infection of COVID–19. Furthermore, severe and long-term COVID-19 and VITT are characterized by the formation of micro- and macrothrombi, resulting in severe thrombocytopenia [73], of which PLE components, such as alkaloids, cardiac glycosides, flavonoids, phenols, saponins, tannins, and terpenes, exhibit antiviral and immunomodulatory properties and can play a significant role in alleviating this event, according to several studies [54],[74],[75].
To date, no in vitro or in vivo studies have been carried out on the potential of PLE as a weapon to combat COVID–19 in a holistic approach. In April 2020, we suggested that an in silico study be first carried out on bioactive compounds from PLE as potential SARS–CoV–2 viral inhibitors [76]. One study has taken up the challenge and targeted SARS–CoV–2 spike protein and RNA–dependent RNA polymerase (PDB code 7BW4). A relatively good docking result was obtained with significant binding energies from –4.20 to –8.90 Kcal/mol and was obtained with the best compound being lutein –8.90 (Kcal/mol) for spike protein and lycopene (–8.71 Kcal/mol) for RNA–dependent RNA polymerase [77] indicating their potential as SARS–CoV–2 inhibitors. However, as no molecular dynamics simulations were carried out, it is difficult to prove the stability of the obtained static poses [78]. With this in mind, we omit any current research using PLE or any other similar phytocompounds that do not carry out MDS works in our discussion and only concentrate on those that did.
In this investigation, significant components from PLE and numerous previously reported SARS–CoV–2 and SARS–CoV inhibitors were evaluated in silico against four SARS–CoV–2 targets; RNA binding domain of nucleocapsid, MPro, RdRp, and spike protein (Wuhan, Delta, and Omicron variants) that are routinely utilized to evaluate potential inhibitors of the virus [37],[79]–[82] and two human targets; TNF–alpha and thrombin, that may lessen or battle the severity of COVID–19 disease. Molecular dynamics simulations using RMSD, RMSF, and a few more analyses were then utilized to validate the interactions between these compounds and their targets. Our findings demonstrate that PLE has remarkable potential in combating the entire spectrum of COVID-19 disease, from viral inhibition to the alleviation of COVID-19 complications in humans.
2.
Materials and methods
2.1. Datasets and ligand preparation
The crystal structures of SARS–CoV–2 protein targets were taken from the protein data bank. The PLE ligands (twenty) were taken from the National Library of Medicine (https://pubchem.ncbi.nlm.nih.gov) and the Laboratories in the Department of Pharmaceutical Chemistry at the University of California, San Francisco (UCSF) or Zinc15 (https://www.zinc.docking.org). The protein receptor is edited by using Pymol (https://www.pymol.org) where water molecules and ions were removed from the protein molecule. Further editing of protein receptor and ligand by adding several parameters before docking was done using the AutoDock tool (http://mgltools.scripps.edu/). The tool was utilized for the addition of hydrogen atoms, merging of nonpolar hydrogens and addition of polar hydrogens to the protein. Subsequently, the protein was saved into a dockable pdbqt format [30].
Ligands were converted into 3D structures and geometrically optimized through energy minimization using the Density Function Theory (DFT) energy minimization protocols of Gaussian16 [83]. By using the AutoDock tool, water was deleted from the protein receptor, and then the hydrogen atoms and the Kollman charges were added. In the ligand preparation, the non–polar hydrogen atoms were merged and the Gasteiger charges were added automatically by the AutoDock tool [34].
2.2. Molecular docking
We utilized Autodock4, which uses a semi–empirical free energy force field to estimate the free binding energy between the binding of two molecules or more in the water [84],[85]. The simulation uses the pair–wise atomic terms during the evaluation of the interaction of the molecules. The force field comprises six pair wise evaluations (V) and the estimated conformation entropy lost during binding is ΔSconf. The estimated free energy of the binding relationship is given as follows:
where: L represents Ligand and P represents Protein.
The estimated loss of torsional free energy upon binding is proportional to the number of rotational bonds (Ntors). The equation is given as follows:
The pair–wise evaluation (V) is given as follows
where W is the weighting constant, E(t) is the angle t dependent hydrogen bond interaction and q is the atomic charge.
The nucleocapsid phosphoprotein crystal structure of SARS–CoV–2 was taken from the protein data bank (6VYO.pdb). For N–protein, the grid box was set to 40 × 60 × 56 Å grid points with 1.0Å spacing. The RdRp's pdb file(7OZU.pdb) was further processed by removing the ligand Molnupiravir/NHC. For RdRp, the grid box was set to the coordinates 82 × 64 × 114 (xyz) (grid spacing 0.375Å) with grid center of 78.711, 96.577, 112.936 (xyz). The main protease MPro/3CLpro (pdb file 7AEH.pdb) was further processed by removing the ligand Fab 2B04 and DMS. The grid box was set to the coordinates of 80 × 108 × 92 (xyz) (grid spacing 0.375Å) with grid center of 11.946, –7.105, 20.021 (xyz). The spike protein was downloaded from the RCSB protein data bank. The spike protein's receptor–binding domain (RBD) pdb file for the Wuhan (7K9I.pdb), Delta (7W9F.pdb) and Omicron variants (7WBP.pdb) were further processed by removing the neutralizing Fab 2B04 for Wuhan variant, the 8D3 monoclonal antibody for the Delta variant and the ACE2 human receptor for the Omicron variant, respectively. The grid box (Omicron) configuration was 62 × 82 × 72 (xyz) (grid spacing 0.375Å) with grid center of –35.572, 26.541, 5.459 (xyz). The TNF–Alpha's pdb file(2AZ5.pdb) was further processed by removing the inhibitor CHEMBL255489. The grid box was set to the coordinates of 78×60×66 (xyz) (grid spacing 0.375Å) with grid center of –13.687, 71.607, 27.002 (xyz). The Alpha Thrombin's pdb file (3VXF.pdb) (heavy chain) was further processed by removing bivalirudin. The grid box was set to the coordinates of 78 × 82 × 86 (xyz) (grid spacing 0.405Å) with grid center of 0.725, 25.026, 24.793 (xyz) (Table 1). Prior to docking, Genetic Algorithm was set using the AutoDock tool. The number of GA Runs was set to 50 with a population size of 300. The docking was done using the Autodock4 and was executed in several batches and the lowest free energy of binding was selected. The ∆G value was calculated by summing the intermolecular energies (van der Waals+hydrogen bond+dissolving+electrostatic energies) minus the torsional free energy of the substrate [86]. The methodology was validated by redocking the inhibitor D–Phe–Pro–Arg chloromethylketone (PPACK) into the human alpha–thrombin active site of (PDB: 1PPB) (RMSD = 1.11 Å) with the ligand superimposed to the native co–crystallized enzyme [87]. Docking and Gromacs were carried out on an Intel 11th Gen Core i7 11700 Processor, with two Nvidia GeForce RTX3070 GPUs using an OS ubuntu desktop 20.04 LTS.
2.3. Molecular dynamic simulation
Molecular dynamics simulations (MDS) were performed using GROMACS (2021 release) in accordance with established protocols. Using 100 ns MD simulations, we compared targeted proteins with or without ligands. The CGenFF service and the 'pdb2gmx' program were used to generate the topology files for the ligand and protein, respectively, for use with the CHARMM36 force field, which is officially supported by GROMACS [88]. Reconstitution of ligand topologies to the refined protein structure completed the system's assembly [29]. The TIP3P water model was then employed to construct a water-solvated system with periodic boundary conditions based on a dodecahedron, which is closest to the ideal spherical shape. A water solvated system was constructed using the TIP3P water model with periodic boundary conditions based on the dodecahedron. The solutes are centered in the simulation box, no more than 10 from the border edges (1.0 nm). Na+ counter-ions were added via the ‘gmx genion’ script to achieve solution neutrality. Thermodynamic parameters such as pressure, temperature and density must be maintained by the barrier to contain molecules in simulation [30]. The complexes were energy minimized at 10 kJ/moL utilizing the steepest descent algorithm for a maximum step of 50,000 based on the Particle Mesh Edward (PME) Coulombic long-range electrostatic interactions using the Verlet cut off-scheme. The system equilibrium was attained in two stages. In the first step, NVT equilibration was performed at 300 K and steps of 5000 ps, and NPT equilibration was performed at 1 bar reference pressure using the Parrinello-Rahman (pressure coupling) and steps of 5000 ps in the second step. Finally, a 100-ns MD run was performed with a time step interval of 2-fs on the protein-ligand and protein complexes involved. Following the completion of MDS, the MD trajectories were analyzed using the GROMACS' g_rmsf and g_rms tools to determine the root-mean-square fluctuation (RMSF) and root-mean-square deviation (RMSD), respectively [20],[29],[89].
2.4. Binding free energy calculation using MM-PBSA
To understand the stability of a given protein-ligand complex, a free energy calculation provides a quantitative estimate of the interactions between the protein and ligand. Using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) method [90], we determined the binding free energy of each protein-ligand complex, which incorporates both the potential energy (Van der Waals and electrostatic interactions) and the free solvation energy (nonpolar and polar solvation energies).
The binding energy is calculated by using the following equation:
where, ΔGbinding = the total binding free energy of the complex, Gcomplex = the total binding free energy of the complex, Greceptor = the binding free energy of free receptor and Gligand = the binding free energy of unbounded ligand.
Furthermore, the free energy for each individual entity can be given by
where, x is the protein−ligand or protein or ligand complexes, <EMM> is the average molecular mechanics potential energy in a vacuum, T and S denote the temperature and entropy, respectively, TS is the entropic contribution to the free energy in a vacuum, Gsolvation is the free energy of solvation.
The MM-PBSA binding free energy calculation was done with ‘g_mmpbsa’ script. The input files used are the mds trajectory file *.xtc, topology-parameter file *.tpr and atom-index file *.ndx [90]. The MM/PBSA calculations provide additional insight into the binding-free energy estimations, including the degree of ligand-protein affinity, the type of interaction, and the residue-wise contributions [35].
3.
Result and discussion
The identification of novel ligands is necessary for medicines that target particular receptors. As part of the journey to drug targets and interactions in the body, pharmacokinetic events, such as Absorption, Distribution, Metabolism and Excretion or ADME, are involved when creating an orally active drug, Lipinski's rule of 5 or Ro5 is another typical discriminating step. Passive diffusion is the only plausible route for drugs to enter cells, according to the rule of five, which ignores the role of transporters. Only around half of the new chemical compounds that are taken orally adhere to the so–called “rule of five”, according to O'Hagan and his colleagues. Many of the 303 compounds being investigated in clinical trials have molecular weights of more than 500 Dalton, including 182 approved drugs prior to 2015 [97]. The most prevalent indications for the use of these drugs were cancer, infection, and cardiovascular disease. As the rule does not apply to natural products [98] we skipped the Ro5 from our screening.
3.1. Molecular docking studies
3.1.1. SARS–CoV–2 RNA binding domain of nucleocapsid molecular docking studies
The docking results show that PLE bioactive compounds are dominant in the top ten best compounds with protodioscin showing the strongest binding with the nucleocapsid protein residues based on the free energy of binding value of –13.83 kcal/mol. The calculated inhibition constant was 73.28 pM. The residues interacting with N protein were Chain B: Leu167, Pro168, Lys169, Gly170, Phe171, Tyr172, Ala173 and Chain C: Pro73, Ile74, Thr76, Gln83, and Thr135 (H bonding residues are in bold) (Table 2) and visualized more comprehensively using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figures 1, clitorin and glycyrrhizic acid, Figures S1 and S2 respectively). In terms of H-bond formations to the main compound protodioscin, this compound forms H bond to chains B and C of the SARS-CoV-2 N-protein with Thr135 forming a H bond at chain C to the alkyl glucoside or more specifically at the o-glucopyranosyl moiety while it forms H bonds with chain B's Leu167 and Phe171 at the 6-deoxy mannopyranosyl moiety of protodioscin (Figure 1(d)). Some of the N–protein residues interacting with protodioscin were also reported to interact with the FDA–approved drug paritaprevir, with the latter residues interacting with SARS–CoV–2 N–protein are Gly170, Ala173, Leu161, Gln163, Thr165, Pro168, Lys169, Phe171 and Tyr172 (H bonding residues are in bold) [44]. The residues Gly170, Phe171 and Tyr172 are reported to significantly affect the dynamics and binding of the RNA to N–protein [99]. Clitorin, the second-best compound binds to N protein and form the same interaction with protodioscin and glycyrrhizic acid (third best compound) at the amino acid residues Thr76, Leu167, Gly170 and Try172. In addition, glycyrrhizic acid, the third–best compound interacted with all these residues and has been reported to also bind to SARS–CoV–2 nucleocapsid [100]. A recent docking exercise among the repurposed drugs such as hydroxychloroquine, chloroquine, ivermectin, remdesivir and favipiravir shows that ivermectin exhibited the highest binding affinity to the predicted active site of SARS–CoV–2 N protein forming four H–bonds (Gln160, Leu161, Gly164, and Thr166) [101], which are very similar to the results found in this study (Table 2).
3.1.2. SARS–CoV–2 RdRp Molecular docking studies
Protodioscin, a PLE compound from the steroidal saponins group that exhibited the highest docking score interacted with the residues Asn496, Asn497, Lys500 and Lys577 (H bond in bold), which include many of the similar residues to streptolydigin and other published RdRP potential inhibitors [109]–[111]. More comprehensive interactions were visualized using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figure 2, manghaslin and kaempferol–3–(2G–glucosylrutinoside), Figures S3 and S4, respectively). In terms of H-bond formations for the main compound, protodioscin, Lys577, Asn496, template C11, U12, G13 formed H bonds with the alkyl glucoside or more specifically at the o-glucopyranosyl moiety of protodioscin, while Asn497, Lys509, Lys545, templates G7, G8, A10, G13, 70K9, products G10, C11 and C11, formed H bonds with the 6-deoxy mannopyranosyl moiety of protodioscin (Figure 2d). All the three top ranking compound shares a similar interaction to Product G10, whilst manghaslin and kaempferol–3–(2G–glucosylrutinoside) interacted with the same amino acid residues Leu854 and Glu857 and the template G15 (Table 3). Interestingly, Vitamin B12, which came out 2nd highest from a screening exercise out of 8,285 compounds to find molecules that can bind to the active site of SARS–CoV–2 RdRp, interacts with all of the residues reported in this study among others [36]. In another docking study, streptolydigin, the antibiotic that works by inhibiting nucleic acid chain elongation by binding to RNA polymerase was evaluated for its potential RdRp inhibition. Streptolydigin was found to bind to nsp–12 active center cleft with free binding energy (Autodock) of –8.11 kcal/mol. Streptolydigin makes four hydrogen bonds with Lys577, Tyr689, Thr591, Ser592 and fifteen hydro–phobic interactions with Asn496, Gln573, Leu576, Ala580, Val588, Ile589, Lys593, Phe594, Trp598, Met601, Ala688, Leu758, Phe812, Cys813 and Gln815 [112]. The ellagitannin punicalin shows the best score in a docking study against SARS–CoV–2 RdRp with the interacting residues are Asn497, Gly590, Asp684, Tyr689, Ile494, Lys577, Asp684 and Ala685 [113].
The SARS–CoV–2 RdRp contains a nidovirus RdRp–associated nucleotidyltransferase (NiRAN) domain that assists in the transfer of nucleoside monophosphate to the RNA and catalytic site. Remdesivir triphosphate and other analogues were discovered to bind not just to the catalytic site but also to the interior of NiRAN, inducing allosteric changes in nsp–12. It has also been claimed that remdesivir triphosphate binds near the catalytic site in the RdRp of SARS–CoV–2, inhibiting the activity of enzymes [114]. Remdesivir is amongst the few FDA–approved SARS–CoV–2 inhibitors that shows appreciable free energy of binding (–8.62 kcal/mol) (Table 3) and the value is within the range of several studies that reported values (Autodock) such as −4.68 [115], –6.1 kcal/mol [116] and −7.4 kcal/mol [115]. The difference in docking score using the same software might be attributed to the use of different version and in some instance the 2D structure of compounds not converted to 3D and not energy minimized before docking (Uttam Pal and Justin Laemkul, pers communications). In comparison to nucleoside inhibitors, non–nucleoside inhibitors including PLE compounds are less toxic and have fewer side effects.
3.1.3. SARS–CoV–2 main protease MPro/3CLpr molecular docking studies
Protodioscin from PLE was again the best compound for potentially inhibiting SARS-CoV-2's MPro. The residues interacting with protodioscin (H bond in bold) were Thr24, Met165, Glu166, Leu167, Pro168, and Gln189 (Table 4) indicating the ability of protodioscin to potentially inhibit the action of SARS–CoV–2 MPro. Of these interacting residues, N–finger of each of the two MPro protomers interacts with Glu166 of the other protomer and thereby helps shape the S1 pocket of the substrate–binding site [80]. More comprehensive interactions were visualized using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figure 3, rutin and kaempferol–3–(2G–glucosylrutinoside), Figure S5 and S6, respectively).
Rutin, the second-best compound binds to MPro and form the same interaction with protodioscin and kaempferol–3–(2G–glucosylrutinoside) (third-best compound) at the amino acid residues His164, Glu166 and Gln189 (Table 4). In terms of H-bond formations for the main compound, protodioscin, Asn142, His164 and Glu166 of MPro interacted with the tetrahydrofuran moiety of the furostane skeleton of protodioscin while Gln189 and Thr190 formed H bonds with the alkyl glucoside or more specifically at the o-glucopyranosyl moiety of protodioscin (Figure 3d). Ivermectin, which has also been demonstrated to bind to MPro in docking studies and in vitro studies with values of free energy of binding (Autodock) of –9.3 kcal/mol [101],[119],[120] shows similar free energy of binding of –8.96 kcal/mol in this study (Table 4). Although the results from Autodock indicate that three hydrogen bonds were forming between protodioscin and MPro (Table 4), Ligplot+ result indicates the presence of more H bonds (albeit weaker interaction or >3 Å). Protodioscin appears to potentially form hydrogen bonds with Thr24, His164, Glu166, Leu167 and Gly170. Telaprevir, a HIV repurposed drug shows in vitro inhibition of SARS–CoV–2 growth in Vero–E6 cells [121]. MD simulation of this drug shows that the main binding residues inside the MPro pocket for telaprevir are His41, Glu166 and Gln189 which occurred between 80 and 99% of the simulation time [121].
The 3–chymotrypsin–like protease enzyme (3CLpro), also known as the “major protease” (MPro), is an appealing target in drug discovery for the treatment of coronavirus infection since MPro lacks a corresponding human enzyme, limiting the likelihood of fashioned inhibitors binding to other human proteases [122],[123]. Crystal structure of SARS–CoV–2 main protease shows residues such as Thr24, Thr25, Thr26, Leu27, His41, Met49, Tyr54, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, Leu167, Pro168, Asp187, Arg188 Gln189, Thr190, Ala191, and Gln192 are a key target for potential inhibitors [122]–[125] among which, His41 and Cys145 are two key hydrophobic catalytic residues [80]. A docking exercise using repurposed drugs such as four HIV protease inhibitors darunavir, lopinavir, nelfinavir, and ritonavir shows that all of these compounds interact hydrophobically with one of the key catalytic residues, which is His41 [125]. Protodioscin interacts with all of these residues, forming an H bond with Glu166 and hydrophobic interactions with His41, Cys145 and Gln189 (Figure 3(d)).
3.1.4. Spike protein
The docking result shows that the protodioscin formed the strongest binding affinity with Wuhan's spike protein followed by Delta and then Omicron (Table 5). Spike protein residues that are less than 3Å from the ACE2 contribute to the firm binding [79] and interacted with protodioscin (Table 5) and more comprehensive interactions were visualized using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figure 4, clitorin and manghaslin, Figures S7 and S8, respectively). In terms of H-bond formations of the main compound protodioscin to chainA of spike protein (Omicron variant), protodioscin forms H bonds with Arg403, Asp405, Val503, Gly502, Gly504 and His505 at the alkyl glucoside or more specifically at the o-glucopyranosyl moiety of protodioscin while it forms H bonds with Tyr453, Arg493 and Ser494 specifically at the 6-deoxy mannopyranosyl moiety of protodioscin (Figure 4(d)).
Clitorin, the second-best compound, binds to spike protein (omicron) and forms the same interaction with protodioscin and manghaslin (third-best compound) at the amino acid residues Ser496, Arg498, Tyr501 and Hs505 (Table 5). Potential residues for inhibitory activity were Tyr453 (Wuhan and Omicron variants), Leu455 (Wuhan and Delta variants), Phe456 (Wuhan and Delta), Arg457 (Wuhan), Gly476, Phe486 (Wuhan), Tyr489 (Wuhan), Gln493 (Wuhan and Delta variants), Gly496 (Omicron variant, Gly496 replaced with Ser496), Gln498 (Omicron, changed to Arg498), Asn501 (Omicron, changed to Tyr501), Gly502 (Omicron variant) and Tyr505 (Omicron, changed to His505) [79],[128]. The results suggest protodioscin can strongly interact with spike protein for all variants at residues that are vital or hotspots for its binding to ACE2. Blocking the interaction of ACE2 to spike protein hotspots may prevent the entry and fusion of SARS–CoV–2 [23].
Any docking works on the spike protein are particularly interesting and important for a number of reasons, including the fact that the spike protein produced by viral infection has been shown to bind to pericytes, the cells that line the heart's tiny arteries, and is also present in the brain. This binding triggers a cascade of events that can impair normal cellular function and produce inflammatory chemicals. This happened even after the virus was no longer linked to the protein [8],[9]. Furthermore, nearly all mRNA–based vaccines now use the spike protein as the primary immunogen, and side effects from spike protein–based vaccine complications are beginning to pile up in the literature [129]–[132]. More evidence is needed to understand the cause of vaccine–induced immune thrombotic thrombocytopenia (VITT), with some researchers pointing to the spike protein as a significant contributor, while others disagree. Furthermore, targeting spike protein is one of the most important aspects of combating COVID–19, either the infection itself or the increasingly reported probable toxicity of the spike protein, which has been proposed as one of the mechanisms that can explain “long COVID” [133]–[136].
The spike protein of WT SARS–CoV–2 has 1273 amino acids and its RBD is composed of residues 319–541 and RBM is of residues 437–507. There are twenty residues in the spike protein that are firmly bound to the ACE2 and they are Val445, Gly446, Leu455, Thr478, Glu484, Gly485, Phe486, Tyr489, Gly496, Gln498, Thr500, Asn501, Tyr505, Gln493, Gly502, Phe456, Tyr449, Phe490, Asn487 and Ser494 [137]. Residues that are less than 3Å from the ACE2 contributing to firm binding are Leu455, Phe456, Arg457, Gly476, Phe486, Tyr489, Gln493, Gly496, Gln498, Thr500, Asn501, Gly502 and Tyr505 [79]. According to Jun et al, there are 14 shared amino acid positions shared by both RBMs (Receptor Binding Motif of Sars Cov–2 and Sars Cov) in the interaction with the ACE2, eight of them having similar residues between the two RBDs, they are Tyr449–Tyr436, Tyr453–Tyr440, Asn487–Asn473, Tyr489–Tyr475, Gly496–Gly482, Thr500–Thr436, Gly502–Gly488 and Tyr505–Tyr491 The receptor–binding motif or RBM that directly interacts with the ACE2 includes the pair of Cys480–Cys488 that joins the loops in the extremity of the RBM [43].
The Delta variant has the mutations L452R and T478K where Lys452 and Thr478 are replaced to Arg452 and Lys478, respectively. The Omicron variant (B.1.1.529) is the current variant of major concern to the world because of the numerous alterations that could affect transmissibility and immune evasion. It has produced many lineages including Arcturus (sublineage XBB.1.16) [138]. On its spike protein, Omicron has up to 30 single point alterations, three deletion mutations, and one insertion mutation as compared to the wild type (WT). Surprisingly, 15 mutations have been found in the Omicron receptor–binding domain (RBD), ten of which are in the receptor–binding motif (RBM) at which human angiotensin–converting enzyme 2 (ACE2) and the majority of monoclonal antibodies (mAbs) engage directly [139]. The 15 mutations of RBDOmicron are not evenly distributed in RBD, but rather crowed in its RBM with 10 residues, viz., N440K, G446S, S477N, T478K, E484A, Q493K, G496S, Q498R, N501Y and Y505H or Asn440, Gly446, Ser477, Thr478, Glu484, Gln493, Gly496, Gln498, Asn501, Tyr505 are replaced with Lys440, Ser446, Asn477, Lys478, Ala484, Lys493, Ser496, Arg498, Tyr501 and His505, respectively. The result analysis shows that the protodioscin formed the strongest binding affinity with 6 hydrogen bonding with the spike protein. The highest rank for the protodioscin shows that all residues involved were reported to interact with the ACE2, implying that protodioscin would be a good inhibitor for the spike protein in preventing the merging of spike protein with the ACE2.
3.1.5. Human TNF–Alpha molecular docking
There are sixteen residues involved in the inhibition of the TNF–alpha, they are six tyrosine residues: Leu57, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121, Gly122 and Tyr151 From Chain A And Leu57, Tyr59, Ser60, Tyr119, Leu120, Gly121 and Tyr151 from chain B [141]. The most crucial residue is the Tyr119 which plays a crucial role in the inhibition process of the TNF–alpha [142]. It can be seen that protodioscin and many PLE compounds, which incidentally form the top–ranking compounds, interacted with most of the residues, with protodioscin showing the best free binding of energy (i.e., –13.64 kcal/mol) (Table 6) with a more comprehensive interactions were visualized using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figure 5, manghaslin and glycyrrhizic acid, Figure S9 and S10, respectively). In terms of H-bond formations of the main compound protodioscin to TNF-apha, Ser60, Gln61, Leu120, Tyr151 binds to Chain A of TNF-alpha at the alkyl glucoside or more specifically at the o-glucopyranosyl moiety of protodioscin. At the at the 6-deoxy mannopyranosyl moiety of protodioscin, several amino acid residues from other chains form H bonds with this protodioscin moiety. This include Val123 from chain B, Ser60 and Leu120 from chain C and Ile58, Ser60, Tyr119, Leu120 and Gly122 from chain D (Figure 5(d)).
Manghaslin, the second-best compound, binds to TNF-alpha and form the same interaction with protodioscin and glycyrrhizic acid (third-best compound) at amino acid residue Tyr119. On the other hand, manghaslin and glycyrrhizic acid did not interacted with chain A residues (Table 6). It is interesting that glycyrrhizin was found to inhibit TNF–alpha–induced apoptosis in the human hepatoblastoma line HepG2 [143] and reduces the level of TNF–alpha in methotrexate–induced enteritis in rats [144], but its role in directly inhibiting TNF–alpha needs further studies in light of the result in this work. In addition, as TNF–alpha levels are routinely assayed using ELISA and glycyrrhizic acid is a large molecule with a molecular weight of 822.9, the possibility of steric hindrance via glycyrrhizic acid interfering with antibody binding to TNF–alpha needs to be explored. Compounds from PLE that interacts with the key Tyr119 residues with values from –9.0 to –10.67 kcal/mol (good binding) are manghaslin, clitorin, kaempferol–3–(2G–glucosylrutinoside), cynaroside (luteolin–7–glucoside), acacic acid, lupeol, betulinic acid, dehydrocarpaine, isoquercetrin, nicotiflorin, benzyl glucosinolate, carpaine, vitexin (in order of decreasing strength of binding) (Tables 6 and S1). Another TNF–Alpha–inhibiting PLE compound; rutin, with a free energy of binding of –8.96 kcal/mol (Table S1) has also been found via a docking study to inhibit TNF–alpha where this compound was found to interact with the residues Gln61, Tyr119, Gly121 and Tyr59 [145], of which two residues (Tyr119 and Tyr59) interaction was seen in this study.
There have been calls for the use of anti–TNF therapy as a possible treatment for COVID–19 [68],[69]. The possible use of PLE in alleviating cytokine storm has been demonstrated in dengue infection in mice models [70]. The role of PLE as an inhibitor to TNF–alpha may play a bigger role than its ability to inhibit IL–6 as it is more upstream than IL–6 in the inflammation pathway cascade [72]. The current biomarkers for predicting the severity of COVID–19 are D–dimer and C–reactive protein [146] and have indicated that pulmonary embolism (PE) is a major cause of death [147]–[151] with thrombocytopenic events a major contributing factor [73],[152],[153]. More recent works indicated that the biomarkers C–X–C motif chemokine ligand 13 (CXCL13) and hepatocyte growth factor (HGF) as better indicators to predict COVID–19 severity and mortality [154]. Despite this, the search for TNF–alpha inhibitors in COVID–19 setting is still ongoing [145],[155],[156]. A recent study that demonstrates the first role of FcγR–mediated monocyte infection, an example of antibody–mediated enhancement (ADE) in SARS–CoV–2 found increased levels of cytokines including TNF–alpha that can potentiate the formation of cytokine storm [157]. This event may provide a bigger role for TNF–alpha inhibitors including PLE compounds in the future.
3.1.6. Human alpha thrombin molecular docking
Docking results show that all of the top compounds inhibiting alpha thrombin were from PLE; protodioscin, manghaslin, apixaban, kaempferol–3–(2G–glucosylrutinoside), cynaroside (luteolin–7–glucoside), vitexin, clitorin and rutin (descending order of binding strength) and they all interacted with the same residues of PPACK, with protodioscin being the strongest inhibitor, having a calculated inhibition constant of 27.07 pM (Table 7). Other interacting compounds are shown in Table S2. More comprehensive interactions were visualized using using the ribbon diagram and three–dimensional representation and also pymol and ligplot+ program (protodioscin, Figure 6, manghaslin and apixaban, Figures S11 and S12, respectively). In terms of H-bond formations of variable strength, the main compound protodioscin forms H bond to the chain H of alpha thrombin with Glu97A, Arg175 and Arg175 at the alkyl glucoside or more specifically at the o-glucopyranosyl moiety while it forms H bonds with Tyr60A, Glu97A, Arg97, Trp96, Leu99 at the tetrahydrofuran moiety. H bonds with Glu39, Leu40, Leu41, Asn143, Thr147, Gly193 and Asp194 were also formed at the 6-deoxy mannopyranosyl moiety of protodioscin (Figure 6(d)).
The catalytic triad of thrombin consists of Ser195, His57, and Asp102 [158] forming a cleft. Four extending polypeptide segments bordered the cleft and gave shape to the cleft. The segments Arg173, Ile174, Asp95, Leu99, Tyr60A, Phe60H and Phe34–Cys42 form on one side and on the opposite side there are two loops which are Gly216–Cys220 and Asn143–Gln151. Protruding into the active–site cleft are the side chains of Ile174, Tyr60A, Trp60D, Lys60F, Phe60H, and Thr147, Trp148. The letters preceding the amino acid numbers (e.g., Tyr60A, Trp60D etc) denote the “chymotryspsin numbering system” (ChyNS) first initiated by Hartley in the 1970s [159]. These segments resulted in the rim to be lined primarily by hydrophobic groups. On the other hand, polar or charged side chains consisting of Glu217, Glu192, Asn143, Gin151, Arg73 form the base [158]. One of the most potent direct thrombin inhibitors is PPACK. Thrombin–PPACK crystal structure complex shows that PPACK interacts with Leu41, His57, Trp60d, Asn143, Gln151, Glu192, Gly193, Ser195, Ser214, Trp215, Gly216 [87].
Manghaslin, the second-best compound binds to alpha thrombin and form the same interaction with protodioscin and apixaban (third-best compound) at the amino acid residues His57, Tyr60A, Trp60D, Lys60F, Glu97A, Glu192, Gly193, Trp215 and Gly216 (Table 7). The direct thrombin inhibitors pachydictyol A, isopachydictyol A, dichotomanol, argatroban and dabigatran all bind to the catalytic side residues either with hydrogen bonding or Van der Waals or hydrophobic interactions [87]. Glycyrrhizic acid has long been known for its thrombin inhibitory activity [160],[161]. A docking study has shown its thrombin binding activity [162] but no residues were identified. In this study, we identify the key residues involved in this compound's interaction with thrombin which shows that it interacts with key residues similar to that of other potent thrombin inhibitors interact (Table 7). In another study, obtusifolin and aurantioobtusin could strongly interact with His–57 and Gly–216 via hydrogen bonding, and with Lys–60 via attractive charging, as well as with Trp–60 via π–π T–shaped interactions and with Cys–220 via alkyl interactions [163]. Thus, as far as PLE compound's ability to inhibit thrombin by binding to key residues that have been demonstrated in other potent thrombin inhibitors suggests the potential use of PLE to treat not only thrombocytopenia in normal COVID–19 progression but also thrombocytopenia was seen in a number of vaccination adverse effects cases [18],[164]–[173].
3.2. Binding-free energy calculations
The binding-free energy calculation was carried out in an effort to learn more about the ligand-protein interaction, the relative affinities of the binding sites, and the roles played by each ligand and residue. Here, the MM/PBSA method [90] was used to estimate the binding free energy, where a larger negative binding energy implies a stronger affinity of the ligand for its specific target pocket. The MM/PBSA is thought to be as precise as the free-energy perturbation techniques but at a significantly lower computational cost [176]. The total binding-free energies and individual energy terms (DGtotal binding ± SD) for top three compounds binding to corresponding SARS-CoV-2 targets (Tables S1–S6) indicates that protodioscin shows the strongest binding to all protein targets with the exception of TNF-alpha, of which glycyrrhizic acid shows a slightly greater binding energy. The results for the MPro protein for protodioscin is comparable to a published result [177] with total binding-free energies ranging from −58.63 to −170.20 kJ/mol for the best marine natural polyketides candidate compounds. As expected, the contribution of the van der Waal energy in all cases is predominant followed by electrostatic energy and solvent-accessible surface area (SASA) energy as reported in many cases [110],[112],[176]–[178].
3.3. Molecular dynamics simulation
The root mean squared deviation (RMSD) is one of the most frequent methods for verifying the validity of docking findings (RMSD). Following receptor superimposition, RMSD measures the average distance between model–ligand atoms and ligand–reference atoms. Near–native" solutions are defined as ligand positions within 2 RMSD heavy–atom distances of the crystallographic pose. 5 While RMSD is widely accepted, it does have flaws that cause postures to be incorrectly classified as either proper or incorrect [179]–[181]. The molecular dynamics simulation (MDS) is performed by using Gromacs. The best substance is used to confirm its binding affinities. The best substance throughout the Molecular Docking using Autodock 4.0 is protodioscin. Although a ligand pose at ≤2 Å RMSD heavy–atom distance with respect to the crystallographic pose is accepted as a “near–native” solution [180], the generally acceptable norm is anything less than 3.00Å [182]–[185], while a cutoff point is > than 4.00Å [186]. Despite this, RMSD is not without its pitfalls, with RMSD still may not be indicative of how well the critical interactions are conserved. For instance, large RMSD values can come from a compound displaying symmetry elements totally flipped over or from flexible atoms not interacting with the protein, like the solvent–exposed parts of the ligand [180].
3.3.1. Nucleocapsid (N) protein
The nucleocapsid (N) protein of coronaviruses aggregates viral genomic RNA into a ribonucleoprotein that shields the viral genome from nucleases in the host. The N protein is made up of two domains: the N–terminal domain (NTD), which is responsible for RNA binding, and the C–terminal domain (CTD), which is also responsible for RNA binding, protein dimerization, and nucleocapsid stability. The structure evaluated in this study is the tetrameric structure consisting of four monomers stabilized with Zn2+ cations that have been evaluated in silico for the potential antiviral binding site [42],[42],[81],[82],[101],[187]–[191]. The backbones RMSD of both apoenzyme and ligand complex fluctuated between 0.1 to 0.5 nm and 0.15 to 0.6 nm, respectively during the first 50 ns of MD simulation. The backbone of apoenzyme stabilizes between 0.5 and 0.55 nm onwards with means and standard deviations of 0.55 ± 0.04 while the backbone of the ligand complex stabilizes between 0.5 and 0.7 nm onwards at 70 ns and above with means and standard deviations of 0.65 ± 0.04 nm. The largest difference between the apoenzyme and its complex occurred at approximately 16.45 ns with a difference of approximately 0.15 nm (Figure 7(a)). The RMSD value for the MDS is quite high (~0.7 nm) and one possibility is that the source of the protein used in this study (6VYO) is a tetramer, the same crystal structure utilized in several SARS–CoV–2 druggability studies including curcumin [82] with an RMSD value of ~0.4 nm, two antiviral moieties from the Asinex databases (5817 and 6799) and Zidovudine with RMSD values ranging from 0.30 to 0.48 nm [191] and the docking of the natural substrate guanosine monophosphate (GMP) to the tetrameric protein [42] with an RMSD value of ~0.4 nm. Docking analysis of the tetramer ryanodine receptor (RyR) shows that at equilibrium, the tetrameric structure has an RMSD value of about 0.4 nm, which is 0.1 nm greater than the monomeric form [192]. It is probable that the tetrameric protein shows a higher RMSD value than the monomeric form, as far as the RNA binding region of the NTD protein is concerned. For example, in one study, the RMSD equilibrium value for HMG–CoA reductase differs by about 1 nm for dimer and tetramer [193] and is also seen in the MDS of tetrameric protofilament and octameric filament of tau R3–R4 [194]. On the other hand, a docking study on the monomeric NTD (6M3M) shows that glycyrrhizin (or glycyrrhizic acid or glycyrrhizic acid), which exhibited a strong inhibitory potential to the protein through in vitro and docking studies show that the NTD–glycyrrhizin complex exhibits high RMSD values (0.4 nm) similar to our findings in this study. The difference in the RMSD profile of monomeric and tetrameric structures upon the same ligand binding, therefore, deserves its own study.
The RMSF method estimates the average change of each amino acid from its reference frame in a protein over time, providing a valuable assessment of the target residues' nonlinear response, which can be thought of as fluctuation and flexibility. There is a dearth of universally accepted cut-off value for indicating support to ligand binding in RMSF trajectories. We adopted the approach of [177] using the difference of root-mean-square fluctuation (ΔRMSF = apoRMSF−holoRMSF) cut-off value of 0.30 Å for main and vicinal interacting residues for a better estimation of the protein local flexibility with amino acids with ΔRMSF above this cutoff value are considered of limited mobility or strong binding to ligand [177]. The RMSF plot for SARS-CoV-2 nucleocapsid protein (Figure 7(b)) shows fluctuations of residues that incidentally occur at predominantly loop regions, which is a characteristic of this region [195] whilst alpha–helix (residues 79–82) and the majority of beta sheets (84–90, 107–112 and 130–134) are more stable and fluctuate less. The root–mean–square–fluctuation (RMSF), which reflects the magnitude of thermal motions of amino acid residues, of the RNA binding domain of nucleocapsid complexes is between 0.1 and 0.9 nm. Residues with ΔRMSF above or approximating the cut-off value (between 0.025 and 0.030Å) for chain C were Pro73, Ile74, Thr76, Gln83 and Thr13 while for chain B they were Leu167, Pro168, Lys169, Gly170 and Phe171 (H bond in bold), which corresponded well to the residues implicated for protodioscin binding region to N protein.
It appears that the binding of a ligand causes an increase in fluctuations of the nucleocapsid chains. Interestingly, this increase in residues' fluctuation upon binding of the ligand in chains is also seen in the binding of curcumin to SARS–CoV–2 nucleocapsid [82]. The observed increase in RMSF is probably due to the partial collapse of the secondary structures upon ligand binding while the overall effect of ligand binding was a decrease in RMSF. It is probable that ligand binding–induced changes in the structural fluctuation of the residues at the ligand binding sites which are then transmitted to and amplified at a distant region translating into an increase in fluctuation compared to the apoenzyme [196] or the high fluctuations of the RMSF might be due to the relatively volatile properties of tetramers compared to monomer, an avenue that needs to be studied as far as docking studies on SARS–CoV–2 NTD (tetramer) is concerned. Despite this issue, the knowledge of the dynamical aspect of proteins thus provides important information to understand the detailed mechanism of the interaction between protein and ligand binding.
3.3.2. RdRp
For the protodioscin–RdRp complex, the backbones RMSD of both apoenzyme and ligand complex fluctuated between 0.15 to 0.25 nm during the first 2 ns of MD simulation. The backbone of apoenzyme stabilizes between 0.16 and 0.22 nm onwards with means and standard deviations of 0.19 ± 0.03 while the backbone of the ligand–complex stabilizes between 0.19 and 0.25 nm at 50 ns onwards with means and standard deviations of 0.20 ± 0.04 nm. The largest difference between the RdRp apoenzyme and its complex occurred at approximately 70 ns with a difference of approximately 0.08 nm (Figure 8).
As before, the RMSF plots of fluctuating regions occurs at predominantly loop regions whilst alpha–helix regions (residues 170–200, 234–248, 303–319, 465–479, 627–640, 686–710) are more stable and fluctuate less. the RMSF values of the RdRp complexes bound to remdesivir exhibit high peaks at the residues Lys91, Thr262, Glu431, Thr853 and Asp893, respectively, which are located outside of the enzyme binding site [189]. In our study, we observed that residues exhibiting high fluctuations (> 0.2 nm) were 220, 228, 262, 302, 431, 496, 823, 853 and 893 of which the residues Glu431, Thr853 exhibit similar fluctuations to remdesivir. A trio of aspartates, Asp618, Asp760, and Asp761, coordinate the catalytic metal ions in SARS–CoV–2 nsp–12 (part of apo RdRp), and Arg555 stabilizes the substrate phosphate, while the residues Asp623, Ser682, and Asn691 are implicated in the 2′–OH interaction with the incoming nucleotide [197]. In our study we observed three residues Asn496 Asn497 and Lys577, exhibiting a stabilization of fluctuation (ΔRMSF) of near or more than 0.3 Å (Figure 8(b)), which was implicated in protodioscin binding to RdRp mentioned previously. These three residues are important as they are the same residues that pralatrexate, one of the 1906 FDA–approved drugs screened for SARS–CoV–2 drugs, binds with a strong inhibitory activity (EC50 value of 8 nM) [37].
3.3.3. MPro/3CLpro
For protodioscin–MPro/3CLpro complex, the backbones RMSD fluctuate between 0.1 to 0.2 nm throughout the MD simulation period with means and standard deviations of 0.22 ± 0.05 and 0.15 ± 0.04 nm from the mean for MPro/3CLpro apoenzyme and MPro/3CLpro–protodioscin complex, respectively, during the last 50 ns. The largest difference between the MPro/3CLpro apoenzyme and its complex occur at approximately 60 ns with a difference of approximately 0.1 nm (Figure 9(a)). The binding of protodioscin to MPro/3CLpro changes the fluctuation dynamics landscape, with some residues depicting a reduction in fluctuation while others reflect an increase in fluctuations, a tell–tale indication that binding has occurred. In the future, statistical support for these fluctuations can be done by repeating the simulations several times.
The highly fluctuating residues seen in this study are also seen in the RMSF of MPro/3CLpro SARS–CoV–2 apoenzyme and with the ligands glycyrrhizin and rhodiolin [189]. The RMSF plot of highly fluctuating regions occurs at predominantly loop regions (especially at residues 71–73, 135–144, 153–155, 167–171, 188–197 and 276–280) while beta–sheets/strands (28–31, 37–41, 86–89, 114–118, 144–150 and 157–166) and alpha–helix regions (residues 202–211, 262–269, 290–296) are more stable and fluctuate less.
The interaction between protodioscin and MPro/3CLpro SARS–CoV–2 increases the fluctuations of many of the residues mentioned above similar to several published results [24],[189],[198]. Ser1, His41, Met49, Gly143, Phe140, Ser144, Cys145, His163, His164, Glu166, Pro168 and Gln189 have been determined in previous studies as lining the active sites of MPro [125] while another study has indicated that the residues Thr24, Thr25, Thr26, Leu27, His41, Met49, Tyr54, Phe140, Leu141, Asn142, Gly143, Ser144, Cys145, His163, His164, Met165, Glu166, Leu167, Pro168, Asp187, Arg188 Gln189, Thr190, Ala191, and Gln192 are a key target for potential inhibitors among which, His41 and Cys145 are two key hydrophobic catalytic residues [80]. In our study we observed the residues Thr24, Glu166, Leu167 and Pro168, exhibiting a stabilization of fluctuation (ΔRMSF) (Figure 9(b)) of near or more than 0.3 Å, which was implicated in protodioscin binding to MPro/3CLpro mentioned previously.
3.3.4. SARS-CoV-2 variants
For the protodioscin complex with the three variants of SARS-CoV-2 spike protein, the backbones RMSD stably fluctuated between 0.1 to 0.25 nm during the end of the MD simulation period (> 60 ns) with a mean and standard deviation (from 60 to 120 ns) for variant Wuhan, Delta, Omicron and apoenzyme of 0.139 ± 0.016 nm, 0.154 ± 0.022 nm, 0.175 ± 0.015 nm and 0.124 ± 0.010 nm, respectively. ANOVA analysis showed that all variants show significantly different (p < 0.05) RMSD values to apoenzyme (Figure 10(a)). The RMSF plot of highly fluctuating regions occurs at predominantly loop regions (especially at residues 340–341, 366–368, 369–371, 475–488) whilst beta sheets (395–403 and 507–515) and alpha–helix regions (residues 349–353,405–410 and 416–422) are more stable and fluctuate less.
The fluctuations were seen at amino acid residues for spike protein RMSF plot centered at the amino acid clusters from residue 338 to 343, 361 to 374, 380 to 387, 475 to 488 and 519 to 277 which are similar to previous findings which reported high variability of the RBD regions centering at 346, 369, 386, 473 to 490 and 519 to 527 [199], where binding of ligand reduces the fluctuation at these clusters especially at clusters 338 to 343, 361 to 374, 475 to 488 while there is a small increase in fluctuation at residue 502 (Figure 10(b)). The residues 475–488 exhibited the most flexibility compared to other regions of the protein. The binding of ACE2 to RBD causes a change >0.2 nm in this region similar to a previous report [200]. The region 438 to 506 represents the receptor–binding motif (RBM) within the RBD [43]. In our study we observed the residues Gly496, Gln498, Thr500, Asn501, Gly502, Val503 and Gly504 exhibiting a stabilization of fluctuation (ΔRMSF) (Figure 10(b)) of near or more than 0.3 Å, which was implicated in protodioscin binding mentioned previously. Incidentally, these residues are part of the spike protein region that are less than 3Å from the ACE2 contributing to firm binding [79].
3.3.5. TNF alpha
TNF alpha suppression is an important therapeutic strategy for combating inflammation in conditions such as rheumatoid arthritis, insulin resistance, and type 2 diabetes. TNF–alpha is a homotrimeric protein with 148 amino acid residues that span the length of each monomer. Only two of the three subunits are involved in the active site's formation. Both subunits contribute 16 amino acid residues to the architectural foundation of the site. Nine of these residues may be found in chain A, whereas the remaining seven can be found in chain B [155],[201]–[203]. These residues are hotspots for efficient TNF inhibition. For protodioscin–TNF–alpha complex, the backbones RMSD fluctuate between 0.1 to 0.4 nm for the first 10 ns of the MD simulation period before stabilizing and fluctuating between 0.45 and 0.48 nm between 50 and 120 ns with means and standard deviations of 0.36 ± 0.05 and 0.49 ± 0.05 nm for TNF–alpha apoenzyme and TNF–alpha–protodioscin complex, respectively. The largest difference between the TNF–alpha apoenzyme and its complex occurs at approximately 112 ns with a difference of approximately 0.24 nm (Figure 11(a)). In our study we observed the residues Tyr59, Ser60, Tyr119, Leu120, Gly121 from chain A exhibiting a stabilization of fluctuation (ΔRMSF) (Figure 11(b)) of near or more than 0.3 Å, which was implicated in protodioscin binding mentioned previously and among the sixteen residues involved in the inhibition of the TNF–alpha [142]. The sixteen residues involved in the inhibition of the TNF–alpha [142], are Leu57, Tyr59, Ser60, Gln61, Tyr119, Leu120, Gly121, Gly122 and Tyr151 from chain A and Leu57, Tyr59, Ser60, Tyr119, Leu120, Gly121 and Tyr151 from chain B [141] exhibited amongst the lowest fluctuations in the RMSF plot for reasons we have discussed previously and was also observed in the RMSF results for the inhibition of TNF–alpha by rutin, hesperidin, and schisantherin A [145]. The RMSF plot of highly fluctuating regions occurs at predominantly loop regions (residues 19–28, 65–75, 84–90, 100–115) whilst beta–sheets/strands (28–29, 54–63, 79–83, 91–97, 118–127, 131–135 and 151–155) and alpha–helix regions (residues 137–142 and 148–152) are more stable and fluctuate less.
3.3.6. Alpha thrombin
For the protodioscin–alpha thrombin complex, the backbones RMSD fluctuation was very minimal and was stable before stabilizing and fluctuating between 0.19 and 0.22 nm with means and standard deviations of 0.209 ± 0.01 and 0.178 ± 0.01 nm for alpha thrombin apoenzyme and protodioscin–alpha thrombin complex, respectively. The largest difference between the alpha thrombin apoenzyme and its complex occurs at approximately 96 ns with a difference of approximately 0.1 nm (Figure 12(a)). The RMSF plot of highly fluctuating regions occurs at predominantly loop regions (especially at residues 21–25, 32–38, 142–154, 190–199 and 213–218) whilst beta–sheets (317–20, 24–29, 40–48, 58–64, 82–87, 137–140, 184–189, 205–210, 224–229 and 237–242) and alpha–helix regions (residues 122–130, 169–176 and 247–255 with the exception of residues 123–128) are more stable and fluctuate less. The RMSF result also shows that the residues constituting the thrombin catalytic triad of His57, Asp102, and Ser195 have relatively low fluctuations, consistent with their required stable positions for thrombin activity [158]. In one study, a general reduction in fluctuations was observed after protodioscin binding at locations similar to the non–sulfated arabinan octasaccharide–thrombin complex [204]. In our study, we observed the residues Tyr60A, Trp60D, Lys60F, Asn143, Gln151, Glu192, Gly193, Trp215, Gly216 from chain H exhibiting a stabilization of fluctuation (ΔRMSF) (Figure 12(b)) of near or more than 0.3 Å, which were implicated in protodioscin binding mentioned previously. The fact that one of thrombin strongest inhibitor, PPACK also binds at many of these residues, [87] indicates a strong candidate for protodioscin as a potentially good inhibitor to alpha thrombin.
The RMSF profile, after considering the ChyNS system, is very similar to the published RMSF profile for alpha thrombin. In our study, we did not include the radius of gyration (Rg), a parameter that is associated with the structural compactness and overall conformational shape of proteins, which was calculated to evaluate the compactness of protein–drug complexes. In one study, the Rg plot did not display remarkable changes which indicated that there was no structural deviation in all three proteins [199]. It is difficult to infer a useful conclusion based on a statistically significant change in the compactness of the binding of small molecules to a large protein, a feature reported in one study where the binding of boceprevir (Boc) and dexamethasone (Dex) drugs to the SARS–CoV–2 protein targets; SARS–CoV–2 main protease (MPro), spike protein C–terminal domain (spike–CTD), and interleukin–6 (IL–6) where the authors concluded that no significant changes were observed and hence no conclusion can be inferred [199]. To assess the change more meaningfully, a repeated simulation of a protein with its ligand should be done many times and the usual robust t–test carried out to test for the existence of a significant change [205]. Future works include redocking the ligands with better resolution crystals (< 2 Å).
Although we show that protodioscin is the best inhibitor for many COVID–19 targets covering the virus and human sides (Figure 13), it does not mean that protodioscin alone can be a “cure–all” strategy for alleviating COVID–19 severity. Despite this, the therapeutic potential of protodioscin in many in vivo and in in vitro experiments has been reported [206]–[211]. The docking studies show numerous PLE compounds at the top of the “food chain” of COVID–19 target inhibitors and coupled with the bioavailability landscape of these numerous compounds, which might be better or worse than the top compound, may mean that taking PLE extract may be better than taking the top compound alone. A total of more than 1000 patients have participated in multiple clinical trials conducted by PLE for dengue fever, with overwhelmingly positive efficacy and safety results. Furthermore, the papaya plant is one of the most widely recognised plants on the planet. Aside from the fruit, the leaves have been used as medicine, food, and drink for hundreds of years, demonstrating the safety of PLE. The leaves are also edible and have been eaten as a routine staple for hundreds of years in several countries in the form of salad, stir fry preparation or tea infusion, chiefly in India, Thailand, Indonesia and Malaysia [212].
Recently, a guideline for the diagnosis, management, and prevention of VITT appearing as cerebral venous sinus thrombosis, which was developed by the Spanish Federation of Medical and Scientific Associations (FACME) (CVST) includes a 50% decrease in baseline platelet count as an alarm for COVID–19 vaccine–induced thrombotic thrombocytopenia (VITT) [63]. As PLE has been clinically shown to return decreasing platelet count to normal in dengue patients, PLE may be of help to combat vaccine–induced complications. The most common way for preparing PLE extract in the traditional setting to treat dengue is to blend one papaya leaves between 20 to 30 cm in length with two cups of water or about half a litre, filtering the juice using a cloth strainer and taking half a cup daily mixed with orange juice or any acidic juices to mask the bitter flavour [213].
PLE may help diabetic and non–diabetic COVID–19 patients with hyperglycemia. PLE can treat the hypercoagulation problem associated with venous and arterial thrombotic issues in severe COVID–19 patients, and PLE may also restore platelet levels in COVID–19 patients who have severe thrombocytopenia due to infection or vaccine complications. The extract of papaya leaves (PLE) is a safe supplemental alternative treatment that has a solid track record in multiple clinical studies involving over 1000 patients. Papaya leaf extract may aid in the battle against COVID–19 infection, but its effectiveness must be evaluated using clinical endpoints such as mortality, time to clinical improvement, and the number of days spent in the intensive care unit.
4.
Conclusions
Major components from PLE and several previously reported SARS–CoV–2 and SARS–CoV inhibitors were screened in silico against four SARS–CoV–2 targets and two human targets that, when combined, may reduce or combat the severity of COVID–19 disease. Docking studies against the virus targets, namely the RNA binding domains of nucleocapsid, MPro, RdRp, and spike protein (Wuhan, Delta, and Omicron variants), as well as the human targets TNF–alpha and thrombin, reveal protodioscin as the top potentially inhibiting agent, with calculated inhibition constants from pM to nM. Other PLE compounds were also in the top rankings for all virus and human targets. The interaction between protodioscin and these targets is confirmed by molecular dynamics simulation using RMSD and RMSF analyses. As a result, PLE may be capable of overcoming the obstacles to successful COVID–19 therapy. These barriers include the inhibition of TNF–alpha and thrombin activities as well as the prevention of viral replication.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
DownLoad:
