Elucidating biophysical basis of binding of inhibitors to SARS-CoV-2 main protease by using molecular dynamics simulations and free energy calculations.

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 (Mpro or 3CLpro) by using all-atom molecular dynamics simulations and free energy calculations. We calculated the total binding free energy (ΔGbind) of both inhibitors and further decomposed ΔGbind 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 (ΔGbind= - 9.05 kcal/mol) compared to Z31792168 (ΔGbind= - 3.25 kcal/mol) against COVID-19 3CLpro. The increase in ΔGbind for α-ketoamide relative to Z31792168 arises due to an increase in the favorable electrostatic and van der Waals interactions between the inhibitor and 3CLpro. Further, we have identified important residues controlling the 3CLpro-ligand binding from per-residue based decomposition of the binding free energy. Finally, we have compared ΔGbind 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.

Introduction

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-to-human 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 2020; Elfiky 2020; 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 (3CLpro) is an essential protein for maintaining the viral life cycle by cleaving the essential polyprotein PP1A to individual functional components. 3CLpro contains 9 α-helices and 13 β-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. 3CLpro plays a crucial role in mediating viral replication and transcription making it an attractive drug target for the treatment of COVID-19.

Figure 1.

(A) The ribbon representation of SARS-CoV-2 3CLpro. The catalytic dyad residues (H41 and C145) are shown in balls and sticks. The inhibitors, α-ketoamide (B) and Z31792168 (C) are represented as sticks.

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 3CLpro 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 3CLpro (Alamri et al., 2020; Abdelli et al., 2020; Khan et al., 2020; Khan 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, α-ketoamide [tert-butyl (1-((S)-1-(((S)-4-(benzylamino)-3,4-dioxo-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 3CLpro and resolved the crystal structure of 3CLpro/α-ketoamide (Zhang et al., 2020). This inhibitor shows a promising binding as well as a pronounced lung tropism and can be administered by the inhalation route (Zhang et al., 2020). In the present work, we have studied the mechanism of binding of this inhibitor to SARS-CoV-2 3CLpro 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-yl-ethanamide) (Fearon et al., 2020) in our study (see Figure 1(C)) and compared with α-ketoamide. Finally, the binding free energy of these two inhibitors is compared with the anti-HIV retroviral drugs, such as lopinavir and darunavir.

Materials and methods

Simulation protocol

The initial coordinates for our molecular dynamics simulations were obtained from the X-ray crystallographic structure of the SARS-CoV-2 3CLpro complexed with the inhibitors α-ketoamide (PDB: 6Y2G) and Z31792168 (PDB: 5R84) (Berman et al., 2002; Zhang et al., 2020). The missing residues in 3CLpro were added by using the Modeller (Fiser & Šali, 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äutler 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.

Dynamic cross-correlation map (DCCM)

The correlated and non-correlated atomic motions of complex protein residues were computed with the help of DCCM (McCammon, 1984; Hünenberger 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α atomic coordinates of each residue. The cross-correlations between the residues were described by the covariance matrix, and the element C of the covariance matrix was calculated using the following equation: where Δr and Δr represent the displacements from the average position of i and j 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.

Principal component analysis (PCA)

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α atomic coordinates were used for this analysis. The atomic fluctuations of Cα-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.

Free energy landscape

The free energy landscape (FEL) calculations were performed by AmberTools19 Cpptraj module of AMBER18 using the below Equation (2) (Frauenfelder et al., 1991); where represents the Boltzmann constant, T is the temperature of each simulated system. N is the population of the ith bin, and N 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.

Energy calculations

In order to estimate the stabilization of binding systems, the binding free energies or affinity of all inhibitors against SARS-CoV-2 3CLpro 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., 2007; Kar et al., 2007; 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 (ΔGbind) of an inhibitor is calculated by utilizing the following equations:

Where ΔEMM, ΔGsolv and -TΔS are the changes in molecular mechanics (MM) energy at the gas phase, desolvation free energy and the conformational entropy upon complexation, respectively. Further, ΔEMM is sum of ΔEinternal (bond, dihedral and angle energies), ΔEelec (electrostatic), and ΔEvdW (van der Waals) energies while ΔGsolv includes ΔGpol (polar solvation energy) and ΔGnp (nonpolar solvation energy). The normal mode analysis (Xu et al., 2011) method was used to calculate the configurational entropy (SMM). 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.

Hydrogen bond criterion

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.

Results and discussion

Structural stability and flexibility analysis

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 3CLpro reached a stable equilibrium after 30 ns while the 3CLpro/α-ketoamide system reached equilibrium after 40 ns. In the case of 3CLpro/Z31792168, the equilibrium is reached after 60 ns.

Figure 2.

Time evolution of root-mean-square deviations (RMSD) of (A) backbone atoms of 3CLpro relative to their respective minimized structure and (B) heavy atoms of inhibitor.

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 3CLpro/α-ketoamide, while the lowest RMSD (1.95 Å) was obtained for 3CLpro/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 α-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).

Table 1.

The average RMSD, solvent accessible surface area (SASA), and radius of gyration for the simulated systems.

System	Backbone RMSD (Å)	SASA (Å2)	Radius of gyration (Rg) (Å)
apo	2.10 (0.01)	14650 (3.5)	22.13 (0.01)
3CLpro/α-ketoamide	2.25 (0.01)	13969 (3.2)	22.03 (0.01)
3CLpro/Z31792168	1.95 (0.01)	14208 (3.6)	21.96 (0.01)

Next, we investigated the flexibility of different parts of 3CLpro by calculating the root mean square fluctuations (RMSFs) of Cα-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-3CLpro, domain II (residues 102-184) shows slightly higher fluctuations, which gets diminished after the ligand binding. In the case of 3CLpro/α-ketoamide, a relatively higher value of RMSF is obtained around residue 50 (domain I) compared to apo or 3CLpro/Z31792168.

Figure 3.

The root-mean-square fluctuations (RMSFs) of Cα atoms for apo (red), 3CLpro/α-ketoamide (green), and 3CLpro/Z31792168 (blue).

The structural compactness of each system was analyzed by estimating the radius of gyration (Rg) from their respective MD trajectories, and the calculated average values are reported in Table 1. A similar Rg 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 3CLpro/α-ketoamide. An intermediate SASA value (14208 Å2) was obtained for 3CLpro/Z31792168.

Dynamic cross-correlation (DCC) analysis

To elucidate the effect of inhibitor-binding on the internal dynamics of 3CLpro, the cross-correlation matrix was calculated by using the coordinates of Cα atoms from MD trajectories, and the dynamic cross-correlation map (DCCM) is displayed in Figure 4. Overall, after the inhibitor binding, a slight increase in the apparent anti-correlated motions is observed for both complexes.

Figure 4.

Dynamic cross-correlation map (DCCM) for (A) apo, (B) 3CLpro/α-ketoamide, and (C) 3CLpro/Z31792168. The extent of correlated and anti-correlated motions are color-coded. Red color shows more anti-correlated motions, and grey color shows correlated motions.

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 3CLpro/α-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 3CLpro/α-ketoamide, more anti-correlated motions are observed in R1 for 3CLpro/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 3CLpro/α-ketoamide compared to 3CLpro/Z31792168. Overall, the binding of α-ketoamide with 3CLpro creates a stable environment near the binding cavity.

Principal component analysis

To further characterize the effect of inhibitor binding on the dynamics of 3CLpro 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 3CLpro 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, 3CLpro/α-ketoamide, and 3CLpro/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 3CLpro. 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 3CLpro/α-ketoamide is more restricted, whereas a wider conformational space is sampled in cases of apo and 3CLpro/Z31792168 (see Figure 5 and Figure S3). This observation suggests that both apo and 3CLpro/Z31792168 systems are more flexible compared to 3CLpro/α-ketoamide, which agrees with the RMSF plot (Figure 3). FELs of apo (Figure 5(A)) and 3CLpro/Z31792168 (Figure 5(C)) are characterized by two dispersed free energy basins, whereas, in the case of 3CLpro/α-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 Supporting Information Figure S4, which was generated via the first mode that shows the largest collective motions of Cα 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 3CLpro/Z31792168 compared to the 3CLpro/α-ketoamide complex.

Figure 5.

Two-dimensional free energy landscapes (FELs) generated by projecting the principal components, PC1 and PC2 for (A) apo, (B) 3CLpro/α-ketoamide, and (C) 3CLpro/Z31792168. The representative structures are shown on the right panel.

Dihedral principal component (dPCA) analysis of loop

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 (φi, ψ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.

Figure 6.

The two-dimensional free energy landscapes (FELs) of (A) apo, (B) 3CLpro/α-ketoamide, and (C) 3CLpro/Z31792168. FEL obtained from dPCA and their representative structures are shown on the right-hand panel.

The apo 3CLpro (Figure 6(A)) has one global minimum (70%) and one local minimum (30%). Similarly, the 3CLpro/α-ketoamide complex (Figure 6(B)) is characterized by two almost equiprobable free energy basins (51% and 49%). The conformational space of 3CLpro/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 3CLpro. Finally, the inhibitor-induced fluctuation of this loop in 3CLpro demands more attention to the design of next-generation drugs.

Binding free energy analysis

To understand the biophysical basis of the molecular recognition of inhibitors by the receptor (3CLpro), the MD/MM-PBSA scheme was employed. It provides different individual components, such as ΔEvdW, ΔEelec, ΔGpol, ΔGnp, and TΔS, to the total binding free energy (ΔGBindsim). 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 3CLpro 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 (ΔEvdW), electrostatic interactions (ΔEelec) and non-polar solvation energy (ΔGnp) favored the binding, while polar solvation free energy (ΔGpol) and the configurational entropy (-TΔS) disfavoured the complexation. For both cases, the net polar contribution (ΔEelec + ΔGpol) is positive, which means that the van der Waals interactions mainly drive the binding.

Figure 7.

Energy components (kcal/mol) for the binding of inhibitors to the 3CLpro. ΔEvdW, van der Waals energy; ΔEelect, electrostatics energy in the gas phase; ΔGpolar, polar solvation energy; ΔGnonpolar, nonpolar solvation energy; TΔSMM, configurational entropy contribution and ΔGbind, total binding energy.

Table 2.

Energetic components of the binding free energy for SARS-CoV-2-inhibitors complexes calculated using MM-PBSA (kcal/mol). Standard errors of the mean (SEM) are provided in parentheses.

Components	α-ketoamide	Z31792168	Lopinavir (Wang, 2020)	Darunavir (Sang et al., 2020)
ΔEvdW	−61.63 (0.10)	−31.61 (0.06)	−20.09 (0.63)	−41.32
ΔEelec	−35.23 (0.13)	−13.64 (0.09)	−52.46 (0.33)	−5.80
ΔGpol	62.94 (0.13)	27.08 (0.09)	46.58 (0.56)	29.01
ΔGnp	−5.31 (0.00)	−2.99 (0.00)	−4.59 (0.02)	−4.75
ΔGsolva	57.63 (0.13)	24.09 (0.09)	41.99 (0.56)	24.26
ΔGpol + elecb	27.71 (0.18)	13.44 (0.13)	−5.88 (0.65)	23.21
ΔEMMc	−96.86 (0.16)	−45.25 (0.11)	−30.56 (0.71)	−47.12
-TΔSd	30.18 (0.81)	17.92 (0.88)	−23.93 (0.01)	NA
ΔGTotale	−39.23 (0.11)	−21.16 (0.01)	−30.56	−22.86
ΔGBindsim	−9.05 (0.82)	−3.25 (0.88)	−6.63 (0.28)	NA
IC50Exp	0.67 ± 0.18 μM	NA	NA	NA

ΔGsolv = ΔGnp + ΔGpol;

ΔGpol + elec = ΔEelec + ΔGpol;

ΔEMM = ΔEvdW + ΔEelec;

ΔS = configuration entropy;

ΔG = EvdW + ΔEelec + ΔGnp + ΔGpol, ΔGBind sim = EvdW + ΔEelec + ΔGnp + ΔGpol - (TΔS).

Further, Table 2 reveals that the estimated binding free energy (ΔGBindsim) of α-ketoamide is higher (-9.05 kcal/mol) than Z31792168 (-3.25 kcal/mol) although in the case of 3CLpro/α-ketoamide, unfavourable contributions from both ΔGpol (62.94 kcal/mol) and -TΔS (30.18 kcal/mol) are relatively higher compared to 3CLpro/Z31792168 (ΔGpol = 27.08 kcal/mol, TΔS = 17.92 kcal/mol). This is because ΔEvdW and ΔEelec are more favorable in 3CLpro/α-ketoamide (ΔEvdW = −61.63 kcal/mol, ΔEelec = −35.23 kcal/mol) than 3CLpro/Z31792168 (ΔEvdW = −31.61 kcal/mol, ΔEelec = −13.64 kcal/mol). Overall, our study suggests that α-ketoamide is more potent against COVID-19 main protease compared to Z31792168.

Next, in our study, the binding affinity of α-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 3CLpro of SARS-CoV-2. Recently, the molecular recognition of lopinavir by the COVID-19 3CLpro 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 α-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 α-ketoamide. In the case of α-ketoamide, the van der Waals interactions energy is more favorable compared to the intermolecular electrostatic interactions. Similarly, in the case of darunavir/3CLpro, 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 3CLpro: α-ketoamide > lopinavir > darunavir > Z31792168 (See Table 2). Therefore, the inhibitor, α-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 receptor-ligand 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 α-ketoamide in comparison to Z31792168 (Met165, His163, and His164). This might be one of the reasons for the higher affinity of α-ketoamide relative to Z31792168 against 3CLpro. 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 α-ketoamide than Z31792168. Overall, the identification of hotspot residues from our analysis can facilitate the discovery of new selective inhibitor against COVID-19 3CLpro.

Figure 8.

Decomposition of ΔGbindsim into contributions from individual residues for (A) 3CLpro/α-ketoamide and (B) 3CLpro/Z31792168.

Table 3.

Decomposition of binding free energy into contributions from individual residues.

Residue	EvdW	Eelec	Gpol	Gnp	Gside_chain	Gbackbone	Gtotal
3CLpro/α-ketoamide
Met165	−2.95	−1.82	1.30	−0.23	−2.03	−1.67	−3.70
Leu27	−2.05	−0.47	0.23	−0.18	−1.79	−0.68	−2.47
His164	−1.80	−4.15	3.70	−0.09	−0.31	−2.03	−2.34
His41	−2.04	−2.49	2.43	−0.22	−2.18	−0.14	−2.32
Glu166	−3.10	−4.14	5.99	−0.51	−0.37	−1.39	−1.76
Cys145	−1.86	−0.56	1.26	−0.16	−1.08	−0.24	−1.32
Pro168	−1.07	−0.02	0.13	−0.23	−0.79	−0.40	−1.19
His163	−0.54	−2.98	2.56	−0.06	−0.96	−0.06	−1.02
3CLpro/Z31792168
Met165	−2.11	−1.07	0.82	−0.20	−1.47	−1.09	−2.56
His163	−0.37	−2.86	1.74	−0.03	−1.43	−0.09	−1.52
His164	−0.95	−2.24	2.12	−0.08	−0.23	−0.92	−1.15

Energetic contributions from the van der Waals (EvdW) and electrostatic interactions (Eelec) as well as polar (Gpol) and nonpolar solvation energy (Gnp) and the total contribution of given residue (Gtotal) for SARS-CoV-2-inhibitor complexes are listed. Gside_chain and Gbackbone represent the side chain and backbone contributions. Only residues with | ΔG | ≥ 1.0 kcal/mol are shown. All values are given in kcal/mol.

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 H-bond occupancies are reported in Table 4. It is observed from Table 4 that in the case of 3CLpro/α-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 3CLpro/Z31792168, the highest H-bond occupancy is obtained for only Glu166 (∼ 35%). This explains why the intermolecular electrostatic interaction is stronger in 3CLpro/α-ketoamide compared to 3CLpro/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 H-bond persists for 3CLpro/α-ketoamide compared to 3CLpro/Z31792168, suggesting the strong bonding of the inhibitor α-ketoamide with 3CLpro.

Table 4.

The hydrogen bonds formed between SARS-CoV-2 and inhibitor and the corresponding average distance and percent determined using the MD trajectories in the MD simulations.

Binding couples		MD
Acceptor	Donor	Avg. Distance (Å)a	Occupancy (%)b
3CLpro/α-ketoamide
Lig@O40	His41@NE2	2.88	43.12
Lig@O22	Glu166@N	2.90	42.53
His164@O	Lig@N38	2.88	41.76
Lig@O48	His163@NE2	2.83	29.97
Glu166@O	Lig@N23	2.87	24.43
Glu166@OE2	Lig@N49	2.84	23.06
Glu166@OE1	Lig@N49	2.84	13.17
Leu27@HD23	Lig@N36	2.82	10.33
Leu27@HD21	Lig@N36	2.83	10.25
Leu27@HD22	Lig@N36	2.84	09.78
Asn142@OD1	Lig@N36	2.82	09.56
Lig@O41	Gly143@N	2.84	07.98
Lig@H20	Gln189@N	2.85	05.10
3CLpro/Z31792168
Lig@O	Glu166@N	2.87	34.55
His164@O	Lig@N	2.84	25.34
Lig@N1	His163@NE2	2.91	09.19

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 α-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 contacts remain in the case of 3CLpro/α-ketoamide, which is in agreement with the finding that the net nonpolar interaction (ΔEvdW + ΔGnp) is more favorable in 3CLpro/α-ketoamide compared to 3CLpro/Z31792168. Overall, our investigations reveal that α-ketoamide is a more potent lead molecule against 2019-nCoV 3CLpro than the FDA approved anti-HIV1 protease inhibitors, such as lopinavir and darunavir which is being recently focused in the treatment of COVID-19.

Figure 9.

The receptor-ligand interaction profile for (A) α-ketoamide and (B) Z31792168. The plots are generated by using Ligplot+. Hydrogen bonds are shown as green dotted lines. Red semicircles show the residues involved in the hydrophobic contacts, and residues involved in hydrogen bonds are represented in green. The inhibitors are described as balls and sticks.

Conclusion

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 α-ketoamide and Z31792168 to COVID-19 main protease (3CLpro) 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 3CLpro-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 α-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 α-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 α-ketoamide has a better binding affinity compared to anti-HIV retroviral drugs, such as darunavir and lopinavir. Overall, α-ketoamide can be used as a lead compound in the development of drugs targeting SARS-CoV-2 3CLpro, and our study might play a useful role for the same.

Disclosure statement

No potential conflict of interest was reported by the authors.