key: cord-337681-579cz2tc
authors: Sk, Md Fulbabu; Roy, Rajarshi; Jonniya, Nisha Amarnath; Poddar, Sayan; Kar, Parimal
title: Elucidating biophysical basis of binding of inhibitors to SARS-CoV-2 main protease by using molecular dynamics simulations and free energy calculations
date: 2020-06-01
journal: J Biomol Struct Dyn
DOI: 10.1080/07391102.2020.1768149
sha: 
doc_id: 337681
cord_uid: 579cz2tc

The recent outbreak of novel “coronavirus disease 2019” (COVID-19) has spread rapidly worldwide, causing a global pandemic. In the present work, we have elucidated the mechanism of binding of two inhibitors, namely α-ketoamide and Z31792168, to SARS-CoV-2 main protease (M(pro) or 3CL(pro)) by using all-atom molecular dynamics simulations and free energy calculations. We calculated the total binding free energy (ΔG(bind)) of both inhibitors and further decomposed ΔG(bind) into various forces governing the complex formation using the Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) method. Our calculations reveal that α-ketoamide is more potent (ΔG(bind)= − 9.05 kcal/mol) compared to Z31792168 (ΔG(bind)= − 3.25 kcal/mol) against COVID-19 3CL(pro). The increase in ΔG(bind) for α-ketoamide relative to Z31792168 arises due to an increase in the favorable electrostatic and van der Waals interactions between the inhibitor and 3CL(pro). Further, we have identified important residues controlling the 3CL(pro)-ligand binding from per-residue based decomposition of the binding free energy. Finally, we have compared ΔG(bind) of these two inhibitors with the anti-HIV retroviral drugs, such as lopinavir and darunavir. It is observed that α-ketoamide is more potent compared to lopinavir and darunavir. In the case of lopinavir, a decrease in van der Waals interactions is responsible for the lower binding affinity compared to α-ketoamide. On the other hand, in the case of darunavir, a decrease in the favorable intermolecular electrostatic and van der Waals interactions contributes to lower affinity compared to α-ketoamide. Our study might help in designing rational anti-coronaviral drugs targeting the SARS-CoV-2 main protease. Communicated by Ramaswamy H. Sarma

The recent outbreak of "coronavirus disease 2019" (COVID-19) poses a severe public health risk. COVID-19 is caused by the novel coronavirus (nCoV) or SARS-CoV-2, which is a highly contagious and pathogenic virus (Aanouz et al., 2020; Pant et al., 2020) . The disease has spread worldwide since its outbreak in the city of Wuhan, China, in 2019 (Boopathi et al., 2020; Hasan et al., 2020; Hui et al. 2020; Sinha et al., 2020; Wu et al., 2020) . So far, the total number of people diagnosed with the COVID-19 viral infection exceeds over 3.7 million, with more than 0.26 million fatalities worldwide. In India, more than 50000 people have been infected, including 1700 deaths (https://www.mohfw.gov.in). The number of positive cases is increasing exponentially due to human-tohuman transmission. The World Health Organisation (WHO) has declared COVID-19 a global pandemic.nCoV is a non-segmented, single-stranded, positive-sense RNA virus that belongs to the beta coronavirus family (Woo et al., 2010; Elmezayen et al., 2020; Enayatkhani et al., 2020; Gupta et al., 2020) . This gene encodes four structural proteins, namely spike glycoprotein (S), a small envelope protein (E), matrix glycoprotein (M), and nucleocapsid protein (N) (Rota et al. 2003; Elfiky & Azzam, 2020) . The S protein attaches with the host receptor by the specific receptor binding domain (RBD) (Wan et al., 2020) , while the N protein binds to RNA in multiple sites to make a helical nucleocapsid structure (Chang et al., 2014; Jin et al., 2020) . Along with these proteins, chymotrypsin-like cysteine protease (3CL pro ) is an essential protein for maintaining the viral life cycle by cleaving the essential polyprotein PP1A to individual functional components. 3CL pro contains 9 a-helices and 13 b-strands (see Figure 1 (A)), creating three domains, domain I (residues 8-101), domain II (residues 102-184), and domain III (residues 201-306). The domains II and III are connected by a long loop (residues 185-200) . The active site of this protease is located between domains I and II, and each domain contributes a single residue to the catalytic dyad, His41 and Cys145, respectively. 3CL pro plays a crucial role in mediating viral replication and transcription making it an attractive drug target for the treatment of COVID-19.

Recently, Jin and co-workers have designed a modified peptide inhibitor for this protease and solved the X-ray crystal structure of the complex (Jin et al., 2020) . The availability of the crystal structure of 3CL pro opens up a lot of opportunities to develop new antiviral drugs. This structure is employed to computationally evaluate the repurposing of different FDA approved drugs against COVID-19 (Das, Sarmah, Lyndem & Singha Roy 2020; Islam et al., 2020; Joshi et al., 2020; Wang, 2020) . Several studies have predicted the binding strength of darunavir, indinavir, like other HIV protease drugs and some other small molecules against 3CL pro (Alamri et al., 2020; Abdelli et al., 2020; Lin et al., 2020; Muralidharan et al., 2020; Sang et al., 2020; Sarma et al., 2020; Umesh et al., 2020; Wahedi et al., 2020) .

Recently, Zhang and co-workers have designed a potent inhibitor, a-ketoamide [tert-butyl (1-((S)-1-(((S)-4-(benzylamino)-3,4dioxo-1-((S)-2-oxopyrrolidin-3-yl)butan-2-yl)amino)-3-cyclopropyl-1-oxopropan-2-yl)-2-oxo-1,2-dihydropyridin-3-yl)carbamate] (see Figure 1 (B)) targeting 3CL pro and resolved the crystal structure of 3CL pro /a-ketoamide . This inhibitor shows a promising binding as well as a pronounced lung tropism and can be administered by the inhalation route . In the present work, we have studied the mechanism of binding of this inhibitor to SARS-CoV-2 3CL pro using molecular dynamics simulations in conjunction with the popular molecular mechanics/Poisson-Boltzmann surface area scheme. We have also considered another inhibitor, Z31792168 (2-cyclohexyl-N-pyridin-3-ylethanamide) (Fearon et al., 2020) in our study (see Figure 1 (C)) and compared with a-ketoamide. Finally, the binding free energy of these two inhibitors is compared with the anti-HIV retroviral drugs, such as lopinavir and darunavir.

The initial coordinates for our molecular dynamics simulations were obtained from the X-ray crystallographic structure of the SARS-CoV-2 3CL pro complexed with the inhibitors a-ketoamide (PDB: 6Y2G) and Z31792168 (PDB: 5R84) (Berman et al., 2002; Zhang et al., 2020) . The missing residues in 3CL pro were added by using the Modeller (Fiser & Sali, 2003; Webb & Sali, 2014) plugin in UCSF Chimera software (Pettersen et al., 2004) . The protonation states of the charged residues were determined using the Propka 3.1 webserver (Olsson et al., 2011) . All simulations were executed by the pmemd.cuda module of AMBER18 (Case et al., 2018) package and analyses were done by the Cpptraj module (Roe & Cheatham, 2013) . We used the latest AMBER ff14SB force field (Maier et al., 2015) to describe the protein structure and the updated generalized Amber force field (GAFF2) (Wang et al., 2004) is used to assign parameters to small molecules. All the missing hydrogen atoms were added by the Leap module of AMBER (Salomon-Ferrer et al., 2013) . The inhibitors were assigned AM1-BCC (Jakalian et al., 2002) charge, which was calculated by utilizing the Antechamber (Wang et al., 2006) module of AMBER18. The systems were solvated in a truncated octahedron periodic box with an explicit TIP3P (Price & Brooks, 2004) water model and a 10 Å buffer distance was considered from the complex along each side.

A suitable integer number of counterions (Na þ ) were added for neutralizing the whole system. The temperature was kept at 300 K and controlled by the Langevin thermostat (Loncharich et al., 1992) . The system pressure was monitored using a Berendsen Barostat (Berendsen et al., 1984) and kept at 1.0 bar. All bond lengths involving hydrogen atoms were constrained by the SHAKE algorithm (Kr€ autler et al., 2001) . We used a time-step of 2.0 fs for the simulation. The particle mesh Ewald summation (PME) (Darden et al., 1993) approach was used to compute the long-range electrostatic interactions. For all cases, the nonbonded cut-off was fixed at 10.0 Å.

Firstly, each complex was optimized using 500 steps of the steepest descent algorithm followed by another 500 cycles of the conjugate gradient scheme. During the minimization, the receptor-inhibitor complexes were restrained to their respective coordinates with a force constant of 2.0 kcal mol À1 Ð À2 . Next, we carried out the minimization without applying any harmonic restraint on the solutes to remove any residual steric clashes. After the minimization, the temperature of each solvated system was gradually increased from 0 K to 300 K at the NVT (canonical) ensemble with a restraint force constant of 2.0 kcal mol À1 Ð À2 acting on the solute atoms. Subsequently, a 50 ps MD simulation with a restraint force constant of 2.0 kcal mol À1 Ð À2 acting on all solute atoms at a constant pressure of 1.0 atm was conducted using Berendsen Barostat (Berendsen et al., 1984 ) at a fixed temperature of 300 K. After 1.0 ns of an equilibration phase, the production simulation was carried for 100 ns at the NPT ensemble. The Cartesian coordinates were stored every 10 ps. Overall, we accumulated 10000 snapshots corresponding to each production simulation.

The correlated and non-correlated atomic motions of complex protein residues were computed with the help of DCCM (McCammon, 1984; H€ unenberger et al., 1995) analysis using Cpptraj (Roe & Cheatham, 2013) module of AmberTools19. Herein, we avoid the apparent correlations between slow side-chain fluctuations and backbone movements. To reduce the statistical noise, we considered only C a atomic coordinates of each residue. The cross-correlations between the residues were described by the covariance matrix, and the element C ij of the covariance matrix was calculated using the following equation:

where Dr i and Dr j represent the displacements from the average position of i th and j th atoms with respect to time, respectively. The angular bracket indicates the time average over the entire trajectory. The cross-correlated values vary between -1 and 1. The positive value represents the positively correlated motions, while negative values represent the anti-correlated motions. We have considered the final 90 ns production simulation trajectories for this analysis.

PCA or principal component analysis (Ichiye & Karplus, 1991) gives us detailed information about residual movements and functional significance of each residue. Similar to DCCM analysis, only C a atomic coordinates were used for this analysis. The atomic fluctuations of C a -atoms of each residue form a covariance matrix, as described in the DCC analysis. The diagonalization of the covariance matrix gives us the othrogonal eigenvectors and the corresponding eigenvalues. The eigenvectors indicate the directions of the concerted motions, and the eigenvalues describe the amplitude of the motions. These eigenvectors and associate eigenvalues represent the set of principal components (PCs), which may be used to describe the movement characteristics. We also calculated the cosine content via GROMACS (g_covar, g_anaeig, and g_analyze modules) (Hess et al., 2008) of first few PCs to check the statistical convergence significance of each trajectory. The higher value of conformational sampling convergence gives a low value of cosine contents. Our first few PC's cosine contents values lie between 0 to 0.6 for each case, which indicates a high conformational sampling convergence.

The free energy landscape (FEL) calculations were performed by AmberTools19 Cpptraj module of AMBER18 using the below Equation (2) (Frauenfelder et al., 1991) ;

where k B represents the Boltzmann constant, T is the temperature of each simulated system. N i is the population of the i th bin, and N m is the population of the most populated bin. Therefore, bins with no population are an artificial barrier as an indication of the lowest probability. To measure the conformational variability of each system in terms of FELs, we considered the first two principal components (PC1 and PC2) as reaction coordinates.

In order to estimate the stabilization of binding systems, the binding free energies or affinity of all inhibitors against SARS-CoV-2 3CL pro were calculated by the widely used molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method (Kollman et al., 2000; Wang et al., 2001; Kar et al., 2011; Kar & Knecht, 2012a; Kar & Knecht, 2012b; Kar & Knecht, 2012c; Kar & Knecht, 2012d; Kar et al., 2013) . We discussed the MMPBSA scheme in detail in our previous studies (Jonniya et al., 2019; Jonniya & Kar, 2020; Roy et al., 2020; Sk et al., 2020) and here we used the same protocol. Briefly, in the MMPBSA scheme, the binding free energy (DG bind ) of an inhibitor is calculated by utilizing the following equations:

Where DE MM , DG solv and -TDS are the changes in molecular mechanics (MM) energy at the gas phase, desolvation free energy and the conformational entropy upon complexation, respectively. Further, DE MM is sum of DE internal (bond, dihedral and angle energies), DE elec (electrostatic), and DE vdW (van der Waals) energies while DG solv includes DG pol (polar solvation energy) and DG np (nonpolar solvation energy). The normal mode analysis (Xu et al., 2011) method was used to calculate the configurational entropy (S MM ). Due to the high computational cost, we considered only 50 configurations from the last 20 ns simulations for entropy calculations. We were also performed the per-residue decomposition of binding energy by molecular mechanics generalize-Born surface area (MM-GBSA) scheme to know the individual residual contributions of each residue. All the parameters used in this calculation were developed by Onufreiv and Bashford (Onufriev et al., 2004) . In this study, 2000 structural frames were applied to run the MM-PBSA calculation.

The hydrogen bonds between the inhibitor and receptor were analyzed using the Cpptraj module in AmberTools19 (Roe & Cheatham, 2013) . The formation of hydrogen bond defined by a distance and an angle cut-off of 3.5 Å and !120 , respectively, between the donor (D) and acceptor (A) atom.

To verify the convergence of our simulations, we estimated the root-mean-square deviations (RMSDs) of backbone atoms relative to the initial structure, and the temporal RMSD of all systems is shown in Figure 2 . From Figure 2 (A), it is clear that apo 3CL pro reached a stable equilibrium after 30 ns while the 3CL pro /a-ketoamide system reached equilibrium after 40 ns. In the case of 3CL pro /Z31792168, the equilibrium is reached after 60 ns.

The average RMSDs for all simulated systems were calculated and reported in Table 1 . The average RMSD was found to vary between 1.95 Å and 2.25 Å for all cases. The highest deviation (2.25 Å) was obtained for 3CL pro /a-ketoamide, while the lowest RMSD (1.95 Å) was obtained for 3CL pro /Z31792168. We also calculated the RMSD of each ligand, and the time evolution of ligand RMSD is shown in Figure 2 (B). As can be seen from Figure 2 (B), two inhibitors reached equilibrium after 50 ns of MD simulations. In the case of a-ketoamide, the frequency distribution of RMSD is characterized by a sharp peak at 2.4 Å while two peaks (0.23 Å and 1.1 Å) are obtained for Z31792168. The K-mean cluster analysis of the ligand suggests two dominant conformations for both ligands, and these conformations correspond to dominant poses in the binding cavity (see Supporting Information Figure S1 ).

Next, we investigated the flexibility of different parts of 3CL pro by calculating the root mean square fluctuations (RMSFs) of C a -atoms for all systems and is shown in Figure 3 . It is evident from Figure 3 that all three systems display more or less similar fluctuations. In the case of apo-3CL pro , domain II (residues 102-184) shows slightly higher fluctuations, which gets diminished after the ligand binding. In the case of 3CL pro /a-ketoamide, a relatively higher value of RMSF is obtained around residue 50 (domain I) compared to apo or 3CL pro /Z31792168.

The structural compactness of each system was analyzed by estimating the radius of gyration (R g ) from their respective MD trajectories, and the calculated average values are reported in Table 1 . A similar R g is obtained for all systems. Finally, the solvent-accessible surface area (SASA) that indicates the degree of solvent exposure was also calculated and reported in Table 1 . It is evident from Table 1 that SASA values vary between 13969 Å 2 and 14650 Å 2 for all systems. The highest value (14650 Å 2 ) was obtained for apo, while the lowest SASA value (13969 Å 2 ) was reported for 3CL pro / a-ketoamide. An intermediate SASA value (14208 Å 2 ) was obtained for 3CL pro /Z31792168.

To elucidate the effect of inhibitor-binding on the internal dynamics of 3CL pro , the cross-correlation matrix was calculated by using the coordinates of C a atoms from MD trajectories, and the dynamic cross-correlation map (DCCM) is displayed in Figure 4 . Overall, after the inhibitor binding, a Figure 3 . The root-mean-square fluctuations (RMSFs) of C a atoms for apo (red), 3CL pro /a-ketoamide (green), and 3CL pro /Z31792168 (blue).

slight increase in the apparent anti-correlated motions is observed for both complexes.

In the case of the apo system, domain I (residues 8-101) displays correlated motions while domain II (or R1 in Figure 4 (A)) (residues 102-184) shows highly anti-correlated movements. However, domain III (residues 201-306) has a lower strength of anti-correlation motions relative to domain I. In the case of 3CL pro /a-ketoamide (Figure 4(B) ), the region R1, which represents the binding cavity, displays a significant reduction in anti-correlated and an increase in correlated motions compared to apo. In contrast to 3CL pro /a-ketoamide, more anti-correlated motions are observed in R1 for 3CL pro /Z31792168 (Figure 4(C) ). Further, it is evident from Figure 4 (A-C) that the inhibitor binding increases the strength of anti-correlated motions in the region R2, which indicates the residual motion of domain III. In region R3, the anti-correlated movements are diminished in 3CL pro /a-ketoamide compared to 3CL pro / Z31792168. Overall, the binding of a-ketoamide with 3CL pro creates a stable environment near the binding cavity.

To further characterize the effect of inhibitor binding on the dynamics of 3CL pro and extract the structural variations in detail, PCA was applied to the coordinates of the backbone atoms of all three systems. The PCA results in terms of eigenvectors (3 N, 3 Â 306 ¼ 918) versus eigenvalues obtained by diagonalization of the atomic fluctuations covariance matrix for complexes and apo 3CL pro are shown in Figure S2 (Supporting Information). The collective motions of the localized fluctuations can be defined by the first few principal components. The first 10 eigenvectors capture ⁓87%, ⁓88%, and ⁓ 89% of the total motion in apo, 3CL pro /a-ketoamide, and 3CL pro /Z31792168.

In general, the two largest principal components (PC1 and PC2) account for more than 70% of the overall fluctuations for apo and the inhibitor-bound 3CL pro . The two-dimensional free energy landscapes (FELs) of all three systems are shown in Figure 5 (A-C). It is notable that the conformational space sampled in the case of 3CL pro /a-ketoamide is more restricted, whereas a wider conformational space is sampled in cases of apo and 3CL pro /Z31792168 (see Figure 5 and Figure S3 ). This observation suggests that both apo and 3CL pro /Z31792168 systems are more flexible compared to 3CL pro /a-ketoamide, which agrees with the RMSF plot ( Figure 3 ). FELs of apo ( Figure 5(A) ) and 3CL pro /Z31792168 (Figure 5(C) ) are characterized by two dispersed free energy basins, whereas, in the case of 3CL pro /a-ketoamide ( Figure 5(B) ), we have obtained a single and stable global free energy minimum confined within a basin on the FEL. Porcupine plots are shown in Figure S4 , which was generated via the first mode that shows the largest collective motions of C a atoms, which has been mapped onto the average structure using the Normal Modes Analysis plugin available in VMD (Humphrey et al., 1996) . To visualize the detailed movements of each protein, corresponding structures were selected based on K-means (Hartigan & Wong, 1979) clusters of PCA distributions. The structures were obtained from cluster centers and superimposed with each other. Overall, it is observed that the mobility is higher for apo and 3CL pro / Z31792168 compared to the 3CL pro /a-ketoamide complex.

Cartesian coordinates based PCA does not provide the fully correct internal and overall dynamics. However, dPCA gives us information about internal and overall conformational dynamics. To further discover the conformational alterations of the loop (residues 185-200) connecting domains II and III, dihedral PCA (dPCA) was performed based on the dihedral angles (u i , w i ) of the peptide group in the loop (see Figure 6 ). The representative structures are obtained from the K-means clustering analysis and shown in Figure 6 . Overall, Figure 6 indicates diverse loop conformations for each system. The apo 3CL pro (Figure 6 (A)) has one global minimum (70%) and one local minimum (30%). Similarly, the 3CL pro / a-ketoamide complex (Figure 6(B) ) is characterized by two almost equiprobable free energy basins (51% and 49%). The conformational space of 3CL pro /Z31792168 is quite wider, as is evident from Figure 6 (C). It is characterized by a global free energy basin (46%) and two local free energy minima (41% and 13%). In this analysis, we confirmed that this loop plays an important role in the inhibitor binding as well as controls the dynamical nature of 3CL pro . Finally, the inhibitorinduced fluctuation of this loop in 3CL pro demands more attention to the design of next-generation drugs.

To understand the biophysical basis of the molecular recognition of inhibitors by the receptor (3CL pro ), the MD/MM-PBSA scheme was employed. It provides different individual components, such as DE vdW , DE elec , DG pol , DG np , and TDS, to the total binding free energy (DG Bind sim ). For binding free energy calculations of both complexes, 2000 structures were taken from the stable region, and 50 structures were considered for the entropy calculations.

A summary of binding components in the binding free energy of inhibitors with 3CL pro is shown graphically in Figure  7 , and the data are listed in Table 2 . It is evident from Table 2 that the intermolecular van der Waals (DE vdW ), electrostatic interactions (DE elec ) and non-polar solvation energy (DG np ) favored the binding, while polar solvation free energy (DG pol ) and the configurational entropy (-TDS) disfavoured the complexation. For both cases, the net polar contribution (DE elec þ DG pol ) is positive, which means that the van der Waals interactions mainly drive the binding. Further, Table 2 reveals that the estimated binding free energy (DG Bind sim ) of a-ketoamide is higher (-9.05 kcal/mol) than Z31792168 (-3.25 kcal/mol) although in the case of 3CL pro /a-ketoamide, unfavourable contributions from both DG pol (62.94 kcal/mol) and -TDS (30.18 kcal/mol) are relatively higher compared to 3CL pro /Z31792168 (DG pol ¼ 27.08 kcal/ mol, TDS ¼ 17.92 kcal/mol). This is because DE vdW and DE elec are more favorable in 3CL pro /a-ketoamide (DE vdW ¼ À61.63 kcal/mol, DE elec ¼ À35.23 kcal/mol) than 3CL pro / Z31792168 (DE vdW ¼ À31.61 kcal/mol, DE elec ¼ À13.64 kcal/mol). Overall, our study suggests that a-ketoamide is more potent against COVID-19 main protease compared to Z31792168.

Next, in our study, the binding affinity of a-ketoamide was further evaluated and compared with the FDA approved anti-HIV protease inhibitors, such as lopinavir and darunavir, which has been reported as potent drugs against 3CL pro of SARS-CoV-2. Recently, the molecular recognition of lopinavir by the COVID-19 3CL pro has been investigated using the MMPBSA scheme (Wang, 2020) , and the binding free energy of lopinavir was found to be lesser (-6.63 kcal/mol) than a-ketoamide (-9.05 kcal/mol) (see Table 2 ). It is further revealed that the electrostatic interaction (-52.46 kcal/mol) favored the complex formation more compared to the van der Waals interactions (-20.09 kcal/mol). This is in contrast to what has been observed for a-ketoamide. In the case of a-ketoamide, the van der Waals interactions energy is more favorable compared to the intermolecular electrostatic interactions. Similarly, in the case of darunavir/3CL pro , the van der Waals force is more favorable (-41.32 kcal/mol) than the electrostatic force (-5.80 kcal/mol) (Sang et al., 2020) . Our study reports that the binding affinity decreases in the following order against COVID-19 3CL pro : a-ketoamide > lopinavir > darunavir > Z31792168 (See Table 2 ). Therefore, the inhibitor, a-ketoamide may be considered as a lead compound in the discovery of rational drugs against COVID-19.

Subsequently, we explored the critical residues involved in the receptor-ligand binding by performing the per-residue decomposition of free energy using MMGBSA. The receptorligand interaction spectra were shown in Figure 8 . A hotspot residue is considered in the MM-GBSA interaction energy when it is higher than -1.0 kcal/mol and reported in Table 3 . It is evident from Table 2 that a more significant number of hotspot residues (Met165, Leu27, His164, His41, Glu166, Cys145, Pro168, and His163) contributes favorably to the binding of a-ketoamide in comparison to Z31792168 (Met165, His163, and His164) . This might be one of the reasons for the higher affinity of a-ketoamide relative to Z31792168 against 3CL pro . It can further be observed from Table 3 that the catalytic dyad, His41 (-2.32 kcal/mol), and Cys145 (-1.32 kcal/mol) contribute more significantly to the binding of a-ketoamide than Z31792168. Overall, the identification of hotspot residues from our analysis can facilitate the discovery of new selective inhibitor against COVID-19 3CL pro .

Subsequently, to compliment the above results, we have performed the hydrogen bond (H-bond) analysis for both the complexes on their respective MD trajectories and the Hbond occupancies are reported in Table 4 . It is observed from Table 4 that in the case of 3CL pro /a-ketoamide, residues His41, Glu166, and His164 ($ 42%) form H-bond with the inhibitor with an occupancy of > 40% during the simulation. However, in the case of 3CL pro /Z31792168, the highest Hbond occupancy is obtained for only Glu166 ($ 35%). This explains why the intermolecular electrostatic interaction is stronger in 3CL pro /a-ketoamide compared to 3CL pro / Z31792168. Besides, the number of H-bond as a function of simulations time is shown in Figure S5 (Supporting Information), which also ensures that a higher number of Hbond persists for 3CL pro /a-ketoamide compared to 3CL pro / Z31792168, suggesting the strong bonding of the inhibitor a-ketoamide with 3CL pro .

A detailed interaction profile of residues involving H-bond and hydrophobic interactions is also computed using LigPlot (Wallace et al., 1995) and shown in Figure 9 . It suggests that one of the catalytic dyad residues, His41 plays a significant role in forming a strong H-bond with a-ketoamide. Also, the hydrophobic contacts resulting from the hydrophobic residues for both the complexes throughout the simulations were computed and shown in Supporting Information Figure  S6 . It also suggests the highest number of hydrophobic .15 a Energetic contributions from the van der Waals (E vdW ) and electrostatic interactions (E elec ) as well as polar (G pol ) and nonpolar solvation energy (G np ) and the total contribution of given residue (G total ) for SARS-CoV-2-inhibitor complexes are listed. G side_chain and G backbone represent the side chain and backbone contributions. Only residues with j DG j ! 1.0 kcal/mol are shown. All values are given in kcal/mol. contacts remain in the case of 3CL pro /a-ketoamide, which is in agreement with the finding that the net nonpolar interaction (DE vdW þ DG np ) is more favorable in 3CL pro /a-ketoamide compared to 3CL pro /Z31792168. Overall, our investigations reveal that a-ketoamide is a more potent lead molecule against 2019-nCoV 3CL pro than the FDA approved anti-HIV1 protease inhibitors, such as lopinavir and darunavir which is being recently focused in the treatment of COVID-19.

The recent outbreak of COVID-19 has caused a severe strain in the public health system in many countries. COVID-19 can cause mild to severe illness. The current situation demands an efficacious vaccine or novel antiviral drugs targeting COVID-19. Herein, we have studied the mechanism of binding of two potential inhibitors, namely a-ketoamide and Z31792168 to COVID-19 main protease (3CL pro ) by using an atomistic molecular dynamics simulation of 100 ns in conjunction with the widely used molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) scheme. Our study shows that the 3CL pro -inhibitor complexation is favored by the intermolecular van der Waals and electrostatic interactions as well as nonpolar solvation free energy. We have also demonstrated that the inhibitor a-ketoamide is more potent compared to Z31792168 due to an increased favorable contribution from the intermolecular van der Waals and electrostatic interactions relative to Z31792168. Furthermore, in the case of a-ketoamide, the nonpolar component of the solvation free energy is also slightly more favorable compared to Z31792168. We have also identified the hotspot residues controlling the receptor-ligand binding. Finally, our study also reveals that a-ketoamide has a better binding affinity compared to anti-HIV retroviral drugs, such as darunavir and lopinavir. Overall, a-ketoamide can be used as a lead compound in the development of drugs targeting SARS-CoV-2 3CL pro , and our study might play a useful role for the same.

No potential conflict of interest was reported by the authors. 

Moroccan Medicinal plants as inhibitors of COVID-19: Computational investigations

In silico study the inhibition of Angiotensin converting enzyme 2 receptor of COVID-19 by Ammoides verticillata components harvested from western Algeria

Pharmacoinformatics and molecular dynamic simulation studies reveal potential inhibitors of SARS-CoV-2 main protease 3CLpro

Molecular dynamics with coupling to an external bath

The protein data bank

Mechanism of Action, Antiviral drug promises and rule out against its treatment

The SARS coronavirus nucleocapsid protein-forms and functions

Particle mesh Ewald: An NÁ log (N) method for Ewald sums in large systems

An investigation into the identification of potential inhibitors of SARS-CoV-2 main protease using molecular docking study

Natural products may interfere with SARS-CoV-2 attachment to the host cell

SARS-CoV-2 RNA dependent RNA polymerase (RdRp) targeting: An in silico perspective

Novel Guanosine Derivatives against MERS CoV polymerase: An in silico perspective

Drug repurposing for coronavirus (COVID-19): In silico screening of known drugs against coronavirus 3CL hydrolase and protease enzymes

Reverse vaccinology approach to design a novel multi-epitope vaccine candidate against COVID-19: An in silico study

PanDDA analysis of COVID-19 main protease against the DSI-poised Fragment Library

Modeller: Generation and refinement of homology-based protein structure models

The energy landscapes and motions of proteins

In-silico approaches to detect inhibitors of the human severe acute respiratory syndrome coronavirus envelope protein ion channel

Algorithm AS 136: A k-means clustering algorithm

A review on the cleavage priming of the spike protein on coronavirus by angiotensin-converting enzyme-2 and furin

GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation

The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health-The latest 2019 novel coronavirus outbreak in Wuhan

VMD: Visual molecular dynamics

Fluctuation and cross-correlation analysis of protein motions observed in nanosecond molecular dynamics simulations

Collective motions in proteins: A covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations

A Molecular Modeling Approach to Identify Effective Antiviral Phytochemicals against the Main Protease of SARS-CoV-2

Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

Investigating specificity of the anti-hypertensive inhibitor WNK463 against With-No-Lysine kinase family isoforms via multiscale simulations

Investigating Phosphorylation-Induced Conformational Changes in WNK1 Kinase by Molecular Dynamics Simulations

Discovery of potential multi-target-directed ligands by targeting host-specific SARS-CoV-2 structurally conserved main protease$

Energetic basis for drug resistance of HIV-1 protease mutants against amprenavir

Energetics of mutation-induced changes in potency of lersivirine against HIV-1 reverse transcriptase

Mutation-induced loop opening and energetics for binding of tamiflu to influenza N8 neuraminidase

Origin of decrease in potency of darunavir and two related antiviral inhibitors against HIV-2 compared to HIV-1 protease

Importance of polar solvation for cross-reactivity of antibody and its variants with steroids

Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease

Dispersion Terms and Analysis of size-and charge dependence in an enhanced Poisson-Boltzmann approach

Systematic study of the boundary composition in Poisson Boltzmann calculations

Targeting SARS-CoV-2: A systematic drug repurposing approach to identify promising inhibitors against 3C-like Proteinase and 2'-O-RiboseMethyltransferase

Identification of Chymotrypsin-like Protease Inhibitors of SARS-CoV-2 Via Integrated Computational Approach

Calculating structures and free energies of complex molecules: Combining molecular mechanics and continuum models

A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations

Molecular modeling evaluation of the binding effect of ritonavir, lopinavir and darunavir to severe acute respiratory syndrome coronavirus 2 proteases

Langevin dynamics of peptides: The frictional dependence of isomerization rates of Nacetylalanyl-N'-methylamide

ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB

Protein dynamics

Computational studies of drug repurposing and synergism of lopinavir, oseltamivir and ritonavir binding with SARS-CoV-2 Protease against COVID-19

PROPKA3: Consistent treatment of internal and surface residues in empirical pKa Predictions

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Peptide-like and small-molecule inhibitors against Covid-19

UCSF Chimera-a visualization system for exploratory research and analysis

A modified TIP3P water potential for simulation with Ewald summation

PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data

Characterization of a novel coronavirus associated with severe acute respiratory syndrome

Investigating conformational dynamics of Lewis Y Oligosaccharides and elucidating blood group dependency of cholera using molecular dynamics

An overview of the Amber biomolecular simulation package

Insight derived from molecular docking and molecular dynamics simulations into the binding interactions between HIV-1 protease inhibitors and SARS-CoV-2 3CLpro

In-silico homology assisted identification of inhibitor of RNA binding against 2019-nCoV N-protein (N terminal domain)

An in-silico evaluation of different Saikosaponins for their potency against SARS-CoV-2 using NSP15 and fusion spike glycoprotein as targets

Exploring the potency of currently used drugs against HIV-1 protease of subtype D variant by using multiscale simulations

Identification of new anti-nCoV drug chemical compounds from Indian spices exploiting SARS-CoV-2 main protease as target

Stilbene-based natural compounds as promising drug candidates against COVID-19

LIGPLOT: A program to generate schematic diagrams of protein-ligand interactions

Receptor recognition by the novel coronavirus from wuhan: An analysis based on decade-long structural studies of SARS coronavirus

Fast identification of possible drug treatment of coronavirus disease-19 (COVID-19) through computational drug repurposing study

Antechamber: An accessory software package for molecular mechanical calculations

Automatic atom type and bond type perception in molecular mechanical calculations

Development and testing of a general amber force field

Comparative protein structure modeling using MODELLER

Coronavirus genomics and bioinformatics analysis. Viruses

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: A modelling study

Fast and accurate computation schemes for evaluating vibrational entropy of proteins

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved a-ketoamide inhibitors