In silico prediction of natural compounds as potential multi-target inhibitors of structural proteins of SARS-CoV-2.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused a colossal loss to human health and lives and has deeply impacted socio-economic growth. Remarkable efforts have been made by the scientific community in containing the virus by successful development of vaccines and diagnostic kits. Initiatives towards drug repurposing and discovery have also been undertaken. In this study, we compiled the known natural anti-viral compounds using text mining of the literature and examined them against four major structural proteins of SARS-CoV-2, namely, spike (S) protein, nucleocapsid (N) protein, membrane (M) protein and envelope (E) protein. Following computational approaches, we identified fangchinoline and versicolactone C as the compounds to exhibit strong binding to the target proteins and causing structural deformation of three structural proteins (N, S and M). We recommend the inhibitory effects of these compounds from our study should be experimentally validated against SARS-CoV-2.Communicated by Ramaswamy H. Sarma.

Introduction

Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) from the Coronaviridae family has created an unprecedented healthcare crisis (Khan et al., 2020). It has disseminated around the world with a mortality rate of 2.15% (22 July, 2021) (WHO). Further, the emergence of new variants of SARS-CoV2 and their high transmission rate has created an unparalleled disruption in human life and still poses a big threat to human health and global economy. In the past decades, several epidemics have been caused by coronavirus, which include SARS-CoV in 2003 with a mortality rate of 10% (Cheng et al., 2007; Lee et al., 2003) and Middle East Respiratory Syndrome (MERS) in 2012 with 35% mortality rate (de Groot et al. 2013; Zaki et al., 2012). The pathogenicity of SARS-CoV-2 is relatively less than SARS-CoV and MERS, but it is infectivity is higher than other harmless human coronaviruses (Dömling & Gao, 2020). The recent launch of several vaccines has opened the possibility of containing the infection and providing protection to a large populace. However, effective drugs for the treatment of infections caused by coronaviruses (Wu et al., 2020) are still not available and a hunt for drugs against SARS-CoV2 remains critical for the therapeutic interventions.

For discovery of new drugs, natural compounds serve as rich resources (Lin et al., 2014) and nearly a quarter of approved drugs have plant origins (Thomford et al., 2018). Several well-characterized natural compounds have been documented for their antimicrobial, anti-inflammatory and other beneficial effects on humans and animals thereby offering promising sources for drug development (Clark, 1996; El Sayed, 2000). Advancement in the synthetic biology also has accelerated the mass production of natural compounds (Chen et al., 2021). The ethno-medical literature has listed various herbal plants for their appealing antiviral activity (Ganjhu et al., 2015; Mukhtar et al., 2008; Potroz & Cho, 2015). Natural products include extracts from herbal plants, phytoconstituents, and precise extracts of seed, root, fruit, stem, and flower of plants (Boukhatem & Setzer, 2020; Rates, 2001) Several derivatives of medicinal herbs have been widely considered for antiviral activity (Ganjhu et al., 2015). Antiherpetic acyclovir’, amoebicide ‘emetine’ and antimalarial quinine are good examples of drugs modeled on natural products (El Sayed, 2000; Ganjhu et al., 2015). Some of the medicinal plants have also shown antiviral activity against SARS-CoV-2 (Benarba & Pandiella, 2020; Mesli et al., 2021; Mouffouk et al., 2021)

SARS-CoV-2 genome contains approximately 30,000 bases but encodes only a few proteins (Dömling & Gao, 2020). The structural proteins, namely, spike protein, envelope protein, membrane protein and nucleocapsid protein are essential to complete the viral structure. Spike (S) protein mediates the entry of virus into host cells (Wrapp et al., 2020); membrane (M) protein assembles the virus via protein interactions between membrane-membrane (M-M), membrane-spike (M-S), and membrane-nucleocapsid (N) (Arndt et al., 2010; Kuo et al., 2016). The envelope (E) protein is a short membrane protein (Pervushin et al., 2009), which interacts with the cellular adaptor proteins and further leads to the viral pathophysiology (Jimenez-Guardeño et al., 2014). N protein is the most abundant structural protein, situated inside the virus and plays an essential role in viral genome replication, transcription, packaging and in virion assembly (Sarma et al., 2020). Non-structural proteins such as proteases, helicase and polymerase play an essential role in the life cycle of the virus. These major structural and non-structural proteins are hence considered as important targets for drug discovery (Prajapat et al., 2020). Several studies on drug discovery have been carried out since the pandemic began (Joshi et al., 2020; Nitulescu et al. 2020; Y. Zhou, Hou, et al., 2020). Earlier, screening of the available chemical and natural compounds from PubChem (https://pubchem.ncbi.nlm.nih.gov/) for antiviral activity has also been studied (Wang et al., 2009) among others.

Our aim was to search and screen the antiviral compounds published in literature but are not yet warehoused in the PubChem libraries (https://pubchem.ncbi.nlm.nih.gov/). We have used a multi-target screening approach as pathogens are less likely to develop resistance against inhibitors that affect more than one target (Xie & Xie, 2019). Since the structural proteins play essential role in the entry, assembly, packaging and release of the virus (Yadav et al., 2021) they serve as promising drug targets. Among all the structural proteins, spike protein has been widely studied (Bojadzic et al., 2020; Unni et al., 2020; Xia et al., 2020). In this study, we have considered all four structural proteins in the search of multi-target inhibitors. Hence, same sets of compounds were screened against each of the 4 structural proteins of SARS-CoV-2 virus.

Material and methods

Data mining

Natural compounds with reported antiviral activity were extracted from the literature using pubmed.mineR (https://cran.r-project.org/web/packages/pubmed.mineR/index.html) (Rani et al., 2015). Total 13,421 PubMed abstracts with query term ‘antiviral natural compounds’ were downloaded and processed.

Molecular docking

The antiviral natural compounds were screened against all four structural proteins of SARS-CoV-2. We examined compounds, which could bind to more than one target protein. To obtain and analyze all possible ligand-receptor complexes, blind docking was carried out using AutoDock Vina 4.0 (Trott & Olson, 2010).

Protein structures

Structure of S protein (PDB ID: 6VSB) with resolution 3.46 Å bound with 2-acetamido-2-deoxy-beta-D-glucopyranose (NAG) (Wrapp et al., 2020) and C-terminal domain (CTD) (PDB ID: 7C22) with resolution 2.0 Å (R. Zhou, Zeng, et al., 2020) of nucleocapsid protein (N) bound with di(hydroxyethyl)ether (PEG) and acetate ion (ACT) were obtained from PDB (Rose et al., 2017). We worked on the CTD of N protein as fewer studies were reported on it than its N-terminal domain. While selecting these structures, we chose the structures showing least number of missing residues. Because the structures of E and M proteins have not yet been solved and published in PDB, I-Tasser, a hierarchical protocol for the prediction of structure and function of the protein (Yang & Zhang, 2015) was used to download the structure of E and M proteins (https://zhanglab.ccmb.med.umich.edu/COVID-19/).

Protein structure preparation

Each of the four proteins (S, M, E and N (CTD)) were individually prepared for docking using Autodock tools (ADT) (El-Hachem et al., 2017), which included addition of hydrogen and removal of water. Minimization of structures was performed using Chimera software (Pettersen et al., 2004). The grid was set large enough to include the whole structure for all possible ligand-receptor interactions. A shell script was used to implement docking of each protein with all natural compounds.

Visual rescoring

Docking was carried out in triplicate and the top hits were recorded using visual inspection of binding residues using PyMol software (Rigsby & Parker, 2016). Antiviral drugs Remdesivir, Favipiravir, Ritonavir, Lopinavir and Hydroxychloroquine were taken as a positive control (Bhatnagar et al., 2020; Cai et al., 2020; Vanden Eynde, 2020, p. 19). LigPlot (Wallace et al., 1995) was used to explore the hydrophobic interactions between proteins and compounds.

Molecular dynamic simulation

Analysis of internal motions of a protein can be helpful for the identification of multiple biological functions along with their profound dynamic mechanisms. Compounds, which showed high binding affinity and also targeted more than one SARS-CoV2 protein (Figure 1) were further analyzed by molecular dynamic simulation to study the consistency and flexibility of the binding mode of selected inhibitor-protein complexes. MD simulation was carried out by using GROMACS (version 5.0) (Van Der Spoel et al., 2005). The topology and coordinates files were prepared by using All-atom GROMOS53a6.ff force field, and the inhibitor topology file was generated from PRODRG server (Schüttelkopf & van Aalten, 2004). Each complex was neutralized by adding a varied number of Na+ ions and enclosed in a periodic box TIP3P water and extended to 10 Å from the solvent. The water model Simple Point Charge (spc216) was used for the solvation. The Particle Mesh Ewald method (Darden et al., 1999) and LINCS algorithm was used to restrain the long-range electrostatics and bond length, respectively. Prior to simulation the protein was minimized by using the steepest descent method. The entire system was equilibrated in two stages as temperature coupling step (NVT) and pressure coupling step (NPT). Once the system was equilibrated, the conformational sampling production was started at constant temperature 300 K and pressure 1 atm MD simulation trajectory was collected for every 2 ps time steps.

Figure 1.

Summary of strategy followed to identify consensus compounds.

Trajectory analysis

GROMACS toolkit was used for the analysis of MD simulation trajectory. Root mean square deviation (RMSD) of initial structure to the completion of simulation was calculated. Root mean square fluctuation (RMSF) determining per residue fluctuation and hydrogen bonds between the compound(s) and the target protein were analyzed to measure the stability of the protein-ligand complexes.

Gibbs free energy landscape

The Gibbs free energy landscape provides the structural features and computational profiles of a protein by using conformational sampling methods (Hsieh et al., 2010). It was projected onto the first principal component (PC1) and second principal component (PC2) with highest eigenvalues calculated from Principal Component Analysis (PCA). It can be defined as

kB  =  Boltzmann constant

T = Temperature

P (PCS1, PC2) = normalized joint probability distribution

Results and discussion

SARS-CoV-2 genome is made of 29,891 nucleotides with 9860 amino acids (Chan et al., 2020), packaged in a circular nucleocapsid protein and encapsulated by envelope proteins (Li, 2016). All four structural proteins included in this study, namely, Spike protein, Membrane protein, Envelope protein, and Nucleocapsid protein play significant roles in the life cycle of the virus and hence are promising therapeutic targets to control COVID-19 (Satarker & Nampoothiri, 2020). In the case of N protein, since its N terminal domain has been widely studied, we have considered only the C terminal domain (CTD) in our study.

We used algorithmic text mining approach to extract literature and retrieve the compounds with known anti-viral activity. Significantly, several of these compounds are reported to be effective, but are not yet included in chemical databases. Text mining also enabled access to more potential variants of the promising antiviral compounds. Our aim was to identify multi-target inhibitor(s) from among these natural compounds, which could potentially inhibit more than one protein in SARS-CoV2. Multi-target inhibitors can not only be more toxic, probability of the virus gaining resistance against such compounds is expected to be low. Keeping this significance in perspective, molecular interaction studies were done between each of the identified natural compounds with the four structural viral proteins. Out of the identified compounds, 7 showed multi-target binding in docking analysis, which we have mentioned as the ‘consesus compounds’. Stability of interaction between the compounds and the target proteins was further evaluated using molecular dynamics simulation.

On text mining of the PubMed abstracts, a total of 106 natural compounds (Supplemental Table 1) with antiviral activity were identified. The available variants of these compounds are also included in the study. Each compound was screened against the four structural proteins of SARS-CoV-2.

Table 1.

Functional and structural details of the seven consensus compounds.

Docking of compounds with each protein showed uneven binding affinity. Hence, from the docking output (Supplemental Table 2), first top hits for all four proteins were considered (Figure 1). Next, the top hit compounds for each of the four proteins were searched to find a consensus compound which could target more than one protein (multi-target effect). Based on binding affinity, 7 compounds were identified as consensus compounds (Table 1), and were further considered for validation by simulation using GROMACS.

Table 2.

Illustration of binding affinity of consensus compounds with and binding residues in the respective proteins.

	Envelope protein		Membrane protein		Nucleocapsid protein(CTD)		Spike protein
Compounds	Binding affinity (kcal/mol)	Binding residues	Binding affinity (kcal/mol)	Binding residues	Binding affinity (kcal/mol)	Binding residues	Binding affinity (kcal/mol)	Binding residues
Fangchinoline	–6.8	LYS53, SER50, ASN48	–7.7	SER197	–9.8	GLN260	–8.3	TYR904
Cepharanthine	–6.9	None	–8.4	None	–10.5	None	–8.1	TYR313
Compound 55	–6.3	ALA32	–8.3	ALA218, SER99	–10.3	GLN260, TYR333	–7.4	SER591
Tetrandrine	–6.2	None	–7.5	ARG107	–9.4	None	–8.1	SER50, GLN853
Compound 3	–6.1	SER55	–7.7	SER99, ASP209	–9.1	ARG277, ALA264, THR296	–7	CYS538, PRO589, SER591
Compound 6	–6.8	None	–8.2	SER4. ARG107	–9.6	AG277	–7.4	ALA288, SER297
Variant of Harmine	–7.6	None	–8	GLY78, TYR196	–9.9	None	–7.5	None
Remdesivir	–5	PHE4, SER6	–6.9	LEU181, SER197, ARG198, VAL139, GLU137, THR77	–8.3	TYR333, GLY335, ILE337	None	None
Favipiravir	–3.6	ASP72, VAL75	–5	GLU137, VAL139, GLY78	–4.3	GLN283, ASN285	–4.5	THR1077
Ritonavir	–5.8	SER55, TYR57, TYR59, ARG61	–6.9	ARG44, ARG101, ASP209	–7.7	ARG277	None	None
Lopinavir	–5.1	ASP72, VAL75	–7.4	HIS125, TRP20	–8.5	GLN260, ARG259, ARG277	–6.8	SER172, GLN173, LEU226, PRO225
Hydroxychloroquine	–4.9	None	–5.7	ALA188	–6	ILE304	–4.8	ASP294, SER297, ASP287
ACT (ACETATE ION)	None	None	None	None	–2.9	Thr271, ARG276, ARG277	None	None
PEG DI(HYDROXYETHYL)ETHER	None	None	None	None	–2.8	PHE286	None	None
NAG2-acetamido-2-deoxy- beta-D-glucopyranose	None	None	None	None	None	None	–4.6	None

Interaction between compounds and proteins

The visual rescoring was performed for deep understanding of binding between ligands and targets (Table 2). Consensus compounds Table 1 were then subjected to PyMol (Rigsby & Parker, 2016). Compounds that did not show hydrogen bonding with the respective target proteins were further searched for hydrophobic interactions using the LigPlot tool (Wallace et al., 1995) (Figure 2).

Figure 2.

Hydrophobic interactions between proteins (E, M, S and N) and consensus compounds. A, B, C and D represent hydrophobic interactions between Envelope protein (E), Membrane protein (M), Spike protein (S) and Nucleocapsid protein (N-CTD), respectively. The red colored semi circles represent the protein residues forming hydrophobic interaction with the protein.

Hydrophobic interactions

After analyzing the hydrogen bonding between the compounds and proteins, compounds lacking hydrogen bonding were screened for hydrophobic interactions using LigPlot as shown in Figure 2.

Molecular dynamics simulation

All four proteins were simulated with each of the consensus compounds regardless of the number of binding residues (except for the compound 6 for E protein, which neither showed hydrogen bonds nor hydrophobic interactions). The preliminary simulation filtered out the least stable and weakly bound compounds and for the remaining compounds, simulation was extended for 100 ns.

Nucleocapsid protein (C terminal domain)

Structure deviations and fluctuations (RMSD and RMSF and RG)

Root mean square deviation, root mean square fluctuations and the radius of gyration were analyzed from the stable region of trajectories to find the structural dynamics of C-terminal domain (CTD) of the N protein complexed with identified compounds, namely, fangchinoline, cepharanthine and compound 3. The average RMSD of apo-CTD protein was 0.60 nm and the average RMSD of CTD-fangchinoline, CTD-cepharanthine and CTD-Compound 3 was 0.54 nm, 0.65 nm, and 0.67 nm, respectively (Figure 3A). It was observed that CTD tends to show stable and low fluctuations when bound to fangchinoline, whereas with cepharanthine and compound 3, fluctuations were observed to be high.

Figure 3.

Structural dynamics of N-CTD. (A) Root mean square deviation (RMSD) plots for CTD as a function of time. (B) RMSD of compounds as a function of time. (C) Root mean square fluctuations (RMSF) plot for CTD. (D) Time evolution of radius of gyration (Rg). Black, blue, green and red color represent values obtained for apo-CTD, CTD-fangchinoline, CTD-cepharanthine and CTD-compound 3 during 100 ns MD simulations, respectively.

The average RMSD values of fangchinoline, cepharanthine and compound 3 were found to be 0.05 nm, 0.04 nm, and 0.11 nm, respectively. Cepharanthine showed least RMSD fluctuation as compared to fangchinoline and compound 3 (Figure 3B).

CTD-fangchinoline showed least residual fluctuations (RMSF) as compared to CTD-cepharanthine and CTD-compound 3, which could be due to strong hydrogen bonding between the respective ligand-receptor complexes (Figure 3C). Radius of gyration (Rg) is linked to the tertiary structural volume of fangchinoline, cepharanthine and compound 3. Proteins with higher Rg are assumed to have less tight packing (Lobanov et al., 2008). The average Rg values for apo-CTD, fangchinoline, cepharanthine and compound 3 were found to be 1.32, 1.39 nm, 1.43 nm, and 1.40 nm, respectively (Figure 3D). Fangchinoline showed tight packing in comparison to the other two compounds. Based on these structural dynamics, fangchinoline appears to attain the highest stable conformation.

SASA and hydrogen bonds

The solvent accessible surface area (SASA) is defined as the surface area of the networks resulting from the complex of a protein with its solvent molecules (Durham et al., 2009). The average SASA values with respect to apo-CTD protein, fangchinoline, cepharanthine and compound 3 were found to be 71.74 nm2, 79.30 nm2, 80.78 nm2, and 77.15 nm2, respectively (Figure 4A). A lower SASA value in the case of compound 3 suggests that its internal residues are not exposed to the solvent.

Figure 4.

Solvent Accessible Surface Area (SASA) and Hydrogen bonds. (A) Solvent Accessible Surface Area as a function of time and (B) the free energy of solvation with reference to time. Black, blue, green and red color represent values obtained for apo-CTD, CTD-fangchinoline, CTD-cepharanthine and CTD-compound 3 during 100 ns MD simulations, respectively. (C) Hydrogen bonds between protein and ligand obtained for fangchinoline, cepharanthine, and compound 3 with CTD.

Free solvation for apo-CTD, fangchinoline, cepharanthine and compound 3 was found to be 145.89, 159.56 KJ/mol/nm2, 170.40 KJ/mol/nm2, and 146.06 KJ/mol/nm2, respectively (Figure 4B). Additionally, the SASA plots were also resolved into hydrophobic and hydrophilic regions. Hydrophobic SASA with respect to backbone for apo-CTD was 26.21 nm2 and SASA with respect to backbone for CTD in complex with fangchinoline, cepharanthine and compound 3 was found to be 27.15 nm2, 28.63 nm2, and 26.57 nm2, respectively. Hydrophilic SASA with respect to backbone for apo-CTD, CTD-fangchinoline, CTD-cepharanthine and CTD-compound 3 was 35.98 nm2, 38.88 nm2, 37.16 nm2, and 36.86 nm2, respectively.

Hydrogen bonds play a vital role in stabilizing the protein conformations. Hydrogen bonding between CTD and compounds, namely, fangchinoline, cepharanthine and compound 3 were analyzed by estimating the hydrogen bond pair within 0.35 nm. Here, compound 3 was found to establish more hydrogen bond contacts with CTD, followed by fangchinoline. Cepharanthine showed weakest interaction (Figure 4C).

Secondary structure analysis

To spot the structural features of CTD, secondary structure analysis was performed. The secondary structure assignments, such as α-helix, β-sheet and turn were split into individual residues and the average number of residues participating in the formation of secondary structure were compared. The percentage of residues participating in structure formation during MD simulation is given in Table 3.

Table 3.

Percentage of residues participating in average structure formation during MD simulations.

Percentage of secondary structure (SS %)
Protein	Structurea	Coil	β-sheet	β-bridge	Bend	Turn	α-helix	310-helix
CTD	42	38	8	2	2	10	22	1
CTD-fangchinoline	41	40	5	2	15	11	23	2
CTD-cepharanthine	44	34	9	2	21	9	25	1
CTD-compound 3	50	31	11	4	18	12	23	1

Structure = α-helix + β-sheet + β-bridge + Turn.

This analysis revealed the average residues participating in structure formation in CTD to be 42%, and after binding with fangchinoline, cepharanthine and compound 3, the percentage of average number of residues participating in structure formation were 41%, 44%, and 50%, respectively (Figure 5). The CTD-fangchinoline complex showed the least values of average secondary structure during MD simulation. Binding of fangchinoline resulted in increased coils as compared to other two compounds, which indicates it might have a denaturing effect on the protein.

Figure 5.

The secondary structure plot as a function of time. Each structural element has been represented by different colors: black (overall structure), red (coil), green (β-sheet), blue (β-bridge), yellow (bend), brown (turn), gray (α-helix), and cyan (3-helix). A) Apo-CTD, B) CTD-fangchinoline, C) CTD-cepharanthine D) CTD-compound 3.

Further, the volume and density of above complexes were calculated. The volumes for apo-CTD, CTD-fangchinoline, CTD-cepharanthine and CTD-compound 3, were 23.07 nm3, 23.63 nm3, 27.75 nm3, and 23.46 nm3, respectively. The average densities were 890.59 g/l, 869.56 g/l, 865.06 g/l, and 875.53 g/l, respectively (Figure 6). The results suggest minimal variations in the volume and density of CTD in complexes with all the three compounds.

Figure 6.

The structural volume and density plot as a function of time. (A) Structural volume and (B) Density of Apo-CTD (black), CTD-fangchinoline (Blue), CTD-cepharanthine (Green), and compound 3 (red) calculated during MD simulations, respectively.

Membrane protein (M protein)

Structure deviations and fluctuations (RMSD and RMSF and Rg)

The structural dynamics of M protein complexed with compound 3 was calculated by root mean square deviation, root mean square fluctuations and the radius of gyration, analyzed from the stable region of their trajectories and was found to show low fluctuations in complex with compound 3. Average RMSD of apo-M protein and complex with compound 3 was 0.55 nm (Figure 7A) and the average RMSD of compound 3 was 0.23 nm (Figure 7B). Average RMSF fluctuations in apo-M protein and M Protein-compound 3 was 0.19 and 0.23 nm, respectively (Figure 7C) and average Rg values for apo-M protein and M Protein-compound 3 complex was found to be 1.63 nm and 1.68 nm, respectively (Figure 7D).

Figure 7.

Structural Dynamics of M Protein. (A) Root mean square deviation (RMSD) plots for M Protein as a function of time. (B) RMSD of compound 3 as a function of time. (C) Root mean square fluctuations (RMSF) plot for M Protein. (D) Time evolution of radius of gyration (Rg). Black, and red color represent values obtained for apo-M Protein, and M Protein-compound 3, respectively during 100 ns MD simulations.

The average SASA values for apo-protein and for Compound 3 were found to be 119.62 nm2, and 127.02 nm2, respectively (Figure 8A). High SASA value suggested internal residues in compound 3 to be exposed to the solvent.

Figure 8.

Solvent Accessible Surface Area (SASA) and hydrogen bonds. (A) SASA as a function of time and (B) the free energy of solvation with reference to time. (C) Hydrogen bonds between compound 3 and M Protein. Black and red color represent values obtained for apo-M protein, and M protein-compound 3, respectively during 100 ns MD simulations.

The free solvation for apo-M protein and compound 3 was 237.24 KJ/mol/nm2 and 240.32 KJ/mol/nm2, respectively (Figure 8B). Additionally, the SASA plots were also resolved into hydrophobic and hydrophilic regions. The hydrophobic SASA with respect to backbone for M protein was 68.58 nm2, whereas SASA with respect to backbone for M protein in complex with compound 3 was 71.09 nm2. The hydrophilic SASA with respect to backbone for apo-M protein, and M protein-compound 3 were found to be 58.43 nm2 and 57.42 nm2, respectively.

Hydrogen bond formations between M protein and compound 3 were analyzed by estimating the hydrogen bond paired within 0.35 nm and compound 3 was observed to establish a high number of hydrogen bond contacts with M protein (Figure 8C).

The secondary structure assignments, namely, α-helix, β-sheet and turns were split into individual residues and the average number of residues participating in the formation of secondary structure were compared. The percentage of residues participated in structure formation during MD simulation are given in Table 4.

Table 4.

Percentage of residues participating in average structure formation during MD simulations.

Percentage of secondary structure (SS %)
Protein	Structurea	Coil	β-sheet	β-bridge	Bend	Turn	α-helix	310-helix
M protein	46	34	10	2	19	12	22	1
M protein-compound 3	41	37	12	3	21	9	17	1

Structure = α-helix + β-sheet + β-bridge + Turn.

In this analysis we found that average residues participating in the structure formation in M protein when free was 46% and on binding with Compound 3 was 41%, (Figure 9). Additionally, an increase in the coils was also observed, suggesting the denaturing effect of the compound on the protein.

Figure 9.

The secondary structure plot as a function of time. Each structural element A) M protein, B) M protein-compound 3 complex, has been represented by different colors. Black (overall structure), red (coil), green (β-sheet), blue (β-bridge), yellow (bend), brown (turn), gray (α-helix), and cyan (3-helix).

Further, on calculating the volume and density of the above shown complexes, less variation was noticed in these two parameters. Volumes of apo-M protein and M protein-compound 3 were 45.70 nm3 and 46.37 nm3, respectively (Figure 10) and the average density were 914.04 g/l and 900.86 g/l, respectively.

Figure 10.

Structural volume and density plot as a function of time. (A) Structural volume and (B) Density of Apo-M protein (black), and compound 3 (red) calculated during MD simulations, respectively.

Spike protein (S protein)

Structure deviations and fluctuations (RMSD, RMSF and RG)

The structural dynamics of S protein complexed with compound 3 was calculated by root mean square deviation, root mean fluctuations and the radius of gyration, analyzed from the stable region of their trajectories and were observed to be low fluctuations in complex with compound 3. The average RMSDs of apo-S protein and its complex with compound 3 was 1.90 nm and 2.21 nm, respectively (Figure 11A). The RMSDs of compound 3 was observed to be 0.22 nm (Figure 11B). We observed low fluctuations in S protein complexed with compound 3. The average RMSF fluctuations in apo-S protein and S Protein-compound 3 complex was found to be 0.72 nm and 0.84 nm, respectively (Figure 11C) and the average Rg values for apo-S protein and S Protein-compound 3 complex were found to be 4.17 nm, and 4.09 nm, respectively (Figure 11D).

Figure 11.

Structural Dynamics of S Protein. (A) Root mean square deviation (RMSD) plots for S Protein as a function of time. (B) RMSD of compounds 3 as a function of time. (C) Root mean square fluctuations (RMSF) plot for S Protein. (D) Time evolution of radius of gyration (Rg). Black, and red color represent values obtained for apo-S Protein, and S Protein-Compound 3, respectively, during 100 ns MD simulations.

The average SASA values for apo-S protein and for compound 3 were found to be 486.76 nm2, and 492.41 nm2, respectively (Figure 12A). High SASA values suggested internal residues in compound 3 to be exposed to the solvent.

Figure 12.

Solvent Accessible Surface Area (SASA) and hydrogen bonds. (A) SASA as a function of time and (B) the free energy of solvation with reference to time (C) Hydrogen bonds between protein and ligand obtained for compound 3 with S protein. Black and red color represent values obtained for apo-S protein and S protein-compound 3 complex during 100 ns MD simulations, respectively.

Free solvation for apo-S protein and compound 3 was 829.97 KJ/mol/nm2, and 863.18 KJ/mol/nm2, respectively (Figure 12B). Hydrogen bonds between S protein and compound 3 within 0.35 nm were observed to be significant in number (Figure 12C).

The secondary structure assignments for α-helix, β-sheet and turn were split into individual residues and the average number of residues participating in the formation of secondary structure were compared. The percentage of residues participated in structure formation during MD simulation are given in Table 5.

Table 5.

Percentage of residues participated in average structure formation during MD simulations.

Percentage of secondary structure (SS %)
Protein	Structurea	Coil	β-sheet	β-bridge	Bend	Turn	α-helix	310-helix
S protein	59	25	28	1	14	8	22	1
S protein-compound 3	59	26	29	2	14	8	21	1

Structure = α-helix + β-sheet + β-bridge + Turn.

In this analysis we obtained that average residues participating in structure formation in S protein are 59%, which remain the same after binding with compound 3 (Figure 13). No alteration in the secondary structure of S protein was observed.

Figure 13.

The secondary structure plot as a function of time. Each structural element A) S protein, B) S protein-compound 3 complex, has been represented by different color such as black (overall structure), red (coil), green (β-sheet), blue (β-bridge), yellow (bend), brown (turn), gray (α-helix), and cyan (3-helix), respectively.

Further, on calculating the volume and density of the above shown complexes, less variation was noticed in these two parameters. Volumes of apo-S protein and S protein-compound 3 were 186.76 nm3 and 187.31 nm3, respectively (Figure 14A) and the average density was 944.20 g/l and 941.39 g/l, respectively (Figure 14B). These results depict less variation in the volume and density of each complex.

Figure 14.

Structural volume and density plot as a function of time. (A) Structural volume and (B) Density of Apo-M protein (black), and Compound 3 (red) calculated during MD simulations, respectively.

Conclusions

Development of new drugs is a time taking and expensive process. Comprehensive efforts are being made globally towards the search of therapeutics against SARS-CoV-2 (Liu et al., 2020; Ojha et al., 2021; Rabaan et al. 2020). Given the urgent requirement for treatment options for Covid19 disease, great thrust has been given on drug repurposing also. Several drugs such as, remdesivir, favipiravir, ritonavir, and lopinavir (Bhatnagar et al., 2020; Cai et al., 2020; Vanden Eynde, 2020, p. 19) have been included in the treatment regimen and shown effective results in several cases. However, to mitigate the current crisis and prepare for future pandemics, there is a continuous need to explore the diversity of natural compounds for their antiviral activity. In this study, we extracted information available in the literature on antiviral phytochemicals and evaluated their potential to interact and potentially inhibit the structural proteins of SARS-CoV-2.

Seven consensus compounds which showed binding to all the four structural proteins on docking are fangchinoline, cepharanthine, compound 55, tetrandrine, compound 3, compound 6 and harmine. While determining the binding residues in each of the four proteins with these 7 compounds using PyMol, we did not observe H-bonding between the compounds tetrandrine, cepharanthine, harmine, and compound 6 with any of the target proteins so we also examined hydrophobic interactions using LigPlot (Figure 4). Compound 6 did not show hydrophobic interaction either (with E protein) and hence was excluded from the simulation studies. Further, four compounds which showed weak binding and low stability during the preliminary simulation were also excluded from this study. Finally, fangchinoline, cepharanthine and compound 3 were studied further by RMSD, RMSF and hydrogen bonding analysis.

Fangchinoline is a bisbenzylisoquinoline alkaloid sourced from Stephania tetrandra (Kim et al., 2019; Weber & Opatz, 2019), known for its anticancer (Mérarchi et al., 2018) and anti-inflammatory (Choi et al., 2000) activities. Fangchinoline has earlier been shown as an inhibitor of the replication of human immunodeficiency virus (HIV). It interfers with proteolytic processing of gp160 in HIV (Wan et al., 2012) and recently its potential in the prevention and treatment of human coronavirus (HCoV-OC43) infection (Kim et al., 2019) has also been reported. In our analysis, not only does fangchinoline bind stably to CTD of N protein, the CTD-fangchinoline complex showed least values of average secondary structure during MD simulation. Further, binding of fangchinoline resulted in increased coils in CTD, indicating it might also have a denaturing effect on the protein.

Aspergillus spp are widely considered as abundant source of secondary metabolites, namely, peptides, alkaloids, poly-ketides, lignans and terpenes with biological properties (Ge et al., 2010; H.-W. Zhang et al., 2010; M. Zhou et al., 2015, 2014). Compound 3, a butyrolactone (versicolactone C) is a subclass of lignans produced by fungi and plants (Matsumoto et al., 2006; Mohagheghzadeh et al., 2006). Butyrolactones are largely known for their anti-bacterial, anti-viral, anti-inflammatory activities ( Haritakun et al., 2010; Qin et al., 2011; H.-J. Zhang et al., 2005). Similar to fangchinoline, compound 3 also causes denaturation of secondary structure in M protein by increasing the number of coils. In S protein, though its secondary structure was not found to be affected, the RMSD and hydrogen bond analysis revealed adequate stability of interaction between the compound.

Taken together, from among the reported bis-benzylisoquinoline alkaloids (tetrandrine, fangchinoline, and cepharanthine), our in silico analysis shows fangchinoline to be a potential inhibitor of NCTD and Compound 3 (versicolactone C) exhibits multi-target binding on the key structural proteins of the coronavirus. Based on these encouraging results, we believe these compounds could be tested experimentally for inhibitory effect on SARS-CoV2 virus for their prospects in antiviral drug discovery.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.