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.
E-mail addresses: nabajyotibaildya@gmail.com (N. Baildya), asoke@klyuniv.ac.in (A.P. Chattopadhyay). 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,0 0 0 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  this plant for in-silico screening against 3CLpro of SARS-CoV-2 by docking and molecular dynamics simulation.

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 −1 nm −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 ( G polar ), non-polar free energy ( G non-polar ), average MM-energy ( E MM ), 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   [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 isosurfaces 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.
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 .

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). 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.

Toxicity prediction
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 .
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.
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. 4 b 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. 4 c 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 .
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 .
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. 7 a 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. 7 b. The figure clearly shows that the docked structure has lesser surface area compared to the undocked structure, which reveals destabilization of the former. 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. 4 a and in Fig. 8 . The structure continues to change as reflected in Fig. 4 a 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. 4 a). 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.

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.

Declaration of Competing Interest
The authors declare no conflicting interest in the present work.