A computational prediction of SARS-CoV-2 structural protein inhibitors from Azadirachta indica (Neem).

The rapid global spread of the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) has created an unprecedented healthcare crisis. The treatment for the severe respiratory illness caused by this virus is primarily symptomatic at this point, although the usage of a broad antiviral drug Remdesivir has been allowed on emergency basis by the Food and Drug Administration (FDA). The ever-increasing death toll highlights an urgent need for development of specific antivirals. In this work, we have utilized docking and simulation methods to identify small molecule inhibitors of SARS-CoV-2 Membrane (M) and Envelope (E) proteins, which are essential for virus assembly and budding. A total of 70 compounds from an Indian medicinal plant source (Azadirachta indica or Neem) were virtually screened against these two proteins and further analyzed with molecular dynamics simulations, which resulted in the identification of a few common compounds with strong binding to both structural proteins. The compounds bind to biologically critical regions of M and E, indicating their potential to inhibit the functionality of these components. We hope that our computational approach may result in the identification of effective inhibitors of SARS-CoV-2 assembly.Communicated by Ramaswamy H. Sarma.

Introduction

Natural compounds like plant products constitute a rich resource for drug discovery (Patridge et al., 2016; Thomford et al., 2018; Wani et al., 1971). Anticancer drugs like taxol, and antimalarial compounds like quinine, were derived from plant sources like Taxus brevifolia and Cinchona spp. (Patridge et al., 2016; Thomford et al., 2018; Wani et al., 1971). Approximately a quarter of FDA or European Medical Agency (EMA) approved drugs are plant based (Thomford et al., 2018), which highlights the importance of plant-inspired compounds in the biomedical arena. Medicinal plants are thought to be a good source for antiviral compounds against multiple viruses including SARS-CoV-2 (Aanouz et al., 2020; Abdelli et al., 2020; Akram et al., 2018; Enmozhi et al., 2020). However, there are major obstacles in generating synthetic drug candidates by mimicking naturally existing compounds. Some common problems include difficulties in extraction of components from original sources and characterization, identifying potential targets, and setting up effective assays for measuring drug efficacy, safety and pharmacokinetics. Recent advances in computational methods can be effective in screening potential targets for newly identified molecules or repurposing licensed drugs.

Neem (Azadirachta indica), a member of the Meliaceae family, is a well-known medicinal plant, especially in the Indian subcontinent. The chemical constituents of Neem include azadirachtin, 7-desacetyl-7-benzoylazadiradione, 17-hydroxyazadiradione, 7-desacetyl-7-benzoylgedunin, nimbin, nimbiol, polyphenolic flavinoids, etc (Alzohairy, 2016). Ethanol extracts of Neem leaves have been shown to exhibit anti-microbial properties, and Neem components have demonstrated free radical scavenging and anti-inflammatory activities (Alzohairy, 2016). Neem bark (NBE), when used at concentrations ranging from 50 to 100 μg/mL, has inhibitory effect on Herpes Simplex Virus (HSV) type-1 propagation (Tiwari et al., 2010). Cells treated with NBE inhibited HSV-1 glycoprotein mediated cell-to-cell fusion and polykaryocyte formation, suggesting a potential role of NBE at the viral fusion step. Leaf extract of Neem (Azadirachta indica A. Juss.) (NCL-11) has virucidal activity against Coxsackievirus B-4, and is thought to interfere at an early stage of the virus replication cycle (Badam et al., 1999). The proven antimicrobial property, and low toxicity of Neem extracts, makes this plant an excellent choice for harvesting and designing of potential antiviral components. A recent study has shown the beneficial effect of Neem bark on neuroinflammation caused by Mouse Hepatitis Virus (MHV), a murine coronavirus, by preventing cell-to-cell spread (Sarkar et al., 2020).

The novel coronavirus/SARS-CoV-2, like other members of the β-coronavirus family, utilizes several structural and non-structural proteins for receptor binding and cellular entry, replication, assembly and cell-to-cell spread. The Spike (S) glycoprotein, RNA-dependent RNA polymerase (RdRp) and 3C-like Main Protease (3CLpro) are utilized for cellular entry and membrane fusion, viral RNA replication and viral polyprotein processing respectively. The 3D structures of these proteins have been resolved by X-ray crystallography or cryoelectron microscopy (Gao et al., 2020; Walls et al., 2020; Zhang et al., 2020). So far, various in vitro and computational studies (Bhardwaj et al., 2020; Caly et al., 2020; Elfiky, 2020; Wu et al., 2020; Zhang et al., 2020) have utilized these components as the main targets for drug screening and repurposing from various sources. Several other less explored protein components, for which 3D structures have not yet been resolved, also have essential roles to play in the viral life cycle and pathophysiology, and can be important drug targets. The membrane (M) protein is involved in virus assembly through M-M, M-S, and M-nucleocapsid (N) protein interactions (Arndt et al., 2010; Kuo et al., 2016). The Envelope (E) protein is a homopentameric, short, integral membrane protein of 76–109 amino acids (8.4 to 12 kDa) (Pervushin et al., 2009; Torres et al., 2006). E is a viroporin that appears to assist virus budding through an unknown mechanism (DeDiego et al., 2007). It also interacts with cellular adapter proteins through its C-terminal (PDZ)-binding motif (PBM), which contributes to cell-cell spread and viral pathophysiology (Jimenez-Guardeno et al., 2014; Schoeman & Fielding, 2019; Teoh et al., 2010; Yang et al., 2005). Although assembly of the viral envelope is coordinated by M, both M and E are required for the production and release of particles (Mortola & Roy, 2004). The removal of E protein from SARS-CoV leads to formation of immature or infection-incompetent progeny, indicating the critical role of E in Coronavirus biology (DeDiego et al., 2007)

We attempted docking of Neem compounds with modeled 3D structures of SARS-CoV-2 proteins M and E, which are essential for formation of virus particles. The top scoring compounds with the highest binding affinity were subjected to 100 ns Molecular Dynamics (MD) simulations to detect stability of binding. The results from these in silico binding studies are reported. It is possible that these compounds may prevent assembly of SARS-CoV-2particles, thus reducing viral propagation. A combination of viral replication and assembly inhibitors may be a more effective regimen for therapeutic intervention.

Materials & methods

Target identification and homology modeling

The protein sequences for SARS-CoV-2 E and M proteins were obtained from the NCBI deposits of Wuhan-Hu-1 isolate (YP_009724392.1 and QHD43419.1 respectively). A 3D model for SARS-CoV-2 E was generated from the ITASSER server (Roy et al., 2010) by using the NMR structure of the SARS-CoV E protein as template (PDB ID: 5 × 29) (Surya et al., 2018), which has 88.71% sequence identity with its SARS-CoV-2 analogue. A pentameric form of the model was generated using MODELLER (Eswar et al., 2006). A 3D model for SARS-CoV-2 M was also generated from the ITASSER server. The models were further refined and energy minimized using the 3D refine server (Bhattacharya et al., 2016). 3D models were validated by Ramachandran plot analysis using the RAMPAGE server (http://mordred.bioc.cam.ac.uk/∼rapper/rampage.php).

Virtual screening and molecular docking

A natural compound ligand library from the Indian Medicinal Plants, Phytochemistry and Therapeutics (IMPPAT) database (Mohanraj et al., 2018) was used for docking. A set of 70 compounds from Neem (Azadirachta indica) was downloaded from the database (Supplementary material Table S1). Virtual screening of the compounds with the 3D models of M and E was carried out using AutoDock Vina (Trott & Olson, 2010) with default configuration parameters. The size of the grid box was chosen to encompass all possible binding sites of E and M. The top-ranked five poses of Azadirachta indica compounds bound to each structural protein were used for further analysis.

Molecular dynamics simulation

The top scoring complexes were subjected to MD simulations for further analysis. Simulations were performed using GROMACS 5 with Gromos force field (Oostenbrink et al., 2004). The parameters and topologies of the ligands were calculated by PRODRG server prior to MD simulation (Schüttelkopf & van Aalten, 2004). Each complex was immersed in a cubic box of 1.2 nm containing SPC water molecules. NA+ and CL− ions were added to neutralize the charge of the systems. Energy minimization was performed using the steepest descent method for 50,000 steps for all systems with a tolerance of 1000 kJ mol−1nm−1. For long-range interactions, the PME method was used with a 1.2 nm cut-off and a Fourier spacing of 0.16 nm. Equilibrations were carried out for 1 ns for each system with constant number of particles, volume, and temperature (NVT; with modified Berendsen thermostat with velocity rescaling (Bussi et al., 2007) at 310 K and a 0.1 ps time step, Particle Mesh Ewald coulomb type (Kawata & Nagashima, 2001) for long-range electrostatics with Fourier spacing 0.16); and constant number of particles, pressure, and temperature (NPT; Parrinello–Rahman pressure coupling (Martoňák et al., 2003) at 1 bar with a compressibility of 4.5 × 105 bar1 and a 2 ps time constant). Finally, the equilibrated systems were subjected to 100 ns MD simulation with time-steps of 2 fs. Bond-lengths were constrained using the Linear Constraint Solver (LINCS) algorithm (Hess et al., 1997). Quality of receptor-ligand complexes were analyzed using in-built Gromacs utilities.

MM-PBSA binding free energy calculations

The g_mmpbsa tool (Baker et al., 2001; Kumari et al., 2014) was employed for calculating the binding free energy of protein–ligand complexes from MD trajectories. Post simulation MM-PBSA (Molecular Mechanics Poisson − Boltzmann Surface Area) binding free energy calculations have been used extensively for screening inhibitors (Khan et al., 2020) and was found to correlate reasonably well with experimental results (Ekhteiari Salmas et al., 2017; Rajkumari et al., 2018). This module estimates Gibb’s free energy of binding using the MM-PBSA method as described by the equations below:

Where, ΔGcomplex, ΔGprotein and ΔGligand signify the total free energy of the protein − ligand complex and total free energies of the isolated protein and ligand in solvent, respectively.

ΔGx is the free energy for each individual entity. TΔS refers to the entropic contribution to the free energy in vacuum where T and S denote the temperature and entropy, respectively. ΔGsolv is the free energy of solvation where Gpolar and Gnonpolar are the electrostatic and non-electrostatic contributions to the solvation free energy. The vacuum potential energy, ΔEMM, includes the energy of both bonded as well as non-bonded interactions. The g_mmpbsa tool also allows decomposing the total binding free energy into the contribution made by each residue.

Pharmacokinetic parameters prediction

Pharmacokinetic properties assist in the early stages of drug discovery by identifying safety and effectiveness of the compounds. Predicted values of toxicity and some pharmacokinetic parameters were obtained from IMPPAT database which were calculated with admetSAR server (Cheng et al., 2012).

Results

Modeling of structural proteins of SARS-CoV-2

The 3D model for SARS-CoV-2 E-protein, generated from I-TASSER server with a C-score of −0.75, was further refined and energy minimized using 3D refine. The best model was chosen based on the lowest 3Drefine (4938.24) and RWplus (-10970.54) scores, which indicated a better and more physically realistic model. The model, whichcontained 65.9% helix, 0% sheet and 34.1% other secondary structure elements, was found to have 92.9% residues in favored region, 5.4% in the allowed region and 1.8% in the outliner region of Ramachandran plot (Supplementary material Figure S2). The3D model for M-protein obtained from I-TASSER with a C-score=-3.35 and was similarly refined and energy minimized, and the best model chosen based on the lowest 3Drefine (14559.9) and RWplus (-43870.59) scores. The 3D model for SARS-CoV-2 M protein (Figure 1B) contained 43.7% helix, 9% sheet and 47.3% other secondary structure elements, and was found to have 86.4% residues in favored, 10.9% in allowed and 2.7% residues in the outlier regions of Ramachandran plot respectively (Supplementary material Figure S2).

Figure 1.

Modeling of SARS-CoV-2 E and M proteins: A) 3D representation of a pentameric form of E-protein, with the N-terminal, transmembrane, and C-terminal domains of each monomer highlighted in blue, green and orange respectively and secondary structure representation is shown on right panel. B) 3D representation of the M-protein, with the N-terminal, transmembrane and C-terminal domains highlighted in blue, green and orange respectively. The C-terminal conserved region is highlighted in red and secondary structure representation is shown on right panel.

Molecular docking

For both E and M proteins, the top five best scoring compounds in docking studies were selected for further analysis. Details of the docking scores and interactions of each compound with the respective proteins are discussed below.

Interaction with SARS-CoV-2 E

All five best scoring compounds in terms of binding energy (Table 1) displayed strong interaction with the transmembrane and C-terminal domains of the E-protein. Almost all compounds showed hydrophobic interaction with residues from all five monomeric chains in the E-protein pentamer. Nimbolin A, which had the highest binding energy of −11.2 kcal/mol, formed 16 hydrophobic interactions. The other compounds: Nimocin, 7-Deacetyl-7-benzoylgedunin, 24-Methylenecycloartanol and Cycloeucalenone -interacted with the E-protein pentamer with binding energies of −11.0 kcal/mol, −10.9 kcal/mol, −10.8 kcal/mol and −10.5 kcal/mol respectively, with 10 or more hydrophobic interactions in each case (Table 1, Figures 2A and 3A).

Table 1.

Docking results of top five ligands with the E-protein of SARS-CoV-2.

Compound	Chemical Structure	Binding Energy (kcal/mol)	H-bonds	Hydrophobic Interactions (Chain ID)
Nimbolin A		−11.2	NA	Leu 18 (A), Leu 19 (A), Ala 22 (A), Leu 19 (B), Leu 18 (B), Ala 22 (C), Phe 23 (C), Phe 26 (C), Asn 15 (D), Leu 18 (D), Leu 19 (D), Val 25 (D), Phe 26 (D), Leu 18 (E), Leu 19 (E), Ala 22 (E)
Nimocin		−11.0	NA	Leu 18 (B), Leu 19 (B), Leu 19 (C), Ala 22 (C), Phe 23 (C), Phe 26 (C), Leu 19 (D), Ala 22 (D),Val 25 (D), Phe 26 (D), Ala 22 (E), Leu 18 (E), Phe 26 (E)
7-Deacetyl-7-Benzoylgedunin		−10.9	NA	Leu 19 (A), Leu 18 (B), Leu 19 (B), Ala 22 (B), Leu 19 (C), Ala 22 (C), Phe 26 (C), Ala 15 (D), Leu 18 (D), Leu 19 (D), Ala 22 (D), Phe 26 (D), Leu 18 (E), Ala 22 (E)
24-Methylenecycloartanol		−10.8	NA	Leu 18 (B), Leu 19 (B), Ala 22 (C), Phe 26 (C),Leu 18 (D), Ala 22 (D), Val 25 (D), Phe 26 (D),Val 29 (D), Phe 26 (E)
Cycloeucalenone		−10.5	NA	Ala 22 (A), Phe 23 (A), Phe 26 (A), Ala 22 (B), Ala 22 (C), Phe 26 (C), Ala15 (D), Leu 18 (D), Leu 19 (D), Ala 22 (D), Leu 18 (E), Ala22 (E), Phe 26 (E)

Figure 2.

2D interaction diagram of five compounds with the highest binding energies against SARS-CoV-2 E and M proteins. The diagrams were generated with LigPlot++.

Figure 3.

3D interaction diagram of five compounds with the highest binding energies against SARS-CoV-2 E and M proteins. The ligands are shown in violet surface representation.

Interaction with SARS-CoV-2 M

The best scoring compounds in terms of binding energy (Table 2) formed primarily hydrophobic interactions with the transmembrane and C-terminal domains. Nimocin had the highest binding energy of −10.2 kcal/mol, with 12 hydrophobic interactions. 24-Methylenecycloartan-3-one, Phytosterol, Beta-Amyrin and Nimbolin A had binding energies of −10 kcal/mol, −9.5 kcal/mol,-9.4 kcal/mol and −9.4 kcal/mol respectively with 10, 9, 8 and 17 hydrophobic interactions. Three of the compounds - 24-Methylenecycloartan-3-one, Phytosterol and Beta-Amyrin - formed hydrophobic interactions with Ser 111 and Phe 112 from the conserved region in the C-terminal domain (Table 2, Figures 2B and 3B).

Table 2.

Docking results of top five ligands with the M-protein of SARS-CoV-2.

Compound	Chemical Structure	Binding Energy (kcal/mol)	H-bonds	Hydrophobic Interactions
Nimocin		−10.2	NA	Gln 36, Ala 38, Asn 41, Arg 42, Phe 45, Gly 78, Ala 81, Ile 201, Gly 202, Asn 203, Tyr 204, Asp 209
24-Methylenecycloartan-3-one		−10	NA	Ile 97, Phe 100,Arg 101,Ala 104,Ser 111, Phe 112, Thr 127, Ile 128, Thr 130, Arg 158
Phytosterol		−9.5	NA	Ile 97, Phe 100, Arg 101, Ala 104, Ser 111, Thr 127, Ile 128, Thr 130, Arg 158.
Beta-Amyrin		−9.4	NA	Ile 97, Phe 100, Arg 101, Ala 104, Ser 111, Phe 112, Ile 128, Thr 130
Nimbolin A		−9.4	NA	Trp 31, Ile 32, Leu 35, Gln 36, Phe 37, Ala 38, Asn 41, Arg 42, Phe 45, Ala 81, Ile 82, Ala 85, Asn 203, Tyr 204, Leu 206, Thr 208, Asp 209.

Molecular dynamics simulation analysis

Molecular dynamics simulations upto 100 ns were performed to analyze the stability of the protein-ligand complexes as described before. The stability of the complexes was assessed by Root Mean Square Deviation (RMSD) of protein backbone and ligand. In case of E-protein complexes, the RMSD of the protein backbone was consistent after 50 ns for all complexes. The overall ligand RMSD was also stable in the final 50 ns of simulations for all ligands except Cycloeucalenone, which showed comparatively higher fluctuations (Figure 4A). In case of M-protein complexes, the protein backbone RMSD stabilized after 80 ns. The overall ligand RMSD was stable throughout simulation time in case of all ligands except Phytosterol, where comparatively higher fluctuations were observed initially, which was stabilized after 80 ns (Figure 4B).

Figure 4.

Molecular Dynamics (MD) simulation results of Protein-ligand complexes: RMSD of protein backbone over time and RMSD of ligands over time for A) E-protein complexes and B) M-protein complexes.

MM/PBSA

In order to further quantify the binding affinity between the structural proteins and Neem compounds, the binding free energy and energy contribution of residues were calculated. The energy components EMM, Gpolar, and Gnonpolar of each complex were calculated from 201 snapshots that were extracted at every 0.1 ns from the production trajectories from the final 20 ns. The binding free energy between the protein and inhibitors were decomposed into the contribution of each residue using MM-PBSA approach. Only the residues with contributions greater than −1.5 kJ/mol were shown.

E-protein

Among the 5 complexes, 7-Deacetyl-7-Benzoylgedunin and Nimbolin A showed the lowest binding free energy values of-297.87 kJ/mol and −297.49 kJ/mol respectively. Nimocin, 24-Methylenecycloartanol and Cycloeucalenone displayed binding free energies of −277.20 kJ/mol, −199.84 kJ/mol and −193.81 kJ/mol respectively (Table 3). Among all observed interactions, the contribution of van der Waals and polar solvation interactions to the binding free energies were more than that of others. Leu 18, Leu 19, Leu 21, Ala 22, Val 25 and Phe 26, were found to be common among all interacting residues in the complexes, and contributed strongly to the binding energy values (Figure 5).

Table 3.

Binding free energy (MM/PBSA) and their components (kJ/mol) of E-protein complexes.

Ligands	Binding Energy kJ/mol	Van der Waal Energy	Electrostatic Energy	Polar Solvation Energy	SASA Energy
Nimbolin A	−297.49 ± 14.05	−329.32 ± 12.37	−5.24 ± 3.53	67.80 ± 8.19	−30.73 ± 1.42
Nimocin	−277.20 ± 9.90	−294.12 ± 10.33	−1.25 ± 1.27	43.28 ± 3.64	−25.10 ± 1.10
7-Deacetyl-7-Benzoylgedunin	−297.87 ± 9.72	−323.11 ± 9.94	−4.18 ± 1.89	54.94 ± 3.76	−25.51 ± 1.10
24-Methylenecycloartanol	−199.84 ± 9.70	−215.18 ± 11.18	−15.12 ± 3.90	52.86 ± 5.11	−22.40 ± 1.38
Cycloeucalenone	−193.81 ± 9.16	−212.18 ± 10.06	0.35 ± 3.03	39.60 ± 8.44	−21.57 ± 1.28

Figure 5.

MM/PBSA per residue energy contribution to the binding energy (only residues scoring ≥1.5 kJ/mol are shown) for individual compounds A) Nimbolin A, B) Nimocin, C) 7-Deacetyl-7-Benzoylgedunin, D) 24-Methylenecycloartanol and E) Cycloeucalenone with E-protein.

M-Protein

The Nimbolin A- M protein complex showed the lowest binding free energy value (-188.99 kJ/mol), followed by 24-Methylenecycloartan-3-one (-173.33 kJ/mol), Beta-Amyrin (-161.55 kJ/mol), Nimocin (-155.77 kJ/mol), and Phytosterol (-123.70 kJ/mol) (Table 4). In this case also, the contribution of van der Waals and polar solvation energies were the highest. Common residues contributing strongly to the binding energy in the complexes of M with Nimbolin A and Nimocin included Met 84, Ala 85, Arg 200 and Tyr 204. Interaction of 24-Methylenecycloartan-3-one, Beta-Amyrin and Phytosterol was with an alternate site of M protein. Per residue contribution (≥1.5 kJ/mol) to the binding energy is shown in Figure 6. Residues in the conserved C-terminal region (SMWSFNPETNIL) of M protein (Arndt et al., 2010), like Met 109, Trp 110, Phe 112, Pro 114 and Ile 118 were some of the highest contributors to the binding energy values.

Table 4.

Binding free energy (MM/PBSA) and their components (kJ/mol) of M-protein complexes.

Ligands	Binding Energy kJ/mol	Van der Waal Energy	Electrostatic Energy	Polar Solvation Energy	SASA Energy
Nimocin	−155.77 ± 15.08	−207.78 ± 17.56	−0.53 ± 4.04	71.52 ± 13.23	−18.98 ± 1.65
24-Methylenecycloartan-3-one	−173.33 ± 11.69	−230.28 ± 12.36	−0.05 ± 2.94	79.45 ± 10.46	−22.45 ± 1.21
Phytosterol	−123.70 ± 13.22	−138.88 ± 13.47	−7.00 ± 6.05	38.32 ± 10.64	−16.14 ± 2.09
Beta-Amyrin	−161.55 ± 11.83	−192.44 ± 12.05	−1.98 ± 6.45	50.70 ± 12.87	−17.83 ± 1.56
Nimbolin A	−188.99 ± 19.16	−235.41 ± 15.09	−10.55 ± 9.22	79.76 ± 16.65	−22.78 ± 1.86

Figure 6.

MM/PBSA per residue energy contribution to the binding energy (only residues scoring ≥1.5 kJ/mol are shown) for individual compounds A) Nimocin A, B) 24-Methylenecycloartan-3-one, C) Phytosterol, D) Beta-Amyrin and E) Nimbolin A with M-protein.

Pharmacokinetic properties of the top compounds

The pharmacokinetic properties from filtering analyses (Supplementary material Table S3) suggested that all top compounds have high probability of absorption and blood brain barrier permeability. All compounds displayed low CYP450 inhibitory promiscuity, and high probability of being non-carcinogenic and non-toxic.

Discussion

The novel coronavirus SARS-CoV-2 constitutes an unprecedented physical, physiological and financial threat to the global population. Till now, this virus has caused ∼3.94 million infections and ∼275,000 deaths all over the world. All movement and activities around the world have come to a grinding halt in an effort to contain the virus. As of now, the treatment is primarily symptomatic; but clinical trials with an investigational antiviral drug Remdesivir has shown reduced recovery time in COVID-19 patients (Hendaus, 2020). Although FDA has issued an emergency use authorization for Remdesivir (Barrett, 2020; Wang et al., 2020), the requirement of more specific and effective antiviral drugs against SARS-CoV-2 is evident. Multiple vaccine candidates are being tested, however the licensing and mass production of an effective vaccine is expected to require a minimum period of 1–1.5 years. Prior to the emergence of SARS-CoV-2, the world has witnessed the emergence of SARS-CoV and MERS-CoV. This indicates that the chance of future emergence of another novel human pathogen from the coronavirus family, with the capability to cause a pandemic, cannot be dismissed (Menachery et al., 2015). Thus, identification of crucial steps in the virus life cycle, and developing specific inhibitors targeted towards these processes, appears to be a viable strategy.

Computational methods tender a fast and cost-efficient approach to design potential inhibitors. We have utilized virtual screening and molecular dynamics simulation to identify small molecules that bind to the structural proteins M and E of SARS-CoV-2 with. M and E are necessary and sufficient to form particles of other β-coronaviruses (Mortola & Roy, 2004), therefore they are thought to have a central role in virus assembly. Interaction of the M-protein with all other major structural proteins like N, S and E is central to virus assembly. The contribution of the ion-channel protein E in membrane interaction and bending is deemed necessary for the progression of virus budding (Schoeman & Fielding, 2019). Previous studies have suggested that the ion-channel activity of E is required for virus propagation; thus inhibitors targeting the E-protein can help to prevent virus production (Boopathi et al., 2020; Nieto-Torres et al., 2014). The binding of the anti-Influenza drug amantadine to ion-channels formed by E-protein has also been reported, although the effect was found to be moderate (Torres et al., 2007). Our docking analysis with E showed that the top five compounds interacted within the central channel formed by pentamer. The two critical interactions made by these compounds that can be linked to the biological activity of the protein, are with transmembrane (TM) domain residues Asp 15 (N15) and Val 25 (V25). Previous studies have shown that mutations at these residues (N15A & V25F) inhibit the ion channel activity, and hamper the oligomerization of E-protein of SARS CoV (Schoeman & Fielding, 2019; Torres et al., 2006, 2007; Verdia-Baguena et al., 2012). These two residues are conserved in multiple isolates of SARS-CoV-2, as well as in major β-coronavirus strains (Supplementary material Figure S4). In our study, V25 was found to be a common residue in all complexes that contributed highly to the MM/PBSA binding free energy. N15 was involved in interactions with 7-Deacetyl-7-Benzoylgedunin and 24-Methylenecycloartanol, and also contributed highly to the binding free energy. Post simulation MM/PBSA analysis showed that residues within the pentameric channel such as Leu 18, Leu 19, Leu 21, Ala 22, Val 25 and Phe 26 were also strong contributors to binding free energy, indicating robust interaction with the screened compounds. Our studies indicate that these compounds may have the ability to inhibit ion channel activity of SARS-CoV-2 Envelope protein, thus hampering virus assembly.

In case of M-protein, docking analysis with Neem compounds revealed two possible binding sites. Nimocin and Nimbolin A interacted primarily with the transmembrane domain, with contributions from the last 20 residues at the C-terminus. The C-terminal tail, including residues Lys 199, Gly 201, Tyr 203 and Arg 204, are required for localization of the M-protein of MERS-CoV in the Trans-Golgi Network (TGN) (Perrier et al., 2019). Localization of M to this site is crucial for virion assembly. Our post simulation MM-PBSA analysis showed that Gly 202 and Tyr 204 are among the strongest contributors to the binding free energy from the C-terminal end of M-protein. The residues GNY (202–204) are conserved in members of the β-coronavirus family (Supplementary material Figure S5).

24-Methylenecycloartan-3-one, Phytosterol and Beta-Amyrin were found to interact with the C-terminal region, primarily with the conserved domain CD (SMWSFNPETNIL),which is essential for virus assembly (Arndt et al., 2010). Thus, the compounds identified in this study may hamper virus assembly by interacting with crucial regions of the M-protein involved in TGN localization and protein-protein interaction.

In conclusion, we have identified a few possible inhibitors of the E and M proteins of SARS-CoV-2 using molecular docking, MD simulation and binding free energy calculations. These compounds, derived from Neem, displayed stable binding and interactions with crucial regions of E and M required for assembly; and were predicted to have good pharmacokinetic properties. Nimbolin A showed the strongest binding free energy with both E and M proteins. Other compounds: Nimocin and Cycloartanols (24-Methylenecycloartanol and 24-Methylenecycloartan-3-one) were also common ligands, binding strongly to both proteins. A recent study has shown ameliorating effect of Neem extract on propagation and pathophysiology of another member of the coronavirus family (Sarkar et al., 2020), reinforcing the potential of these compounds as potential therapeutic options. Experimental validation and optimization of these natural compounds might add value to the development of specific therapeutics against SARS-CoV-2.