A QM/MM Evaluation of the Missing Step in the Reduction Mechanism of HMG-CoA by Human HMG-CoA Reductase

Statins are important drugs in the regulation of cholesterol levels in the human body that have as a primary target the enzyme β-hydroxy-β-methylglutaryl-CoA reductase (HMGR). This enzyme plays a crucial role in the mevalonate pathway, catalyzing the four-electron reduction of HMG-CoA to mevalonate. A second reduction step of this reaction mechanism has been the subject of much speculation in the literature, with different conflicting theories persisting to the present day. In this study, the different mechanistic hypotheses were evaluated with atomic-level detail through a combination of molecular dynamics simulations (MD) and quantum mechanics/molecular mechanics (QM/MM) calculations. The obtained Gibbs free activation and Gibbs free reaction energy (15 kcal mol−1 and −40 kcal mol−1) show that this hydride step takes place with the involvement of a cationic His405 and Lys639, and a neutral Glu98, while Asp715 remains in an anionic state. The results provide an atomic-level portrait of this step, clearly demonstrating the nature and protonation state of the amino acid residues involved, the energetics associated, and the structure and charge of the key participating atoms in the several intermediate states, finally elucidating this missing step.


Introduction
β-hydroxy-β-methylglutaryl-CoA reductase (HMGR, EC 1.1.1.34) has been studied for decades through X-ray crystallography [1][2][3]. HMGR plays an important role in the biosynthesis of various lipids such as cholesterol and compounds such as isoprenoids and steroids [4,5]. The reduction of HMG-CoA to mevalonate is the rate-limiting step, and this reaction is one of the rare examples of an oxidoreductase that uses two molecules of the cofactor, in this case, NADP(H) (Nicotinamide adenine dinucleotide phosphate) [6][7][8][9][10]. The HMGR is a tetramer or "dimer of dimers" with an active site on each chain [11].
There are two classes of HMGRs, Class I (eukaryotes) and Class II (prokaryotes). The Class I HMGRs are mainly membrane-bound and consist of transmembrane and a catalytic domain carboxyl-terminal and an amino-terminal membrane anchor domain, whereas the Class II HMGR enzymes are all cytosolic and lack the membrane anchor domain [12,13]. Due to dissimilarities between the two classes, the existence of significant differences in the reaction mechanism is possible. The reaction mechanism of this enzyme has been studied both by computational models and experimentally.
Despite this, many questions have remained open to dispute. Earlier studies have shown that the reaction mechanism involves two hydride transfers, a cofactor exchange and a hemithioacetal decomposition step; thus, the hypothesis of a hemithioacetal intermediate was developed but never proven due to its non-reactivity toward NADP(H) [6][7][8][9][10]12].
However, the aldehyde intermediate, mevaldehyde, has never been detected experimentally [6][7][8][9][10]12]. Qureshi et al. [9] described a hydride transfer mechanism and hemithioacetal breakdown, whereas Veloso et al. [14] used pH-dependent kinetic analysis and showed that the protonation state of the histidine affects the redox form of the NADP(H). Organisms that can survive on mevalonate as a source of carbon are known as Pseudomonas mevalonii (Pm), and it has been discovered that large quantities of PmHMGR are involved in mevalonate catabolism, which is a reversed mevalonate pathway [15]. Later, it was observed that Glu83 and His381 (analogous to Glu98 and His405) were essential for PmHMGR activity, as well as Asp766 (analogous to Asp715) in the mammalian enzyme [16]. Hence, Frimpong et al. [17] proposed a new reduction mechanism with the aspartate as the general acid. The PmHMGR structure confirmed that Asp283, Lys267, and Asn271 (analogous to Asp715, Lys639, and Asn643) were conserved in the active site and that the Lys267/Lys639 is located close to Glu83/Glu98 and Asp283/Asp715; thus, Tabernero et al. offered a new mechanism where Lys267/Lys639 protonated the carbonyl of the substrate [18].
Istvan et al. [19] studied the human HMGR mechanism and concluded that the acidic residue Glu559 is protonated, while Asp767 (analogous to Glu98 and Asp715) is negatively charged. Moreover, Glu559/Glu98 and Lys735 (analogous to Lys639) form an oxyanion hole that stabilizes the intermediates after the hydride transfer [3,20]. Asp767/Asp715 supports the position of Lys691 (analogous to Lys639), which controls the protonation state of the glutamate, whereas the possible role of His866 (analogous to His405) was previously described by Frimpong et al. in a study of Syrian-Hamster and other studies [3,17,19].
Haines et al. [3] conducted a computational study focused on the protonation state of Glu83/Glu98 in the PmHMGR mechanism, which suggested that the first hydride transfer was the rate-limiting step and that the acid/base for this reaction was Glu83/Glu98, which would be neutral at the beginning of the reaction. Another Haines et al. [12] study proposed a revised mechanism where the thioester of HMG-CoA is reduced and the oxyanion is stabilized through Glu83/Glu98 and Lys267/Lys639. The rate-limiting step has not been determined experimentally, but it has been determined computationally that the energy barrier for the first hydride transfer is 21.8 kcal/mol, which corresponds to a rate constant of 1 s −1 /min [3,12]. A turnover rate of 1 s −1 to 1 min −1 (18.2-20.7 kcal/mol) is mentioned in the literature [20].
Finally, Oliveira et al. [20] studied the first hydride transfer in hHMGR and showed that the most favorable reaction mechanism is when Glu559/Glu98 is protonated and Lys691/Lys639 acts as a proton donor with the Gibbs free activation energy of 17.8 kcal/mol. Oliveira et al. [20] explored the possibility of mevaldehyde formation and mevaldyl-CoA intermediate, which led to the conclusion that the protonation of the hemithioacetal by Lys691/Lys639 was the more stable and favorable intermediate. In the reaction mechanism, several residues in the active site are ionizable, and therefore, various studies have proposed different protonation states for Glu98, His405, Lys639, and Asp715 [2,3,12,14,16,17,19,20]. A recent review on HMG-CoA reductase and its mechanism was done by Gesto et al. [21]. Moreover, a detailed reaction mechanism based on the studies of Oliveira et al. [20] and Haines et al. [12] is shown in Figure 1.
It was also noted that the cofactor exchange happens before aldehyde formation, which is consistent with the studies of Qureshi et al. [9] and Dokainish et al. [22]. Additionally, flap domain movements are involved in the binding of a substrate and cofactor, suggesting that flap domain motions have some kind of an impact on the formation of the mevaldehyde [12].
To design advanced medications with fewer side effects and improve current cholesterollowering drugs such as statins, which are a key molecule in the competitive inhibition of hHMGR in a liver [23], it is essential to understand the reduction mechanism of the hHMG-CoA. With that in mind, we studied the second reduction step of the hHMGR catalytic mechanism to finally understand the second hydride transfer reaction at the atomic level. of the hHMG-CoA. With that in mind, we studied the second reduction step of the hHMGR catalytic mechanism to finally understand the second hydride transfer reaction at the atomic level. Figure 1. The reaction mechanism catalyzed by the HMG-CoA reductase based on the previous studies of Oliveira et al. [20] and Haines et al. [12], as described in detail in [21]. His866/His405 protonates -SCoA, whereas Lys691\Lys639 (-) or Glu559/Glu98 (-) is a proton donor after the hydride transfer. Asp767 corresponds to Asp715, and Glu559 corresponds to Glu98 in the current study.

Model Preparation
All systems were modeled from a human wild-type structure with Protein Data Bank (PDB) ID: 1DQA (resolution 2.0 Å), which represents the human enzyme in two dimeric structures (dimers AB and CD) with bound NADP+ cofactor, CoA, and a mevalonate analogue. Dimer AB was chosen to model the system, as chains A and B have the smallest distances between the mevalonate analogue and the cofactor, and thus, they are more likely prone to react. We used the mevalonate analogue structure to model a mevaldehyde intermediate for the second reduction step. The side-chain protonation states of The Neutral Lys639 System were taken from the most favorable model of Oliveira et al. [20], which was built for the study of the first reduction step, whereas The Neutral His405 System and The Cationic His405 System (Figure 1) are based on previous studies [2,3,12,14-17,19,20]. We modeled in total three different systems, all starting from the usually Figure 1. The reaction mechanism catalyzed by the HMG-CoA reductase based on the previous studies of Oliveira et al. [20] and Haines et al. [12], as described in detail in [21]. His866/His405 protonates -SCoA, whereas Lys691\Lys639 (-) or Glu559/Glu98 (-) is a proton donor after the hydride transfer. Asp767 corresponds to Asp715, and Glu559 corresponds to Glu98 in the current study.

Model Preparation
All systems were modeled from a human wild-type structure with Protein Data Bank (PDB) ID: 1DQA (resolution 2.0 Å), which represents the human enzyme in two dimeric structures (dimers AB and CD) with bound NADP+ cofactor, CoA, and a mevalonate analogue. Dimer AB was chosen to model the system, as chains A and B have the smallest distances between the mevalonate analogue and the cofactor, and thus, they are more likely prone to react. We used the mevalonate analogue structure to model a mevaldehyde intermediate for the second reduction step. The side-chain protonation states of The Neutral Lys639 System were taken from the most favorable model of Oliveira et al. [20], which was built for the study of the first reduction step, whereas The Neutral His405 System and The Cationic His405 System ( Figure 1) are based on previous studies [2,3,12,[14][15][16][17]19,20]. We modeled in total three different systems, all starting from the usually accepted assumption that -SCoA leaves the active site in the thiolate form and is protonated in the solvent. The difference between the three systems was the protonation state of the active site, which is crucial to determine the reaction mechanism. In the first system, both Glu98 and Lys639 were neutral, while His405 was cationic. In the second system, Glu98 and His405 were neutral, while Lys639 was cationic, whereas in the third system, Glu98 and Lys639 were neutral, and His405 was cationic, as described in Table 1. The active site is represented in Figure 2. Table 1. Description of the differences among each modeled system. All systems included a mevaldehyde intermediate and NADPH.

Residues in the Active Site The Neutral Lys639 System The Neutral His405 System The Cationic His405 System
Glu98 Anionic Processes 2021, 9, x FOR PEER REVIEW 4 of 12 accepted assumption that -SCoA leaves the active site in the thiolate form and is protonated in the solvent. The difference between the three systems was the protonation state of the active site, which is crucial to determine the reaction mechanism. In the first system, both Glu98 and Lys639 were neutral, while His405 was cationic. In the second system, Glu98 and His405 were neutral, while Lys639 was cationic, whereas in the third system, Glu98 and Lys639 were neutral, and His405 was cationic, as described in Table 1 The active site is represented in Figure 2.

Residues in the Active Site The Neutral Lys639 System The Neutral His405 System The Cationic His405 System
Glu98

Molecular Dynamics Simulations
The protein, substrate, and cofactor were solvated using TIP3P [24] water molecules with a minimum distance of 10 Å from the protein and each face of the box. The system was neutralized with 8 Na + . The total size of the system used in the MD simulations was of ca. 87,000 atoms treated with periodic boundary conditions. Mevaldehyde and NADPH were parameterized using the Hartree-Fock (HF/6-31G(d)) method to calculate RESP atomic charges, whereas the antechamber from the AMBER12 package [25] was used to assign the atom types. We chose one representative structure of the production phase of the MD simulation, which was later used in all QM/MM calculations. This structure had

Molecular Dynamics Simulations
The protein, substrate, and cofactor were solvated using TIP3P [24] water molecules with a minimum distance of 10 Å from the protein and each face of the box. The system was neutralized with 8 Na + . The total size of the system used in the MD simulations was of ca. 87,000 atoms treated with periodic boundary conditions. Mevaldehyde and NADPH were parameterized using the Hartree-Fock (HF/6-31G(d)) method to calculate RESP atomic charges, whereas the antechamber from the AMBER12 package [25] was used to assign the atom types. We chose one representative structure of the production phase of the MD simulation, which was later used in all QM/MM calculations. This structure had the most favorable reaction distances and most appropriate geometry between the substrate and the cofactor, as well as catalytic residues in the active site (Glu98, His405, Lys639, and Asp715).
We used a four-stage refinement protocol in which the constraints on the enzyme were gradually removed. In the first stage, harmonic constraints have been applied to all atoms except the ones from water molecules with a force constant of 50 kcal/mol/Å 2 for a total of 2000 steps. In the second stage, these constraints were applied to all heavy atoms for 4000 steps, and in the third stage, they were applied to the Ca and N-type atoms (backbone α-carbon and nitrogen atoms) for 8000 steps. The fourth stage was a full-energy minimization, which ran for a maximum of 50,000 steps until the rms gradient was smaller than 0.0002 kcal/mol. The MD simulations were performed using AMBER12 software [25], where the temperature was increased from 0 to 310.15 K in small increments using the Langevin thermostat during the first 20 ps (equilibration phase) with 2 ps −1 collision frequency and a 2 fs timestep. After that, a production simulation considering an isothermal-isobaric (NPT) ensemble was run for 20 ns. The General Amber Force Field (GAFF) [26] and ff10 force fields [27] were used, as well as SHAKE algorithm which was used to constrain the bonds involving hydrogen atoms during the MD simulations.
The QM/MM model included A and B chains, mevaldehyde, and NADPH in each active site of the dimer. All the water molecules 8 Å away from the protein were cut from the system, and all the water molecules outside a 5 Å radius from the active site were fixed (frozen) during the QM/MM calculations. The high layer was described at the B3LYP/6-31G(d) level and consisted of 91-92 atoms, which included atoms from mevaldehyde, NADPH, Glu98, His405, Lys639, and Asp715 (a detailed list can be found in Table 1). The low layer was described with the AMBER force field implemented in Gaussian09 [42] and consisted of more than 20,000 atoms. The B3LYP functional was chosen because of its reliability in the study of biological systems [43].
All the structures were fully optimized to the minima or transition state (the latter starting from the structure of the higher energy point of linear transit scans) and confirmed through vibrational frequency calculations. The transition states were confirmed through the existence of a single imaginary frequency whose transition vector confirmed the associated chemical transformation, whereas the corresponding reactants and products were confirmed through vibrational frequency calculations with no imaginary frequencies.
The reaction coordinate was between NADPH's hydride (-:H) and mevaldehyde's carbon. Energies of the minima and the TS were further refined through single-point calculations using the more complete basis set 6-311++G(2d,2p) in the QM layer with added thermal corrections and entropy, the latter estimated through the rigid rotor/harmonic oscillator formalism. Consequently, the activation and reaction energies were calculated through the difference between the corresponding Gibbs free energy of TS or the product and reactant. The correction of energies by the addition of a larger basis set through a single point energy calculation done in geometries determined by smaller basis sets or less demanding methods is an obvious approximation that in some cases can yield problems. However, it is the most common approach in computational enzymatic catalysis by ONIOM QM/MM [44]. In this case, as seen in the results, it was done to confirm that the final energies and conclusions did not depend on the specific combination of method/basis set chosen.

Results
The goal of this study was to determine the mechanism of the second reduction step of the catalytic reaction of hHMGR. There are two different main assumptions, where -SCoA can either leave the active site or stay in it before the second reduction step takes place. There is not only one channel from which -SCoA can leave the active site, but two, and this can be observed using any visualization software that can draw and show the surface of the enzyme. However, whether the -SCoA stays in the active site or not will probably not affect the hydride transfer and reaction itself, as it is approximately 4 Å far from the active site and atoms that are involved in the reaction. Our starting point was that the -SCoA leaves the active site.
The reaction itself involves two separate hydride transfers and a proton exchange that could be concerted, as shown in previous studies with some differences in the proton donor for mevaldehyde, depending on the class of the HMGR. According to the previous studies, the proton donor for the prokaryotic enzyme should be Glu83/Glu98 [3,12], whereas in hHMGR, it should be Lys691/Lys639 [20]. All the results obtained through QM/MM studies are described below. We explored three different hypotheses for the mechanism, each of them involving a different protonation state of the His405 and Lys639. Therefore, three systems were modeled: The Neutral His405 System, The Neutral Lys639 System, and The Cationic His405 System.
We tested the hypothesis where Lys639 was neutral, Asp715 was anionic, Glu98 was neutral, and His405 was cationic (The Neutral Lys639 System). We expected this model to be the least favorable, as there are two negative charges, those of Asp715 and mevaldehyde, which are close to each other without cationic Lys639 to stabilize them. In this case, the system included 91 atoms in the high layer and the only difference from other systems is the protonation state of the Lys639. Oliveira et al. [20] showed that Lys691/Lys639 protonates a mevaldyl-CoA intermediate and thus becomes neutral. Hence, this system is a direct follow up from the first hydride transfer fully characterized by Oliveira et al. [20]. The reaction itself is quite similar to the proposed mechanism for P. Mevalonii by Tabernero et al. [18]. Despite the linear transit scan, which indicated that the activation and reaction energies are 69.2 kcal/mol and 47.1 kcal/mol, respectively, the TS structure was not reasonable and therefore not taken into account.

The Neutral His405 System
In this system, His405 is neutral, whereas all the remaining residues in the active site are in the same protonation states as in The Cationic His405 System; thus, Glu98 was neutral, Lys639 was cationic, and Asp715 was anionic. Assuming that the reaction mechanism follows the study of Oliveira et al. [20], after the first hydride transfer, His405 is protonated by another residue or solvent and protonates Glu98, which protonated Lys639 before the second hydride transfer. Proton exchange between Glu98 and Lys639 would be immediate due to the pKa coupling. After the hydride transfer from NADPH to mevaldehyde, which stabilizes its negative charge through proton exchange with Glu98, we obtain mevalonate. Indeed, the Glu98 stays deprotonated close to the negatively charged Asp715, as well as close to His405 and Lys639. While modeling this system, the protonation state of His405 was chosen to maintain the strong interaction between His405's hydrogen and Glu98's oxygen and a possible hydrogen bond.
In more detail, in both reactants and the transition state structure, Glu98 and mevaldehyde are 3.20 Å and 2.41 Å apart, whereas the newly formed bond between Glu98-H and O-mevaldehyde is 1.01 Å. In addition, Glu98 stays close to His405, with the distance ranging from 2.03 Å in the reactants to 1.79 Å in the products, whereas the distance between Glu98 and Lys639 changes from 2.70 to 1.80 Å. Moreover, Glu98 is quite far from the negatively charged Asp715, and the distance does not change greatly in both reactants and products, 3.82 Å and 3.78 Å, respectively. In contrast, Lys639 is cationic and thus stays extremely close to the anionic Asp715 during the entire reaction with distances of 1.53 Å and 1.63 Å, respectively. Mevaldehyde is relatively close to the Lys639, with distances of Processes 2021, 9, 1085 7 of 12 2.13 Å and 2.41 Å in reactants and products. Interestingly, we observed Glu559/Glu98 as a proton donor and not Lys691/Lys639 as in the previous study by Oliveira et al. [20] The Lys691/Lys639 as a proton donor was explained through the difference in the class of reductase [20], but it can be further explained in the absence of strong interaction and hydrogen bond between His405 and Glu98 as we chose the protonation state of the His405, which would maintain this interaction throughout the simulation.
In addition, we observed a change in the Mulliken charges of all the atoms involved in the reaction mechanism ( Figure S1). The charge of O-mevaldehyde changes from −0.4814 in the reactants to −0.2935 in the products, whereas in C-NADPH, it goes from −0.0904 in reactants to 0.0571 in products. Hence, we expected this reaction mechanism to be the most favorable one.
The obtained activation and reaction Gibbs free energy were calculated as 18.0 kcal/mol and -16.0 kcal/mol, as shown in Figure 3. There is no proton exchange between His405 and Glu98, although they remain rather close during the entire simulation.
between Glu98 and Lys639 changes from 2.70 to 1.80 Å. Moreover, Glu98 is quite far from the negatively charged Asp715, and the distance does not change greatly in both reactants and products, 3.82 Å and 3.78 Å, respectively. In contrast, Lys639 is cationic and thus stays extremely close to the anionic Asp715 during the entire reaction with distances of 1.53 Å and 1.63 Å, respectively. Mevaldehyde is relatively close to the Lys639, with distances of 2.13 Å and 2.41 Å in reactants and products. Interestingly, we observed Glu559/Glu98 as a proton donor and not Lys691/Lys639 as in the previous study by Oliveira et al. [20] The Lys691/Lys639 as a proton donor was explained through the difference in the class of reductase [20], but it can be further explained in the absence of strong interaction and hydrogen bond between His405 and Glu98 as we chose the protonation state of the His405, which would maintain this interaction throughout the simulation.
In addition, we observed a change in the Mulliken charges of all the atoms involved in the reaction mechanism ( Figure S1). The charge of O-mevaldehyde changes from -0.4814 in the reactants to −0.2935 in the products, whereas in C-NADPH, it goes from -0.0904 in reactants to 0.0571 in products. Hence, we expected this reaction mechanism to be the most favorable one.
The obtained activation and reaction Gibbs free energy were calculated as 18.0 kcal/mol and -16.0 kcal/mol, as shown in Figure 3. There is no proton exchange between His405 and Glu98, although they remain rather close during the entire simulation. Reaction mechanism for a system with neutral His405 with corresponding imaginary frequency and distances between atoms important for reaction. All three structures represent optimized structures and the imaginary frequency corresponds to the TS. The reaction coordinate followed was marked in red, whereas newly formed bonds are marked in orange. After the hydride transfer and the protonation by Glu98, mevalonate is formed. Figure 3. Reaction mechanism for a system with neutral His405 with corresponding imaginary frequency and distances between atoms important for reaction. All three structures represent optimized structures and the imaginary frequency corresponds to the TS. The reaction coordinate followed was marked in red, whereas newly formed bonds are marked in orange. After the hydride transfer and the protonation by Glu98, mevalonate is formed.

The Cationic His405 System
In the last system studied, both His405 and Lys639 were cationic, whereas Glu98 was neutral and Asp715 was deprotonated and thus anionic. Although the Gibbs free activation energy we obtained corresponds and follows the previous data, again, we observed neutral Glu559/Glu98 as a proton donor for mevaldehyde, which is opposed to Lys691/Lys639 in the study of Oliveira et al. [20].

The O-H bond between Glu98-H and O-mevaldehyde is
Processes 2021, 9, 1085 8 of 12 0.98 Å, whereas the bond length between hydride (H:-) of an NADPH and mevaldehyde; thus, the newly obtained mevalonate is 1.09 Å. As to that, we observed the stabilization of Glu98 through proton exchange with His405 (1.02 Å), which then becomes neutral. Mevaldehyde is in the proximity of the Lys639 and Glu98 during the entire simulation around 2.13 Å and 2.17 Å in reactants and 0.98 Å and 1.98 Å in products, respectively. Lys639 remains next to Asp715 with 1.51 Å and 1.98 Å distances in both reactants and products. In addition, it seems to keep the mevaldehyde in the right conformation to facilitate the proton exchange between mevaldehyde and Glu98 and therefore the hydride transfer. As expected, the negatively charged Asp715 stays rather far from the Glu98 (more than 3.5 Å). The most important interactions and distances can be seen in Figure 4. and mevaldehyde; thus, the newly obtained mevalonate is 1.09 Å. As to that, we observed the stabilization of Glu98 through proton exchange with His405 (1.02 Å), which then becomes neutral. Mevaldehyde is in the proximity of the Lys639 and Glu98 during the entire simulation around 2.13 Å and 2.17 Å in reactants and 0.98 Å and 1.98 Å in products, respectively. Lys639 remains next to Asp715 with 1.51 Å and 1.98 Å distances in both reactants and products. In addition, it seems to keep the mevaldehyde in the right conformation to facilitate the proton exchange between mevaldehyde and Glu98 and therefore the hydride transfer. As expected, the negatively charged Asp715 stays rather far from the Glu98 (more than 3.5 Å). The most important interactions and distances can be seen in Figure 4.
The Gibbs free activation energy obtained with the B3LYP/6-311++G(2d,2p) functional and with added thermal correction is 15.0 kcal/mol, whereas the Gibbs free reaction energy is −40.0 kcal/mol. Correspondingly, we superimposed the crystallographic structure PDB ID 1DQA with the fully optimized structure of the products, where we obtained the RMSd of 2.95 Å and 1.48 Å for the entire system and the backbone of the protein, respectively. There is a notable change in the Mulliken charges of each atom involved in the reaction mechanism, as seen in Figure S2. Reaction mechanism for a system with cationic His405 with corresponding imaginary frequency and distances between atoms important for reaction. All three structures represent optimized structures, and the imaginary frequency corresponds to the TS. The reaction coordinate followed was marked in red, whereas newly formed bonds are marked in Figure 4. Reaction mechanism for a system with cationic His405 with corresponding imaginary frequency and distances between atoms important for reaction. All three structures represent optimized structures, and the imaginary frequency corresponds to the TS. The reaction coordinate followed was marked in red, whereas newly formed bonds are marked in orange. After the hydride transfer and the protonation by Glu98, mevalonate is formed with the proton exchange between His405 and Glu98.
The Gibbs free activation energy obtained with the B3LYP/6-311++G(2d,2p) functional and with added thermal correction is 15.0 kcal/mol, whereas the Gibbs free reaction energy is −40.0 kcal/mol. Correspondingly, we superimposed the crystallographic structure PDB ID 1DQA with the fully optimized structure of the products, where we obtained the RMSd of 2.95 Å and 1.48 Å for the entire system and the backbone of the protein, respectively. There is a notable change in the Mulliken charges of each atom involved in the reaction mechanism, as seen in Figure S2.

Discussion
In this extensive study of the second hydride transfer of the HMGR's reduction mechanism, we tested several different hypotheses which included different protonation states of His405 and Lys639. The obtained results suggest that the system with cationic His405 is the most favorable. There is a possibility of a small error in the obtained energies of 18.0 kcal/mol and −16.0 kcal/mol for activation and reaction Gibbs free energies (The Neutral His405 System), but it seems that the energies of the last system, The Cationic His405 System, would still be more favorable-the activation Gibbs free energy was 15.0 kcal/mol, whereas the Gibbs free reaction energy was −40.0 kcal/mol. In both systems, we observed Glu98 as a proton donor after a hydride transfer, which does not follow the study of Oliveira et al. [20] where Lys691/Lys639 would be a proton donor, but instead follows the previous study of Haines et al. [12]. We chose a protonation state of His405 in The Neutral His405 System to maintain strong interaction between His405 and Glu98, which facilitated the proton exchange between Glu98 and mevaldehyde. As to that, in the most favorable system, His405 is cationic, which again supports Glu98 as a proton donor, as His405 exchanges a proton with Glu98, and there is additional stabilization by Lys639. A choice of starting structure or the lack of -SCoA could impact the result. Oliveira et al. [20] stated that the first hydride transfer occurs with the activation barrier of 17.8 kcal/mol, whereas Haines et al. [3,12] confirmed that the second hydride transfer in PmHMGR occurs with the energy barrier of 2.9 kcal/mol lower than in the first hydride transfer, which confirms this hypothesis as a second reduction step.
Haines et al. [12] showed that two hydrogen bonds or a highly acidic group for proton transfer were necessary for the stabilization of the intermediate, although Oliveira et al. [20] concluded that the reaction is highly asynchronous. The change in the atomic charge at Glu98's carbon is minimal as it goes from 0.6033 in reactants to 0.6301 in both transition state and products. This would perhaps imply that the reaction proceeded in two steps rather than one. In the first step, NADPH donates a hydride to mevaldehyde, which then stabilizes its negative charge through proton exchange with Glu98 in the following step. A comparison of the energy profiles of the two hypotheses can be seen in Figure 5. In the end, the question that remains unsolved is whether -SCoA stays in the active site or leaves the active site before the proton exchange. This theory should be further tested as we observed two channels while looking at the surface of the HMG-CoA reductase, which suggests that the -SCoA could stay inside the active site. Presumably, it should not affect the reaction itself, as it is too far from the active site and the atoms involved in the reaction mechanism. The Cationic His405 System (−)).

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Figure S1: The mulliken charges for each atom involved in the reaction mechanism for The Neutral His405 System (left) with a corresponding optimized structure of the products (right). Figure S2: The mulliken charges for each atom involved in the reaction mechanism for The Cationic His405 System (left) with a corresponding optimized structure of the products (right). Figure S3: