Next Article in Journal
Facile Synthesis of MXene/MnO2 Nanocomposites for Efficient Removal of Radionuclide Uranium
Next Article in Special Issue
Self-Assembly and Conformational Change in the Oligomeric Structure of the Ectodomain of the TBEV E Protein Studied via X-ray, Small-Angle X-ray Scattering, and Molecular Dynamics
Previous Article in Journal
Thermal Decomposition and Solidification Characteristics of BFFO
Previous Article in Special Issue
Probing the Role of a Conserved Phenylalanine in the Active Site of Thiocyanate Dehydrogenase
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

QM/MM Study of a Nucleophilic Substitution Reaction Catalyzed by Uridine Phosphorylase from Vibrio cholerae

by
Alexander A. Lashkov
1,*,
Polina A. Eistrich-Geller
1,
Valeriya R. Samygina
1,2 and
Sergey V. Rubinsky
1
1
A.V. Shubnikov Institute of Crystallography, Federal Scientific Research Centre “Crystallography and Photonics”, Russian Academy of Sciences, 59, Leninskii Prospect, 119333 Moscow, Russia
2
National Research Centre “Kurchatov Institute”, 1, Akademika Kurchatova pl., 123182 Moscow, Russia
*
Author to whom correspondence should be addressed.
Crystals 2023, 13(5), 803; https://doi.org/10.3390/cryst13050803
Submission received: 6 April 2023 / Revised: 2 May 2023 / Accepted: 5 May 2023 / Published: 11 May 2023
(This article belongs to the Special Issue Protein Crystallography: Achievements and Challenges (Volume II))

Abstract

:
Uridine phosphorylases are used for biotechnological synthesis of pyrimidine derivatives and, moreover, their substrates and inhibitors are used in medicine. Therefore, studies of the mechanisms of the chemical reaction catalyzed by the enzyme and its specificity for various substrates are relevant. The research into the enzymatic reaction main stage—nucleophilic substitution of the nitrogenous base in uridine with an orthophosphate or orthovanadate group by hybrid QM/MM methods—was carried out. A comparison of various levels of theory and calculation schemes showed that preliminary optimization of the reactants’s geometry, as well as calculation of the initial trajectory of the minimum energy path, can be achieved by semi-empirical methods. At the same time, for the minimum energy path clarification, transition state geometry optimization, and calculation of the thermochemical parameters, it is preferable to use density functional theory in combination with modern ab initio methods. In comparison with the calculations of the activation barrier carried out in a solvent without an enzyme, differences in the kinetics of the enzymatic reaction due to the orientation and concentration actions of amino acid residues of the enzyme were revealed. This led to lowering the activation barrier by 20 kcal/mol and contributed to the reaction under physiologically acceptable conditions. It was shown that the free activation energy during the nucleophilic attack for uridine with hydrovanadate ion is 2 kcal/mol lower than for the hydrophosphate ion and this is consistent with the literature data.

1. Introduction

Uridine phosphorylase (UP) is an enzyme of the pyrimidine nucleoside phosphorylase family found in eukaryotic and prokaryotic cells. These enzymes catalyze the reversible phosphorolytic cleavage of pyrimidine nucleosides and their derivatives. UP can be used as a target for antitumor and antiparasitic drugs [1,2] and also as a biocatalyst for the production of pharmacologically active compounds, pyrimidine nucleoside derivatives. An understanding of the mechanism of enzyme action will help to optimize the biocatalysis and develop specific regulators of UP activity.
The study of catalytic processes in enzymes at the structural and functional levels using experimental methods requires further development of time-resolved protein crystallography, in particular, free-electron lasers (FELs) and serial data collection. Nevertheless, it is currently possible to preliminarily describe the enzymatic reaction pathway and calculate the energy parameters of the process and/or its individual steps using hybrid quantum-mechanical/molecular-mechanical (QM/MM) methods. The choice of methods and levels of theory for QM calculations is a difficult task, because of their wide variety; besides, the accuracy of calculations using different methods and their performance may differ by several orders of magnitude.
We studied the nucleophilic substitution as a key step of the UP-catalyzed enzymatic reaction using hybrid QM/MM calculations. The accuracy and calculation speed at different levels of theory for the description of quantum subsystems to calculate the reaction pathway, electronic energy, thermochemical parameters of the reaction, and geometry optimization were compared.
In the work [3], the experimental study of kinetics and, in particular, the kinetic isotope effect of uridine phosphorolysis reaction, catalyzed Trypanosoma cruzi UP was accompanied by quantum–chemical calculations of this reaction in the gas phase. The transition state was searched by scanning the potential energy surface and was analyzed. It was found that the nucleophilic substitution reaction proceeds via the SN2 mechanism. In the study [4], the structure of the conformer bounded by UP from Trypanosoma cruzi was calculated and its differences to the free uridine conformer were analyzed. The charge distribution obtained from the natural bond orbital populations was analyzed. However, these works did not take into account the effect of charged groups of amino acid residues (aa), as well as the solvent effect. In addition, the scanning of the potential energy surface along at least two coordinates is not the best method with which to find the transition state in the case when high-resolution structural data are available, as is the case for homologous bacterial UP from V. cholerae [5,6] and E. coli [7]. We believe that our work overcomes the above shortcomings.
In the recent work [8], it was experimentally shown that the orthovanadate ion can serve as a substrate for pyrimidine phosphorylases along with the phosphate ion. Up to a certain concentration, the vanadate ion is a more efficient substrate than the phosphate ion for some pyrimidine phosphorylases. Our paper compares the nucleophilic substitution reactions catalyzed by bacterial UP with these two substrates based on QM/MM calculations. There are differences in both the free energy of activation and the reaction energy.

2. Materials and Methods

2.1. Preparation and Parameterization

The structure of bacterial UP from Vibrio cholerae with a sulfate ion and uridine previously solved by X-ray diffraction analysis (RCSB PDB ID: 5LHV) was used as the starting model of the UP + uridine + phosphate ion double complex. All dimers of the hexameric quaternary structure were removed from the structure except for one, and water molecules except for two, which were constantly present in the active site of the bacterial UPs. The sulfate ion was replaced by an orthophosphate ion. Hydrogen atoms were added to the model of the complex using the CHARMM-GUI WEB server [9]. Two protonated forms were calculated for the phosphate ion: phosphate ion and monohydrogen phosphate ion. Considering the literature data, the form of the monohydrogen phosphate ion was taken as the main form [2,3]. The model of the complex was parameterized using the CHARMM-GUI server in the CHARMM36m force field [10] for amino acid residues and CgenFF [11] for ligands. Additional water molecules were described using an explicit three-point water model—TIP3P and were added to a cubic box (with an edge of 88 Å), involving the structure of the complex. Additionally, Na+ and Cl ions were added to the box at a physiological concentration of 0.1 M to neutralize the total charge of the system. Preliminary geometry optimization of the system was performed using the NAMD program 2.1.14 [12]. Equilibration at a temperature of 303.15 K was carried out for 2500 picoseconds with an integration step of 2 femtoseconds. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME, [13]) method. The LJ-potential with forces smoothly switched to zero in the range 10–12 Å was used to describe van der Waals interactions.
Further calculations were carried out using the Orca [14] quantum chemical calculation software package version 4.2.1. [15]. The studied model was divided into quantum and classical parts. The quantum part included the atoms of the ligands and the chemical groups of the active site residues (His7*, Gly25, Arg29, Arg 47*, Arg90, Arg92, Thr93, Gln165, Arg167, Glu197, Arg222, in which * indicates amino acid residues from another subunit, Figure 1a) of the enzyme bonded through donor–acceptor (hydrogen) bonds. The quantum part also included two water molecules that are constantly present in the active site of bacterial UPs and are possibly directly involved in the enzymatic reaction. The quantum part consisted of 150 atoms. The classical part included all other atoms of the system, but the region to be optimized contained only the rest of the atoms of amino acid residues involved in the quantum part and covalently bounded to them in the main chain (6*, 8, 24, 26, 28, 30, 46*, 48*, 89, 91, 164, 166, 168, 196, 198, 221, 223), as well as the active site fragment (94–98 aa). The shifted Coulomb’s potential was used at a cut-off of 12 Å to calculate electrostatic interactions. Hydrogen atoms were used as link-atoms for treating the QM/MM boundary. Electrostatic embedding was used. The total electronic energy of the system was calculated using the additive scheme as the sum of the energies of quantum and classical parts.
All models for the protein complexes with substrates and products of the enzymatic reaction were obtained from one system—UP + uridine + monohydrogen phosphate ion—by editing the coordinates and types of atoms and then optimizing the geometry using the QM/MM method. The structure of the E. coli UP complex with the corresponding ligands (RSCB PDB ID: 1TGY) was used as the conformational template and the mutual arrangement of the reaction products (ribose-1’-phosphate and uracil).

2.2. Levels of Theory Used

Several theory levels were used to compare the accuracy and timings of calculations. The fastest one was the semi-empirical extended tight-binding method (xTB) with the GFN2 parameterization [16,17]. This method was recently proposed and is intended for calculating the atomic geometry and vibrational frequencies in large (about 1000 atoms) organic and organometallic systems. The slowest one is based on the density functional theory (DFT) with hybrid functional B3LYP, the 6-31**G basic set for carbon atoms and the 6-311++G** basic set for the rest of the atoms of the quantum part. DFT integration grids were set to GRID4 and FINALGRID5. The D3(BJ) empirical dispersion correction [18] was applied. A combined method based on DFT-B97 with the optimized basis set def2-mTZVP: B97-3c [19] was intermediate in speed. To accurately estimate the electronic energy, the coupled cluster method CCSD(T) with the domain-based local pair natural orbital (DLPNO) computational scheme [20], the ma-def2-TZVPP basis set, and the auxiliary def2-TZVPP/C basis set [21] was used. The DPLNO computational scheme allows precision calculations of the electronic energy in supercomputers for organic systems consisting of hundreds of atoms, while the standard CCSD(T) scheme is applicable only to systems containing several tens of atoms. The disadvantage of the latter scheme is the difficulty in calculating the energy gradient and Hessian, which greatly limits the applicability of this scheme. The TightOpt convergence criteria for the self-consistent field (SCF) was used at all levels of theory.

2.3. Calculations of the Reaction Path and the Transition State

Calculations of the minimum energy path (MEP) for the nucleophilic substitution reaction were performed using the NEB-TS method [22]. The latter is a combination of the climbing image-nudged elastic band method (NEB-Cl, [23]) and the transition state geometry optimization using the eigenvector following (EF) method. Sixteen intermediate images were used to describe the reaction pathway. The initial trajectory was generated by the image-dependent-pair-potential (IDPP, [24]) or a combination of IDPP with the NEB method at the xTB-GFN2 level of theory. Before calculating the MEP, the models of the complexes were optimized at different levels of theory. The correspondence of the structure to a local energy minimum was validated by calculating the partial Hessian, as described in the next section. The partial Hessian was also calculated for the transition state obtained using the NEB-TS method.

2.4. Calculations of Thermochemical Parameters

To calculate the zero point energy and the oscillatory entropy, as well as the correction for non-zero absolute temperature, a partial Hessian (PHVA, [25]) was calculated numerically for the quantum part and the optimized region of the classical part. In partial Hessian calculations, atoms located at distances shorter than the sum of the van der Waals radii of atoms + 1 Å from the atoms of the optimized region were recorded as immobile atoms. The quasi-rigid rotor harmonic oscillator approximation (qRRHO, [26]) was used to calculate the oscillatory entropy of the system. The Hessian calculations were carried out at the B3LYP/6-311++G**/D3 level of theory only for one model, UPhase + uridine + phosphate ion, since calculations at this level of theory are difficult to perform because of their resource intensity and duration.

2.5. Calculations for the Nucleophilic Attack Reaction without an Enzyme

The reaction path calculations, the geometry optimization of the transition state, and the evaluation of the thermochemical parameters of the nucleophilic substitution reaction in water without the enzyme were performed at the B3LYP/6-311++G** (with D3(BJ) empirical correction) and B97-3c. The electronic energy was calculated at the DPLNO-CCSD(T)/ma-def2-TZVPP/def2-TZVPP/C level of theory. The Hessian was calculated analytically. The solvent effect was taken into account using the implicit C-PCM model [27] with a set of parameters for water.

3. Results and Discussion

3.1. Optimized Structures of Substrate Complexes and Products with UP

Geometry optimization of the QM/MM systems of the protein complexes with products and substrates was carried out at three different levels of theory (GFN-XTB2, B3LYP/6-311++G**/D3, B97-3c) for the QM part using the CHARM36m force field for the MM part. The RMSD for the atoms of the quantum part did not exceed 0.1 Å at different levels of theory.
Figure 1a shows the structure of the QM-part of the V. cholerae UP complex with uridine and monohydrogen phosphate; Figure 1b—the structure of the putative pre-reaction state in a solution without an enzyme (B97-3c level of theory). It is worth noting that, in water without an enzyme, the hydrogen phosphate ion tends to form hydrogen bonds with hydroxyl groups of the ribose moiety of uridine, which lowers the energy of the pre-reaction state compared to the simple sum of the electronic energies of the reactants in solution. In combination with the protein, there is only one such hydrogen bond between the substrates due to the orienting effect of the active site residues—Glu197 on the ribose moiety of uridine and the arginines on the hydrogen phosphate ion (Figure 1b).
There are differences in the mutual arrangement of the attacking oxygen atom of the hydrogen phosphate ion and the glycosidic carbon of uridine (Figure 1)—in the structure of the putative pre-reaction complex in solution, they are closer to each other by 0.6 Å. On the one hand, this fact is not evidence for the enzymatic process. On the other hand, in the pre-reaction complex in water, the angular position is less favorable for the attack of the hydrogen phosphate ion on a N-C-glycoside bond compared to the enzyme active site.
Differences are also observed in the  θ ( O 4 C 1 N 1 C 6 )  dihedral angle of uridine (−105.4° in an aqueous solution and −78.6° in combination with the enzyme, Figure 1). This change in the angle is apparently due to the joint effect of uridine binding by the conserved residues—Gln165 and Arg167 of the uracil-binding site and Glu197 of the ribosyl binding site. According to the  θ ( O 4 C 1 N 1 C 6 )  dihedral angle [28], the uridine molecule in the enzyme adopts a higher energy conformation compared to that in an aqueous solution without the enzyme. The  θ ( O 4 C 1 N 1 C 6 )  dihedral angle of uridine in complex with the enzyme is close to the  θ ( O 4 C 1 N 1 C 6 ) = 75  observed for the Michaelis complex [4]. Meanwhile, the dihedral angle in an aqueous solution differs from −117°, which can be attributed to the effect of the hydrogen phosphate ion on the conformation of uridine in the putative pre-reaction complex.
To estimate the difference in the electronic energies of the conformations of uridine, a calculation was carried out at the level of the DPLNO-CCSD(T) theory using the implicit water model C-PCM of bounded and free conformers. The energy of the bounded conformer was −909.7026 Ha and that of the free conformer was −909.7215 Ha ( Δ E  = 11.86 kcal/mol). This calculation further confirms that the uridine molecule in the enzyme adopts a higher energy conformation compared to that in an aqueous solution without the enzyme.
The N-C glycosidic bond lengths are nearly equal in both structures, but these bonds have different polarity. According to the calculated Mulliken charges, the glycosidic bond is polarized in the presence of the active site residues of UP, while the glycosidic carbon acquires a larger positive charge. Although different types of charges obtained at different levels of theory cannot be quantitatively compared, there are differences in the qualitative charge distribution on uridine atoms in solution and the enzyme (Figure 1) compared to the data presented in [4]. In this study, the charge distribution and polarity of the N-C glycosidic bond of uridine are the same in the associated and unbound state of UP. However, the authors of the study [4] did not take into account the polarizing effect of the protein environment of the substrate, in particular the Gln165 and Arg167 side-chain groups, and the resulting redistribution of charges along the aromatic ring.
On the contrary, the oxygen atom of the hydrogen phosphate ion involved in the nucleophilic attack becomes more nucleophilic under the influence of guanidine groups of arginines 29, 47, 90 of the phosphate-binding site of UP (Figure 1). This leads to an increase in the rate of the nucleophilic substitution reaction. Additionally, the higher rate of the direct reaction in the presence of the enzyme is due to the higher electronic energy of the uridine conformer in a complex with the enzyme, compared to the minimum energy of the uridine conformer in solution (Table 1).

3.2. Reaction Pathway and Transition State for the Nucleophilic Attack of Hydrogen Phosphate Ion on Uridine

The Hessian calculations showed the absence of imaginary frequencies in the structures under consideration. Therefore, they were used to construct the MEP for the transition of reactants into products and for the preliminary determination of the saddle point corresponding to the transition state. The energy profiles along the reaction coordinate calculated by the DFT method are given at two levels of theory (B3LYP/6-311++G**/D3 and B97-3c) for the nucleophilic substitution reaction in the aqueous environment (Figure 2, left) and in the presence of an enzyme (Figure 2, right). The change in the reaction coordinates during the substrate–product transition in this system corresponds to an increase in the distance between C1-atom of ribozyl moiety and N1-atom of uracil and a decrease in the distance from C1-atom of ribozyl and O-atom of hydrophosphate ion. As expected, the nucleophilic substitution of uridine proceeds through the SN2 mechanism without the formation of an intermediate with one peak in the plot corresponding to the transition state. Similar energy profiles were obtained by DFT calculations at two levels of theory, and their maxima are shifted relative to each other by 4 kcal/mol. The energy profile of the reaction in an aqueous solution significantly differs from that in the protein environment. The maximum energy at the same level of theory was found to be three times lower in the presence of the enzyme, confirming and supplementing the results and conclusions made earlier when analyzing the reactant complexes.
As opposed to MEP calculations using DFT, the profiles obtained using the semi-empirical XTB-GFN2 method were found to be different with a significantly (three times) larger maximum of the electronic energy. However, a comparison of the positions of atoms in the images at the maximum point of the profile showed a fairly good consistency of the XTB-GFN2 method with DFT (RMSD = ∼0.2 Å) in the calculations of the geometry of transition states with a significant difference in the calculation speed (∼70×).
The transition state obtained using the NEB-CI method was not optimized; in addition to the imaginary frequency corresponding to the transition state, several smaller imaginary frequencies were determined. The molecular system in the transition state was subjected to further geometry optimization using the EF method. The vibrational frequency calculations for the optimized transition states of both ligand complexes in water and in the enzyme showed the presence of one imaginary frequency corresponding to the saddle point of the system in the transition state.
It is worth noting that the vibrational frequency calculations should be carried out, if possible, at the level of theory, at which the geometry was optimized. However, for approximate QM/MM calculations of enzymes without heavy element co-factors, one can use semi-empirical methods parameterized based on thermochemical data, such as xTB with the parameterization of GFN or GFN2, to estimate the vibrational frequencies.
A comparison of various levels of theory and methods for calculating the electronic reaction energies ( Δ E R ) and activation barriers ( Δ E A ) was carried out. In the presence of the enzyme, different values were obtained (10 kcal/mol for DFT and 30 kcal/mol for xTB-GFN2, see Table 1). Meanwhile, the value calculated by the coupled cluster method was 14 kcal/mol (Table 1). In all cases, the semi-empirical xTB method overestimated the energy of the transition state compared to DPLNO-CCSD(T)/ma-def2-TZVPP, while DFT underestimated this energy. With the use of modern computing equipment and computational schemes, it has become possible to apply high-precision ab initio methods, such as the coupled cluster method, to calculate the difference in the electronic energy of large organic systems, which should be used when calculating the activation barriers. In the case transition states, as well as atypical bonding situations, these methods give a smaller error than DFT or semi-empirical methods. This is explained by the fact that these states are usually not included in training sets for functionals’ construction or selection of parameters in semi-empirical methods. As for the changes in the substrate–product energy, it is possible to use DFT with acceptable accuracy (1–2 kcal/mol) and even modern semi-empirical methods (Table 1).
The calculated vibrational frequencies were used to determine the entropy part of the free energy of the system and the contributions to the enthalpy part due to the movements of the nuclei of atoms for the determined transition states, reactants, and products. The free energy of activation and the reaction energies ( Δ G A Δ G R ), respectively, were calculated (Table 1).
The structures of the transition states are shown in Figure 3. It is interesting to compare the distances from the glycosidic carbon atom of the ribose moiety to the nearest oxygen atom of the phosphate ion (r(C-O)) and from the glycosidic carbon atom to the N1 atom of the heterocyclic base (r(C-N)) for the transition state in water solution and enzyme. The distances r(C-N) for the transition state in water and the enzyme are almost identical Figure 3, whereas the r(C-O) values were different. A larger value of r(C-O) for the transition state in the enzyme is apparently due to the strong electrostatic attraction between the hydrogen phosphate ion and the positively charged guanidine groups of arginine residues. In contrast, for the transition state in water, the hydrogen bonds between substrates, which were previously described for the pre-reaction complex, are retained. Due to these bonds, r(C-O) in water is reduced.
The distances r(C-O) and r(C-N) differ from those reported in [3], which can be attributed both to the influence of the environment (atoms of the residues and the solvent model) and different methods for calculating the transition state. It is interesting to note that r(C-O) is close to the value obtained in this work for the transition state in the enzyme, but not in the solution without the enzyme. This fact can be explained not only by the influence of the solvation model, but also by the fact that, in our work, the transition state in water was modeled taking into account the formation of a pre-reaction complex.
A solvent without UP was also calculated using methods of quantum chemistry, in order to take into account the dependence of the activation energy on the molecular translational entropy and the conformational changes, assuming that a stable pre-reaction complex is not formed [29]. For this purpose, the geometry of each substrate and product was optimized separately at the level of theory B3LYP/6-311++G**/D3. Then, the electronic and free energies were calculated for each molecule. However,  Δ G A  was almost identical (Table 1) to that calculated involving the pre-reaction complex. The  Δ G R  were different, but this is not of decisive importance, since  Δ G A  in solution without the enzyme is very high, making the reaction impossible under physiologically acceptable conditions.

3.3. Nucleophilic Substitution Reaction Involving Hydrogen Orthovanadat Anion

Figure 4 shows MEP for the nucleophilic substitution of the pyrimidine moiety of uridine with a hydrogen orthovanadate ion in comparison with the hydrogen phosphate ion, calculated at the B97-3c level of theory. At the beginning of the reaction pathway, the energy profiles are well consistent with each other, but once the activation barrier is overcome, significant differences are observed. The electronic energy difference between the substrates and products is negative and more modulo compared to that for the hydrogen phosphate ion. This is indicative of a shift of the chemical equilibrium towards the direct conversion of uridine into uracil. After the geometry optimization of the transition state and free energy calculations,  Δ G A  for the reaction with the hydrogen ortovanadate ion was reduced by 2 kcal/mol, mainly due to a change in the electronic energy of the optimized transition state (Table 1). A decrease in the activation barrier leads to a slight increase in the reaction rate of the nucleophilic attack involving the hydrogen orthovanadate ion compared to the hydrogen phosphate ion. However, it should be noted that, in the above calculations, the inhibitory effect of high concentrations of the vanadate ion on pyrimidine phosphorylases [8] has not appeared in any way. Apparently, this process is associated with the binding of vanadat outside the QM-simulated region of the enzyme molecule and requires further study using experimental structural biology methods.
The quantum subsystem of the optimized structure of the UP complex with uridine and hydrogen orthorvanadat-ion is presented in Figure 5a. No significant differences were found in the conformation of uridine and, accordingly, in the energy of the conformation when comparing the pre-reaction enzyme-substrate complexes of uridine with vanadate and with phosphate. Figure 5b shows the structure of the corresponding transition state. The distances r(C-O) and r(C-N) and the bonds orders remained unchanged in the case of the transition state involving the vanadate ion as compared with the phosphate ion.

4. Conclusions

The nucleophilic substitution as a key step of the UP-catalyzed enzymatic reaction was studied using hybrid QM/MM calculations. Different levels of the theory and schemes for calculating the electronic energy, MEP, thermochemical parameters, and geometry optimization of reactants, products, and the transition state was compared. A comparison of various levels of theory and calculation schemes showed that the preliminary optimization of the reactants’/products’ geometry as well as the calculation of the initial trajectory of MEP and obtained an approximate MEP/transition state that can be made by using high-throughput semi-empirical methods (e.g., xTB-GFN2). At the same time, for the MEP refinement, transition state geometry optimization, and calculation of the thermochemical parameters, it is preferable to use DFT in combination with modern ab initio methods.
A comparison of the results of QM/MM calculations carried out in a solvent without the enzyme and with the enzyme revealed the differences in the energy profile of the nucleophilic substitution reaction of uridine with a hydrogen phosphate ion and the structure of the transition state due to the orienting and concentration effects of the enzyme residues. These effects lead to a decrease in the activation barrier by 20 kcal/mol, allowing the reaction to occur under physiologically acceptable conditions. It was shown that a pre-reaction complex, with an energy reduced both relative to the sum of the energies of the reactants in water and the complex of reactants in the active site of UP, can be formed in an aqueous solution without the enzyme.
It should be noted that all calculations were carried out in the harmonic approximation, which does not always correspond to the real dynamics of protein molecules due to their high conformational mobility, and this issue needs further research and clarification.
The free energy of activation for the enzymatic reaction with hydrogen vanadate was reduced by 2 kcal/mol in comparison to hydrogen phosphate ion. A decrease in the activation barrier leads to a slight increase in the reaction rate of the nucleophilic attack involving the hydrogen vanadate ion compared to the hydrogen phosphate ion. It should be noted that, in the above calculations, the inhibitory effect of high concentrations of vanadate ions on pyrimidine phosphorylase is not manifested in any way and requires further study using experimental structural biology methods.

Author Contributions

Conceptualization, A.A.L. and V.R.S.; methodology, A.A.L.; software, A.A.L. and S.V.R.; validation, A.A.L. and S.V.R.; investigation, A.A.L. and P.A.E.-G.; writing—original draft preparation, A.A.L.; writing—review and editing, S.V.R. and V.R.S.; visualization, P.A.E.-G.; supervision, A.A.L. and V.R.S.; project administration, A.A.L. and V.R.S.; funding acquisition, V.R.S. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financially supported by the Russian Foundation for Basic Research (RFBR, project No. 19-29-12054; the assessment of the accuracy and performance of different QM/MM methods) and was performed within the State Assignment of FSRC “Crystallography and Photonics” RAS (in part of modeling the nucleophilic substitution reactions with various substrates).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the lack of the publication’s technical feasibility.

Acknowledgments

The research was carried out using the infrastructure of the Shared Research Facilities «High Performance Computing and Big Data» (CKP «Informatics») of the FRC CSC RAS (Moscow).

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
QM/MMQuantum Mechanics/Molecular Modelling
UPUridine Phosphorylase
aaamino acid residues
xTBextended Tight-binding method
DFTDensity Functional Theory
MEPMinimum Energy Path
NEB-CIClimbing Image Nudged Elastic Band method
EFEigenvector Following method
IDPPImage-Dependent-Pair-Potential method
PHVAPartial Hessian
qRRHOquasi-Rigid Rotor Harmonic Oscillator
URAUracil moiety
RIBRibosyl moiety
PMEParticle-Mesh Evald method
DPLNODomain-based Local Pair Natural Orbital computational scheme
CCSD(T)Coupled Cluster Single-Double-Triple method
SCFSelf-Consistent Field

References

  1. Lashkov, A.A.; Zhukhlistova, N.E.; Seregina, T.A.; Gabdulkhakov, A.G.; Mikhailov, A.M. Uridine phosphorylase in biomedical, structural, and functional aspects: A review. Crystallogr. Rep. 2011, 56, 560–589. [Google Scholar] [CrossRef]
  2. Paul, D.; Oleary, S.E.; Rajashankar, K.; Bu, W.; Toms, A.; Settembre, E.C.; Sanders, J.M.; Begley, T.P.; Ealick, S.E. Glycal formation in crystals of uridine phosphorylase. Biochemistry 2010, 49, 3499–3509. [Google Scholar] [CrossRef] [PubMed]
  3. Silva, R.G.; Vetticatt, M.J.; Merino, E.F.; Cassera, M.B.; Schramm, V.L. Transition-state analysis of Trypanosoma cruzi uridine phosphorylase- catalyzed arsenolysis of uridine. J. Am. Chem. Soc. 2011, 133, 9923–9931. [Google Scholar] [CrossRef]
  4. Silva, R.G.; Kipp, D.R.; Schramm, V.L. Constrained Bonding Environment in the Michaelis Complex of Trypanosoma cruzi Uridine Phosphorylase. Biochemistry 2012, 51, 6715–6717. [Google Scholar] [CrossRef] [PubMed]
  5. Prokofev, I.I.; Lashkov, A.A.; Gabdulkhakov, A.G.; Balaev, V.V.; Seregina, T.A.; Mironov, A.S.; Betzel, C.; Mikhailov, A.M. X-ray structures of uridine phosphorylase from Vibrio cholerae in complexes with uridine, thymidine, uracil, thymine, and phosphate anion: Substrate specificity of bacterial uridine phosphorylases. Crystallogr. Rep. 2016, 61, 954–973. [Google Scholar] [CrossRef]
  6. Eistrikh-Heller, P.A.; Rubinsky, S.V.; Samygina, V.R.; Gabdulkhakov, A.G.; Kovalchuk, M.V.; Mironov, A.S.; Lashkov, A.A. Crystallization in Microgravity and the Atomic-Resolution Structure of Uridine Phosphorylase from Vibrio cholerae. Crystallogr. Rep. 2021, 66, 777–785. [Google Scholar] [CrossRef]
  7. Caradoc-Davies, T.T.; Cutfield, S.M.; Lamont, I.L.; Cutfield, J.F. Crystal structures of Escherichia coli uridine phosphorylase in two native and three complexed forms reveal basis of substrate specificity, induced conformational changes and influence of potassium. J. Mol. Biol. 2004, 337, 337–354. [Google Scholar] [CrossRef]
  8. Antipov Alexey, N.; Okorokova Natalya, A.; Safonova Tatyana, N.; Veiko Vladimir, P. Vanadate as a new substrate for nucleoside phosphorylases. JBIC J. Biol. Inorg. Chem. 2022, 27, 221–227. [Google Scholar] [CrossRef]
  9. Lee, J.; Cheng, X.; Swails, J.M.; Yeom, M.S.; Eastman, P.K.; Lemkul, J.A.; Wei, S.; Buckner, J.; Jeong, J.C.; Qi, Y.; et al. CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field. J. Chem. Theory Comput. 2016, 12, 405–413. [Google Scholar] [CrossRef]
  10. Huang, J.; Rauscher, S.; Nawrocki, G.; Ran, T.; Feig, M.; De Groot, B.L.; Grubmüller, H.; MacKerell, A.D. CHARMM36m: An improved force field for folded and intrinsically disordered proteins. Nat. Methods 2016, 14, 71–73. [Google Scholar] [CrossRef]
  11. Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671–690. [Google Scholar] [CrossRef] [PubMed]
  12. Phillips, J.C.; Hardy, D.J.; Maia, J.D.; Stone, J.E.; Ribeiro, J.V.; Bernardi, R.C.; Buch, R.; Fiorin, G.; Hénin, J.; Jiang, W.; et al. Scalable molecular dynamics on CPU and GPU architectures with NAMD. J. Chem. Phys. 2020, 153, 044130. [Google Scholar] [CrossRef]
  13. Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089–10092. [Google Scholar] [CrossRef]
  14. Neese, F. The ORCA program system. WIREs Comput. Mol. Sci. 2012, 2, 73–78. [Google Scholar] [CrossRef]
  15. Neese, F. Software update: The ORCA program system, version 4.0. WIREs Comput. Mol. Sci. 2018, 8, e1327. [Google Scholar] [CrossRef]
  16. Bannwarth, C.; Ehlert, S.; Grimme, S. GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions. J. Chem. Theory Comput. 2019, 15, 1652–1671. [Google Scholar] [CrossRef]
  17. Bannwarth, C.; Caldeweyher, E.; Ehlert, S.; Hansen, A.; Pracht, P.; Seibert, J.; Spicher, S.; Grimme, S. Extended tight-binding quantum chemistry methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2021, 11, e1493. [Google Scholar] [CrossRef]
  18. Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. [Google Scholar] [CrossRef]
  19. Brandenburg, J.G.; Bannwarth, C.; Hansen, A.; Grimme, S. B97-3c: A revised low-cost variant of the B97-D density functional method. J. Chem. Phys. 2018, 148, 064104. [Google Scholar] [CrossRef]
  20. Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E.F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory. J. Chem. Phys. 2016, 144, 024109. [Google Scholar] [CrossRef]
  21. Zheng, J.; Xu, X.; Truhlar, D.G. Minimally augmented Karlsruhe basis sets. Theor. Chem. Acc. 2011, 128, 295–305. [Google Scholar] [CrossRef]
  22. Ásgeirsson, V.; Birgisson, B.O.; Bjornsson, R.; Becker, U.; Neese, F.; Riplinger, C.; Jónsson, H. Nudged Elastic Band Method for Molecular Reactions Using Energy-Weighted Springs Combined with Eigenvector following. J. Chem. Theory Comput. 2021, 17, 4929–4945. [Google Scholar] [CrossRef] [PubMed]
  23. Henkelman, G.; Uberuaga, B.P.; Jónsson, H. Climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 2000, 113, 9901–9904. [Google Scholar] [CrossRef]
  24. Smidstrup, S.; Pedersen, A.; Stokbro, K.; Jónsson, H. Improved initial guess for minimum energy path calculations. J. Chem. Phys. 2014, 140, 214106. [Google Scholar] [CrossRef] [PubMed]
  25. Li, H.; Jensen, J.H. Partial Hessian vibrational analysis: The localization of the molecular vibrational energy and entropy. Theor. Chem. Acc. 2002, 107, 211–219. [Google Scholar] [CrossRef]
  26. Grimme, S. Supramolecular binding thermodynamics by dispersion-corrected density functional theory. Chem.—Eur. J. 2012, 18, 9955–9964. [Google Scholar] [CrossRef] [PubMed]
  27. Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995–2001. [Google Scholar] [CrossRef]
  28. Yildirim, I.; Stern, H.A.; Kennedy, S.D.; Tubbs, J.D.; Turner, D.H. Reparameterization of RNA ξ Torsion Parameters for the AMBER Force Field and Comparison to NMR Spectra for Cytidine and Uridine. J. Chem. Theory Comput. 2010, 6, 1520–1531. [Google Scholar] [CrossRef]
  29. Ryu, H.; Park, J.; Kim, H.K.; Park, J.Y.; Kim, S.T.; Baik, M.H. Pitfalls in Computational Modeling of Chemical Reactions and How To Avoid Them. Organometallics 2018, 37, 3228–3239. [Google Scholar] [CrossRef]
Figure 1. Structure of the quantum part of the UPhase complex from V. cholerae with uridine (URI) and hydrogen phosphate ion (Pi) (a); the structure of a possible pre-reaction state in an aqueous environment without the enzyme (b). The Mulliken charges of atoms are shown in blue (for negative) or red (for positive). The value of the dihedral angle ( θ ( O 4 C 1 N 1 C 6 ) ) is indicated with an arc. Amino acid residues from another subunit is indicated by asterix.
Figure 1. Structure of the quantum part of the UPhase complex from V. cholerae with uridine (URI) and hydrogen phosphate ion (Pi) (a); the structure of a possible pre-reaction state in an aqueous environment without the enzyme (b). The Mulliken charges of atoms are shown in blue (for negative) or red (for positive). The value of the dihedral angle ( θ ( O 4 C 1 N 1 C 6 ) ) is indicated with an arc. Amino acid residues from another subunit is indicated by asterix.
Crystals 13 00803 g001
Figure 2. Electronic energy difference for the QM/MM system versus the nucleophilic substitution reaction coordinate in an aqueous solution without the enzyme (on the left) and in the presence of V. cholerae UP (on the right), involving uridine and a monohydrogen phosphate ion. The calculations were performed using the NEB-CI method at the B3LYP/6-311++G**/D3 (blue plot) and B97-3c (red plot) levels of theory.
Figure 2. Electronic energy difference for the QM/MM system versus the nucleophilic substitution reaction coordinate in an aqueous solution without the enzyme (on the left) and in the presence of V. cholerae UP (on the right), involving uridine and a monohydrogen phosphate ion. The calculations were performed using the NEB-CI method at the B3LYP/6-311++G**/D3 (blue plot) and B97-3c (red plot) levels of theory.
Crystals 13 00803 g002
Figure 3. Structure of the quantum part in the transition state for the V. cholerae UP complex with uridine (URA—uracil moiety, RIB—ribosyl moiety) and monohydrogen phosphate anion Pi (a); structure of the transition state for the nucleophilic substitution reaction of uridine with the monohydrogen phosphate anion in an aqueous environment without the enzyme (b). Amino acid residues from another subunit is indicated by asterix. The calculations were performed using the NEB-TS method at the B97-3c level of theory.
Figure 3. Structure of the quantum part in the transition state for the V. cholerae UP complex with uridine (URA—uracil moiety, RIB—ribosyl moiety) and monohydrogen phosphate anion Pi (a); structure of the transition state for the nucleophilic substitution reaction of uridine with the monohydrogen phosphate anion in an aqueous environment without the enzyme (b). Amino acid residues from another subunit is indicated by asterix. The calculations were performed using the NEB-TS method at the B97-3c level of theory.
Crystals 13 00803 g003
Figure 4. Electronic energy difference for the QM/MM system versus the nucleophilic substitution reaction coordinate in the presence of uridine phosphorylase from V. cholerae involving uridine and a monohydrogen phosphate anion (blue plot); involving uridine and a monohydrogen orthovanadate anion (red plot). The calculations were performed using the NEB-CI method at the B97-3c level of theory.
Figure 4. Electronic energy difference for the QM/MM system versus the nucleophilic substitution reaction coordinate in the presence of uridine phosphorylase from V. cholerae involving uridine and a monohydrogen phosphate anion (blue plot); involving uridine and a monohydrogen orthovanadate anion (red plot). The calculations were performed using the NEB-CI method at the B97-3c level of theory.
Crystals 13 00803 g004
Figure 5. Structure of the quantum part in the UP complex from V. cholerae with uridine (URA—uracil moiety, RIB—ribosyl moiety) and monohydrogen vanadate anion Vi (a) and its transition state (b). The value of the dihedral angle ( θ ( O 4 C 1 N 1 C 6 ) ) is indicated with an arc. Amino acid residues from another subunit is indicated by asterix. The calculations of the transition state were performed using the NEB-TS method at the B97-3c level of theory.
Figure 5. Structure of the quantum part in the UP complex from V. cholerae with uridine (URA—uracil moiety, RIB—ribosyl moiety) and monohydrogen vanadate anion Vi (a) and its transition state (b). The value of the dihedral angle ( θ ( O 4 C 1 N 1 C 6 ) ) is indicated with an arc. Amino acid residues from another subunit is indicated by asterix. The calculations of the transition state were performed using the NEB-TS method at the B97-3c level of theory.
Crystals 13 00803 g005
Table 1. Electronic ( Δ E ) and Gibbs free ( Δ G ) energies (in kcal/mol) of the reaction (R) and activation (A) in the nucleophilic substitution reactions. Values obtained at different levels of theory used by energy evaluation are shown in the columns’ headers (s/e: semi-empirical xTB-GFN2, DFT: density functional theory (B3LYP/6-311++G**/D3 or B97-3c), CC: coupled cluster method CCSD(T) with DPLNO computational scheme).
Table 1. Electronic ( Δ E ) and Gibbs free ( Δ G ) energies (in kcal/mol) of the reaction (R) and activation (A) in the nucleophilic substitution reactions. Values obtained at different levels of theory used by energy evaluation are shown in the columns’ headers (s/e: semi-empirical xTB-GFN2, DFT: density functional theory (B3LYP/6-311++G**/D3 or B97-3c), CC: coupled cluster method CCSD(T) with DPLNO computational scheme).
Calculation Methods Δ ER (s/e) Δ EA (s/e) Δ ER (DFT) Δ EA (DFT) Δ ER (CC) Δ EA (CC) Δ GR (CC) Δ GA (CC)
XTB-GFN2 (MEP 1, OPT 2, FREQ 3)−2.631.23.313.41.920.23.617.7
XTB-GFN2 (MEP, FREQ) + B3LYP 4 (OPT)−1.032.01.510.50.216.60.614.0
B3LYP (MEP, OPT) + XTB-GFN2(FREQ)−1.032.51.510.60.216.70.614.0
B97-3c--−1.29.10.216.40.114.7
Solvent (water, C-PCM) w/o enzime
B3LYP--10.733.69.038.46.835.3
B3LYP, each reactant/product separately--−3.413.8−3.421.5−2.635.6
B97-3c--7.924.88.634.48.031.6
B97-3c, each reactant/product separately----−8.313.0−7.526.5
Nucleophilic substitution reaction involving hydroorthovanadat-ion
B97-3c--−6.17.3−3.814.4−3.712.6
1 MEP—calculations of the minimum energy path; 2 OPT—geometry optimization of reactant, product and transition state; 3 FREQ—calculations of vibrational frequencies; 4 Calculations using the B3LYP density functional were carried out with the 6-311++G** basis set using the D3(BJ) correction for the dispersion interaction energy.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Lashkov, A.A.; Eistrich-Geller, P.A.; Samygina, V.R.; Rubinsky, S.V. QM/MM Study of a Nucleophilic Substitution Reaction Catalyzed by Uridine Phosphorylase from Vibrio cholerae. Crystals 2023, 13, 803. https://doi.org/10.3390/cryst13050803

AMA Style

Lashkov AA, Eistrich-Geller PA, Samygina VR, Rubinsky SV. QM/MM Study of a Nucleophilic Substitution Reaction Catalyzed by Uridine Phosphorylase from Vibrio cholerae. Crystals. 2023; 13(5):803. https://doi.org/10.3390/cryst13050803

Chicago/Turabian Style

Lashkov, Alexander A., Polina A. Eistrich-Geller, Valeriya R. Samygina, and Sergey V. Rubinsky. 2023. "QM/MM Study of a Nucleophilic Substitution Reaction Catalyzed by Uridine Phosphorylase from Vibrio cholerae" Crystals 13, no. 5: 803. https://doi.org/10.3390/cryst13050803

APA Style

Lashkov, A. A., Eistrich-Geller, P. A., Samygina, V. R., & Rubinsky, S. V. (2023). QM/MM Study of a Nucleophilic Substitution Reaction Catalyzed by Uridine Phosphorylase from Vibrio cholerae. Crystals, 13(5), 803. https://doi.org/10.3390/cryst13050803

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop