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 (Cys 115 ) 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 M pro .
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-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-J o u r n a l P r e -p r o o f Journal Pre-proof 
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 Figure 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 The QM/SQM energies were calculated using the following equation as the separate energies of a QM calculation on the QM region (E model , QM) and an SQM calculation on the whole system (E real , SQM) minus an SQM calculation on the QM region (E model ,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:
J o u r n a l P r e -p r o o f Journal Pre-proof 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-J o u r n a l P r e -p r o o f
Lewis NBO (j). For each electron-donor NBO (i) and electron-acceptor NBO (j), the stabilization energy E 2 ij associated with delocalization ij is calculated as (2) where is the donor orbital occupancy, and are orbital energies, and F(i,j) is the offdiagonal 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 M pro 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.
As stated above, we started our investigation with the crystal structure of ligand-free M pro (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 His 41 , Gly 143 , Ser 144 , Cys 145 , His 164, and His 165 . Our objective is to understand the characteristics and relative strength of hydrogen bonds and other noncovalent interactions occurring in the active part of the enzyme in the absence of the warhead of inhibitor. Figure 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 Figure S1b 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 Figure S2 . Aiming to have a complete overview of the ketoamide inhibitor, we optimized a simple -ketoamide N-methyl-2-oxopropanamide. We find that J o u r n a l P r e -p r o o f Journal Pre-proof there are three possible isomers of -ketoamide C 1 , C 2 , and C 3 ( Figure 1b) . Obviously, the C 3 isomer is more stable than the other two isomers C 1 and C 2 , 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 Figure 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 Figure S3 ).
Herein, the His41 protonated in ε-position instead of δ was used. The protonation of His 41 in εposition prevents this residue from being directly involved in the proton transfer from Cys 145 , the role that was demonstrated in studies on proteolysis in M pro . 37
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 1 is formed.
Then, this complex proceeds to form a tetrahedral complex (Int 1 ) via a transition state (TS 1 ), by nucleophilic attack of Cys 145 S onto C atom, with a simultaneous transfer of a proton from SH to O.
We optimized all these geometries ( Figure 3 ) and for only TS 1, we found one imaginary frequency of However, in the protein environment at  = 4, the energy these values are lowered, for example, E act and  r E are ca. 19.1 and 13.1 kcal/mol, respectively ( Figure 5 and Table 1 ). 
To model this approach, we propose that the reaction started from the initial inhibitor-enzyme complex EI 2 , 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 J o u r n a l P r e -p r o o f Journal Pre-proof system. In this case, a water molecule serves as a shuttle to transfer proton of Cys 145 SH to -keto oxygen. As the reaction proceeds, EI 2 transforms to tetrahedral-intermediate Int 2 via the 6-membered transition state TS 2 . The TS 2 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, His 145 functions as an H-bond donor and not as a proton
shuttler. The optimized geometries along with important geometrical parameters of EI 2 , TS 2, and Int 2 are shown in Figure 6 . In EI 2 , the critical S-C distance is 2.07 Å and becomes shorter to reach 1.96 Å at indicating the participation of water in the transition states (See Figure 4 and Table S1 ).
Furthermore, the NBO charge analysis shows that, in the conversion of EI 2 to TS 2 , 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 2 and crystallographic structure (See Figure 6 and J o u r n a l P r e -p r o o f Journal Pre-proof relative to EI 2 ( Figure 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.
Optimized geometries of the intermediates and transition state of the concerted water-assisted reaction (as in Pathway 2) between M pro 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.
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 His 41 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 3 is shown in Figure 7 . In EI 3 formed by C. We consider two paths for such proton transfer: direct and water-mediated, and consequently, we obtained two corresponding transition states TS 3_pt and TS 3_ptw as shown in Figure 7 .
TS 3_pt and TS 3_ptw are characterized by the imaginary frequencies of 416i and 313i cm -1 , respectively.
The transition state TS 3_pt shows that the transferring proton is exactly shared by the imidazole N of 
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 
Journal Pre-proof an explicit water molecule in the proton transfer. Therefore, herein, a water-mediated proton transfer is assumed to occur. Figure 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 His 41 or CS of Cys 145 (as shown in Figure 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 (His 41 , Cys 145, and His 164 ) along with Met 165 and two water molecules (W 1 and W 2 ). We define the reaction proceeds from the reactant complex (R z ) to the ion-pair product complex E z via the transition state (TS z ). As shown in Figure 8 (Figure 4 and Table S1 )
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 Figure 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 M pro with Michael acceptors. However, our study shows that this step is exotheric for -ketoamide as an inhibitor of M pro . 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, Cys 145 and His 41 residues. For other residues, we have applied PM6 method. In
Journal Pre-proof this step, the Cys 145 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 ( Figure 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 4 , 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 4 . From the SC surface scan, we assumed EI 4 that resides very close to TS 4 in the PES, and thus, we optimized EI 4 with a constant distance of 2.33 Å. The optimized structures of the stationary points are shown in Figure 8 . In EI 4 , the -keto oxygen is stabilized by a strong hydrogen bond of 1.69 Å from His 41 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 Cys 145 and transfer of a proton from His 41 N to -keto oxygen occurs synchronously. Additionally, we find two main NBO pairs  and  to support the donoracceptor properties ( Figure 4 and (Figures 8 and 10 ). The energy barrier of the step is low 
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 5 (Figure 9 ). Compare to the Si face reaction with CoV-2 M pro , here, S…C distance is significantly longer, ca. Table S1 ). With the progress of the reaction, the main distance SC decreases from 3.09 Å in EI 5 to attain 2.59 Å in TS 5 , and 1.89 Å in Int 5 , and thereby, the change of hybridization at C is from sp 2 to sp 3 to form an oxyanion (which we will discuss later in detail).
Consequently, the CO  increases from 2.21 (in EI 5 ) to 2.22 (in TS 5 ) to 2.30 Å (in Int 5 ), however, this is not shown in Figure 9 for clarity. The optimized structure of Int 5 shows that the oxyanion is stabilized Figure 9 and 10. The energy barrier this reaction is 4.6 kcal/mol in gas phase and 6.8 at  = 4 ( Figure 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. 
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 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.
J o u r n a l P r e -p r o o f Journal Pre-proof Table 2 The values of electron density  and its Laplacian  2  at the important BCPs along with Hbond energy in the interaction between inhibitor-warhead and active enzyme site of the crystal structures of R-and S-thiohemiketal. See Figure 10 for the labelling of BCP. 
In the present paper, we report a theoretical investigation of the reaction mechanism of -ketoamide based inhibitor with coronavirus main protease M pro 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 E act of the reverse reaction, that is, Int 3 EI 3 , 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  r E 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.
J o u r n a l P r e -p r o o f Journal Pre-proof 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.
Figures S1-S5 and Tables S1 along the Cartesian coordinates of all the relevant species are presented in Supporting Information.
Corresponding Author * S. Banerjee. Email: sbanchem@gmail.com 