key: cord-0888564-cm4jc2vk authors: Sobhia, M. Elizabeth; Ghosh, Ketan; Sivangula, Srikanth; Kumar, Siva; Singh, Harmanpreet title: Identification of potential SARS-CoV-2 M(pro) inhibitors integrating molecular docking and water thermodynamics date: 2021-01-08 journal: Journal of biomolecular structure & dynamics DOI: 10.1080/07391102.2020.1867642 sha: cc1c03a8f361b6c4405dde255d187d0d79be58f0 doc_id: 888564 cord_uid: cm4jc2vk The COVID-19 pandemic is an ongoing global health emergency caused by a newly discovered coronavirus SARS-CoV-2. The entire scientific community across the globe is working diligently to tackle this unprecedented challenge. In silico studies have played a crucial role in the current situation by expediting the process of identification of novel potential chemotypes targeting the viral receptors. In this study, we have made efforts to identify molecules that can potentially inhibit the SARS-CoV-2 main protease (M(pro)) using the high-resolution crystal structure of SARS-CoV-2 M(pro). The SARS-CoV-2 M(pro) has a large flexible binding pocket that can accommodate various chemically diverse ligands but a complete occupation of the binding cavity is necessary for efficient inhibition and stability. We augmented glide three-tier molecular docking protocol with water thermodynamics to screen molecules obtained from three different compound libraries. The diverse hits obtained through docking studies were scored against generated WaterMap to enrich the quality of results. Five molecules were selected from each compound library on the basis of scores and protein-ligand complementarity. Further MD simulations on the proposed molecules affirm the stability of these molecules in the complex. MM-GBSA results and intermolecular hydrogen bond analysis also confirm the thermodynamic stability of proposed molecules. This study also presumably steers the structure determination of many ligand-main protease complexes using x-ray diffraction methods. Communicated by Ramaswamy H. Sarma Newly emerging pathogens are becoming worldwide threats to public health in recent times. In December 2019, a cluster of patients with 'pneumonia of unknown origin' was reported in Wuhan city of Hubei province, China. The etiology of this disease was later associated with a virus of the coronavirus family. The virus was later termed as Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) by experts of the International Committee on Taxonomy of Viruses (ICTV) due to its apparent similarity to SARS-CoV from 2002. The World Health Organization (WHO) labelled this infectious disease as Coronavirus Disease-2019 (COVID-19) and declared it as a global health emergency. Though the earlier outbreaks of severe acute respiratory syndrome coronavirus (SARS-CoV) in China and middle east respiratory syndrome coronavirus (MERS-CoV) in the Middle east, posed a significant threat, with the COVID-19 numbers increasing exponentially, the global public health has been under serious stress (Paules et al., 2020; Wu et al., 2020) . As of 6 September 2020, WHO reports more than 27 million confirmed cases and 900,000 deaths worldwide (WHO, 2020) . An enormous amount of research is underway to develop a vaccine or identify drugs against this virus (Abd El-Aziz & Stockand, 2020; Al-Kassmy et al., 2020; Elfiky, 2020; Huggins, 2020) . Considerable work on the development of drugs for other coronaviruses such as SARS-CoV and MERS-CoV has enabled rapid identification of drug targets for SARS-CoV-2 (Anand et al., 2003; Ghosh et al., 2005; Kumar et al., 2016; ) . The RNA genome of SARS-CoV-2 encodes for various structural, non-structural and accessory proteins. The non-structural proteins like 3-chymotrypsin-like protease (3CL pro ), papain-like protease, helicase, and RNA-dependent RNA polymerase as well as structural proteins such as spike glycoprotein are attractive targets for developing antiviral drugs (Li & De Clercq, 2020) . The main reading frame of RNA, ORF 1ab encodes for two overlapping polyproteins pp1a and pp1ab which undergo proteolysis through two viral proteases, the papain-like protease along with 3CL pro also referred to as main protease (M pro ) to produce 16 non-structural proteins. The remaining portion of the RNA codes for accessory proteins and structural proteins viz., spike (S) protein, membrane (M) protein, envelop (E) protein, nucleocapsid (N) protein, hemagglutinin esterase (HE) glycoprotein (Ullrich & Nitsche, 2020; Ziebuhr et al., 2000) . The cleavage of viral polypeptides pp1a and pp1ab is an important step for viral replication. Therefore, inhibiting viral proteases can be an attractive option to inhibit viral replication. Such an approach has been successfully applied earlier to develop several antiviral drugs that are used to treat Hepatitis C (HCV) and HIV infections (Agbowuro et al., 2018) . Since the advent of COVID-19 antivirals including ritonavir and lopinavir were suggested for the treatment of the disease however they were found ineffective or of little benefit in the treatment of disease Smolders et al., 2020) . Several other investigations are underway testing the efficacy of broad-spectrum antivirals like interferon-a/c (Felgenhauer et al., 2020) and other HCV or HIV protease inhibitors like danoprevir and darunavir Nicolini et al., 2020) to treat the infection. The structure of M pro is rather conserved among coronaviruses. Moreover, till now there is no known human homolog for viral M pro which makes it an ideal drug target (Ullrich & Nitsche, 2020; . Intensive efforts are being made to identify options to treat SARS-CoV-2 infection focusing on enzymes involved in various viral pathways such as TMPRSS2 (Huggins, 2020) , RNA polymerase (Ivashchenko et al. 2020; Yin et al., 2020) , Spike protein (Lei et al., 2020; Liu et al., 2020) , M pro Ma et al., 2020; and helicase (Shu et al., 2020; Yuen et al., 2020) . Previous studies based on SARS-CoV and MERS-CoV recognised a-ketoamides as broadspectrum inhibitors of coronaviruses . This had provided leads to the scientific community to identify new inhibitors for SARS-CoV-2. Several covalent inhibitors have already been identified for SARS-CoV-2 M pro and their co-crystallized structures successfully solved Ma et al., 2020; . Two covalent protease inhibitors of ketoamide category boceprevir and GC376 exhibited strong in vitro activity against SARS-CoV-2 M pro ). An organoselenium compound ebselen has recently emerged as a promising drug lead candidate in cell-based assays which targets this crucial enzyme Menendez et al., 2020) . In this study, we focused on viral M pro to identify molecules with potential to bind the active site of the SARS-CoV-2 M pro structure (PDB ID -5R82). Despite similar studies have published earlier, we followed a three-tiered docking protocol using commercially available ligand libraries and WM/MM-GBSA calculations to meet our objective. Our study focuses on small molecule candidates of natural or synthetic origin which can accelerate drug discovery efforts to combat the COVID-19 pandemic. Table S1 contains the details of all the experimentally solved crystal structures of M pro . The resolution of these structures having 306 residues ranges from 1.31 Å to 2.20 Å. Amino acid sequence alignment SARS-CoV-2 with SARS-CoV M pro reveals a high degree of sequence identity ($96 percent). Superimposition of their crystal structures also revealed high structural similarity with conserved active site region (Ullrich & Nitsche, 2020) . This protein structure forms a homodimer, crystallized with a small fragment molecule bound to the substrate-binding site. Each monomer forms three subdomains with substrate binding site harbouring in a cleft between domain II and domain III. The 3D-structural analysis of 5R82 reveals that the binding cavity is larger than the space occupied by the co-crystallized ligand. The substrate-binding pocket is further divided into five sub-pockets S1, S2, S3, S4 and S1' which accommodate different amino acid residues of the substrate polypeptide. The two monomers of M pro are connected through non-covalent linkages where SER1 of each monomer interacts with GLU166 of the other monomer which stabilizes the S1 sub-pocket of the binding site in M pro protein. The proteolytic activity of M pro is believed to follow a multi-step mechanism using CYS145 and HIS41 residues (catalytic dyad) with help of a water molecule, present in S1' sub-pocket of the binding cavity Ullrich & Nitsche, 2020) . A ligand that can form interactions with these catalytic residues along with hydrophilic and hydrophobic interactions with residues of other sub-pockets can act as a potent SARS-CoV-2 M pro inhibitor. Figure 1 shows the biological assembly of SARS-CoV-2 homodimer with substrate-binding pocket highlighting the various sub-pockets and catalytic residues. The residues surrounding the binding cavity are shown in Table 1 . Among these residues THR26, HIE41, SER46, PHE140, GLY143, CYS145, GLU166 and GLN 189 are important for main protease inhibition. Three databases namely Enamine (Shivanyuk et al., 2007) , Specs (www.specs.net), Natural product (Sterling & Irwin, 2015) were refined by eliminating duplicates and adding missing hydrogen atoms. All the ligand libraries were prepared using Ligprep module of Schr€ odinger suite (Release 2019-1: LigPrep, Schr€ odinger, LLC, New York, NY, 2020). The 3D structures were generated using the new improved Optimized Potential for Liquid Simulations-3e (OPLS-3e) force field (Harder et al., 2016) . Ligands were desalted and possible tautomers and stereoisomers were generated while preserving the specified chiralities. Epik (Shelley et al., 2007) utility was used to generate possible ionisation states at pH 7.0 ± 2.0. We selected PDB 5R82 from the list of available structures for our study as it has the best resolution (1.31 Å) amongst all the available structures. The protein structure file was prepared using the Protein Preparation Wizard of Schr€ odinger's Maestro (Release 2019-1: Protein Preparation Wizard, Schr€ odinger, LLC, New York, NY, 2020). This involved addition and optimization of hydrogen bonds, disulphide bonds creation, atomic clashes removal, addition of formal charges to the hetero groups and then optimizing at neutral pH. Finally, the structure was minimized using the OPLS-3e force field. The crystal water molecules from the crystal structure near the catalytic site and substrate binding site were retained for the molecular docking studies. The receptor grid was generated at the centroid of the bound ligand in the crystal structure of SARS-CoV-2 M pro . The grid box was extended up to 15 Å for the inner box and 26 Å for the outer box covering the entire binding site cavity completely. We used a series of hierarchical filters to search for probable locations of the ligand in the binding site of a receptor. Docking studies were performed using the Glide module of Schr€ odinger suite (Release 2019-1: Glide, Schr€ odinger, LLC, New York, NY, 2020). In a three-tiered docking strategy, the prepared ligands from the 3D databases were docked in three stages, following the docking protocols, which are High Throughput Virtual Screening (HTVS), Standard Precision (SP) and Extra Precision (XP). The top $5 percent of the screened ligands from the first stage of HTVS docking were passed on to the second stage of SP docking. 500 screened molecules from SP docking were then docked using a more accurate and computationally intensive XP mode. The HTVS and SP docking use the same scoring function, however, HTVS reduces the number of intermediate conformers that reduce the thoroughness of the final torsional refinement. On the other hand, XP mode uses a more extensive sampling than SP and employs an improved scoring function that rewards or penalises certain interactions (Friesner et al., 2006) . No functional group constraints or torsional constraints were applied during the docking process. Post docking minimization was performed for the ligands and all other parameters were set to default when following the three-tiered docking strategy. WaterMap calculations were used to assess the thermodynamic parameters of waters in the binding site. These parameters can be used for ligand design by examining how ligands displace/replace these hydration sites. It is based on a small molecular dynamics simulation of explicit solvent followed by clustering and generation thermodynamic properties of hydration sites. It classifies water hydration sites into two categories, namely stable water and unstable waters based on the thermodynamic properties like entropy (-TDS), enthalpy (-DH), free energy (DG) (Abel et al., 2011; Riniker et al., 2012) . WaterMap calculations for the prepared protein structure were performed using default simulation parameters (TIP4P solvent model at 300 K, 1 atmospheric pressure and 2 ns of simulation time) and treating existing waters as a part of explicit solvent (Release 2019-1: WaterMap, Schr€ odinger, LLC, New York, NY, 2020). The binding site was defined using the coordinates of the co-crystallized ligand. The truncate protein option was left unchecked for this study. The generated hydration sites were examined according to their enthalpy, entropy and free energy values. The WM/MM scoring approach was employed to score the ligands which uses a combination of MMGBSA and WaterMap data to assess the best binding pose. The scoring was performed using VSGB solvation model and defining residues within the vicinity of the ligand as non-flexible. Two complexes of protein with ligands ZINC31167921 and ZINC67912395 were considered for Molecular dynamics (MD) simulations to understand the dynamic changes in the complex and binding interactions. These two complexes were submitted for MD simulations of $100ns using AMBER18 (Case et al., 2005) . The protein topologies were prepared using ff14SB force field (Maier et al., 2015) parameters with the help of AMBER LEaP module while antechamber module and GAFF force field with am1bcc charges (Wang et al., 2004) were utilized to generate the parameters for ligands. The system was solvated in TIP3P water molecules and neutralized by adding the necessary amount of Clions in electrostatically preferred positions as the protein is positively charged at neutral pH (Harrach & Drossel, 2014) . The electrostatic interactions were estimated using Particle Mesh Ewald (PME) method with non-bonded cutoff kept at 10 Å. The solvated complexes were minimized three times for 1000, 500 and 250000 steps through a short steepest descent minimization for 50% steps followed by conjugate gradient minimization for the remaining steps. Heating was performed using NVT ensemble for 50 ps where the protein-ligand complex was restrained with force constant of 2 kcal/mol/Å 2 and gradually heated from 0 to 300 K. Three density equilibration steps were carried out; two steps of 50 ps each with the restrained weight of 2 kcal/mol/Å 2 and 1 kcal/mol/Å 2 respectively followed by one step of 100 ps unrestrained density equilibration. The system was equilibrated to free simulation for 1 ns at 300 K temperature and 1 atm pressure. Finally, a long-range simulation for 100 ns production run was performed at 300 K temperature and 1 atm pressure. Coordinate trajectories were recorded after every 10 ns for complete MD run. MD was performed to calculate a single trajectory using pmemd program of AMBER18 running on NVIDIA Tesla K20Xm GPU workstation (Salomon-Ferrer et al., 2013) . The outline of the strategy followed for this study is mentioned in a flowchart ( Figure 2 ). Main protease crystal structure (PDB:5R82) offers the opportunity to employ a structure-based drug design strategy in identifying novel COVID-19 inhibitors. In the experimental structure, ligand RVZ forms a hydrogen bond with GLN189 and a weak pi-pi stacking interaction with HIS41 as shown in Figure 3 . It is clear that the binding site is surrounded by various loop regions that owe to the flexible nature of the binding cavity. Cysteine and serine proteases usually possess three catalytic residues while M pro has only two catalytic residues in its active site (Ullrich & Nitsche, 2020) . It can be presumed that a water near the catalytic HIS41 residue plays the role of the third residue by making a strong hydrogen bond with it. A water molecule near the HIS41-CYS145 catalytic dyad was preserved for docking studies. Another crystal water was preserved near the catalytic site stabilized by three hydrogen bonds with HIS41, HIS164 and ASP187 residues. The molecules screened using XP docking were subsequently rescored using WM/MM-GBSA scoring to find the relative binding affinities towards SARS-CoV-2 M pro . The combined WaterMap MM-GBSA scoring protocol generally provides a better correlation with experimental activity than either of them when used alone (Abel et al., 2011) . Approximately 7.7 lakh molecules from three databases were screened by HTVS followed by SP and XP docking. The top docked molecules from each database were chosen for further analysis. The molecules were filtered according to the docking scores and analysed for hydrogen bond interaction with any of the residues in the binding site region especially HIS41, CYS145, GLU166 and GLN189. The data of five selected molecules from each database is presented in Table 2 . These results include parameters like docking scores, XP GScore, MMGBSA DG and WM/MM DG values. Average ligand binding free energies for selected molecules range from around À45 kcal/mol to À60 kcal/mol. The 3D interaction diagram of the top 5 molecules from each database are shown in Supplementary Figure S1. The selection of the molecules was based on the fact that these molecules possess required bonded and non-bonded interactions with residues involved in catalytic mechanism and engulfing the S1, S1', S3 and S4 sub-pockets of 5R82. Most of the selected molecules occupy the catalytic S1 and S1' subsites showing interactions with atleast one of the residues of catalytic dyad and hydrogen bonds with the other residues of the binding site. The docking analysis revealed that a common pi-stacking interaction usually forms by the alignment of aromatic rings of the ligands and HIE41 side chain. Ligands Z2229246315 and ZINC15675289 formed both hydrogen bond and pi-stack interaction with HIE41. The selected ligands ZINC38139713, ZINC31165691, Z3061991917, ZINC49781425 and ZINC44305556 form multiple hydrogen bonds/salt bridge with GLU166 (residue involved in enzyme dimerization along with SER1 of the other monomer) that can interfere with the biological assembly of the enzyme and inhibit its activity. Ligands ZINC38139713, ZINC49781425 and ZINC67912395 from the Specs and natural products database showed a significantly stronger affinity towards the M pro indicated by multiple hydrogen bonds and hydrophobic interactions that may help ligand to stabilize into the binding pocket. Also, a good MMGBSA DG and WM/MM DG scores for these three ligands affirmed their likelihood to act as M pro inhibitors. Molecules like Z3040398264, Z2063273911, ZINC67911229, ZINC5675289 and ZINC44305556 contain an electrophilic warhead group (carbon atom of the carbonyl group) in close proximity to the catalytic thiol-group of CYS145 which may behave as peptidomimetic substrate to form strong covalent interaction with the enzyme. A general comparison of these ligands to active in vitro molecules boceprevir and GC376 displayed similar interaction patterns. The two crystal water molecules conserved for docking simulation didn't show any major water bridged interaction except in ligand Z2063273911 which forms a water bridged interaction with HIE164. Although water molecules are involved in the catalytic mechanism of the protein, they may not have a significant role in ligand binding interactions. A list of the type of ligand interactions with the corresponding residues involved is provided in Table 4 . The WaterMap calculations were used to compute the free energy (DG) of the explicit waters around the binding site of protein. The hydration sites generated for the crystal structure 5R82 are displayed in Figure 4 . Hydration sites are coloured according to their free energy values ranging from green to red. The most stable ones are represented using green while the most unstable or energetically unfavourable ones are shaded in red. There are particularly three hydration sites with high DG values of 5.47, 4.96 and 3.67 kcal/mol located in hydrophobic site S3 and site S4 of the binding cavity. These waters have low entropy due to positional restrictions and a high enthalpy due to lack of polar interactions. High DG values along with high enthalpy (DH ) 0kcal/ mol) indicate that it is favourable to displace them. Any ligand with nonpolar atoms entering these sites will tend to displace these waters and get rewarded with a higher binding affinity. No favourable replacement sites were present in the binding cavity. Except for a few sites, most of the other hydration sites in the binding cavity are highly conserved and can only be utilized to form bridging interactions. One of the conserved water molecules used for our study partially overlaps this kind of hydration site. Simple docking calculations ignore the free energy contributions from desolvation of binding pocket. Data from WaterMap calculations can be used to optimize the ligand pose in the binding site. WM/MM scoring was performed for PlogPo/w: Predicted octanol/water partition coefficient. (-2.0 -6.5). PlogHERG: Predicted IC 50 value for blockage of HERG K þ channels. (concern below -5). Natural products Database a) ZINC38139713 THR25, THR26, SER46, GLY143, GLU166, HIE164 -HIE41 THR25, THR26, LEU27, HIE41, SER46, MET49, LEU141, ASN142, GLY143, CYS145, HIS163, HIE164, MET165, GLU166, ASP187, ARG188, GLN189 b) ZINC31165691 HIE41, ASN142, GLU166, GLN189 --THR25, LEU27, HIE41, THR45, SER46, MET49, LEU141, ASN142, CYS145, HIS163, MET165, GLU166, LEU167, ARG188, GLN189 c) ZINC31167921 ASN142, GLU166, GLN189, THR190 -HIE41 THR25, THR26, LEU27, HIE41, MET49, LEU141, ASN142, GLY143, CYS145, HIE164, MET165, GLU166, LEU167 --THR25, LEU27, HIE41, CYS44, THR45, SER46, MET49, ASN142, GLY143, CYS145, HIE164, MET165, GLU166, LEU167, PRO168, ASP187, ARG188, GLN189, THR190, GLN192 h) Z2229246315 THR25, HIE41, GLN189 -HIE41 THR24, THR25, THR26, LEU27, HIE41, THR45, SER46, MET49, GLY143, CYS145, HIE164, MET165, GLU166, ARG188, GLN189 i) Z3061991917 GLY143, GLU166, GLN189 GLU166 HIE41 THR25, THR26, LEU27, HIE41, CYS44, THR45, SER46, MET49, PHE140, LEU141, ASN142, GLY143, CYS145, HIS163, MET165, GLU166, ARG188, GLN189 j) Z2063273911 LEU141, ASN142, HIE164, GLU166, GLN189 -HIE41 THR25, LEU27, HIE41, SER46, MET49, LEU141, ASN142, CYS145, HIS163, HIE164, MET165, GLU166, ASP187, ARG188, GLN189 Specs Database k) ZINC49781425 THR26, PHE140, LEU141, GLU166, ARG188 -HIE41 THR25, THR26, LEU27, HIE41, THR45, SER46, PHE140, LEU141, ASN142, GLY143, SER144, CYS145, HIS163, MET165, GLU166, LEU167, PRO168, ARG188, GLN189, THR190, GLN192 l) ZINC67912395 GLU166, LEU167, GLN189 -HIE41 THR24, THR25, THR26, LEU27, HIE41, CYS44, THR45, SER46, MET49, LEU141, ASN142, GLY143, CYS145, HIS163, HIE164, MET165, GLU166, LEU167, PRO168, ASP187, ARG188, GLN189 m) ZINC67911229 LEU141, GLY143, GLU166, GLN189 --THR24, THR25, THR26, HIE41, THR45, SER46, MET49, LEU141, ASN142, GLY143, SER144, CYS145, HIS163 Combining the WM/MM scoring with docking protocol changes the ligand rankings significantly which may be related to experimental binding affinities in a much better way (Abel et al., 2011) . The selected molecules had WM/MM scores in a range of around À18 to À26 kcal/mol as depicted in Table 2 . Molecules ZINC38139713, ZINC31165691, ZINC31167921, ZINC49781425, ZINC67912395 and ZINC67911229 had good WM/MM DG bind values along with high docking scores. From the visual analysis, it was observed that many of the screened molecules had significant overlap with the three most unstable hydration sites in the binding cavity and a few of them showed a more significant binding affinity in terms of the calculated parameters. Figure 4 indicates the interaction of one such ligand ZINC31167921 with the WaterMap. The ligand ZINC31167921 effectively overlapped the unfavourable hydration sites identified using WaterMap calculations. The aromatic ring of the ligand fits snuggly into binding pocket S3 with an attached fragment entering the S4 pocket contributing to a high WaterMap DG bind value. In a similar fashion, ligands Z2706637055 and ZINC67912395 a satisfactory overlap with those hydration sites. Overall, WM/MM scoring helped in identification of ligands containing a balance of polar and hydrophobic interactions necessary for viral enzyme inhibition. ADMET analysis ADME properties of the selected compounds were predicted using Qikprop module of Schr€ odinger. Various basic physiochemical properties such as PlogPo/w (Predicted octanol/ water partition coefficient), PlogHERG (Predicted IC50 value for the blockage of HERG K þ channels), PPCaco (Predicted Caco-2 cell permeability for the gut blood barrier), PlogBB (predicted brain/blood partition coefficient) and PlogKhsa (Predicted binding to human serum albumin) of these compounds were calculated. The calculated parameter values of the selected compounds are given in Table 3 . Most of the selected molecules fall in a decent range of values for the given parameters. In silico toxicity parameters namely mutagenicity and hepatoxicity were assessed using pkCSM server (Pires et al., 2015) which makes use of graph-based signatures to calculate properties. The predicted toxicity parameters viz.; Drug induced liver injury (DILI) and AMES toxicity for the selected compounds are listed in Table 3 . Molecular dynamics simulations are essential to confirm the stability of the predicted ligand binding mode (Liu et al., 2017) . We performed a 100 ns MD simulation for two of the best molecules to inspect their mechanical stability in complex and ensure stable binding under dynamic conditions. The degree of variation in RMSD is inversely related to the stability of complex: the smaller the variation, the greater the stability. The RMSD of both complexes showed minimum difference (RMSD < 1 Å) in reference to their initial coordinates which states that the complexes were stable throughout the simulation, as shown in Figure 5 . Later we focused on the local atomic level fluctuations of the protein in dynamic conditions. The RMSF of complex 1 (ZINC31167921) showed a higher degree of fluctuation compared to complex 2 (ZINC67912395). The RMSF analysis of each complex revealed a high fluctuation region which corresponds to domain III (between residues position 201-303) of the protein. This portion is far from the complex's interface or the active site and is associated with protein dimerization. Most of the binding site residues in domain I and domain II have lesser fluctuations with RMSF values below 1 Å (shown as D1 and D2 in Figure 6 ). Although not involved directly in the ligand binding, residues like THR24 form a part of the catalytic site as well. It has slightly higher fluctuation, which may be attributed to its location in the turn region between the two b-sheets. Another residue GLY138 with a little high RMSF value in complex 2 is located in the loop region at the interface of domain II and domain III. The fluctuations in all the three domains were more or less same in both complexes with slightly higher values in complex 1. Our interaction analysis showed that small-molecule ZINC31167921 had formed a more number of intermolecular H-bonds with the protein in complex 1 compared to ZINC67912395 in complex 2. The oxygen atom of THR25 can form H-bond with oxygen of ZINC31167921with an occupancy of 90.8%. Meanwhile oxygen atom as well as nitrogen of GLU166 displayed two hydrogen bonds with ZINC31167921 during the course of simulation with average occupancy of 54.5 and 43.2 percent. Another H-bond with an occupancy of 48 percent was formed with the side chain nitrogen at delta position which was absent earlier. Other residues like GLN189, THR26 and GLN192 displayed some propensity to form H-bonds with ZINC31167921. The gamma oxygen of THR25 as well as oxygen and nitrogen atoms of THR26 displayed hydrogen bonds with ZINC67912395 with an occupancy of 51.6, 65.6 and 46.4 percent. Oxygen atoms of ASP187 and GLU166 displayed H-bonds with ZINC67912395 with occupancy of 43.5 and 34.1 percent. The ligand ZINC67912395 also displayed H-bond interactions with side chains of ASN142 and GLN189 with average occupancy of 34.1 and 23.2 percent. The graphs highlighting the hydrogen bond interaction data are presented in supplementary data Figure S2 . MM-GBSA analysis revealed a high van der Waals contribution in the protein and small molecule binding. Overall binding free energy (DG) for ZINC31167921 was found to be À50.2062 kcal/mol and À24.6658 kcal/mol for ZINC67912395. Residues THR25, THR26, HIE41, MET49, CYS145, GLU166 and GLN189 contributed to ligand binding. Catalytic residues of M pro (HIE41 and CYS145) contributed to the ligand-binding in both cases. Per residue scores indicated that THR25 contributed maximum to the binding free energy, as evident from Figure 7 . Contribution of all the above-mentioned residues was higher for ZINC31167921 than ZINC67912395, which is apparent from the overall binding free energy (DG) values. The COVID-19 pandemic calls for a quick response from the scientific realm and with no specific remedy available, computational approach is the key for identification of potential inhibitors that can be repurposed for effective therapeutics. The main protease involvement in viral replication makes it an important target to inhibit the virus growth. The present study has offered insights into the possible route of drug design against SARS-CoV-2 M pro . An organized approach involving three-tiered molecular docking and water thermodynamics was followed to perform this study. We report 15 candidates in this paper that show optimal binding characteristics to the viral main protease. Potential leads obtained from the binding affinity screening studies can be further utilized for in vitro, in vivo analysis. From the results obtained we can expect that the identified molecules ZINC31167921 and ZINC67912395 may exhibit good antiviral activity against the new coronavirus. Molecular dynamics and MM-GBSA results verified that the complexes are not only thermodynamically stable but structurally robust. Although this data doesn't confirm the antiviral activity of the selected molecules, yet it provides a basis for further in vitro/ in vivo studies. The information derived from water thermodynamic studies offers possible routes to design and optimize these hits to increase their selectivity towards SARS-CoV-2. Our subsequent studies will concentrate on understanding the covalent interaction abilities of identified molecules possessing reactive warhead groups (Z3040398264, Z2063273911, ZINC67911229, ZINC5675289 and ZINC44305556) and establishing their stability in the complex. Recent progress and challenges in drug development against COVID-19 coronavirus (SARS-CoV-2) -an update on the status. Infection, Genetics and Evolution : journal of Molecular Epidemiology and Evolutionary Genetics in Infectious Diseases Contribution of explicit solvent effects to the binding affinity of small-molecule inhibitors in blood coagulation factor serine proteases Proteases and protease inhibitors in infectious diseases Vaccine candidates against coronavirus infections. Where does COVID-19 stand? Viruses Coronavirus main proteinase (3CLpro) structure: Basis for design of anti-SARS drugs The Amber biomolecular simulation programs First Clinical Study Using HCV Protease Inhibitor Danoprevir to Treat Naive and Experienced COVID-19 Patients. medRxiv: 2020 Antiviral activity and safety of darunavir/cobicistat for the treatment of COVID-19 Anti-HCV, nucleotide inhibitors, repurposing against COVID-19 Inhibition of SARS-CoV-2 by type I and type III interferons Extra precision glide: Docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes Both Boceprevir and GC376 efficaciously inhibit SARS-CoV-2 by targeting its main protease Structure-based design, synthesis, and biological evaluation of peptidomimetic SARS-CoV 3CLpro inhibitors Design and synthesis of peptidomimetic severe acute respiratory syndrome chymotrypsin-like protease inhibitors OPLS3: A force field providing broad coverage of drug-like small molecules and proteins Structure and dynamics of TIP3P, TIP4P, and TIP5P water near smooth and atomistic walls of different hydroaffinity Structural analysis of experimental drugs binding to the SARS-CoV-2 target TMPRSS2 AVIFAVIR for Treatment of Patients with Moderate COVID-19: Interim Results of a Phase II/III Multicenter Randomized Clinical Trial Structure of M(pro) from SARS-CoV-2 and discovery of its inhibitors Structural basis for the inhibition of SARS-CoV-2 main protease by antineoplastic drug carmofur Identification and evaluation of potent Middle East respiratory syndrome coronavirus (MERS-CoV) 3CLPro inhibitors Identification, synthesis and evaluation of SARS-CoV and MERS-CoV 3C-like protease inhibitors Neutralization of SARS-CoV-2 spike pseudotyped virus by recombinant ACE2-Ig Therapeutic options for the 2019 novel coronavirus (2019-nCoV) Efficacy and Safety of Lopinavir/ Ritonavir or Arbidol in Adult Patients with Mild/Moderate COVID-19: An Exploratory Randomized Controlled Trial Exploring the stability of ligand binding modes to proteins by molecular dynamics simulations Potent neutralizing antibodies against multiple epitopes on SARS-CoV-2 spike Boceprevir, GC-376, and calpain inhibitors II, XII inhibit SARS-CoV-2 viral replication by targeting the viral main protease ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB Molecular characterization of ebselen binding activity to SARS-CoV-2 main protease Reply to: "Antiviral activity and safety of darunavir/cobicistat for treatment of COVID-19 Coronavirus infectionsmore than just the common cold pkCSM: Predicting Small-Molecule Pharmacokinetic and Toxicity Properties Using Graph-Based Signatures Free enthalpies of replacing water molecules in protein binding pockets routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald Epik: A software program for pK( a ) prediction and protonation state generation for drug-like molecules Enamine real database: Making chemical diversity real SARS-coronavirus-2 Nsp13 Possesses NTPase and RNA helicase activities that can be inhibited by bismuth salts SARS-CoV-2 and HIV protease inhibitors: Why lopinavir/ritonavir will not work for COVID-19 infection ZINC 15-ligand discovery for everyone The SARS-CoV-2 main protease as drug target Development and testing of a general amber force field Coronavirus disease 2019 (COVID-19): weekly epidemiological update A new coronavirus associated with human respiratory disease in China Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by remdesivir SARS-CoV-2 nsp13, nsp14, nsp15 and orf6 function as potent interferon antagonists. Emerging Microbes & Infections a-Ketoamides as broad-spectrum inhibitors of coronavirus and enterovirus replication: Structure-based design, synthesis, and activity assessment Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved a-ketoamide inhibitors Virus-encoded proteinases and proteolytic processing in the Nidovirales The authors thank the National Institute of Pharmaceutical Education and Research (NIPER) SAS Nagar, Department of Pharmaceuticals, Ministry of Chemicals and Fertilizers, New Delhi, Government of India for providing the facility. The authors declare no competing interests. M.E.S designed the study. All the authors contributed to computations, analysis and preparation of the manuscript.