Evaluation of the active constituents of Nilavembu Kudineer for viral replication inhibition against SARS-CoV-2: An approach to targeting RNA-dependent RNA polymerase (RdRp).

The World Health Organization has declared the novel coronavirus (COVID-19) outbreak a global pandemic and emerging threat to people in the 21st century. SARS-CoV-2 constitutes RNA-Dependent RNA Polymerase (RdRp) viral proteins, a critical target in the viral replication process. No FDA-approved drug is currently available, and there is a high demand for therapeutic strategies against COVID-19. In search of the anti-COVID-19 compound from traditional medicine, we evaluated the active moieties from Nilavembu Kudineer (NK), a poly-herbal Siddha formulation recommended by AYUSH against COVID-19. We conducted a preliminary docking analysis of 355 phytochemicals (retrieved from PubChem and IMPPAT databases) present in NK against RdRp viral protein (PDB ID: 7B3B) using COVID-19 Docking Server and further with AutoDockTool-1.5.6. MD simulation studies confirmed that Orientin (L1), Vitexin (L2), and Kasuagamycin (L3) revealed better binding activity against RdRp (PDB ID: 7B3B) in comparison with Remdesivir. The study suggests a potential scaffold for developing drug candidates against COVID-19. PRACTICAL APPLICATIONS: Nilavembu Kudineer is a poly-herbal Siddha formulation effective against various diseases like cough, fever, breathing problems, etc. This study shows that different phytoconstituents identified from Nilavembu Kudineer were subjected to in silico and ADME analyses. Out of the former 355 phytochemical molecules, Orientin (L1), Vitexin (L2), and Kasuagamycin (L3) showed better binding activity against RdRp viral protein (PDB ID: 7B3B) in comparison with the synthetic repurposed drug. Our work explores the search for an anti-COVID-19 compound from traditional medicine like Nilavembu Kudineer, which can be a potential scaffold for developing drug candidates against COVID-19.

INTRODUCTION

At the end of 2019, several cases of pneumonia were reported in Wuhan, Hubei province, with an idiopathic cause. Later, the screening of samples obtained from the lower respiratory tract of patients confirmed the causative agent as severe acute respiratory syndrome coronavirus‐2 or SARS‐CoV‐2. About 153,187,889 confirmed cases of COVID‐19, including 3,209,109 deaths, have been reported by the WHO (Di Gennaro et al., 2020; https://covid19.who.int/). SARS‐CoV‐2 belongs to a large class of RNA virus family, Coronaviridae, and constitutes an enveloped structure with a single‐stranded RNA genome. Among the four different genera (α, β, γ, and δ) of viruses, SARS‐CoV‐2 belongs to the β‐class of coronaviruses. SARS‐CoV‐2 is revealed to have more sensitivity toward lung angiotensin‐converting enzyme (ACE‐2) receptors and thus infects the lower respiratory tract, producing significant infections in the bronchioles and alveoli (Kumar et al., 2021). Coronavirus disease 2019 or COVID‐19 is transmitted from an infected person via droplets produced from talking, sneezing, or coughing (Baby et al., 2021; Singhal, 2020). The SARS‐CoV‐2 enters the host cells via binding to the human ACE‐2 receptor. The viral Spike protein binds with the ACE‐2 receptor and forms a viral envelope that gets fused and allows the envelope to enter the host cells through endocytosis and autophagic pathways (Adithya et al., 2021). The viral replication inside the host cell is aided by RNA‐dependent RNA polymerase (RdRp) and open reading fragments (Baby et al., 2021; Loeffelholz & Tang, 2020). The viral RNA is released into the host cells and undergoes a translation process to form various structural and non‐structural proteins. The viral progeny thus produced will worsen the immune responses and promote the infection to a more aggravated stage (Adithya et al., 2021; Baby et al., 2021) (Figure S1).

Herbal medicines are well‐known for their pharmacological activities; their uses are still in practice in different parts of the globe (Adhikari & Paul, 2018). Traditional herbal medicines symbolize the ancient form of the health‐care system, which helps to prevent and treat various ailments (Adhikari & Paul, 2018; Yuan et al., 2016). Siddha medicine system is implied as the oldest means of medicinal practice in South India (Adhikari & Paul, 2018). Nilavembu Kudineer is a polyherbal formulation used in the Siddha system of medicine against multiple viral infections (Christian et al., 2015). Nilavembu Kudineer formulation contains an equivalent amount of nine ingredients: Nilavembu (Andrographis paniculata), Vilamichai ver (Plectranthus vettiveroides), Vetiver (Vetiveria zizanioides), Santanam (Santalum album), Korai kizhangu (Cyperus rotundus), Parpatakam (Mollugo cerviana), Chukku (Zingiber officinale), Pei Pudal (Trichosanthes dioica), and Milagu (Piper nigrum) (Chitra et al., 2021). Nilavembu Kudineer was found to have immunomodulatory activities, thereby boosting the immune response against viral infections. The ingredients present in Nilavembu Kudineer possess antipyretic, antiviral, cardioprotective, and hepatoprotective activities (Kamalarajan et al., 2019). A recent study by Jain et al. (2020) reports Nilavembu Kudineer as an effective antiviral agent against Dengue and Chikungunya. The present study aims to identify potent bioactive as a viral replication inhibitor from NK against SARS‐CoV‐2 via computational investigations and further validate with an in vitro study.

MATERIAL AND METHODS

Natural product

Nilavembu Kudineer formulation constitutes primarily nine ingredients, namely, Nilavembu (Andrographis paniculata), Vetiver (Vetiveria zizanioides), Santanam (Santalum album), Parpatakam (Mollugo cerviana), Chukku (Zingiber Officinale), Pei Pudal (Trichosanthes dioica), Milagu (Piper nigrum), Vilamichai ver (Plectranthus vettiveroides), Korai kizhangu (Cyperus rotundus) (Chitra et al., 2021). Around 355 phytoconstituents from different ingredient herbal plants have been identified in Nilavembu Kudineer formulation with reported antiviral activity against Dengue and Chikungunya. The structures of all the 355 phytoconstituents present in nine different NK herbal plants were retrieved from the IMPPAT (https://cb.imsc.res.in/imppat/home) and the PubChem (https://pubchem.ncbi.nlm.nih.gov/) databases. The affinity of phytoconstituents from each ingredient of NK formulation with the standard Remdesivir was investigated using the COVID‐19 database (https://ncov.schanglab.org.cn/) and further confirmation with AutoDockTool‐1.5.6. The active phytoconstituents present in each herbal plant are summarized in Table S1.

Protein identification and preparation

The three‐dimensional structures of an RdRp (PDB ID: 7B3B) viral proteins have been identified by the Research Collaboratory of Structural Bioinformatics Protein Data Bank (RCSB PDB), www.rcsb.org.

Protein characterization

The physicochemical properties of the protein were analyzed using the PROTPARAM tool. The parameters include molecular weight, theoretical Pi, instability index, half‐life (hours), aliphatic index, and gravy value (Table S2).DThe secondary structure characterization was predicted using the SOPMA tool.

Ligand identification and preparation

Three hundred and fifty‐five ligand molecules were identified from the Nilavembu Kudineer formulation. The structures of all ligand molecules were downloaded from the PubChem database and saved in the SDF format. Of the 355 phytochemical molecules, Orientin (L1, PubChem CID: 5281675), Vitexin (L2, PubChem CID: 5280441), Kasuagamycin (L3, PubChem CID: 65174) showed significant activity in comparison with Remdesivir against the viral proteins, RdRp (Kamalarajan et al., 2019).

Computational methodology

Retrieval and preparation of protein and ligands and molecular docking

SAR‐CoV‐2 RNA‐dependent RNA polymerase (RdRp) X‐ray crystal structure was obtained from the RCSB Protein Data Bank (accession code: 7B3B) (Kokic et al., 2021). The 3D structure consists of an RdRp monomer, cofactors [a non‐structural protein 7 (nsp7), and a non‐structural protein 8 (nsp8)] structure, and two atoms of zinc metal each coordinated to histamine (HIS) and three cysteines (CYS) residues (Ahmad et al., 2020). The study of RdRp was carried out without removing the cofactors, zinc atoms, nsp7, and nsp8 since they play a crucial biological role in the stability of RdRp protein (Ahmad et al., 2020). The average C‐α atom distance between the active site (SER759, ASP760, ASP 761) of the RdRp protein and proximal three residues of the nsp7 (VAL33, ASN37, LEU40) and nsp8 (ASP78, ASN118, ALA163) were measured using the Chimera‐1.14. Therefore, the average distance from the RdRp binding cavity to cofactors nsp7 and nsp8 was 26.05 and 37.21 Å, respectively. This finding implies that Nsp8 is closer to the RdRp active site than Nsp7; thus, it may impact the protein more stable. In addition, the distance between the active site and each zinc atom was 28.28 and 31.59 Å, respectively. With this established, we can be sure the impact of the proposed inhibitors is accounted for or neutralize the stability impact of the cofactors to the studied protein. The ligands were drawn in the MarvinSketch‐17.21 (http://www.chemaxon.com) (WHO, n.d.), and their geometries were optimized with the General Amber Force Field (GAFF) in Avogadro‐1.2.0 (Hanwell et al., 2012). The ligands (L1, L2, and L3) were the three of the promising compounds we identified in our laboratory against SARS‐CoV‐2 proteins. The receptor and inhibitors were prepared using the Chimera‐1.14 (Pettersen et al., 2004) and the AutoDockTool‐1.5.6 (Goodsell & Olson, 1990). Moreover, the systems were opened and saved on the Molegro Molecular Viewer‐2.5 (Thomsen & Christensen, 2006) for the molecular dynamics simulation study. The two zinc metals coordinated to two CYS and one HIS residues were left unstripped from the RdRp structure during docking since the metals confer stability to the protein. Preparing the unbound and ligands‐bound complexes involves removing nonstandard molecules, including water, NAP, and sodium (Adewumi, Elrashedy, et al., 2020). Figure 1 provides the structures of the RdRp conformation showing some active site residues (A), and the chemical structures of inhibitors L1, L2, and L3. The active site of the SARS‐CoV‐2 RNA‐dependent RNA polymerase was obtained from the binding of Remdesivir by zoning and its surrounding and selected all atoms/bonds and residues within the region of <5.0 angstroms. Additionally, residues ASP760 and 761 were treated as flexible during docking. Docking calculations were carried out with L1, L2, and L3 on RdRp active site. Each compound is assessed in Molegro Molecular Viewer‐7.0. (MMV) (http://molexus.io/molegro‐molecular‐viewer) to ascertain the correct bond angles and hybridization state. Ligands were minimized using the steepest descent method and GAFF force field in Avogadro‐1.2.0 (http://avogadro.cc/) (Hanwell et al., 2012). For this study, we used the holoenzyme of RdRp, which consisted of its cofactors nsp7 and nsp8 since reports showed that the latter enzymes contributed to the RdRp's replication/transcription functions. Each compound was separately docked into the active site of each ACE2 protein. Grid box generation for docking was specified using AutodockTool‐1.5.6 software (Hanwell et al., 2012; Trott & Olson, 2010) and defined the center as X = 94.919, Y = 92.113, and Z = 97.261 with dimensions X = 54, Y = 60, and Z = 52. The grid box was centered around the critical residues binding Remdesivir binding up to <5 angstroms, corresponding to its specific binding site consisting of TYR619, CYS622, ASN695, MET756, ILE757, LEU758, SER759, ASP760, ASP761, ALA762, VAL763, GLU811, CYS813, and SER814 (Adewumi, Ajadi, et al., 2020). The structural and visual analyses of the RdRp and the ligands‐bound systems were carried out using the Discovery studio Visualizer Software‐4.0 (Ahmad et al., 2020) (http: www.accelrys.com) and Chimera Tool.

FIGURE 1

Some active site residues of RdRp protein (a). The molecular docking interaction networks of L1‐bound RdRp (d), L2‐bound RdRp (b), and L3‐bound RdRp (c).

Zinc Parameterization and Molecular dynamics (MD) simulations

The structure files were converted to Amber format to be prepared for parameterization by removing all hydrogen atoms before MD simulations. The zinc coordination site was modeled with the parameters generated using the MCPB package (Kumar et al., 2021). The chosen bonded model allows the reflection of the partial charges on every atom within the coordination complex and prevents ligand substitution during simulations. The bond and angle parameters for the zinc ions were generated using the empirical method during the MCPB.py procedure. We further investigated the proposed inhibitors of RdRp protein compared to the unbound system by evaluating the complexes' dynamic characteristics using the all‐atom molecular dynamics (MD) simulation for 220 ns. The MD simulation was carried out in Amber using the GPU version of the PMEMD engine provided by the AMBER18 package at the Lengau CHPC server (http://www.ambermd.org) (Silva et al., 2014). The FF18SB force field variant was applied to describe the inhibitors‐RdRp protein complexes (Adewumi, Ajadi, et al., 2020; Case et al., 2005; Lee et al., 2018). The antechamber obtained atomic partial charges for the ligands using the restrained electrostatic potential (RESP) and the General Amber Force Field (GAFF) procedures (Zheng et al., 2016). The Leap module of Amber 18 enabled the addition of hydrogen atoms, sodium (Na+), or chloride (Cl−) counter ions to the systems for neutralization (Oluyemi et al., 2022). The complexes, implicitly, were suspended within an orthorhombic box of TIP3P water molecules to contain all the atoms within 8 Å of any box edges (Case et al., 2005; Florova et al., 2010). The systems were initially minimized for 2,500 steps with 500 kcal/mol Å (WHO, n.d.) restraint potential, and a whole minimization step of 5,000 steps was further run without restraint using the conjugate algorithm. The heating step was done in a canonical ensemble condition (NVT) by a gradual heating simulation (0–300 K) executed for 5 ps, while several atoms and volume were fixed. The systems‐containing solutes were imposed with a potential harmonic restraint of 10 kcal/mol Å (WHO, n.d.) and a collision frequency of 1.0 ps−1. An equilibration estimation of 1 ns was carried out while keeping the operating temperature (300 K) constant.

Moreover, the number of atoms and pressure were constant, depicting an isobaric‐isothermal ensemble (NPT). Using the Berendsen barostat, the complex pressure was kept at 1 bar (Kufareva & Abagyan, 2012; Ryckaert et al., 1977). We used the SHAKE algorithm in constructing the hydrogen atoms bonds. Coordinates and trajectories were printed and analyzed every 1 ps using the PTRAJ module available in AMBER18GPU. The stabilities and flexibilities of the free protein and ligands‐bound complexes were investigated statistically by estimating the averages of RMSD (Silva et al., 2014) and RMSF (Kufareva & Abagyan, 2012), respectively. The detailed methods for obtaining RMSD (Dong et al., 2018), RMSF (Lobanov et al., 2008), and RoG have been reported severally in the works of literature, including our previous results.

Thermodynamic calculations and per‐residue energy decomposition analysis

The binding affinities of the bound and unbound systems were obtained by computing the binding free energies using the Molecular Mechanics Generalized‐Born Surface Area method (MMGBSA) (Gapsys et al., 2015). The free energy was calculated based on the average of 40,000 snapshots extracted from a 40 ns trajectory. The estimated free binding energy, ΔG, for each molecular system, including the complex, inhibitor, and protein, can be given as (Hayes & Archontis, 2012):

Moreover, the details of the method can be found in our previous work (Silva et al., 2014).

Furthermore, individual residue energy contributions to the total binding free energy in the complexes L1‐RdRp, L2‐RdRp, and L3‐RdRp were obtained, carrying out the per‐residue energy decomposition analysis using the MMGB(PB)SA method in AMBER14 (Aftab et al., 2020; Kollman et al., 2000).

ADME analysis

ADME prediction was performed using Swiss ADME predictor software. ADME properties like number of heavy atoms, number of heavy aromatic atoms, fraction Csp3, number of rotatable bonds, number of H‐bond acceptors, number of H‐bond donors, and solubility were analyzed. All the molecules possessed ADME properties.

Pseudo virion assay for SARS‐CoV‐2

The assay is based on the lentiviral backbone expressing ZsGreen as a traceable marker. We have utilized stable colon cancer cell HEK293T expressing human ACE‐2 as the SARS permissive cells. The procedure involves transfection of HEK Lenti Cells (Invitrogen) with the expression vector encoding ZsGreen, a plasmid expressing Spike, and plasmids expressing the minimal set of lentiviral proteins necessary to assemble viral particles (Gag/Pol, Rev). The cells were transfected with the expression vectors prepared via Qiagen Midi prep using lipofectamine 2000 as per the manufacturer's instruction. After 6 hr, the cells were replaced with fresh medium containing serum. From the transfected cells, SARS‐CoV‐2‐ Spike‐pseudotyped lentiviral particles were collected at 72 hr, filtered using a 0.45‐micron filter, and infected the HEK293T– hACE2 cells using polybrene as per the standard protocol. The test samples were incubated with pseudovirions containing medium at indicated dilutions. The media diluted pseudovirion sample acts as the control. After 48 hr, the cells were imaged under a fluorescent microscope, and cells expressing ZsGreen fluorescence were counted. Percentage positivity was calculated based on the total number of cells in the field.

RESULTS AND DISCUSSION

Molecular docking

The three hundred and fifty‐five phytoconstituents from the different ingredients of Nilavembu Kudineer formulation with antiviral activity were selected for the study. The phytoconstituents of each ingredient have been retrieved from the two databases IMPPAT (https://cb.imsc.res.in/imppat/home) and PubChem (https://pubchem.ncbi.nlm.nih.gov/). The affinity of these 355 phytoconstituents against Remdesivir (standard) was analyzed using an online docking tool, the COVID‐19 database (https://ncov.schanglab.org.cn/) and further with AutoDockTool‐1.5.6 for the three main compounds. Among these 355 compounds, only three phytoconstituents, Orientin (L1), Vitexin (L2), and Kasuagamycin (L3) were found to have a higher affinity against the standard compound. Vitexin (L2) was found to have a higher docking affinity (−6.84, Schrodinger Database, and −7.90, COVID‐19 Docking Server) toward the RdRp viral protein (PDB ID: 7B3B) in comparison with the positive control, Remdesivir (−6.09, Schrodinger Database, and −7.40, COVID‐19 Docking Server), followed by Orientin and Kasuagamycin (Tables S3 and S4). We have selected L1, L2, and L3 with higher affinity for further confirmatory studies (Tables 1, 2, 3).

TABLE 1

ADME data of selected phytochemical molecules present in Nilavembu Kudineer formulation

Compounds	Molecular weight (g/mol)	No. of heavy atoms	No. of aromatic heavy atoms	Fraction Csp3	Rotatable bonds	HBA	HBD	Solubility
Vitexin	432.38	31	16	0.29	3	10	7	Soluble
Kasuagamycin	379.36	26	0	0.86	4	11	8	Highly soluble
Remdesivir	602.58	42	15	0.48	14	12	4	Moderately soluble

TABLE 2

The molecular docking value of the most active phytoconstituents present in Nilavembu Kudineer against RdRp

Compound	Docking score	Software
Orientin	−7.70
Vitexin	−7.90	COVID‐19 Docking Server
Kasuagamycin	−7.6
Remdesivir	−7.40

TABLE 3

Molecular docking results of the proposed inhibitors of SARS‐CoV‐2 RdRp AutoDockTool‐1.5.6

Ligand name	Estimated binding energy (kcal/mol)	Interacting residues
L1	−8.3	ASN497, SER501, ASN543, VAL557, GLY559, ARG569
L2	−8.1	ASP452, ARG553, ARG555, THR556, ASP618, LYS621, ARG624, ASP760
L3	−7.9	LYS500, GLY559, ARG569, LEU576, GLY683, ALA685

The critical active site residues of RNA‐dependent RNA polymerase (RdRp) are involved in the viral replication and transcription of the SARS‐CoV‐2 genome. The molecular docking technique was used to predict the binding mode of L1, L2, and L3 compounds docked in the active site, including ASP760 and ASP761 and other residues within the RDRp pocket. The best docking calculations from L1, L2, and L3 gave the lowest binding energy in complex with SARSCOV‐2 RNA‐dependent RNA polymerase enzyme. Table 4 shows the binding energies (kcal/mol) of L1 (−8.3), L2 (−8.1), and L3 (−7.9) docked on the RdRp co‐crystallized with its cofactors (nsp7 and nsp8). Compared with the reported binding affinity of the Remdesivir (−14.06) (Singhal, 2020) to RdRp, these energy values are lower. From the chosen docking poses of L1 to L3, the docked poses of L1 and L2 were close to the loop covering the active site, while the L2 pose showed proximity with the active site residues ASP760. The docking results of L1 and L3 displaced no direct interaction with the aspartic residues (ASP760 and ASP761), although the ligands interacted with other residues involved in the RdRp interaction. In addition, the docking of L1, L2, and L3 did not show any evidence of coordination with the zinc in the protein. In L1‐bound RdRp, ASN497, SER501, ASN543, GLY559, and ARG569, strong H bonding with the ligand, while VAL557 formed pi‐alkyl hydrophobic interaction. Similarly, the L2‐bound RNA‐dependent RNA polymerase's interaction network showed that all the complex's residues formed a strong conventional H bond. The L2‐bound RdRp residues involved in all hydrogen bond formation are shown in Table 3. However, molecule L3 interactions with residues ARG553, ARG555, LYS621, and ARG624 did not affect H‐bonding. ARG553, ARG555, and LYS621 for med pi‐alkyl hydrophobic, ARG555 showed an unfavorable electrostatic interaction with L3, and ARG624 formed donor‐donor exchange (Figure 1).

TABLE 4

Thermodynamics analysis: summary of MM/GBSA‐based binding free energy results of the L1‐RdRp, L2‐RdRp, and L3‐RdRp complexes

Complex	Energy components (kcal/mol)
	ΔEvdW	ΔEelec	ΔGgas	EGB	ESA	ΔGsolv	ΔGbind
L1 ‐RdRp	−48.58 (±3.82)	−32.80 (±7.53)	−81.38 (±8.60)	44.20 (±5.83)	–6.19 (±0.40)	38.02 (±5.60)	43.36 (±4.44)
L2 ‐RdRp	−21.11 (±4.28)	−46.83 (±18.97)	−67.93 (±19.42)	53.25 (±15.72)	–3.54 (±0.66)	49.71 (±15.34)	–18.22 (±6.54)
L3 ‐RdRp	−42.42 (±10.46)	−25.52 (±10.82)	−67.94 (±18.26)	40.14 (±10.83)	–5.57 (±1.34)	34.56 (±9.78)	–33.37 (±10.55)

Abbreviations: ΔEelec, electrostatic; ΔEvdw, Van der Waals; ΔGbind, calculated total free binding energy; ΔGgas, gas‐phase energy; ΔGsolv, solvation free energy.

Post‐molecular dynamics simulations analysis

Binding free energy computations

The binding free energy analysis is a critical parameter for validating ligand‐protein interactions (Chai & Jhon, 2000). The Molecular Mechanics Generalized‐Born Surface Area (MM/GBSA) has been extensively used to compute small molecule ligands' binding free energy (BFE) to biomolecules. Tables 1 and 2 shows the individual total binding power of L1, L2, and L3 to the SARSCOV‐2 RdRp protein. The binding free energy (kcal/mol) values of RdRp to L1, L2, and L3 were −43.36 ± 4.44, −18.22 ± 6.54, and −33.37 ± 10.55, respectively. These MD results indicate that the binding energy of the L1‐RdRp complex is the highest, as revealed by the docking analysis, which agrees with docking binding energy (BE). However, the BFE of L3‐RdRp > L2‐RdRp thus disagreed with the values obtained in the molecular docking calculations. This docking result revealed the BE of L1‐RdRp > L2‐RdRp > L3‐RdRp. There is a significant binding affinity difference (11.73 ± 4.01) between L2‐RdRp and L3‐RdRp. The interaction forces Van der Waals and electrostatic in L1‐RdRp and L3‐RdRp systems contribute to their high total binding energies. Moreover, the formation of hydrophobic interactions between the interacting residues and L1 or L3 increases their BFE. In addition, multiple hydrophobic residues, including VAL, ILE, and ALA, might have characterized the high total binding free energy in the L1‐RdRp and L3‐RdRp systems.

Residue–ligand interaction network

Although we investigated the inhibitory mechanisms of the proposed compounds L1, L2, and L3 in the presence of the cofactors nsp7 and nsp8 proteins, we observed that none of the interacting residues showed binds to a residue of the cofactors. The residue interaction network shows the details of all the protein residues and interaction options formed in a protein–ligand system. Per‐energy decomposition (PRED) analysis assesses individual pocket residue's energy contributions to the total binding free energy. PRED Analysis also reveals molecular insights into the protein dynamics effect of the extent of binding forces. A 40 ns MD simulation was run on complexes (L1‐RdRp, L2‐RdRp, L3‐RdRp) to obtain energy contributions of the residues by computing PRED from the MM/GBSA calculations. Figure 2b shows that ASN395 (−2.35 kcal/mol), VAL455 (−1.85 kcal/mol), THR463 (−1.50 kcal/mol), SER399 (−1.35 kcal/mol), GLY455 (−1.10 kcal/mol), ALA583 (−1.00 kcal/mol), and ILE460 (−1.00 kcal/mol) were the highest contributing residues interacted with compound L1 in the L1‐RdRp system. Strong hydrophobic interactions (Pi‐alkyl) were observed between the VAL458 side chain and the ‐CH3 group of L1. Moreover, alkyl and Pi‐alkyl interactions were formed between VAL455, LYS398, and ALA583 side chain alkyl groups and the benzene ring of L1. In contrast, energetically favorable SER399 and THR463 formed residues hydrogen bonds with hydroxyl groups of the compound.

FIGURE 2

(Color online) Structure of L1‐RdRp at 40 ns (a). Visual illustration of the interaction options in the L1‐RdRp system (b). Per‐residue energy contributions of the RdRp interacting residues (c).

In addition, ASN395, TYR414, ASN441, GLY457, and ILE460 formed at least one conventional or non‐classical hydrogen bond with the terminal ‐OH group of L1. Figure 2a illustrates the binding of L1 in the active site of RdRp, while Figure 2c shows the individual energy contributed to total binding energy by each interacting residue.

Figure 3b illustrates the interaction types binding L2 to the interacting residues of RdRp protein. The hydrogen bonds within L2‐RdRp complex included residues ALA695 (−0.82 kcal/mol), TRP698 (−1.24 kcal/mol), THR699 (−0.04 kcal/mol), and GLU700 (−0.06 kcal/mol), LYS705 (−0.06 kcal/mol), while GLU709 (0.4 kcal/mol) formed Pi‐Pi T‐shaped hydrophobic interactions with benzopyran group of L2 compound. Figure 3a shows the binding of L2 to RdRp active site at 40 ns, while the individual residue and their respective energy contribution are presented in Figure 3c.

FIGURE 3

(Color online) Structure of L2‐RdRp at 40 ns (a). Visual illustration of the interaction options in the L2‐RdRp system (b). Per‐residue energy contributions of the RdRp interacting residues (c).

Lastly, Figure 4b shows the residues involved in the binding interactions between L3 and RdRp, and the energy values affected by each residue are illustrated in Figure 4c. The side chains of GLY581 and ALA586 formed conventional hydrogen bonds, while the GLY581 side chain also formed a carbon‐hydrogen bond. Moreover, the side chains of residues LYS398 and TYR587 began strong Pi‐alkyl and Pi‐Pi T‐shaped hydrophobic interactions with L3 benzene rings. Lastly, the VAL455 side chain interacted with a ‐CH2 of L3 to form an alkyl bond. Figure 4b shows the binding of L3 to RdRp active site at 40 ns.

FIGURE 4

(Color online) Structure of L3‐RdRp at 40 ns (a). Visual illustration of the interaction options in the L3‐RdRp system (b). Per‐residue energy contributions of the RdRp interacting residues (c).

System stability and flexibility

An MD simulation of 230 ns was carried out on the unbound and bound RdRp systems to ensure a detailed investigation of well‐relaxed and equilibrated systems (Alberts et al., 2003; Oluyemi et al., 2022). The stability of three‐dimensional structures of RdRp backbone atoms and their potential energy fluctuations were assessed and compared to the starting structures. The structural changes in the SARS‐CoV‐2 RdRp enzyme were characterized by the binding of L1, L2, and L3 separately in the pocket using the all‐atom molecular dynamics simulations (Figure 5a). Subsequently, the MD trajectory analysis was performed by investigating every printed‐out trajectory for systems' stability, flexibility, and compactness by computing the root mean square of deviation (RMSD), root mean square fluctuation (RMSF), and the radius of gyration (RoG), respectively.

FIGURE 5

The RMSD, RMSF, and RoG plot C‐α atoms of the bound and unbound protein systems. (a) The superimposed structures of L1‐RdRp, L2‐RdRp, and L3‐RdRp. (b) RMSD, (c) RMSF, and (d) RoG plots for ApoRdRp (black), L1‐RdRp (red), L2‐RdRp (blue), and L3‐RdRp (green).

Figure 5b presents the C‐α atoms RMSD of the ligand‐bound and unbound RdRp systems, which revealed differential conformational changes in the stability and convergence of the systems. As shown in the figure, the observed unbound and bound systems vary with time. The mean RMSD values of the residues C‐α atoms for the ApoRdRp, L1‐RdRp, L2‐RdRp, and L3‐RdRp were 2.10, 1.81, 1.68, and 4.42 Å, respectively. The plot showed that the RdRp protein is stable and convergent even upon binding the proposed inhibitors L1 and L2. Contrastingly, L3‐RdRp did not attain stability and convergence until about 92 ns but remained stable and converged throughout the later simulations experiment. The reductive deviations in the L1‐ and L2‐bound RdRp systems could indicate functional and structural inactivity. The results suggest that the binding of compound L3 induced instability on the RdRp protein due to high motions in the backbone atoms. This increased atomistic deviation could imply the mechanism of inhibition of the RdRp functions, which indicate structural activity. The nsp7 and nsp8 residues appeared to be more stable and converged than RdRp residues, as observed in plot 5B. Perhaps, the reason was that nsp7 and nsp8 proteins play stability role in the RdRp protein structure.

Amino acids play a critical role in the conformational structures of proteins (Henchman et al., 2005) conformational changes occur when a mechanical event or chemical reaction takes place (Henchman et al., 2005). Hence, binding a ligand to the active sites often induces conformational changes in its structure and impacts its functions. Specifically, ligand‐induced motion causes changes in the protein's conformation (Ahmad et al., 2013). RMSF matrices measure the average atomic mobility of protein backbone atoms (N, C‐α, and C) during MD simulations (Al‐Sehemi et al., 2020; Lobanov et al., 2008). These findings were further corroborated with residue backbone C‐α atoms RMSF calculation, and flexibility prediction.

Figure 5c presents the RMSF plots of the Apo and bound RdRp systems. These results revealed that ApoRdRp, L1‐RdRp, L2‐RdRp, and L3‐RdRp have mean values of 1.07, 0.97, 0.96, and 6.07 Å, respectively. Compared to the unbound protein, the reductive fluctuation differences in the L1‐RdRp and L2‐RdRp systems seemed insignificant. The highest fluctuations were observed generally in the regions of the nsp7 (932–993) and nsp8 (817–931) protein residues of the bound systems. More significant fluctuations were observed in the cofactors residues regions for the L3‐bound complex, which suggests unstable interaction. In comparison, the minimal shift in the fluctuations of RdRp upon the binding of L1 and L2 could mean adopting an inactive conformation state.

Compactness of RdRp systems

The residues' backbone C‐α atoms RoG of RdRp were calculated for the bound and unbound systems to understand the impact of the proposed inhibitors on the protein. RoG measures the compactness of a protein to provide insights into the changes in the molecular structures (Yanga et al., 2012). A high RoG value meant a less tightly packed structure and increased mobility (Oluyemi et al., 2022). The RoG of the ApoRdRp was used to identify the atomic distribution variations in L1‐, L2‐, and L3‐bound systems as observed in Figure 5d. The mean RoG (Å) values of ApoRdRp, L1‐RdRp, L2‐RdRp, and L3‐RdRp were 30.46, 30.52, 30.57, and 31.00, respectively. However, the RoG values of the L3‐RdRp system were high at different simulation times between 0 and 92 Å. The average RoG value indicates almost insignificant distinctions in the compactness compared to other systems. However, the order of increasing compactness will be L3‐RdRp < L2‐RdRp < L1‐RdRp < ApoRdRp, which implies that the unbound system is most compact with high mobility, while L3‐RdRp is least compact with reduced residue mobility. Furthermore, this result agrees with the RMSF findings of diminished flexibilities for L2‐RdRp and L1‐RdRp systems.

Pseudovirion assay

Pseudovirion assay offers a safe‐effective protocol to study highly infectious and pathogenic viruses such as SARS‐CoV‐2. SARS‐CoV‐2‐spike‐pseudotyped lentiviral particles produced in transfected HEK293T‐hACE2 cells have green fluorescence due to ZsGreen traceable marker. After 48 hr treatment, Chloroquine and Vitexin produced a significant reduction in the number of cells with green fluorescent compared with control cells, which indicates potential anti‐SARS‐CoV‐2 activity of the selected candidate, Vitexin as illustrated in Figure 6.

FIGURE 6

demonstrates the lentiviral particle count on HEK293T‐hACE2 cells upon 48 hr treatment with chloroquine (62.51 μM) and Vitexin (46.25 μM).

CONCLUSION

COVID‐19 pandemic has remained challenging to the general global populace, and the search for a potent drug against the pathogen (SARS‐CoV‐2) has been the research focus of several interdisciplinary Scientists. Therefore, the current study was designed to investigate the inhibitory effect of study compounds against SARS‐CoV‐2 RNA‐dependent RNA polymerase. A potent inhibitor should interact with the critical amino acids of the RdRp active site to prevent the binding of the templated positive‐sense RNA. The blocking of RdRp RNA positive‐sense hinders the further synthesis of new viral proteins. In this study, only the L2 molecule interacted with ASP760 and other crucial residues of the active site. However, L1 and L3 interacted with other residues. However, L2 showed stable interaction with the active site residues with higher binding energy than L3. Although L1 did not bind to many critical active site residues, rather L1 bound to the proximity residues, its binding affinity to RdRp is higher than L2. From the ADME analysis, Orientin, Vitexin, and Kasuagamycin represent good ADME properties and can be developed as promising molecules if further validated properly through in vitro and in vivo conditions.

AUTHOR CONTRIBUTIONS

LRN & BM conceptualized and designed the work. LRN, BM, MK, & BL wrote the paper, AK, MAA, ATA, & MES conducted the experiments.

Funding information

This work is done as a part of the Amrita Rapid‐Response‐Research COVID 19 (ARC) Initiative to LRN and Amrita Vishwa Vidyapeetham PG Research Grant to AK.

CONFLICT OF INTEREST

The author declares no conflicts of interest.

Appendix S1

Click here for additional data file.