key: cord-0894761-82yhugm0
authors: Taranto, Alex G.; Carvalho, Paulo; Avery, Mitchell A.
title: QM/QM studies for Michael reaction in coronavirus main protease (3CL(Pro))
date: 2008-05-09
journal: J Mol Graph Model
DOI: 10.1016/j.jmgm.2008.05.002
sha: d55e637f236eb9452e2b248b0e48ba584927c908
doc_id: 894761
cord_uid: 82yhugm0

Severe acute respiratory syndrome (SARS) is an illness caused by a novel corona virus wherein the main proteinase called 3CL(Pro) has been established as a target for drug design. The mechanism of action involves nucleophilic attack by Cys145 present in the active site on the carbonyl carbon of the scissile amide bond of the substrate and the intermediate product is stabilized by hydrogen bonds with the residues of the oxyanion hole. Based on the X-ray structure of 3CL(Pro) co-crystallized with a trans-α,β-unsaturated ethyl ester (Michael acceptor), a set of QM/QM and QM/MM calculations were performed, yielding three models with increasingly higher the number of atoms. A previous validation step was performed using classical theoretical calculation and PROCHECK software. The Michael reaction studies show an exothermic process with −4.5 kcal/mol. During the reaction pathway, an intermediate is formed by hydrogen and water molecule migration from a histidine residue and solvent, respectively. In addition, similar with experimental results, the complex between N3 and 3CL(Pro) is 578 kcal/mol more stable than N1-3CL(Pro) using Own N-layer Integrated molecular Orbital molecular Mechanics (ONIOM) approach. We suggest 3CL(Pro) inhibitors need small polar groups to decrease the energy barrier for alkylation reaction. These results can be useful for the development of new compounds against SARS.

Severe acute respiratory syndrome (SARS) is an illness caused by a novel corona virus wherein the main proteinase called 3CL Pro has been established as a target for drug design. The mechanism of action involves nucleophilic attack by Cys145 present in the active site on the carbonyl carbon of the scissile amide bond of the substrate and the intermediate product is stabilized by hydrogen bonds with the residues of the oxyanion hole. Based on the X-ray structure of 3CL Pro co-crystallized with a trans-a,bunsaturated ethyl ester (Michael acceptor), a set of QM/QM and QM/MM calculations were performed, yielding three models with increasingly higher the number of atoms. A previous validation step was performed using classical theoretical calculation and PROCHECK software. The Michael reaction studies show an exothermic process with À4.5 kcal/mol. During the reaction pathway, an intermediate is formed by hydrogen and water molecule migration from a histidine residue and solvent, respectively. In addition, similar with experimental results, the complex between N3 and 3CL Pro is 578 kcal/mol more stable than N1-3CL Pro using Own N-layer Integrated molecular Orbital molecular Mechanics (ONIOM) approach. We suggest 3CL Pro inhibitors need small polar groups to decrease the energy barrier for alkylation reaction. These results can be useful for the development of new compounds against SARS.

ß 2008 Elsevier Inc. All rights reserved.

the active site. Thus, the cysteine thiol moiety possesses a thiolate form at physiological conditions [9] . After the substrate binds to the active site, the thiolate nucleophile attacks the carbonyl carbon of the scissile amide bond, leading to the formation of an oxyanion tetrahedral intermediate [7, 10] . This intermediate is stabilized by strong hydrogen bonds donated by amide groups of the enzyme. The so-called ''oxyanion hole'', in the case of 3CL pro , is formed by the backbone of Gly143 and Cys145 [7, 9] . Several inhibitors have been developed for HRV and HCoV virus [11] [12] [13] [14] [15] . The most successful of them is AG7088 (rupintrivir), which is undergoing testing for its efficiency in the treatment of common cold by HRV, and has advanced to phase II/III clinical trials [16] [17] [18] . This compound is a peptidomimetic with a trans-a,bunsaturated ethyl ester moiety (Michael acceptor), which is able to form an irreversible covalent bond within the active site of 3CL pro proteinase [16] [17] [18] . Similar strategy was used for the development of four wide-spectrum coronavirus inhibitors (N1, N3, N9 and I2) [19] . These inhibitors were designed based on the P4-P1 substrate for recognition of S1, S2 and S4 regions of the active site of M Pro , and they have a trans-a,b-unsaturated ethyl ester moiety as well. High-resolution crystallographic studies identified 21 different amino acids residues located in the active site of M Pro . These inhibitors have been shown to interact directly with the side chains of Gly143, Cys145, His163, His164, Glu166, Gln189, Thr190, and with the main chain of Ala191, Pro168, Arg188, Gln192, Asp187, Met165, Phe140, His41, Met49, Thr26, Thr25, Leu141, Asn142 and His172. These inhibitors, similar to AG7088, have the same kinetic mechanism. Initially, the enzyme forms a reversible complex with the inhibitor, and then undergoes nucleophilic attack by the thiolate cysteine [16, 19] .

The increasing interest in search of a coronavirus inhibitor with higher biological activity has addressed the need of theoretical studies as suitable tool for the design of new compounds [20, 21] . In addition, the knowledge of inhibitors such as analogous fumarate, vinyl sulfones, and vinyl esters, which contain an active double bound, has shown activity for cysteine proteinase and papain-like enzymes [15, 22] . The active site of these enzymes undergoes an addition leading to an alkylated enzyme, as a classical Michael reaction [15, 23] . Initially, the reactivity of the Michael reaction with biological nucleophiles was studied using the semiempirical method AM1 [24] . In this study, the activation enthalpy of the reagents (methylamine, imidazole, both as Michael donors), acrylic acid (as Michael acceptor), and the products (corresponding anions) were characterized. The authors showed that the reaction occurs through an endothermic process, and the rate of the Michael reaction is a function of the rate constant (derived from energy of activation) and the concentration of acceptor and nucleophile [24] . In another study, AM1 was used to study the reactivity of salicylaldehyde acylhydrazone derivatives, acting as Michael acceptors. Conformational analysis, lowest unoccupied molecular orbital (LUMO), and molecular electrostatic potential (MEP) were performed to localize the most favorable atom for nucleophilic attack by the thiolate ion [25] . On the other hand, a set of semiempirical PM3 models using different nucleophilic agents and Michael acceptors were carried out. These models proposed a mechanism through an exothermic process, through tetrahedral intermediate which is more stable than the reagent. A correlation between second-order rate constant and the value of the exponential term in the Arrhenius equation was found, suggesting the use of semiempirical method PM3 for the calculation of the reaction pathways of this system [26] . Another study, Hartree-Fock with 6-31+G(d) and Monte Carlo with a hybrid method AM1/TIP3P calculation were used to study the nucleophilic addition reaction of methanethiolate to N-methylacetamide both gas phase and aqueous solution [27] . Differently from previous PM3 calculation [26] , a stable tetrahedral intermediate were not formed, but these calculations suggested that this structure is a transition state, implying a concerted mechanism [27] . Finally, the reaction for the alkylation of methyl thiolate by electrophilic three-membered heterocycles was carried out by density-functional theory (DFT), using BLYP/TZV+P and COSMO approaches [28] . From this small model of cysteine proteinase, the authors concluded that improved inhibition potency was not resultant from influences on the alkylation step, but from a favorable ionic interaction between carboxylate and imidazolium ion which can stabilize the noncovalent enzyme inhibitor complex. Furthermore, a water molecule, available in the active site, could act as a proton donor to decrease the energy barrier in the ring-opening reactions [28] .

Although several researches have been devoted to study the Michael reaction using models for cysteine proteinase with considerable potential, none of these calculations considered the complete active site, the oxyanion hole, with their respective amino acids (Cys145, His45 and Gly143), and specially the enzyme environment with reasonable theoretical level. Moreover, the cartesian coordinates of atoms in all models were generated randomically, and they did not take into account the enzymatic geometry [20] [21] [22] [23] [24] [25] [26] [27] [28] . Considering these observations, we carried out computational studies including force-field-based methods (molecular mechanics-MM), density-functional theory, and the hybrid approach (QM/QM and QM/MM) denoted by our own n-layered integrated molecular orbital and molecular mechanics (ONIOM). As M Pro consists of two protomers, A and B [19] , in our studies only protomer A will be evaluated. It has 309 amino acids with a total of 4.738 atoms (without water molecules) [19] . Consequently, it has higher computational cost for quantum calculations. Therefore, the enzyme needs to be described by an adequate model, combining accuracy with computational tractability. Combined quantum mechanical and molecular mechanical (QM/MM) methods have gained widespread use in the computational study of large molecular systems, including the enzyme reaction mechanism. The QM/MM approach was described for the first time by Warshel and Levitt [29] . These methods involve modeling a small region of the system that requires quantum mechanical formalism (QM), where occur break and formations of bonds for example, by semiempirical, ab initio or density-functional theory methods and treating the remainder of protein and solvent environment with a low cost molecular mechanics (MM) method. This hybrid approach has been developed because of higher computational cost of conventional QM methods when applied for large systems [30] . Following this guideline, QM/MM methods were applied to study Michael reactions in M Pro enzyme as detailed below. Our attention was focused on the molecular level of the ligand and in the interaction between the ligand and the active site in the pseudoactive site model and protein. This approach can describe the interactions between the ligand and amino acids present in the active site of M Pro . These calculations can be useful for the design of new inhibitors with non-peptidic structure. Peptidyl inhibitors are very susceptible to other proteases and they have low oral bioavailability [11, 15] . Proteinases are potentially important chemotherapeutic targets and have involvement in many others diseases such as viral and parasitic infections, arthritis, cancer, and osteoporosis [15] . One such example is the successful use of proteinase inhibitors in the treatment of acquired immune deficiency syndrome (AIDS) [31] .

The publication of the X-ray structure of SARS M Pro cyteine proteinase in the Research Collaboratory for Structural Bioinformatics Database (access code 1WOF) [19, 32] permitted that our studies were initiated. Three different models were constructed using the atomic coordinate from 1WOF, which they were systematically increased.

In the first model, only the ligand (N1) was fully optimized with tight self-consistent field (scf) routine implemented in Gaussian 03 software [33] . Becke's three-parameter hybrid functional [34] , along with the nonlocal correlation functional of Lee, Yang and Parr (B3LYP-DFT) [35, 36] with 3-21G, 6-31G, 6-31G(d), 6-31G(d,p) and 6-31+G(d,p) basis set [37] and semiempirical Parametric Method 3 with MM correction for peptide bond (PM3MM) [38] methods were used in this study. Afterward, two models with two-layer hybrid were created by ONIOM approach [39, 40] . The first one, B3LYP/basis set and PM3MM were defined as higher and lower layer respectively. Similarly, the second ONIOM model was constructed using PM3MM and the universal force field (UFF) [41] as higher and lower respectively. The higher layer was defined by the double bond and ester group present in trans-a,bunsaturated system, and whole ligand was defined as lower layer (Fig. 1A) . The main goal is observe the effect in Mulliken net atomic charges and geometry in our ONIOM models. In this step, a validation between classical QM and ONIOM was performed. For this purpose, the tight keyword was used to ensure adequate convergence. Of course, the computational cost was increased. The orbital of N1, N3 and AG7088 ( Fig. 1) were carried out using PM3MM and B3LYP/6-31G(d) as well.

Afterward the validation result, another model was designed to study the additions reaction between N1 and the active site of SARS-CoV. The model was constructed through cartesian coordi- nates of N1 and the active site (including oxyanion pocket) obtained from protomer A of 1WOF, which was previously prepared by AMBER8 [42] . Different models were constructed with ONIOM in two and three-layers. However, only a small model with N1, His41, Cys145, Gly143, and three water molecules using B3LYP/6-31G(d,p):PM3MM approach was successfully minimized. In other words, this model considered the trans-a,b-unsaturated moiety, catalytic system, and oxyanion hole in higher layer with 32 atoms, and others parts of ligand, amino acids, and two water molecules in lower layer with 133 atoms. As the amino acids were dissected from protomer A, hydrogen atoms were added to occupy the free valence. The cartesian coordinate of a-carbons and hydrogens added were keep fix during the geometry optimization. Neutral charge was assumed for the whole model, but thiolateimidazolium ion pair was considered in started geometry (Fig. 2 ). The first model of this set was called Reagent, the second of Inter1, which the C 1 -S distance was fixed by 2.720 Å , and the last one was called by Hprod, which consists of an alkylated model. These structures described the nucleophilic attack from thiolate to transa,b-unsaturated moiety of N1. Mulliken net atomic charges, geometrical parameters and relative energies were evaluated in detail for this Michael reaction.

A third model was created, which consists in the whole protein. This model was constructed to demonstrate the influence of longrange electrostatic interactions in Michael adduct reactions, including all atoms of protein. Initially, the geometry of protomer A was obtained from PDB, as describe previously. This structure was prepared following the recommendation by AMBER software, and the three waters molecules were added using the coordinates of PDB file. The S-C 1 bond was broken, and a geometry optimization was carried out by classical force field UFF using Gaussian (default routines). After, ONIOM approaches was carried to build a model with two-layer (PM3MM:UFF and B3LYP/6-31G(d,p):UFF). The ligand, water molecules, His41, Gly143 and Cys145 were definite in higher layer and the others atoms were definite in lower layer, similar to previous model, but Cb and respective hydrogens of His41 were included in QM part. The same ONIOM model was used to study the complexation between N3 and M Pro as well. A validation of stereochemical quality of optimized structures was performed by Ramachandran plot [43] using PROCKECK software [44] . All calculations were performed using the Redwood machine present in Mississippi Center for Supercomputing Research-MCSR.

In general, the lower layer environment directly influences the QM region in hybrid approaches. Consequently, it may be significant in an enzymatic system. As our system has several atoms, it was not possible to avoid a separation between higher and lower layer by a covalent bond. In the case of ONIOM approach, the separation between the regions introduces a link atom. Link atoms have the theoretical disadvantage of introducing extra degrees of freedom, and it is possible that interactions arising as a result of the presence of the link atom might be overcounted. Also, it has to be taken into consideration that hybrid approaches, in general, demand more time than classical calculations, thus it is necessary to limit the number of atoms in the higher layer. In addition, enlarging the higher layer does not necessarily result in a more accurate model [45] . However, it is suggested that QM part should be sufficiently large to incorporate the atoms where chemical and electronic changes occur, including charged QM group, and should not disrupt conjugated system. Considering all observations above, a validation of our system was carried out to check if ONIOM models can give results similar to those obtained by classical QM calculations.

Thus, selected Mulliken net atomic charges (Table 1 ) and geometrical parameters ( Table 2) of trans-a,b-unsaturated system (atoms 1-5, Fig. 1 ) of N1 inhibitor were compared using classical QM and ONIOM models. Table 1 shows Mulliken net atomic charges for selected atoms (Fig. 1 ) calculated by PM3MM, B3LYP/ basis set, B3LYP/basis set:PM3MM, B3LYP/6-31G(d,p):UFF and PM3MM:UFF. As it can be seen in all cases, the Mulliken net atomic charges changed according to the theoretical level, as expected. Mulliken net atomic population analysis is an arbitrary scheme for assigning charges. Indeed, all such schemes are ultimately arbitrary. Atomic charges, unlike the electron density, are not a quantum mechanical observable, and are not unambiguously predictable from first principle [46] . However, the most positively charged atom is C 3 , except for B3LYP/6-31+G(d,p) that is C 2 , and the most negative atom is O 4 , except for B3LYP/3-21G, B3LYP/6-31G and B3LYP/6-31+G(d,p):PM3MM that is O 5 . In general, a good correlation coefficient (r 2 ) was obtained between classical QM and ONIOM models, ranging from 0.97 to 0.99, except between B3LYP/ 6-31+G(d,p) and correspondent ONIOM model (B3LYP/6-31+G(d,p):PM3MM (r 2 = 0.517)). This difference can be explained. 6-31+G(d,p) basis set has diffuse functions that permit orbitals to occupy a large region of space only for heavy atoms. It is possible that orbitals between heavy atoms of quantum part and link atoms, which are hydrogens, are not interacting, maybe it is necessary a double plus of this basis set, 6-31++(d,p), which adds diffuse functions for hydrogen atoms as well. Although 6-31++(d,p) basis set calculations are necessary to answer this question, the size of our system, with 43 heavy atoms, makes this practically unviable.

The same methodology to validate Mulliken net atomic charges was carried out to study geometrical parameters in ONIOM models. Table 2 shows the trans-a,b-unsaturated geometry of N1. Atomic distances, angles, and dihedral angles of classical QM calculations were compared with ONIOM and X-ray geometry. As shown in this table, a high correlation between classical QM and ONIOM geometry (r 2 = 1.00) was obtained. Only two atomic Fig. 2 . Model of N1 complexed with oxyanion and catalytic system. The colors red, black, and * symbol were used to assign atoms with higher layer, lower layer, and frozen respectively. distances are different from the others, C 1 C 2 and C 2 -C 3 , in X-ray structure. This is a consequence of sp 3 hybridization of C 1 and C 2 atoms that they have in protein structure, which is alkylated by Cys145 [19] . Although these differences exist, a good correlation among X-ray, classical QM and ONIOM geometry was obtained, with r 2 = 0.99 in all cases.

Molecular orbital (MO) HOMO (highest occupied molecular orbital) and LUMO were carried out for N1, N3, and AG7088 by PM3MM (N1 only) and B3LYP/6-31G(d). MO are expected to be responsible for the specificity of the interactions between protein and ligand. As it can be seen in Fig. 3 , HOMO is localized mainly on the g-lactam ring, whereas LUMO is localized mainly on the isoxazole ring. Semiempirical and DFT methods have practically the same localization of molecular orbitals, even though DFT MO are more extended, as a result of inclusion of d orbitals from basis set. Only a small region of HOMO is localized in the trans-a,bunsaturated system. These features suggest the most important part of these inhibitors, and indicating that heteroatoms of lactam and isoxazole ring have a nucleophilic and electrophilic character respectively, responsible for electronic interaction with receptor. In other words, lactam and isoxazole rings are important for the reversible complex formed with the enzyme, before the nucleophilic attack by the thiolate cysteine [16, 19] . In the case of inhibitors of 3CL pro , previous SAR studies described that an amine in P1 and small hydrophobic groups in P2 positions are required to improve antiviral activity, which are the same positions of the MO found in our studies [16, [47] [48] [49] .

In this first validation, semiempirical PM3MM and ONIOM models demonstrated ability to reproduce results obtained by DFT, principally with B3LYP/6-31G(d) and B3LYP/6-31G(d,p). Furthermore, ONIOM models using MM as lower layer gave reasonable results. These ONIOM models can be used to study the mechanism of Michael reaction considering the whole enzyme.

As demonstrated earlier, validation of models is very important for the next steps, which Michael reaction will be performed. It can have a dramatic influence on the results that will be obtained from this. As can be seen, small region in higher layer can reproduce similar charge and geometrical when compared with classical QM 

As described previously, three models were constructed to study this Michael reaction. However, one of them called Oxyprod will be discussed differently from the others. Select geometrical parameters, Mulliken charges, and Reaction coordinates for all models are provided in Tables 3 and 4 , and Fig. 4 respectively. For a better analysis of these results, four groups were created using nine atomic distances and one dihedral angle. In the first group, the nucleophilic attack was described by C 1 -S, C 1 -C 2 , and Ht-C 2 distances. In the second group, the effect of catalytic system will be discussed analyzing the distances among Ht, Nt, and S atoms. The effect of the oxyanion pocket will be discussed comparing the distances among O 4 , H 8 and H 9 , as third group. Finally, conformational changes will be discussed by H 7 (w1)-C 2 , O(w2)-Hp distances, and C-Ca-Cb-Cg(His41) dihedral angle.

As it can be seen by the first group of geometrical parameters (Table 3) , the distance between C 1 and S(Cys145) decreased from 3.8 to 1.87 Å , while the C 1 -C 2 distance (typical double bond) increased from 1.34 to 1.54 Å in Reagent and Hprod models. These results suggest that, while S(Cys145) attack C 1 , a negative charge (À0.23) is formed on C 2 ( Table 4 , Inter1). This negatively charged atom can abstract a proton (Ht) from Nt(His41), as seen in Table 3 by the Ht-C 2 distance that changed from 1.10 Å (Reagent) to 3.95 Å (Hprod). It shows that the hybridization of C 1 and C 2 changed from planar sp 2 to tetrahedral sp 3 as well. The analysis of catalytic system (Ht-Nt and Ht-S distances) showed that Ht is located between S(Cys145) and Nt(His41), 2.11 and 1.36 Å from Nt and S atom respectively. S and Nt atoms have charge density very close between each other (À0.14 and À0.12, respectively (Table 4) ). It is implying that the energy barrier for catalysis process should be very low, and ionic and neutral form of His41 and Cys145 are very close in energy in studies performed in vacuum system. After nucleophilic attack, the distance between Ht and Nt increased from 2.11 to 2.95 Å , while the negative charge on Nt decreased from À0.12 to À0.11, as seen in Tables 3 and 4 respectively. In addition, a negative charge is formed on S and C 2 atoms (À0.07 and À0.17 respectively). The oxyanion hole has been described as formed by backbone hydrogens of Cys145 and Gly143 [19] . In our model, it is formed by H 8 and H 9 . However, Table 3 shows that only O 4 -H 8 distance is a standard hydrogen bond (2.37-2.62 Å ).

The minimum energy reaction path connecting Reagent and Hprod is depicted in the diagram coordinate shown in Fig. 4 . Our calculations show that the product obtained has a relative energy of À4.85 kcal/mol. Although frequency calculation are not available in ONIOM approach, our calculations suggest that a transition state structure was obtained (Inter1) with a distance of 2.72 Å between S and C 1 atom, and it has 19.31 kcal/more less stable than Reagent model.

A search for a transition state found an unexpected saddle point structure, called Oxyprod (Fig. 5) , which has C 1 -S distance less than Hprod (1.81 Å , Table 3 ). Weak hydrogen bonds are formed between Ht(His41) and O 5 (2.89 Å ). In addition, a higher negative charge is formed in C 2 (À0.69), which is stabilized by Ht(His41) and another proton from water molecule (w1) with distances of 2.89 and 2.41 Å respectively. This model has more 50.05 kcal/mol than Reagent model. Similar structure had described for papain inhibitors [28] . In this case, a favorable ionic interaction between imidazolium and ligand was observed, even though the cartesian coordinate of atoms was not fixed and the geometry of enzyme was not considered. As a result, the atoms could find favorable position in the PES.

In this study, a Michael reaction has investigated using QM/QM models of ligand (N1), oxyanion pocket (Cys145 and Gly143), and catalytic system (Cys145 and His41). Theses results indicate an exothermic mechanism, as described previously [26, 28] . In addition, a transition state can be described for the whole protein addressing for a concerned mechanism, because S and Ht atoms are approximating of trans-a,b-unsaturated in the same time. Furthermore, an interaction between hydrogen from water molecule (H 7 w1) and C 2 was obtained (2.29 Å ) due the negative charge formed after nucleophilic attack. It is probably that this water molecule decreases the activation energy in transition state structure as well, very similar with Oxyprod (Fig. 5) . Likewise with our results, QM/MM studies obtained by Byun and Gao showed that a tetrahedral intermediate does not exist for cysteine protease [27] . Another point, it is very probable that small polar groups of Michael acceptor can decrease the energy barrier for this reaction. It can be observed for inhibitors of papain [22] . Consequently, they should be more active than apolar groups. In addition, previous studies have revealed that all proteases have a His-Cys catalytic dyad [6, 9, 19] . However, our calculations have shown an interaction between a water molecule (w2) and Hp(His41). Perhaps, the w2 has important role as a third member, similar with others proteases [15] .

In this part of our studies, the whole protein was considered, which active site, oxyanion hole, and the trans-a,b-unsaturated moiety are in QM part. Initially, the stereochemical quality of PM3:UFF was compared with X-ray structure of 1WOF, AMBER force field and UFF (Fig. 6a-d) . As can be seen by Ramanchandran plots, the quality of models depend of methodology used. In general, Asn88 was found in disallowed region for all models. Amber force field showed more than two amino acids in this region, Thr28 and Leu286 (Fig. 6b ). In addition, the hybrid model showed Asn88 and Val237 in disallowed region, and Ser117, Val118, Asn123, Ser127, Ile217, Pro256, Asn281 in generously allowed region (Fig. 6c) . Although the PROCHECK showed a set of amino acids which have geometrical distortion in Phi and Psi angles, none of them are present in active site. In other words, the hybrid model can be used to study the interaction between the ligand and M Pro [16, 19] .

After this validation, which geometrical aspects were compared among different force field and the ONIOM model, the geometry of optimized active site, oxyanion hole, and the trans-a,b-unsaturated moiety were demonstrated by Fig. 7 . As can be seen, the water molecules (w1 and w2) are very close of trans-a,bunsaturated and His41 respectively, for both N1 and N3. The distance between w1 and Ca is 3.46 and 4.01 Å ; and the distance between w2 and His41 is 2.35 and 2.43 Å for N1 and N3 respectively, characterizing a typical hydrogen bond between w2 and His41. As described before, it is possible that these water molecules can contribute for decrease the active energy in the alkylation process. Before the nucleophilic attack by Cys145, a reversible complex is formed between inhibitor and protease. This complex is characterized by the distance among ligand and enzyme. First, the sulphur atom of Cys145 distance 3.40 and 3.56 Å from Cb for N1 and N3 respectively. Second, the oxyanion hole is well characterized by distance among Gly143, Cys145 and ligand. The hydrogen of Gly143 and Cys145 distance 3.44 and 3.03 Å from oxygen of ligand for N1 and N3 respectively, characterizing a hydrogen bond. In the case of Cys145, this distance increases to 4.27 and 4.49 Å . Furthermore, a hydrogen bond was formatted between His41 and the oxygen of ligand, which distance 2.55 and 3.07 Å for N1 and N3. This result shows that His41 can contribute for complex stage stabilizing the oxyanion hole, and after to nucleophilic attack of sulphur atom of Cys145, which Ht can approximate of Ca.

The binding energy between M Pro and two inhibitors (N1 and N3) [9] is described in Table 5 . As can be seen, both inhibitors can complex with enzyme with favorable binding energy. N1 and N3 can stabilize the M Pro by À340.997 and À918.992 kcal/mol respectively. In other words, N3 can stabilize the enzyme 577,995 kcal/mol more than N1. This difference of binding energy between N1 and N3 can be attributed by intermolecular interactions. In order to observe this with more details, the interactions between ligand and protein, the number of hydro- phobic interactions, hydrogen bond (HB), and bad contact among atoms in active site were carried out by Rasmol and Swiss-PDBViewer softwares [50, 51] . Rasmol showed that ethyl group of N1 distance 6.25 and 5.08 Å , from Met49 and Leu27, respectively, while aromatic ring of N3 distance only 3.40 and 4.54 Å from the same amino acids. Following, Swiss-PDBViewer showed that N1 and N3 have 27 and 29 HB, and 7 and 5 atoms with undesirable contact between ligand and protein, respectively. According with observations, it seems that N3 fits better than N1 in active site, decreasing the binding energy. This result shows the first step of reaction, which the inhibitor initially forms a reversible complex with the protease. The evaluation of the ONIOM model addresses for N3 as more active than N1, similar to experimental results [9] , which enzyme inhibition of 10.7 and 9.0 mM were obtained by N1 and N3 respectively.

In this study, the models constructed were validated with classical quantum-mechanic calculation (DFT) and X-ray structure, which assured trustful models. As the whole system has 4.738 atoms, two strategies were adopted to decrease the complexity of the system. The first one, the number of atoms was decreased resulting in a small model which only the ligand, active site and oxyanion hole were studied by QM/QM, using B3LYP/6-31G(d,p):PM3MM model. On the other hand, the second strategy, the whole system was studied by PM3MM/UFF. The results of both models indicated the main conformational changes during complexation and alkylation reaction between M Pro and N1 (and N3, in the case of QM/MM studies). These approaches allowed us to study the Michael reaction by QM/QM and QM/MM, which showed an exothermic and concerned mechanism. Furthermore, the results suggested that small polar groups and the water molecule (w1) can address the reaction decreasing the energy barrier for the second step of reaction, which it would be a good starting point for new leader compounds. In addition, the water molecule (w2) could have an important participation as a third member of the catalytic system. Another point, although only the binding energy of two compounds was performed, the QM/MM model was able to estimate the theoretical activity of them. Although further calculations are necessary to understand the whole process, these findings can provide information for development of new inhibitors against a future emerging SARS disease, and showed new perspectives about the proteolytic mechanism of Cys proteases.

A major outbreak of severe acute respiratory syndrome in Hong Kong

Identification of a novel coronavirus in patients with severe acute respiratory syndrome

Enteric involvement of severe acute respiratory syndrome-associated coronavirus infection

The piconarviral 3C proteinases: cysteine nucleophiles in serine proteinase folds

Coronavirus genome: prediction of putative functional domains in the non-structural polyprotein by comparative amino acid sequence analysis

Coronavirus main proteinase (3CL PRO )

Inhibitors of 3C cysteine proteinases from picornaviridae

Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra a-helical domain

The Crystal structures of severe acute respiratory syndrome virus main protease and its complex with an inhibitor

Exploring the binding mechanism of the main proteinase in SARS-associated coronavirus and its implication to anti-SARS drug design, Bioorgan

Thiol-dependent enzymes and their inhibitors: a review

High-throughput screening identifies inhibitors of the SARS coronavirus main proteinase

Synthesis and evaluation of ketoglutamine analogues as potent inhibitors of severe acute respiratory syndrome 3CL pro

Design and synthesis of peptidomimetic severe acute respiratory chymotrypsin-like protease inhibitors

Inhibitors of cysteine proteases

Structure-assisted design of mechanism-based irreversible inhibitors of human rhinovirus 3C protease with potent antiviral activity against multiple rhinovirus serotypes

Structure-based design, synthesis, and biological evaluation of irreversible human rhinovirus 3C protease inhibitors. 8: Pharmacological optimization of orally bioavailable 2-pyridone-containing peptidomimetics

Conservation of amino acids in human rhinovirus 3C protease correlates with broad-spectrum antiviral activity of rupintrivir, a novel, human rhinovirus 3C protease inhibitor

Design of wide-spectrum inhibitors targeting coronavirus main proteases

Human rhinovirus 3C protease: generation of pharmacophore models for peptidic and nonpeptidic inhibitors and their application in virtual screening

A new lead nonpeptidic active-site-directed inhibitors of the severe acute respiratory syndrome coronavirus main protease discovered by a combination of screening and docking methods

Structure-activity relationships for inhibition of papain by peptide Michael acceptors

March's Advanced Organic Chemistry: Reaction, Mechanism, and Structure

Modeling the reactivity of acrylic acid and acrylate anion with biological nucleophiles

A possible molecular mechanism for the inhibition of cysteine proteases by salicylaldehyde N-acylhydrazones and related compounds

An approximation to the mechanism of inhibition of cystein proteases: nucleophilic sulphur addition to Michael acceptors type compounds

A combined QM/MM study of the nucleophilic addition reaction of methanethiolate and N-methylacetamide

Theoretical studies about the influence of different ring substituents on the nucleophilic ring opening of three-membered heterocycles and possible implications for the mechanism of cysteine protease inhibitors

Theoretical studies of enzymic reactions: dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme

Combined quantum mechanical/molecular mechanical methodologies applied to biomolecular systems

Goodman & Gilman's The Pharmacological Basis of Therapeutics

Gaussian 03, Revision A.7

Density-functional thermochemistry. 3: The role of exact exchange

Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density

Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis

Special issue: MOPAC-A semiempirical molecular-orbital program

A new ONIOM implementation in Gaussian98. Part I: The calculation of energies, gradients, vibrational frequencies and electric field derivatives

Combining quantum mechanics methods with molecular mechanics methods in ONIOM

UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulation

Conformation of polypeptides and proteins

Manual Procheck v.3.0: Programs to Check the Stereochemical Quality of Protein Structures

Quantum Medicinal Chemistry

Exploring Chemistry with Electronic Structure Methods

Structure-based design, synthesis, and biological evaluation of irreversible human rhinovirus 3C protease inhibitors. 2: Peptide structure-activity studies

Synthesis and evaluation of peptidyl Michael acceptors that inactivate human rinovirus 3C protease and inhibit virus replication

Tripeptide aldehyde inhibitors of human rhinovirus 3C protease: design, synthesis, biological evaluation, and cocrystal structure solution of P1 glutamine isosteric replacements

RasMol biomolecular graphics for all

Deep View: The Swiss-PDBViewer User Guide Version 3.7

The author would like acknowledgments The University of Mississippi, in especially for Mississippi Center for Supercomputing Research-MCSR, Dr. Prashant Desai, and Universidade Estadual de Feira de Santanta-UEFS.