An insight into the interaction between α-ketoamide- based inhibitor and coronavirus main protease: A detailed in silico study.

The search for therapeutic drugs that can neutralize the effects of COVID-2019 (SARS-CoV-2) infection is the main focus of current research. The coronavirus main protease (Mpro) is an attractive target for anti-coronavirus drug design. Further, α-ketoamide is proved to be very effective as a reversible covalent-inhibitor against cysteine proteases. Herein, we report on the non-covalent to the covalent adduct formation mechanism of α-ketoamide-based inhibitor with the enzyme active site amino acids by QM/SQM model (QM = quantum mechanical, SQM = semi-empirical QM). To uncover the mechanism, we focused on two approaches: a concerted and a stepwise fashion. The concerted pathway proceeds via deprotonation of the thiol of cysteine (here, Cys145 SγH) and simultaneous reversible nucleophilic attack of sulfur onto the α-ketoamide warhead. In this work, we propose three plausible concerted pathways. On the contrary, in a traditional two-stage pathway, the first step is proton transfer from Cys145 SγH to His41 Nδ forming an ion pair, and consecutively, in the second step, the thiolate ion attacks the α-keto group to form a thiohemiketal. In this reaction, we find that the stability of the tetrahedral intermediate oxyanion/hydroxyl group plays an important role. Moreover, as the α-keto group has two faces Si or Re for the nucleophilic attack, we considered both possibilities of attack leading to S- and R-thiohemiketal. We computed the structural, electronic, and energetic parameters of all stationary points including transition states via ONIOM and pure DFT method. Additionally, to characterize covalent, weak noncovalent interaction (NCI) and hydrogen-bonds, we applied NCI-reduced density gradient (NCI-RDG) methods along with Bader's Quantum Theory of Atoms-in-Molecules (QTAIM) and natural bonding orbital (NBO) analysis.

Introduction

A new human coronavirus (HCoV), SARS-CoV-2, began spreading in December 2019 in Wuhan City, China causing pneumonia called COVID-19. All the affected countries are trying to prevent the spread of the new mysterious coronavirus. Despite this, the number of infected persons is continually increasing each day. As the viral maturation mainly relies on the Mpro's activity, the coronavirus main protease (Mpro) (also called 3C like protease or 3CLpro) is also an attractive target for anti-CoV drug design. Recently, it is reported that deadlier coronavirus is highly similar with 96% sequence identity to that of the SARS-CoV Mpro [1]. Due to this worldwide pandemic, considerable research is focused to design many inhibitors for the novel coronavirus [[2], [3], [4], [5], [6], [7], [8], [9], [10]]. For drug designing, the reversible covalent inhibitors are generally more preferable than irreversible ones due to its risks in long term treatment and side effects [[11], [12], [13], [14]]. In addition to the development of novel drugs, researchers look for effective substances among drugs that are already approved or developed for other purposes.

Covalent inhibitors, that bind with cysteine proteases, are generally substrate with an electrophilic center, including Michael acceptor (—CHCH—CO—), aldehydes (—CHO), nitriles (—C N), vinyl sulfone (—SO2 —CHCH2), and α-ketoamide (—CO—CO—NH—) [[15], [16], [17], [18], [19], [20], [21], [22]]. Of them, α-ketoamide proved to be very effective against cysteine proteases [23]. The advantage of the α-ketoamides is that their warhead can interact with the catalytic center of the target proteases by two hydrogen-bonding interactions, namely, via the α-keto oxygen and amide oxygen [1]. Whereas, warheads like aldehydes or Michael acceptors interacts with the catalytic center by only one hydrogen bond, the α-ketoamide motif is emerging as the most promising group with a possible reversible inhibition of the enzyme's activity better than the other investigated C-terminus warheads [23,24] The α-keto group of the α-ketoamide forms a reversible tetrahedral adduct with the active cysteine enzyme. For example, calpain I inhibitors form tetrahedral hemiketal adduct by the cysteine residue (Cys115) onto the α-keto group of the α-ketoamide [25]. Chen J-C et al. [26] synthesized a series of small molecules bearing an α-ketoamide warhead and evaluated for their ability to inhibit cathepsin S. In another study, Ota. E. et al. [27] showed that thienyl-substituted α-ketoamide is a less hydrophobic reactive group for photo-affinity labeling. Recently, Hilgenfeld and his coworkers [1,28], designed a series of α-ketoamide as inhibitors of SARS-CoV-2 Mpro.

Whereas, few theoretical studies about the reaction mechanism of a cysteine protease with Michael acceptor, aldehydes, and nitriles and similar many researches are present in literature [[29], [30], [31], [32], [33], [34], [35], [36], [37], [38], [39]], including the interesting work of Zanetti-Polzi et al. [40] and Arafet et al. [41], but precise reaction mechanism of α-ketoamide and cysteine is rarely reported [42]. It is essential to understand the inhibitory action of these drugs to design a more promising drug against coronavirus. Commonly, it is believed that in these mechanisms, the first step is deprotonation of thiol of cysteine with a basic side chain, usually by a histidine residue followed by a nucleophilic attack of the anionic sulfur onto the α-keto group of the inhibitor forming an S—C covalent bond and an oxyanion. In the last step, the oxyanion gets stability by hydrogen bonds and/or by accepting a proton. To date, few studies have been made to understand how the keto group binds with the cysteine proteases. In particular, to the best of our knowledge, the literature lacks the detailed quantum chemical study of α-ketoamide based inhibitors.

In the present paper, for the first time, we will investigate the detailed reaction mechanism for the conversion of noncovalent adduct to covalent adduct by the interaction of Mpro of coronavirus and two α-ketoamide compounds with the help of semi-empirical and density functional (DFT) studies. We extracted these compounds from the pdb file 6y2f and 5n19 (henceforth, we name them as S—In and R-In). We considered two approaches for the reaction: i) a concerted and ii) a stepwise mechanism. In the first case, the nucleophilic attack and such proton transfer occur in a single step. Herein, we considered four pathways, for example, Pathways 1–3, as shown in Scheme 1 . On the other hand, the second case involves a two-step process (Pathway 4). In the first step, a zwitterionic ion pair (CysS−/HisN+) is formed from the fast proton transfer from cysteine (Cys145) to nearby histidine residue (His41) (Scheme 2 ). In the second step, the thiolate anion of cysteine Cys145 S− will then attack the α-keto group of the inhibitor to form a C—S covalent bond and an oxyanion, followed by abstraction of a proton from the nearest base histidine (His41) (Scheme 2). Moreover, as shown in Fig. 1 , the α-keto group has Re and Si faces for the nucleophilic reaction. Thus, addition from Re and Si faces leads to R- and S- stereochemistry adducts. All the adducts, formed by pathways 1–4, caused by the attack from Si face and hence have S stereochemistry. Additionally, using SARS-CoV Mpro, we also considered the attack of Cys S from Re face leading to R-thiohemiketal. However, for the Re face reaction, we only examined the reaction by thiolate nucleophile (Cys S−) (Scheme 3 ). Finally, to find out the stability of R- and S- thiohemiketal and the oxyanion/hydroxyl group, we will apply modern theoretical tools like NCI-reduced density gradient (NCI-RDG) methods along with Bader's Quantum Theory of Atoms-in-Molecules (QTAIM) and natural bonding orbital (NBO) analysis.

Scheme 1

Proposed concerted reaction mechanism between the active site of SARS-CoV-2 Mpro and α-ketoamide inhibitor-warhead from Si face based on the present calculations.

Scheme 2

Proposed two-step reaction mechanism, according to pathway 4, between active site of SARS-CoV-2 Mpro with α-ketoamide inhibitor-warhead from Si face based on the present calculations.

Fig. 1

a) Enantioselective reaction between SARS-CoV Mpro and α-ketoamide. b) Possible confirmation of α-ketoamide. The labels on conformation C3 indicates NBO charges in e unit. c) MEP plot of α-ketoamide d) Isosurface plots of HOMO and LUMO of α-ketoamide.

Scheme 3

Proposed two-step reaction mechanism between the active site of SARS-CoV Mpro with α-ketoamide inhibitor-warhead from the Re face to generate R-thiohemiketal based on the present calculations.

Method and computational studies

Here, we applied ONIOM-QM/SQM model (QM = quantum mechanical; SQM = semi-empirical QM) treating the protease and inhibitor computationally. In this model, we have considered a large number of atoms from the nearest interacting amino acids as shown in Fig. 2 . Here, to fulfill the need, we consider only the active species and key hydrogen bonding interactions. We have used hybrid B3LYP [43,44] and Semi-empirical PM6 [45] level of theory (ONIOM-B3LYP/6-31G(d):PM6) [46] for optimization. To avoid the expensive QM/SQM molecular dynamics method, here we explored the reaction mechanism in static conditions. For the more accurate single-point energy calculation, we used dispersion corrected ωB97X-D/6-31G(d,p):PM6 level of theory on the B3LYP/6-31G(d):PM6 optimized geometries. Hydrogens were used as link atoms. We performed all calculations unless indicated otherwise, by using the Gaussian 09 suite of programs [47]. The effect of the low polarity protein environment, on the active-site model, was estimated by performing single-point calculations using the polarizable continuum model (PCM) with UAKS radii [48] method at the same level of theory as the geometry optimizations. The dielectric constant was set to the standard value ε = 4 to represent this environment [49].

Fig. 2

Schematic diagram indicating the atoms of the S- and R-thiohemiketal adducts taken from pdb code a) 6y2f and b) 5n19. Green lines indicate hydrogen bonds. Not all atoms are shown for the sake of clarity. Red rectangles indicate the α-ketoamide inhibitor-warhead of interest in the present paper. The figure is drawn with ligplot [50]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The QM/SQM energies were calculated using the following equation as the separate energies of a QM calculation on the QM region (Emodel, QM) and an SQM calculation on the whole system (Ereal, SQM) minus an SQM calculation on the QM region (Emodel,MM). The use of such ONIOM model is capable of optimizing and characterizing all critical structures properly along with the energy profile.

The energy is calculated using the following equation:

To explore the activation energy barrier of each step and the energy profile, we optimized all the intermediate and transition state (TS) geometries. We found no imaginary frequencies for the geometries of minima indicating that the obtained stationary points are the local minimum at potential energy surface (PES). On the other hand, only one imaginary frequency at the TS geometry demonstrated that this geometry is a first-order saddle point on the PES. To verify the reliability of the TSs, intrinsic reaction coordinates (IRC) were also computed.

Natural bond orbital (NBO) [51] analysis at B3LYP/6-31G(d) level was carried out to understand the orbital interactions and charge delocalization for the title complexes. In such analysis, NBOs are first defined for each covalent bond, lone pair, and anti-bonding orbital by using the molecular orbitals obtained by quantum-chemical calculations, and the orbital interaction energies are subsequently computed for all possible interactions between electron-donor Lewis-type NBO (i) and acceptor non-Lewis NBO (j). For each electron-donor NBO (i) and electron-acceptor NBO (j), the stabilization energy E2 ij associated with delocalization i → j is calculated aswhere qi is the donor orbital occupancy, εi and εj are orbital energies, and F(i,j) is the off-diagonal NBO Fock matrix element.

Additionally, we used a powerful tool NCI-RDG method to distinguish and visualize weak interactions: stabilized (hydrogen bonding), destabilized (steric repulsion) and delocalized weak (van der Waals). Specifically, the plot of S and sign(λ2)ρ reveals troughs, which signify the different noncovalent interactions. The troughs corresponding to attractive interactions (such as H-bonds) lie in the region of negative sign(λ2)ρ, whereas repulsive interactions lie in the positive region.

The initial structure of coronavirus protease is extracted from the crystallographic data of the ligand-free enzyme of the SARS-CoV-2 (2019-nCoV) main protease (pdb:6y2e). Furthermore, the initial structure of the enzyme-inhibitor complexes are obtained from the crystallographic data of R- and S-adduct of SARS-CoV and SARS-CoV-2, respectively, main protease and α-ketoamide inhibitor using 6y2f and 5n19 pdb files from the protein data bank.

Now in the simplified model, during optimization, there may be possibilities to collapse the rigid enzyme structure leading to unrealistic structure. To prevent this, we kept some atoms fixed based on the parameters found in the crystal structure. The inclusion of such constrain allows obtaining a more realistic rigid structure. However, such constrain optimization sometimes leads to small frequencies, but here, as the number of frozen atoms is very small, in few structures, we obtained imaginary frequencies, that too below 2 cm−1, and thus neglected.

Here, the mechanism determined based on the reduced model of inhibitor that involves only small portion (warhead) of the full molecule bound in the X-ray structure of Mpro can be different from the final mechanism once substituents responsible for the recognition process are attached to the α-ketoamide warhead. Additionally, as recently demonstrated [52] change of the chemical surrounding can meaningfully influence the electrophilic character of the carbon atom and in consequence the reaction barrier.

Results and discussion

Structure of the inhibitor-free enzyme and α-ketoamide

As stated above, we started our investigation with the crystal structure of ligand-free Mpro (pdb:6y2e). Initially, the protonation and crude optimization were performed using the H++ web server at pH 7.3 to mimic the experimental conditions [53]. Then, we optimized a truncated version of the enzyme by using some active site selected residues including His41, Gly143, Ser144, Cys145, His164, and His165. Our objective is to understand the characteristics and relative strength of hydrogen bonds and other non-covalent interactions occurring in the active part of the enzyme in the absence of the warhead of inhibitor.

Fig. S1 (in supporting information) shows that the H-bond network among the residues in the active site creates a pocket large enough to accommodate the α-ketoamide. Also, we detect the other intermolecular interactions as depicted by NCI-RDG isosurface plot (shown in Fig. S1b by red arrows) and important bond critical point (bcp) indicating the inter-residual H-bonds, along with a cage critical point (Fig. S1c, green sphere).

Additionally, to compare the structure of enzyme-inhibitor complex and inhibitor-free enzyme, we have deleted the inhibitor from the enzyme-inhibitor complex (6y2f) and compare the structure with free enzyme active sites as shown in Fig. S2. Aiming to have a complete overview of the α-ketoamide inhibitor, we optimized a simple α-ketoamide N-methyl-2-oxopropanamide. We find that there are three possible isomers of α-ketoamide C1, C2, and C3 (Fig. 1b). Obviously, the C3 isomer is more stable than the other two isomers C1 and C2, in this case, by 9.75 kcal/mol computed at b3lyp/6-31G(d) level for the presence of H-bond of 2.20 Å. The natural population analysis (NPA) shows that the natural charge on amidic carbon is +0.58e while that on αC is +0.52e. However, the molecular electrostatic potential (MEP) plot shows that αC is more electrophilic than that of amidic carbon. The results show that HOMO-LUMO energy gap of this molecule is 4.86 eV and the isosurface plot of HOMO and LUMO are depicted in Fig. 1d. Interestingly, the LUMO of S—In contains the π* orbital for accepting nucleophile, while in R-In, the π* orbital is located on LUMO+1 (See Fig. S3).

Herein, the His41 protonated in ε-position instead of δ was used. The protonation of His41 in ε-position prevents this residue from being directly involved in the proton transfer from Cys145, the role that was demonstrated in studies on proteolysis in Mpro [37].

Reaction mechanism of α-ketoamide with Mpro from Si face to form S-thiohemiketal

Pathway 1: concerted self-catalyzed reaction without water and histidine as a function of the shuttle

Herein, the system contains 68 atoms and as the system is small, we have applied DFT method on the whole system. In this mechanism, first, a non-covalent enzyme-inhibitor complex EI is formed. Then, this complex proceeds to form a tetrahedral complex (Int ) via a transition state (TS ), by nucleophilic attack of Cys145 S onto αC atom, with a simultaneous transfer of a proton from SH to O. We optimized all these geometries (Fig. 3 ) and for only TS we found one imaginary frequency of −1332 cm−1 indicating a sharp TS, and the animation at this mode corresponds to the change in the reaction coordinate as stated above. No imaginary frequencies are observed for EI and Int . In EI , the S—H and CO bonds are oriented in parallel, with Cys145 S atom pointing at αC and thiol H atom pointing at O(αC) atom. The S—αC and H—O(αC) bond distances are 2.23 and 2.55 Å, respectively. We identified a bond critical point bcp between incoming S and electrophilic αC by applying QTAIM analysis. The electron density (ρ), Laplacian of electron density (∇2ρ), and energy density (H(r)) are 0.0479, 0.2949, and − 0.0068 a.u. The positive sign of ∇2ρ and negative sign of H(r) supports a non-covalent nature of the S-αC interaction. We find a donor-acceptor NBO pair nS → πC=O ∗ with delocalization energy of E2 = 50.4 kcal/mol. NBO charge analysis also shows that the natural charge on sulfur increases from 0.213e in EI1 to 0.254e in TS1, and subsequently, on the O(αC), the negative charge increases from −0.332e (in EI1) to −0.357e in TS1 indicating a charge transfer from Cys145 moiety to the α-ketoamide. In TS , these bond distances are reduced to 1.92 and 1.43 Å, respectively. At TS1, two important NBO pairs are detected (Table S1 and Fig. 4 ): the first one is nO(αC) → σS−H ∗ which indicates the S—H bond cleavage O—H bond formation, and the second one nO(αC) → σNδ−H ∗ supports that His41, by forming H-bonds with O(αC), helps to make the αC more electrophilic and thus more susceptible to the nucleophilic attack (Fig. 3). The geometry at αC changes from planar to tetrahedral as the reaction proceeds from EI to TS to Int . In Int , the values of ρ, ∇2ρ, and H(r) at the bcp on S-αC bond are 0.1631, −0.2526, and −0.1336. The negative sign of ∇2ρ and H(r) indicates the covalent nature of C—S bond. The energy barrier (Eact) for this mechanism is 21.0 kcal/mol computed at ωB97X-D/6-31G(d,p)//B3LYP/6-31G(d):PM6 and the reaction are exoergic, ΔrE ≈ 25.1 kcal/mol. However, in the protein environment at ε = 4, the energy these values are lowered, for example, Eact and ΔrE are ca. 19.1 and − 13.1 kcal/mol, respectively (Fig. 5 and Table 1 ).

Fig. 3

Optimized geometries of the intermediates and transition state of the concerted self-catalyzed reaction (as in Pathway 1) between Mpro and α-ketoamide of the inhibitor from Re face. Atoms marked with asterisks are fixed during geometry optimization. For the sake of clarity, only polar hydrogen atoms are shown.

Fig. 4

The donor-acceptor NBOs of transition states involved in the proposed mechanism according to different pathways: i) TS: a and b, ii) TS: c, d, and e, iii) TS: f, iv) TS: g, v) TS: h, vi) TS: i; vii) TS: j, viii) vii) TS: k.

Fig. 5

Calculated energy profiles of all the pathways for R- and S-thiohemiketal according to a) pathway 1, b) pathway 2, c) pathway 3, d) pathway 4, and e) reaction from Re face according to Scheme 3.

Table 1

Activation parameters and the reaction energies (kcal/mol) of investigated paths explored in the present study in the gas phase and PCM media at ε = 4.a

	Gas phase			ε = 4
	ΔEact(forward)	ΔEact(reverse)	ΔrE	ΔEact(forward)	ΔEact(reverse)	ΔrE
TS1	21.0	46.1	−25.1	19.1	32.4	−13.3
TS2	2.0	29.7	−27.7	1.8	6.0	−4.2
TS3	4.9	13.4	−8.5	−1.0	2.8	−3.8
TS3_ptb	2.4 (10.9)	12.3	−1.4	14.2 (18.0)	20.4	1.4
TS3_ptwb	2.5 (11.0)	12.4	−1.4	2.9(6.7)	5.3	1.4
TSz	12.1	13.2	−1.1	2.5	15.6	−13.1
TS4	1.0	20.2	−19.2	1.0	10.9	−9.9
TS5	4.8	6.0	−1.2	6.8	31.6	−24.8

All ΔEact(forward) and ΔrE values are relative to their corresponding enzyme-inhibitor non-covalent complexes, for example, for TS1, the results are relative to EI1. The ΔEact(reverse) values are relative to their corresponding enzyme-inhibitor covalent complexes, for example, for TS1, ΔEact(reverse) value is calculated relative to Int1.

Data within the bracket are relative to Int3.

Pathway 2: concerted water-assisted reaction without histidine as a function of a shuttle

To model this approach, we propose that the reaction started from the initial inhibitor-enzyme complex EI , which is formed by the interaction between active site amino acids and inhibitor. Herein, the system contains 77 atoms and as the system is small, we have applied DFT method on the whole system. In this case, a water molecule serves as a shuttle to transfer proton of Cys145 SH to α-keto oxygen. As the reaction proceeds, EI transforms to tetrahedral-intermediate Int via the 6-membered transition state TS . The TS is characterized by an imaginary frequency of −567 cm−1, and the imaginary mode clearly identified that the TS corresponds to the synchronous process of proton transfer and nucleophilic attack. Importantly, here, His145 functions as an H-bond donor and not as a proton shuttler. The optimized geometries along with important geometrical parameters of EI , TS and Int are shown in Fig. 6 . In EI , the critical S-αC distance is 2.07 Å and becomes shorter to reach 1.96 Å at TS and 1.88 Å in Int . In EI , we find a bcp between S and αC with ρ and ∇2ρ values of 0.0923 and 0.4328 a.u. The non-negative value of ∇2ρ supports the non-covalent interaction of S and αC. On the other hand, in Int , the negative values of ∇2ρ and H(r) are indications of the covalent nature of C—S bond. In TS , the water molecule is assisting the proton-transfer by bridging between O(αC) and Cys145 S, with two H-bonds of 1.64 and 1.50 Å, respectively.

Fig. 6

Optimized geometries of the intermediates and transition state of the concerted water-assisted reaction (as in Pathway 2) between Mpro and α-ketoamide of the inhibitor from Re face. Atoms marked with asterisks are fixed during geometry optimization. For the sake of clarity, only polar hydrogen atoms are shown.

The NBO analysis shows that, in EI , the main donor-acceptor NBOs nS → παC=O ∗ (E2 = 101.4 kcal/mol) clearly indicates the attack of Cys145 S to αC of ketoamide (Fig. 4). Also, in TS , we properly identified three main donors -acceptor NBO pairs nS → σαC=O ∗, nS → σO−H(W) ∗, and nO(αC) → σO−H(W) ∗ indicating the participation of water in the transition states (See Fig. 4 and Table S1). Furthermore, the NBO charge analysis shows that, in the conversion of EI to TS , the charge transfer on S is not significant, whereas, the negative charge on O(αC) changes from −0.339e to −0.440e indicates that the protonation of O(αC) initiates the nucleophilic attack. However, we note that, in this pathway, the overall resemblance between Int and crystallographic structure (See Fig. 6, Fig. 10, later) is not good. For example, S—C distance in Int is 1.88 but the crystallographic distance is 1.81 Å. Similarly, there a moderate disagreement in the H-bond network between the His41 and tetrahedral-intermediate. The calculated energy of TS is 2.0 kcal/mol in the gas phase and 1.8 kcal/mol at ε = 4 relative to EI (Fig. 5 and Table 1). This step is also highly exoergic of −27.5 and − 4.2 kcal/mol in the gas phase and 4, respectively.

Fig. 10

Hydroxyl group of the R- and S-thiohemiketal of the 5n19 and 6y2f crystal structures, respectively. a) H-bonds between the oxygen atoms of hydroxyl group and amide moiety in warhead with the hydrogen atoms of coronavirus enzyme residue. b) Some selected QTAIM bond critical points indicating by orange sphere. (See Table 2 for details) c) Scatter plots of reduced density gradient (S) vs sign(λ2)ρ. d) Color-filled RDG isosurfaces depicting noncovalent interaction region. The color coding scheme used for isosurfaces is as follows: blue for attractive, red for repulsive, and green for intermediate interactions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Pathway 3: concerted reaction with two water and histidine as a function of shuttle

In this approach, we assumed that the concerted proton transfer and nucleophilic attack is assisted by two water molecules and the nearby neutral Nε protonated His41 residue (see Scheme 1). Herein, the system contains 63 atoms and as the system is small, we have applied DFT method on the whole system. The optimized structure of the enzyme-inhibitor complex EI is shown in Fig. 7 . In EI , two water molecules bridges between Cys145 SH and His41 δNH with three H-bonds of 1.75, 1.31, and 1.47 Å. The reactant proceeds from EI to Int via TS . We detect a bcp at EI , between S and α C with ρ, ∇2ρ, and H(r) values of 0.0540, 0.3303, and − 0.0121 a.u., indicating non-covalent nature of the bond. To support from NBO study, we find a donor-acceptor NBO pair nS → σαC=O ∗ of E2 = 57.1 kcal/mol. The optimized TS is characterized by an imaginary frequency of 208i cm−1. The vibrational movement at imaginary mode indicates the simultaneous transfer of a proton from Cys145 SH to His41 δN via two bridging water molecules W1 and W2 together with the nucleophilic attack by the Cys145 SH onto the αC of the inhibitor. Additionally, the donor-acceptor properties, involved in TS , are well supported by the NBO analysis (Fig. 4). We find that, in TS , four important NBO pairs, namely, nS → σO−H(W1), ∗ nO(W1) → σO−H(W2), ∗ nNδ → σO−H(W2), ∗ and nS → σαC−O ∗ well supports the donor and acceptor atoms. The optimized geometries show that the main distance S-αC in reduces from 2.23 Å at EI to 2.18 Å in TS , 1.97 Å in Int , and 1.83 Å in Int ′, and consequently, the αC-O− increases from 2.22 to 2.24 to 2.28 to 2.41 Å. Furthermore, Int is a tetrahedral oxyanion-intermediate stabilized by three H-bonds: His41 εH···O−, W3 H···O−, and intramolecular α-amidic H···O− with 2.16, 1.76 and 2.31 Å, respectively. Similar to the previous pathways, here also, the non-covalent S—C bond becomes covalent in Int (∇2ρ ≈ −0.2321, H(r) ≈ −0.1240 a.u.). The energy barrier for this step shows 4.9 kcal/mol relative to EI in the gas phase. While the reaction is barrierless (Eact = −1 kcal/mol) at the protein environment ε = 4 (Table 1). Finally, the intermediate Int converts to Int ′ by transferring a proton from His41 δN to oxyanion formed by αC. We consider two paths for such proton transfer: direct and water-mediated, and consequently, we obtained two corresponding transition states TS and TS as shown in Fig. 7. TS and TS are characterized by the imaginary frequencies of 416i and 313i cm−1, respectively. The transition state TS shows that the transferring proton is exactly shared by the imidazole Nδ of His41 and the oxyanion. In TS , we detect an NBO pair nNδ → σO−H ∗ of E2 = 230.1 kcal/mol to support the donor-acceptor properties. In TS , the two important NBO pairs are nO(W3) → σNε−H ∗ and nO(αC) → σO−H(W3) ∗ with delocalization enrgy of 149.6 and 263.1 kcal/mol. Energy profile diagram shows that proton transfer for via TS is more preferable (Eact ≈ 2.9 kcal/mol) than via TS (Eact ≈ 14.2 kcal/mol) in protein environment at ε = 4. Furthermore, the results indicate that, in the gas phase, the ion pair Int is less stable, by ca. 1.4 kcal/mol, than its neutral counterpart Int ′ while the reverse is true in continuum solvent model, that is, at ε = 4, Int is more stable than Int ′ by ca. 1.4 kcal/mol. In Int ′ the S-αC distance is 1.83 Å. The computed structure of Int ′ almost resembles the crystallographic structure within the full enzyme environment.

Fig. 7

Optimized geometries of the intermediates and transition state of the concerted reaction according to Pathway 3 between Mpro and α-ketoamide of the inhibitor from Re face. Atoms marked with asterisks are fixed during geometry optimization. For the sake of clarity, only polar hydrogen atoms are shown.

Pathway 4: stepwise process using histidine as proton donor-acceptor

First-step

Herein, the system contains 74 atoms and as the system is small, we have applied DFT method on the whole system. As stated in the introduction, the first step is the proton transfer from Cys145 SH to Nε protonated His41 Nδ. Before starting the calculation, for correct protonation site and fixing any missing atoms at experimental pH 7.3, we run the crystallographic molecule in H++ server. The results show that, in the crystal, the distance between Cys145 sulfur and His41 Nδ (S⋯Nδ) is 4.94 Å, and this is too large to transfer. The results of H++ server also show that the S—H bond is oriented in such a direction from where the proton transfer is not possible. Besides, as the nearest water molecule, H2O637 (in pdb:6y2e), is 3.49 Å away from the Nε atom of His41, to overcome the problem, we include an explicit water molecule in the proton transfer. Therefore, herein, a water-mediated proton transfer is assumed to occur. Fig. S2 shows that there is a rotation of the imidazole ring to accept the inhibitor. Thus, it can be hypothesized that, in the presence of inhibitor, because of the rotamer of either the imidazole ring of His41 or C—S of Cys145 (as shown in Fig. S4) or by both, the distance between the donor atoms Cys S and His Nδ becomes 3.54 Å. This distance favors the proton transfer to form CysS−/HisH+ ion pair. To optimize and characterize all critical structures properly along with the energy profile and to find the transition state as a first-order saddle point, we used the minimum number of proteins as in the catalytic triad (His41, Cys145, and His164) along with Met165 and two water molecules (W1 and W2). We define the reaction proceeds from the reactant complex (R ) to the ion-pair product complex E via the transition state (TS ). As shown in Fig. 8 , in optimized R , the water molecule W1 connects the donor atoms Cys145 S and His Nδ with two H-bonds of 1.69 and 2.04 Å, respectively, by forming a bridge. The transition state TS is characterized by an imaginary frequency of −117 cm−1. In TS , the bridging W1 has two O—H distances of 1.26 and 1.17 Å. The S—H distance is increased from 1.40 Å in R to 1.68 Å in TS to 3.09 Å in E . The animation of imaginary mode surely indicates the transfer of the proton from donor to acceptor. Also, in the TS, the NBO pairs nS → σO−H ∗ and nNδ → σO−H ∗clearly shows that proton transfer from Cys145 SH to His41 δN (Fig. 4 and Table S1).

Fig. 8

Optimized geometries of the intermediates and transition state of the two-step concerted reaction according to Pathway 4 between Mpro and α-ketoamide of the inhibitor from Si face. Atoms marked with asterisks are fixed during geometry optimization. For the sake of clarity, only polar hydrogen atoms are shown.

The activation barrier of the proton transfer is about 2.5 kcal/mol at ε = 4, while this value is high (12.1 kcal/mol) in the gas phase (Table 1 and Fig. 5). As shown in Table 1, this value is comparable with the recent works [41,54] where the authors showed the inhibition of SARS-CoV-2 Mpro with Michael acceptors. However, our study shows that this step is exotheric for α-ketoamide as an inhibitor of Mpro. The process is an exoergic process of about −13.1 kcal/mol at ε = 4. These results are also in agreement with the previously reported data for another similar nitrile-based inhibitor [ 29].

Second-step

Herein, the system contains 115 atoms and as the system slightly is large, we applied DFT method α-ketoamide, Cys145 and His41 residues. For other residues, we have applied PM6 method. In this step, the Cys145 thiolate ion, formed in by the previous step, attacks the α-keto of the inhibitor to form a C—S covalent bond. In this model, to find the reaction profile, we increase S—C distance gradually, as the reaction coordinate, from 1.73 Å to 3.30 Å with an increment of 0.1 Å at each step (Fig. S5). The result shows that the guess transition structure is at about 2.22–2.23 Å and accordingly, we found the correct transition state TS , which is characterized by an imaginary frequency of 109i cm−1. Unfortunately, from the IRC calculation, we could not obtain the initial enzyme-inhibitor complex EI . From the S—C surface scan, we assumed EI that resides very close to TS in the PES, and thus, we optimized EI with a constant distance of 2.33 Å. The optimized structures of the stationary points are shown in Fig. 8. In EI , the α-keto oxygen is stabilized by a strong hydrogen bond of 1.69 Å from His41 Nε. Also, one crystalline water molecule stabilizes the amidic oxygen to make the αC more electrophilic. Notably, the animation of imaginary mode shows that, at the TS, the attacking of Cys145 and transfer of a proton from His41 Nε to α-keto oxygen occurs synchronously. Additionally, we find two main NBO pairs nS → σαC−O ∗ and nNε → σO−H ∗ to support the donor-acceptor properties (Fig. 4 and Table S1). Finally, TS converts to tetrahedral thiohemiketal Int , where the S-αC bond distance becomes 1.83 Å from 2.22 Å in TS . Consequently the αC-O increases from 1.24 to 1.33 to 1.39 Å as the reaction proceeds from EI  → TS  → Int . Here also, the covalent nature of Int is confirmed from their ∇2ρ, and H(r) values, −0.1663 and − 0.0922 a.u. The results show that all the structural parameters of Int resemble well the crystallographic data (Fig. 8, Fig. 10). The energy barrier of the step is low ca. 1 kcal/mol both at ε = 1 and 4. The conversion EI  → Int is exothergic of −9.9 and − 19.2 kcal/mol in the gas phase and ε = 4, respectively. The exoergicity is similar to the recent work [41,54] for the inhibition with Micheal acceptors with a similar kinetically controlled sulfur‑carbon bond formation reaction.

Reaction mechanism of α-ketoamide with Mpro from Re face to form R-thiohemiketal

Herein, the system contains 95 atoms and as the system is small, we have applied DFT method on the whole system. To model this reaction mechanism, we used the crystal structure of SARS coronavirus main protease in complex with the alpha-ketoamide of R-thiohemiketal (pdb:5n19). We assume that the reaction started from the initial enzyme-inhibitor complex EI (Fig. 9 ). Compare to the Si face reaction with CoV-2 Mpro, here, S…C distance is significantly longer, ca. 3.09 Å. To verify whether S and αC are in fact interacting or not, we detect a bcp between them with ρ, ∇2ρ, and H(r) values are 0.0048, 0.0381, and 0.0024 a.u. Besides, we detect a donor-acceptor NBO pair nS → παC−O ∗ with low delocalization energy of 8.3 kcal/mol. Both the QTAIM and NBO analysis indicates that the interaction between Cys S and αC is very weak and non-covalent in nature. Interestingly, as the Cys145 thiolate approaches the α-keto of the inhibitor from Re face, both the attacking nucleophile and the electrophilic-warhead, experience completely different interacting sites from the Si face reaction. For example, in EI , amide oxygen is stabilized by His41 Nε with a H-bond of 2.36 Å. Similarly, the oxygen of α-keto group is stabilized by the H-bonds from the donor atom of Cys145, Gly143, and one water (W1) with bond distances of 2.90, 2.23, and 2.38 Å, respectively. The reaction proceeds from EI to transition state TS to tetrahedral intermediate Int . The transition state is characterized by one imaginary frequency of 81i cm−1 and the animation of the imaginary mode correctly corresponds to the nucleophilic attack of S− on αC. We find an NBO pair nS → παC−O ∗ to support the direction of electron migration at TS (Fig. 4 and Table S1). With the progress of the reaction, the main distance S-αC decreases from 3.09 Å in EI to attain 2.59 Å in TS , and 1.89 Å in Int , and thereby, the change of hybridization at αC is from sp2 to sp3 to form an oxyanion (which we will discuss later in detail). Consequently, the αC-O− increases from 2.21 (in EI ) to 2.22 (in TS ) to 2.30 Å (in Int ), however, this is not shown in Fig. 9 for clarity. The optimized structure of Int shows that the oxyanion is stabilized by strong H-bonds from Cys145, Gly143, and W1 of 1.66, 2.07, and 1.71 Å, respectively. Notably, as stated above, Int has R-stereochemistry. At Int5, the values of ∇2ρ, and H(r) of bcp are negative which again supports the covalent nature of S—C bond. To verify the validity of the computed results, we protonated and optimized the oxyanion of Int (named as Int -H), and compared it with the crystallographic result. The computed results resemble closely that of the experimental data. For example, the main covalent bond distance S-αC is 1.82 in Int -H is good agreement with the experimental data 1.80 Å [1]. Also, all the H-bonds are in close agreement with the crystallographic data as shown in Fig. 9, Fig. 10 . The energy barrier this reaction is 4.6 kcal/mol in gas phase and 6.8 at ε = 4 (Fig. 5 and Table 1). Like the reaction from Si face, here also, the reaction is exoergic of 24.8 kcal/mol in enzyme environment. However, in the gas phase, exoergicity is lowered.

Fig. 9

Optimized geometries of the intermediates and transition state of the reaction between Mpro and α-ketoamide of the inhibitor from Re face. Atoms marked with asterisks are fixed during the geometry optimization. For the sake of clarity, only polar hydrogen atoms are shown.

Table 2

The values of electron density ρ and its Laplacian ∇2ρ at the important BCPs along with H-bond energy in the interaction between inhibitor-warhead and active enzyme site of the crystal structures of R- and S-thiohemiketal. See Fig. 10 for the labeling of BCP.

BCP	Connected atoms	ρ (a.u.)	∇2ρ (a.u.)	H-bond energy(kcal/mol)
R-thiohemiketal
cp1	α-keto O-Cys145 NH	0.0248	0.0849	−4.79
cp2	α-keto O-Gly143 NH	0.0082	0.0331	−1.08
cp3	α-keto OH-water O	0.0277	0.0963	−5.43
cp4	α-ketoamidic O-His41 NεH	0.0850	0.3334	−18.21
cp5	α-keto C-Cys145S	0.1657	−0.1884
S-thiohemiketal
cp1	α-keto O-His41 NεH	0.0513	0.1492	−10.70
cp2	α-keto O-α-keto amidic N	0.0458	0.1975	−9.47
cp3	α-ketoamidic O-Cys145 NH	0.0195	0.0707	−3.60
cp4	α-ketoamidic O-Ser144 NH	0.0127	0.0495	−2.09
cp5	α-ketoamidic O-Gly143 NH	0.0312	0.1122	−6.21
cp6	α-keto C-Cys145S	0.1613	−0.1718

The noncovalent interaction and stability of oxyanion/hydroxyl group of R- and S-thiohemiketal

To explore the noncovalent interactions including H-bonding to the oxyanion, we extracted the crystal structure and added hydrogen atoms in the missing position and reoptimized with respect these added hydrogen. In order to examine the oxyanion, first, we tried to optimize the geometry by deleting the proton on oxyanion, but, in case of S-adduct, all attempt to obtain oxyanion in presence of His41 NεH are failed, resulting in the transfer of proton from Nε H to O−. However, in this adduct, we notice that the hydroxyl group is stabilized by four nearby polar atoms, namely, Nε of His41, amide H of warhead, and slightly by H of accessible H2O (W1) and another amide H of inhibitor (See Fig. 10). The corresponding distances are 1.72, 1.77, 2.96, and 3.00 Å, respectively. Notably, only two of them form hydrogen-bonds with the oxyanion as evident from Bader's QTAIM plot (Fig. 10b), and the orientation of non-covalent interaction is such that a hydroxyl group is formed. We detected two bcps, namely, cp1 and cp2, correspond to two strong H-bonds with the electron densities of 0.0513 and 0.0458 a.u., respectively. Also, we computed the H-bond energy associated with the bcps. To calculate the H-bond energy, we used the following latest formula derived by Emamian et al. [55] at the very accurate CCSD(T)/jul-cc-pVTZ level including basis set superposition error (BSSE) bond correction as well as SAPT2 + (3)δMP2/aug-cc-pVTZ level, the latter allows the bond energy to be decomposed into physically meaningful components to shed light on the electronic nature of the considered hydrogen bond interactions:

Furthermore, the amidic oxygen is stabilized by NH of triad residues Cys145-Ser144- Gly143, forming three H-bonds of 2.06, 2.77, and 1.88 Å, respectively. These H-bonds are characterized by three bcps, namely, cp3, cp4, and cp5, and their respective ρ and ∇2ρ values are shown in Table 2 . The computed H-bond energy of them are −3.60, −2.09, and − 6.21 kcal/mol, respectively. Thus, such a huge stability, by the H-bonds, of the α-ketoamide warhead reminds us of the advantage of the α-ketoamide compared to the other inhibitor-warhead.

Fig. 10c and d show a scatter plots of RDG S vs sign(λ2)ρ and color-filled RDG isosurfaces depicting noncovalent interaction region. Although we obtained several π-π, CH-π, and H-bond interactions from the calculations, we discussed only the interactions related to α-ketoamide motif.

The thiohemiketal of R stereochemistry experiences completely an opposite interacting residues due to the addition of incoming nucleophile from Re face. For example, the hydroxyl group at thiohemiketal forms three hydrogen bonds, with the triad residues Cys145-Ser144- Gly143 at Mpro, of 1.99, 3.09, and 2.52 Å, respectively, and additionally, forms another H-bond with crystallographic water W1 (Fig. 10a). The QTAIM analysis (Table 2) found three are bcps, namely, cp1, cp2, and cp3 correspond to the H-bonds α-keto O-Cys145 NH, α-keto O-Gly143 NH, and α-keto OH-W1 O but not with the Ser144 residue. For these H-bonds bonds, the H-bond energy computed according to Eq. (3) are −4.79, −1.08, and − 5.43 kcal/mol, respectively.

Besides, unlike S-adduct, the amidic oxygen is stabilized by a very strong H-bond with His41 NεH of 1.40 Å. The QTAIM result shows that the ρ and ∇2ρ values at the bcp (cp4) are 0.0850 and 0.3334 a.u., respectively. Obviously, the H-bond energy is also high: −18.21 kcal/mol.

Finally, the covalent nature of S-αC bond in both the thiohremiketal, R- and S- stereochemistry, is characterized from their bcps (cp5 and cp6) and their ρ and ∇2ρ values. The covalent bond is usually characterized by a negative Laplacian and a negative energy density. Here, for the R-adduct, the energy density and ∇2ρ values are −0.1052 and − 0.1884 a.u., respectively, while for the S-adduct, the energy density and ∇2ρ values are −0.0980 and − 0.1718 a.u., respectively.

Conclusion

In the present paper, we report a theoretical investigation of the reaction mechanism of α-ketoamide based inhibitor with coronavirus main protease Mpro by using DFT and semi-empirical quantum mechanical method within ONIOM model. Our calculation was based on the crystal structure of free coronavirus protease (pdb:6y2e) and the coronavirus-inhibitor complexes of both S and R stereochemistry (respectively, pdb:6y2f and pdb:5n19).

We fulfilled main objectives, namely, to characterize of the noncovalent and covalent interactions between the enzyme active site and inhibitor warhead by using NCI-RDG and Bader's QTAIM, and to uncover the mechanism of thiohemiketal formation between cysteine protease and α-ketoamide inhibitor.

Apart from the common two-step pathway via CysS−/His41+ ion pair formation (pathway 4), we propose three more concerted pathways without forming an ion pair. We have identified and characterized all the stationary points including intermediates and transition states on the PES. Since for the reversible covalent-inhibitor, the change of non-covalent to covalent state is crucial [40], we identified each state properly by applying AIM and NBO analyses.

The computed kinetic and thermodynamic parameres show that although Pathway 1 has a high activation barrier, pathway 2 and 3 have comparable activation barriers with the traditional pathway 4. Also, all these steps are highly exoergic due to the high stability of the oxyanion/hydroxyl group. Nevertheless, the exoergicity of pathway 3 is slightly less (ca. -4.2 kcal/mol) making the ΔEact of the reverse reaction, that is, Int  → EI , low (ca. 2.8 kcal/mol). This indicates the reversible binding nature of the inhibitor. Another notable feature, derived from PCM calculation at ε = 4, is that for all except Re face reaction, the nucleophilic-attacking steps have higher ΔrE in the gas phase. This is probably due to greater stabilization of non-covalent enzyme-inhibitor complex (reactant side) than that of covalent S-thiohemiketal complex (product side).

In conclusion, our computed results suggest that in addition to the traditional mechanism (pathway 4) the reaction can also proceed via the other plausible paths with similar activation barriers.

In the α-ketoamide, the LUMO is π* in nature and resides on α-keto group to accept a lone-pair of electrons from Cys S. Interestingly, here, the LUMO of S—In contains the π* orbital for accepting nucleophile, while in R-In, the lowest energy π* orbital is located on LUMO+1 for such acceptance.

Finally, we have found that in R-thiohemiketal, the α-keto oxygen in the oxyanion/hydroxyl group is stabilized by three hydrogen bonds, while for S-thiohemiketal, the α-keto oxygen is stabilized by two H-bonds but slightly stronger than in R-thiohemiketal. As a result, in both adducts, the sulfur forms S—C covalent bond of similar strengths.

Indeed, the application of an extensive QM/MM approach by using full enzyme residues may slightly alter the position of potential energy surfaces from our results, that is, using static calculation; however, we hope that our model is sufficiently accurate for this work.

In the present situation, we believe that our work is providing extremely relevant and useful information for designing new anti-coronavirus drug using α-ketoamide as an inhibitor.

Declaration of Competing Interest

The authors declare no competing financial interest.