The novel coronavirus outbreak was first reported in Wuhan city of China on December 2019 which causes the Coronavirus disease 2019 (COVID-19) (Sohrabi et al., 2020 ) is now spread worldwide and effecting the entire human race. The virus strain responsible for this ongoing pandemic is called Severe Acute Respiratory Syndrome 2 (SARS-CoV-2) which is the successor to SARS-CoV-1 and also related to Middle East coronavirus (MERS-CoV) (Petrosillo et al., 2020) having the highest mortality rate among all human coronaviruses (Elfiky & Azzam, 2020) . All three viruses have a common host and have originated from bats (Hemida & Ba Abduallah, 2020; Hu et al., 2015) . The spread of SARS-CoV-1 (in the year 2002) and MERS-CoV (in the year 2012) to humans is believed to be through civet cats (Wang & Eaton, 2007) and camels (Azhar et al., 2014) unfortunately the mediator for SARS-CoV-2 is still unknown but researchers suspect that ant-eating pangolins may be responsible for this spread as Pangolin-CoV is 91.02% identical to SARS-CoV-2 (T. .
During this pandemic several drug targets have been identified which includes RNA-dependent RNA polymerase (Anwar et al., 2020; Arya & Dwivedi, 2020; Babadaei, Hasan, Vahdani, et al., 2020; Elfiky, 2020b) , Spike glycoprotein (Adeoye et al., 2020; Basit et al., 2020; de Oliveira et al., 2020; Elfiky, 2020a; Sinha et al., 2020; Wahedi et al., 2020) , Angiotensin-converting enzyme 2 (ACE2) (Abdelli et al., 2020; Ahmad et al., 2020; Hasan et al., 2020; V. Kumar, Dhanjal, Bhargava, et al., 2020; Veeramachaneni et al., 2020; Wahedi et al., 2020) , Nucleocapsid protein (Sarma et al., 2020) , Papain-like proteinase (Quimque et al., 2020) , 3 C-like proteinase (Das et al., 2020; Enmozhi et al., 2020; Islam et al., 2020; Umesh et al., 2020) and they may be useful for drug development. However, design and search of multi-epitope vaccine candidate molecules have also been initiated in different laboratories (Abraham Peele et al., 2020; Enayatkhani et al., 2020) .
The SARS-CoV-2 is a positive-sense single-stranded RNA virus (Hendaus, 2020; Lu et al., 2020) which utilizes non-structural protein (nsp) complexes for the process of transcription and replication in human cells. The replicase polyprotein is processed by viral protease known as Papain like protease (PLP) and 3 C-like protease (3CL protease) or also known as the main protease which is attractive drug targets for the present pandemic COVID-19 (Al-Khafaji et al., 2020; Bhardwaj et al., 2020; Gupta et al., 2020; Gyebi et al., 2020; A. Kumar, Choudhir, Shukla, et al., 2020; D. Kumar, Kumari, Jayaraj, et al., 2020; Mittal et al., 2020; Sk et al., 2020) and is most widely studied target for COVID-19 which also remains conserved among several coronaviruses (CoV), so we select it for our study. The X-ray crystal structures of SARS-CoV-2 3CL-protease were surveyed from RCSB Protein Data Bank (Berman et al., 2000; Burley et al., 2019) which showed the protein consists of three domains. The Domains I (consisting of 1-99 amino acids), Domain II (consisting of 100-183 amino acids), and Domain III (consisting of 200-306 amino acids) in which domain I contributed His41 and domain II contributed Cys145 together forming a catalytic dyad found in entire CoV superfamily. However, the third catalytic residue is missing and the position seems to be occupied by an invariant water molecule (L. Zhang, Lin et al., 2020) . The Domain I is connected to Domain III through a connecting loop consisting of 184-199 amino acid residues (Figure 1(a) ). The ligand-binding site is comprised of Thr26, His41, Phe140, Asn142, Gly143, Cys145, His163, His164, Met165, Glu166, Gln189 and Thr190 residues as can be seen in 6lu7 crystal structure (Jin et al., 2020) (Figure 1(b) ).
The Hydroxychloroquine is an FDA approved drug that has been used to treat malaria(Al- Bari, 2015; Tanenbaum & Tuffanelli, 1980) rheumatoid arthritis (Rynes, 1988) , chronic discoid lupus erythematosus (Wahie & Meggitt, 2013) , and systemic lupus erythematosus (Chew et al., 2020) . During this pandemic emergency, HCQ has been showing promising effect in the treatment of COVID-19 (Boopathi et al., 2020; Hendaus & Jomha, 2020) so, several theoretical and experimental studies are conducted to understand the inhibiting mechanism of this malarial drug to SARS-CoV-2 (Amin & Abbas, 2020; Beura & Chetti, 2020; Beura & Prabhakar, 2020; Devaux et al., 2020; Gautret et al., 2020) .
The X-ray crystallographic studies have indicated the inhibitor/drug binding site on 3CL-protease which helps in identification for new potential inhibitor molecules, consequently a large number of theoretical studies are conducted that attempts identification of repurposing drugs for COVID-19 (Aanouz et al., 2020; Choudhury, 2020; Elmezayen, Al-Obaidi, Şahin, & Yelekçi, 2020; Joshi et al., 2020; S. A. Khan, Zia, Ashraf, Uddin, & Ul-Haq, 2020; ; Lobo-Galo, Terrazas-L opez, Mart ınez-Mart ınez, & D ıaz-S anchez, 2020; Mahanta et al., 2020; Muralidharan, Sakthivel, Velmurugan, & Gromiha, 2020; Pant, Singh, Ravichandiran, Murty, & Srivastava, 2020) . In our study we carefully investigated the catalytic and inhibitor binding site of 3CL-protease in the crystal and MD-simulated structures. We also report the binding mechanism and stabilization of Hydroxychloroquine with SARS-CoV-2 main protease through non-covalent interactions as well as attempt has been made to model some new structural analogs (Mod I and I / , Mod II and II / ) and their recognition dynamics to conserved watermediated catalytic triad.
The X-ray structures of some 3CL-protease from SARS-CoV-2 have been surveyed which were selected from RCSB Protein Data Bank (Supplementry material Table S1 ). The position and hydrogen bonding interaction of catalytic residues were examined using the UCSF Chimera program (Pettersen et al., 2004) .
The ligand (HCQ, Mod I, Mod I / , Mod II, Mod II / ) structures ( Figure 2 ) were built by using the Gaussview program (Dennington et al., 2009) . The structures were then subjected to geometry optimization with B3LYP/6-31G level of theory, solvation energies were added using the conductor-like Polarizable Continuum Models (CPCM). The partial atomic charges of the atoms were retained during classical MD-simulation. All the QM calculations were performed by using the Gaussian 09 package (Frisch et al. (2009) Gaussian 09 A.02).
The 3CL-protease has one ligand binding site which has been indicated in the X-ray crystallographic studies. Ligand-receptor docking was performed using Autodock Vina v.1.1.1 (Trott & Olson, 2009 ). The 6LU7 X-ray structure (excluding water and ligand molecules) was considered as a receptor. Two PDBQT files were generated for the receptor protein, one for rigid portion and another for flexible (His41 and Cys145) side chains using AutoDock Tools v.1.5.4 (Morris et al., 2009 ) by assigning Kollman united atom charges (Weiner et al., 1984) . The QM optimized structures of ligand molecules (HCQ and Mod-I, I / , II, II / ) were converted into PDBQT files after including their partial atomic charges using the Gasteiger method (Gasteiger & Marsili, 1980) . Grid point spacing was set at 1 Å and 20 grid points were taken in each direction. As the location of the ligand binding site was already known (PDB id: 6LU7), so grid box was centered at that site. Vina automatically calculated the grid map for searching. All other docking parameters were assigned to their default values. The docked structures were visually inspected and selected for further work. The ligand-receptor docking was also validated with DockThor (A receptor-ligand docking program) server (De Magalhães et al., 2014; Santos et al., 2020) .
Molecular dynamics simulation was performed using NAMD v.2.12 (Kal e et al., 1999; Phillips et al., 2005) with the CHARMM36 force field (Brooks et al., 1983; Huang & MacKerell, 2013; MacKerell et al., 1998) . Then each structure was converted to Protein Structure File (PSF) by Automatic PSF Generation Plug-in within VMD program v. 1.9.2 (Humphrey et al., 1996) . Subsequent energy minimization was performed by the conjugate gradient method. The process was conducted in two successive stages; initial energy minimization was performed for 2000 steps by fixing the backbone atoms, followed by a final minimization for 5000 steps that were carried out for all atoms of the system to ensure the removal of any residual steric clashes. Then the energy minimized structures were simulated at constant temperature (310 K) and pressure (1 atm) by Langevin dynamics (Gullingsrud et al., 2001) using periodic boundary condition. The Particle Mesh Ewald method was applied for full-electrostatics and the Nose-Hoover Langevin piston method was used to control the pressure and dynamical properties of the barostat (Feller et al., 1995) . All-atom molecular dynamics simulation for 100 ns was carried out for HCQ, Mod I, and Mod II docked 3CL-protease structure. The atomic coordinates of simulated structures were recorded at every 2 ps for further analysis.
The ligand-receptor binding energies were calculated by Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) method using Calculation of Free Energy (CaFE) program (Liu & Hou, 2016) . It is capable of handling the simulation trajectories generated by various force field and is powered by VMD and NAMD. Three different energy terms are calculated by the MM-PBSA method which includes energy difference between the complex and separate receptor and ligand molecule in the gas phase from NAMD. Then, the polar solvation free energy was calculated using the APBS program (Baker et al., 2001) and finally, the difference of solvent-accessible surface area (SASA) is measured and the nonpolar solvation free energy was estimated. The binding energies are calculated by using the following equation:
The binding energies of the ligand-receptor complexes were calculated at 20 ns interval over the simulation trajectories.
The analysis of 3CL-protease crystal structures from SARS-CoV-2 showed that the His 41 of Domain I lies close to Cys145 of Domain II which together form the catalytic dyad. An invariant water molecule, HOH445 (W1) is positioned within hydrogen bond distance to His41 ND1 , Arg40 NB , His164 ND1 , and Asp187 OD2 atoms (Figure 3) . 
The optimized structure of HCQ was docked in 6lu7 crystal structure, the ligand molecule binds to a groove on the surface of the protein with binding energy $-5.4 kcal/mol. The catalytic dyad (His41 and Cys145) interacts with the si-face (under the aromatic ring) of the ligand molecule through non-covalent interaction. The His41 endowed parallel displaced p‧ ‧ ‧ p interaction with aromatic ring 1 of HCQ and the Cys145 shown S-H‧‧‧p interaction with aromatic ring 2 of HCQ. The Thr24 OG1 located at the entrance of the inhibitor binding site is stabilized by the O21 hydroxyl group of HCQ through hydrogen bond interaction (Supplementry material Figure S1 ).
During simulation, the hydrogen bond interaction between Cys145 SG and His41 NE2 is ranging from 3.1 to 3.3 Å and the average H-bond interaction of W1 water center to His41 ND1 , Arg40 NB , His164 ND1 , and Asp187 OD2 atoms are observed to be $ 3.0-3.3, 3.1-3.4, 2.8-3.4, 2.7-3.2 Å respectively ( Figure 5) .
A distorted trigonal pyramidal geometry has been observed around the W1 water center where Arg40 NB and His164 at the apex position and His41, His164, and Asp187 are present at the three corners of the plane (Figure 6 ) with $100% residential frequency. In this HCQ complex structure, the Cys145 SG atom is observed to lie near C8 and C9 carbon centers with $3.88 and 3.72 Å distance respectively suggesting these two centers may be important for inhibitor modeling so, these carbon centers were modified by substituting with different functional groups for easy nucleophilic attack by Cys145 catalytic residue.
Observing the structural arrangement and interaction of His41 and Cys145 with HCQ, some ligand molecules were modelled ( Figure 2) and docked into the inhibitor binding cavity of the 6lu7 crystal structure. The docked structures of modified ligand molecules (Mod I, I / , II, and II / ) win the 6lu7 crystal structure are shown in Supplementry material Figure S2 .
The binding energies (B.E.) for Mod I and Mod II are observed to be À6.1 and À5.5 kcal/mol whereas that of Mod I / and Mod II / are À5.4 and À5.3 kcal/mol respectively. The results indicate that the modification of HCQ at the C8 position is energetically favourable compared to the C9-position. Moreover, we may say the formation of carbocation at C8-position of Modified ligand molecules would be more stable than at C9-position due to -I effect of N7-atom in the aromatic ring of the ligand molecule so, Mod I and Mod II structures were selected for further MD-simulation studies. The molecular dynamics simulation of Mod I and Mod II complexed 6lu7 structure shows that the Cys145 SG lies in close proximity to O25 atom of the ligand with average distance ranging from 3.2 to 3.7 Å for Mod I and 3.5 to 4.0 Å for Mod II complexed structure (Figure 7 and 8) . The hydrogen bond interaction between the catalytic dyad Cys145 SG ‧‧ ‧ His41 NE2 varies from 3.2 to 3.5 Å. The W1 water molecule forms strong hydrogen bond interaction with His41 ND1 and Asp187 OD2 thus bridging the gap between these two residues. The His41 ND1 ‧‧‧W1 distance ranges from $3.1 to 3.4 Å and 3.0 to 3.3 Å whereas that of W1‧ ‧ ‧ Asp187 OD2 ranges from 2.7 to 2.9 and 2.7 to 3.0 Å for the respective Mod I and Mod II complex structures (Figure 7 and 8) .
In classical proteases (like chymotrypsin, trypsin or elastase), the catalytic triad is made up of an acid, base and a nucleophile but in coronavirus main proteases (SARS-CoV-2 Mpro (PDB Id 6lu7), SARS-CoV Mpro (PDB Id 2zu5) and MERS-CoV Mpro (PDB Id 5wkk)) presence of water-mediated Aspartic acid residue along with base (His) and nucleophile (Cys) seem to be unique and reflecting the formation of water-mediated catalytic triad so, the occupation of W1 water center adjacent to the base (His41) and its recognition to Asp187 through Asp187 OD2 ‧‧‧W1‧‧‧His41 NE2 interaction could be an evolutionary change in the catalytic site of CoV-Mpro for unknown reasons.
During simulation, the recognition of Asp187 to Mod I/II is observed to made through Asp187 OD2 ‧‧‧W1‧‧‧His41 ND1/ NE2 ‧‧‧Cys145 SG ‧‧‧Mod I/II (O25) path. The carbonyl oxygen (O25) is stabilized by amide backbone atoms of Cys145, Ser144, and Gly143 residues (Figure 9 ).
The binding free energies of the ligand-receptor complexes by the MM-PBSA method (Table 1 ) are calculated. The final binding energy is the cumulative sum of Electrostatic, Van der Waal, Polar solvation, and SASA energies which show that the Mod II structure has better binding energies when compared to Mod I and HCQ molecules. The Mod I structure shows low binding and electrostatic energies compared to the HCQ complex structure. The Mod II complexed structure shows high binding energy of $À213.09 kcal/mol at 80-100ns time interval.
So, analyses of crystal and ligand-bound simulated structures from SARS-CoV-2 Mpro a mechanism has been proposed which may shed some light on the mechanistic role of the above residues in ligand complexed Mpro structures ( Figure 10) .
In SARS-CoV-2 Mpro the catalytic triad is made up of Cys145, His41, and Asp187 residues, a catalytic water (HOH445) which is associated with His41 and Asp187 through hydrogen bond interaction. Initially, His41 abstracts a proton from Cys145 making it a better nucleophile (cysteine thiyl radical, CysS -) for attack on the carbonyl carbon (C24) of the ligand molecule. The Asp187 and HOH445 anchor the His41 molecule to its correct conformation during the formation of transition state complexes. The reaction is 
The Molecular docking followed by MD-simulation of the HCQ complexed structure of SARS-CoV-2 Mpro reveals the binding of this ligand molecule on the surface of the protein structure near to the catalytic dyad (His41 and Cys145). The Cys145 is observed to lie near the C8 and C9 centers of the ligand molecule which encouraged us to modify these centers by addition of -COOH and -COOCH 3 functional groups. The modified ligand molecules (Mod I and Mod II) are found to have higher binding energy towards the inhibitor binding domain compared to their counter partners (Mod I / and Mod II / ). During the simulation of Modified ligand molecules with SARS-CoV-2 Mpro, the Cys145 is observed to be perfectly positioned near the C24 atom of ligand molecules thus suggesting the possibility of a nucleophilic attack to this site during catalysis. The W1 water center hydrogen-bonded to His41 ND1 seems to be crucial for the recognition of ligand molecules to catalytic dyad and Asp187 (located on the connecting loop between Domain I and Domain III) through Mod-I/II(O25)‧‧‧Cys145‧‧‧His41 NE2/ND1 ‧‧‧W1‧‧‧Asp187 OD2 interaction. So, from the results, it may be presumed that the Asp187 may function as the third catalytic residue where W1 water center maintains the structural recognition of this catalytic triad.