Inhibitory effect of compounds extracted from Monochoria hastata (L.) Solms on SARS-CoV-2 main protease: An insight from molecular docking and MD-simulation studies.

Using molecular docking and other studies, 20 compounds extracted from Monochoria hastata (L.) Solms were screened, and their inhibitory efficiency examined against main protease (3CLpro) of SARS CoV-2. All the compounds were found to binding with 3CLpro through van der Waals and electrostatic forces of attractions. Among them, Azelaic dihydrazide (ADZ) was found to have the highest docking score. 3CLpro-ADZ complex was studied by MD simulation. ADZ was found to disrupt the structure of 3CLpro after 2 ns. RMSD and RMSF analysis along with sequence and binding energy analysis suggest that ADZ can be a potential drug against SARS CoV-2.

Introduction

The devastating effect of novel coronavirus disease, COVID-19, progressively affects global health and economy. Such an emergency throws a challenge to the entire scientific community to work together and come up with an effective antidote against this deadly virus (SARS-CoV-2). Numerous attempts have been made to find useful drugs or vaccines to combat this virus [1], [2], [3], [4]. But till date, no specific drug has been found. Drugs like chloroquine, hydroxychloroquine, favipiravir, remdesivir, azithromycine etc [5], [6], [7], [8], [9]., extracts from different plants [10], [11], [12], normally used in other diseases, are under clinical trial against this virus, but they are being used with caution.

Coronavirus was first found in 1960 [13]. It returned in new avatars in 2002 and 2012 to haunt us [14], [15], [16], [17]. However, the novel coronavirus viz. SARS-CoV-2 or COVID-19 has created a pandemic because of its higher transmission rate (R0) compared to previous cases. The SARS-CoV-2 belongs to the β-coronavirus group, and its genome is quite similar to that of previous SARS-CoV. Although exact origin of this virus is not clear [18], it is transmitted rapidly among human beings by community transmission [19]. The outer surface of this virus contains spike glycoprotein, an enveloped membrane and nucleocapsid protein containing a positive sense single strand of RNA. It mainly binds with human angiotensin-converting enzyme-2 (ACE-2) receptors using its spike glycoprotein. Along with spike proteins it contains the cellular serine protease TMPRSS2, cysteine protease, main protease (3CLpro), papain-like protease (PLpro) and different non-structural proteins (nsp) which are basically used for its replication. Therefore, drugs that effectively block one or more of these proteins or the ACE-2 enzyme involved in binding to host cell or those active in replication process are much focused [20], [21], [22].

From ancient times, herbs have been used as effective medicinal resource. There are approximately 25,000 plants in India which are used in traditional and folk medicine [23]. An aquatic medicinal plant of India, Monochoria hastata (L.) Solms, belonging to the family Pontederiaceae, is being used as remedy in Indian villages for various disorders. The extracts obtained from its leaf show antibacterial and anti-viral effects. It contains different types of alkaloids, phenols, terpenoids, flavonoids, glycosides, monoacyl glycerol etc. [24]. 20 different compounds were selected from extracts of this plant for in-silico screening against 3CLpro of SARS-CoV-2 by docking and molecular dynamics simulation.

Methodology

Virtual screening of Monochoria hastata (L.) Solms extracts with 3CLpro

The crystal structure of 3CLpro of SARS-COV-2 (PDB ID: 6LU7) was obtained from the Protein Data Bank. The sdf files of the drugs were downloaded from PubChem (National Library of Medicine) and were converted to pdb format to check ligand-protein binding interactions. UCSF Chimera [25] package was used to clean the structure of 3CLpro. Autodock Vina [26] was used to investigate docking between 3CLpro and the drug molecules and necessary files were prepared by using Autodock Tools.

Molecular dynamics simulation studies of ADZ docked 3CLpro complex

Azelaic dihydrazide (ADZ) was found to have the highest docking score with 3CLpro. The minimum energy docked configuration of 3CLpro-ADZ complex was selected for molecular dynamics (MD) simulation studies. GROMACS (version 5.1) [27] was used to perform the MD simulation using CHARMM36-mar2019 force-field [28] with TIP3P solvation model [29]. Ligand parameters and topology files were generated with the help of CHARMM General Force Field server. Periodic boundary conditions were applied to maintain a cubical box in which the 3CLpro-ADZ complex was at least 1 nm from the edges for successive imaging of the complex. To maintain electroneutrality, four Na+ ions were added to the system. Energy minimization was carried out by the steepest descent algorithm followed by conjugate gradient protocol, until the system reached a minimum force of 10 kJmol−1nm−1. For equilibration of the system, isochoric-isothermal (NVT) ensemble was used at 300 K for 100 ps. The time step was set to 2 fs. This was followed by equilibration in an isothermal-isobaric (NPT) ensemble for 100 ps at 300 K. Modified Berendsen thermostat was used for the NPT ensemble. Here also the time step was set to 2 fs. For both the NVT and NPT equilibration, 1 nm cut-off was kept for van der Waals and electrostatic interactions. For the calculation of long range interactions, smooth particle mesh Ewald (PME) [30] method was used. Finally, MD simulation was performed for 10 ns using the equilibrated ensembles with the same cut off as used before. A modified Berendsen thermostat and a Parinello-Rahman barostat at 300 K temperature and 1 bar pressure were used. Snapshots of the trajectory were taken after each nanosecond of simulation.

Analysis of MD simulations

The trjconv tool was used for analysis of structural trajectories of the 3CLpro-ADZ complex during simulation. Drug, water, counter ions and 3CLpro were re-centered within the cubical box using the same tool. Variation of root mean square deviation (RMSD), root mean square fluctuations (RMSF), radius of gyration (Rg) and solvent accessible surface area (SASA) were plotted using xmgrace plotting tool. g_mmpbsa tool [31,32] was used to calculate average binding free energy, polar free energy (ΔGpolar), non-polar free energy (ΔGnon-polar), average MM-energy (ΔEMM), van der Waals energy, SASA energy and electrostatic energy. The change of binding energy with time was plotted using Origin 8.0.

Results and discussions

20 compounds, extracted from Monochoria hastata (L.) Solms, were screened against SARS-CoV-2 main protease (3CLpro). Binding affinity of these compounds was found to be in the range of −3.6 to −5.8 Kcal/mol as shown in Fig. 1 . The highest binding affinity of −5.8 Kcal/mol was observed for Azelaic dihydrazide (ADZ) against 3CLpro. Previous studies showed that Azelaic Acid and its hydrazide derivation have high binding affinity with DNA Polymeraze I (2KFN) due to enhance H-bonding showing its potential inhibitory activity against 2KFN [33]. This observation is in accordance with our simulated results. It has been also reported that Azalaic Acid showed potential anti-viral inhibitory activity [34]. To analyze the stereochemistry of the drug-protein complex we have plotted contour diagram of drug with neighbouring residue of the 3CLpro as shown in Table S2. van der Waal donor-acceptor iso-surfaces of drug-protein complex clearly indicates electrostatic interaction is exists between them while the 2D contour plot demonstrates van der Waal and the number of H-bonding interaction between drug and protein.

Fig. 1

Binding affinity (Kcal/mol) of Monochoria hastata (L.) Solms extracts.

To predict toxicity of these compounds, ADME toxicity analysis was performed using pkCSM online server [35]. All the estimated pharmacokinetic properties are listed in Table 1 .

Table 1

Toxicity parameters for 20 Monochoria hastata (L.) Solms compounds studied.

Compound	AMES toxicity	Max. tolerated dose (human)	hERG I inhibitor	hERG II inhibitor	Oral Rat Acute Toxicity (LD50) (mol/kg)	Oral Rat Chronic Toxicity (LOAEL)(log mg/kg_bw/day)	Hepato-toxicity	Skin Sensiti-sation	T.Pyriformis toxicity (log ug/L)	Minnow toxicity (log mM)
1-Methyl −2-pyrrolidinone	No	1.077	No	No	2.153	1.255	No	Yes	−0.887	2.861
Azelaic dihydrazide	No	0.907	No	No	1.762	2.87	No	Yes	0.143	2.839
4,5,6,7-Tetrahydro-1-benzothiophene-3-carbonitrile	No	0.438	No	No	2.049	1.882	No	Yes	1.431	1.005
Phytol	No	0.05	No	Yes	1.607	1.043	No	Yes	1.884	−1.504
2(4H)-Benzofuranone, 5,6,7,7a-tetrahydro	No	0.783	No	No	1.94	2.412	No	No	0.037	2.224
2-Pentadecanone, 6,10,14-trimethyl	No	0.244	No	No	1.532	1.094	No	Yes	2.14	−1.478
2H-Pyran-2-one, tetrahydro-6-octyl-	No	0.276	No	No	1.826	2.429	No	Yes	1.491	0.216
1,2,3-Propanetriol, 1-acetate	Yes	1.506	No	No	1.708	2.789	No	No	−0.616	3.366
Z,E-3,13-Octadecadien-1-ol	No	−0.207	No	Yes	1.578	1.106	No	Yes	1.703	−1.156
Diethyl phthalate	No	1.37	No	No	2.09	2.648	No	No	0.656	1.441
Phthalic acid, 6-ethyl-3-octyl isobutyl	No	1.324	No	No	1.712	2.42	No	No	1.468	−1.526
Dodecanoic acid	No	−0.34	No	No	1.511	2.89	No	Yes	0.954	−0.084
Tetradecanoic acid	No	−0.559	No	No	1.477	3.034	No	Yes	0.978	−0.601
Methyl 9- methyltetradecanoate	No	0.337	No	No	1.653	2.832	No	Yes	2.268	−1.044
Pentadecanoic acid	No	−0.642	No	No	1.458	3.107	No	Yes	0.922	−0.842
Dodecanoic acid, 10-methyl-, methyl ester	No	0.421	No	No	1.697	2.679	No	Yes	2.166	−0.558
9-Hexadecenoic acid, methyl ester, (Z)	No	0.125	No	No	1.608	2.926	No	Yes	2.005	−1.245
Tridecanoic acid, methyl ester	No	0.301	No	No	1.645	2.779	No	Yes	2.191	−0.638
2-Hexyldecanoic acid	No	−0.574	No	No	1.564	3.068	No	Yes	0.457	−0.923
Oleic Acid	No	−0.81	No	No	1.417	3.259	No	Yes	0.676	1.438

ADME toxicity prediction

Pharmacokinetic analysis reveals that the compounds have low blood-brain barrier (BBB) permeability, with values ranging between −0.016 and 0.816. The drugs were found to have no effect on CYP2C9, CYP2D6 and CYP3A4 inhibitors. Skin permeability values were found between −1.384 and −4.376 for the drugs. These drugs also do not inhibit renal OCT2 substrate as shown Table S3 (in Supporting Information).

Toxicity prediction

Table 1 represents the toxicity level of the studied compounds. None of the compounds have AMES toxicity and they do not have hERGI inhibitory activity. Estimated oral rat toxicity values were found to lie between 1.417 and 2.153 mol/kg. Estimates of chronic oral rat toxicity (LOAEL) values are in between 1.094 and 3.259 (log mg/kg_bw/day). None of the compounds show hepatotoxicity, though most of the compounds show skin sensitization.

From the above analysis it is clear that all 20 compounds under consideration could serve as drugs for SARS CoV-2. The bound structures of some of the drugs with 3CLpro are shown in Fig. 2 .

Fig. 2

Docked structures of some of the drugs considered with 3CLpro.

From the docking results it is clear that Azelaic dihydrazide (ADZ) showed the highest binding affinity to the main protease (3CLpro) of SAR-CoV-2. The docked structure along with a 2D contour plot of ADZ with 3CLpro are shown in Fig. 3 considering the electrostatic and van der Waals interactions.

Fig. 3

Docked structure (left) and 2D contour plot of 3CLpro-ADZ complex.

From the above figure it is clear that ADZ binds strongly with 3CLpro via formation of five strong H-bonds ranging from 2.15 to 2.97 Å along with donor-acceptor interactions ranging from 2.37 to 2.62 Å which are very short with respect to other systems.

The lower value of inhibition constant (Ki) between ADZ and 3CLpro indicates stronger binding affinity between them. The neighboring residues are shown in both the views of the docked structure. It was decided to carry out molecular dynamics (MD) simulation on the 3CLpro-ADZ complex to investigate the effect of the drug on 3CLpro in greater detail. Fig. 4 a represents the RMSD plot of undocked 3CLpro and ADZ, and of 3CLpro-ADZ complex. Fig. 4b represents the RMS fluctuation plot of docked and undocked 3CLpro. The RMSD plot for the 3CLpro showed much lesser fluctuation as compared to 3CLpro-ADZ complex after 2 ns which suggests decreased stability of the docked structure vis-a-vis the undocked one. It is clear from the above analysis that ADZ creates a huge structural disruption in 3CLpro. This is further supported by the RMSF plot which reveals that ADZ docked 3CLpro showed less fluctuations in the residues compare to the undocked 3CLpro. Moreover, RMSF plots Fig. 4c accounts less fluctuation of the neighboring catalytic residues (HIS41, TYR54, GLY143, SER144, CYS145, ASP187) to the drug molecule in the binding pocket of 3CLpro. The structural changes of ADZ docked 3CLpro is also corroborated by the sequence analysis of the docked and undocked 3CLpro as shown in Fig. 5 .

Fig. 4

RMSD (a) and RMSF (b) plots of docked and undocked 3CLpro, the color codes have been mentioned in the figure.

Fig. 5

Sequence analysis of 3CLpro after MD in undocked and docked structure.

Sequence analysis revealed that the residues after 202 had completely changed their positions indicating huge conformational change within the 3CLpro. The binding free energy of polar and non-polar groups of the docked structure with respect to time are shown in Fig. 6 c and 6d, respectively. With the passage of time, polar binding free energy (H-bonding and electrostatic) decreased indicating much stronger binding of the drug ADZ in the 3CLpro cavity. It is to be noted that the stronger binding of ADZ with 3CLpro explains the huge structural change and protein chain sequence changes in 3CLpro. Magnitudes of some major interactions in the 3CLpro-ADZ complex are listed in Table 2 .

Fig. 6

Variation of different binding energy components: (a) Binding energy, (b)ΔEMM, (c) ΔGpolar and (d) ΔGnon-polar with time.

Table 2

Average values of various binding free energy components of 3CLpro-ADZ complex.

System	Binding free energy (kJ/mol)	van der Waal energy, ΔEvdW(kJ/mol)	Electrostatic energy, ΔEelec (kJ/mol)	Polar solvation energy, ΔGpolar (kJ/mol)	SASA energy (kJ/mol)
3CLpro+ADZ	- 49.1 ± 18.8	- 82.7 ± 17.4	- 91.3 ± 20.0	137.3 ± 26.8	- 12.3 ± 1.5

The compactness of a system is measured by the radius of gyration (Rg) analysis. Variation of Rg with time is shown in Fig. 7 a for undocked and docked 3CLpro with ADZ. Form Fig. 7a it is clear that the docked structure is more compact as compared to the undocked one. The solvent accessible surface area (SASA) of docked and undocked 3CLpro are plotted against time in Fig. 7b. The figure clearly shows that the docked structure has lesser surface area compared to the undocked structure, which reveals destabilization of the former.

Fig. 7

Radius of gyration (a) and SASA (b) plot of docked and undocked 3CLpro.

Fig. 8 represents structural changes during MD simulation. From the above analysis it is clear that ADZ creates a large impact to destabilize the structure of 3CLpro. This is also clearly reflected in Fig. 8. Snapshots are given in the figure after every 2 ns starting from the beginning. The RMSD plot and Fig. 8 are in agreement with each other. At 2 ns, a substential conformational change is indicated in Fig. 4a and in Fig. 8. The structure continues to change as reflected in Fig. 4a and b, and also in Fig. 8. The docked structure does not seem to attain equilibrium up to 10 ns, i.e. till the end of simulation, whereas the undocked structure of the protease reaches equilibrium after 3 ns, as revealed from the RMSD plot (Fig. 4a). It is clear, therefore, that the molecule ADZ from M. hastata destabilized the structure of 3CLpro on binding with it. This is expected to affect the activity of SARS-CoV-2 to noticeable extent. 3CLpro plays an essential role in processing the ploy proteins translated from viral RNA and thus aids in replication of the coronavirus. Inactivation of 3CLpro is therefore expected to lead to inhibiting viral replication of SARS-CoV-2.

Fig. 8

Conformational changes in 3CLpro upon binding with ADZ during MD-simulation.

Conclusions

20 compounds from the extract of the ethnomedical aquatic herb M. hastata (L) Solms were screened against the main protease (3CLpro) of the novel coronavirus SARS-CoV-2. All the compounds were found to have good docking affinity with 3CLpro, as also desirable properties of pharmacokinetics, as obtained from ADME and related analysis. Among all the compounds, Azelaic dihydrazide (ADZ) was found to have the best docking score with 3CLpro. This docked structure of 3CLpro-ADZ was then subject to molecular dynamics (MD) simulation for 10 ns. Analysis of the simulation data viz. RMSD, RMSF, Rg and SASA plots, as well as from the structures of the 3CLpro-ADZ complex during simulation, the following conclusion can be drawn: ADZ creates a large conformational change on 3CLpro after binding or docking with it. This may cause noticeable change of activity of 3CLpro, and thereby inhibiting the replication process SARS-CoV-2. This calls for further investigations on effect of same and other compounds of M. Hastate (L.) Solms on SARS-CoV-2 and ACE2 of host cell.

Data availability

Data is available upon request to the corresponding author.

CRediT authorship contribution statement

Nabajyoti Baildya: Conceptualization, Writing – original draft, Methodology, Supervision. Narendra Nath Ghosh: Conceptualization, Writing – original draft, Methodology. Asoke P. Chattopadhyay: Conceptualization, Writing – original draft, Supervision. Vivekananda Mandal: Conceptualization, Writing – original draft, Methodology. Sourav Majumdar: Writing – original draft, Methodology. Delwar Ansary: Writing – original draft, Methodology, Data curation. Md Muttakin Sarkar: Writing – original draft, Methodology, Data curation.

Declaration of Competing Interest

The authors declare no conflicting interest in the present work.