A Theoretical Approach to Engineering a New Enzyme

Density function theory, a subfield of quantum mechanics (QM), in combination with molecular mechanics (MM) has opened the way to engineer new artificial enzymes. Herein, we report theoretical calculations done using QM/MM to examine whether the regioselectivity and rate of chlorination of the enzyme chloroperoxidase can be improved by replacing the vanadium of this enzyme with niobium through dialysis. Our calculations show that a niobium substituted chloroperoxidase will be able to enter the initial steps of the catalytic cycle for chlorination. Although the protonation state of the niobium substituted enzyme is calculated to be different from than that of the natural vanadium substituted enzyme, our calculations show that the catalytic cycle can still proceed forward. Using natural bond orbitals, we analyse the electronic differences between the niobium substituted enzyme and the natural enzyme. We conclude by briefly examining how good of a model QM/MM provides for understanding the mechanism of catalysis of chloroperoxidase.


Introduction
Density functional theory (DFT), a subfield of quantum mechanics (QM), is successful in calculating the ground state energy of large molecules, but it cannot be applied to macromolecules like enzymes which have thousands of atoms [1]. On the other hand, molecular mechanics (MM), which applies Newtonian mechanics to molecules, has been applied to model the conformational states of enzymes; however, it does poorly in modeling chemical reactions in which covalent bonds are broken or formed. For modeling chemical reactions, DFT is needed. The combination of the DFT and MM (i.e., "QM/MM") has been applied to study a number of enzymes [2]. In QM/MM studies, DFT is used for modeling some portion of the active site, and MM is used to model the remainder of the enzyme.
Herein we utilize QM/MM to study vanadium chloroperoxidase (VCPO), which catalyzes chlorination of electrophilic substrates and compare our results with the literature [3][4]. QM/MM has also been applied to determine the 51 V chemical shifts of VCPO as well as the peroxo form of the enzyme [5]. DFT studies [6] have been carried on just the vanadium cofactor and methyl amine and guanidine to model the active site; however, QM/MM studies showed that long range electrostatic forces have very important effects on the protonation state of the vanadium co-factor [3].
VCPO catalyzes oxidation of chloride in the presence of cellular hydrogen peroxide: In this process, the vanadium co-factor remains in the same oxidation state (i.e., +5) throughout the catalytic cycle, which is outlined in figure 1. An X-Ray crystal structure at a resolution of 2.15 reveals the active site as depicted in figure 2 [7]. This resolution is such that one cannot ascertain the protonation state of the VO4 co-factor. The VO4 cofactor has a trigonal bipyramidal structure and is held in the enzyme by hydrogen bonds from Lys353, Ser402, His404, Arg360 and Arg496 and possibly by a weak bond to His496. One can easily remove the vanadium cofactor by dialysis and reincorporate it [8]. There is no reason in principle why some other metal could not be substituted in the reincorporation step.
Our long range goal is to determine whether we can improve the rate and regioselectivity of VCPO by altering the transition metal in VCPO. We thus began this study by determining whether replacement of vanadium with niobium would even lead to a viable enzyme that could enter the catalytic cycle of VCPO.

Theory
The QM/MM approach to modelling of enzymes involves dividing the enzyme into a quantum mechanical part and molecular mechanical part. Thus, the overall energy of the enzyme is calculated as follows [9]: To find the energy, one assumes that "fictitious" orbitals, i , can be assigned to represent each electron, 1, 2, 3, . . . i, . . . N where N is the total number of electrons in the system. The first two terms are calculated self-consistently in a way parallel to Hartree-Fock theory.
The remaining exchange-correlation energy functional is approximated by the adiabatic connection method. For this study, we used the hybrid B3LYP functional which involves parameterization as follows: The three parameters (a,b,c) were selected to fit a set of reference data (ie, the G2 data base) consisting of atomization energies, ionization energies, proton affinities and some total energies. DFT with B3LYP gives the energies within 1 kcal/mol for a surprisingly wide variety of systems [1]. Each of the terms of equation 4 is a complicated expression. For example, the first term is the following: The probability density is now represented by rs instead of and a different set of empirical constants (A, x0, b and c) are needed.
The third term of equation (1), / QM MM E , includes terms from van der Waals contributions which are calculated by MM. Additionally, one must treat the covalent bonds which are cut by the boundary between the QM and MM regions. We used the frozen orbital approach of Friesner and co-workers [2] in which the coefficients of the basis functions for the orbitals on the boundaries are set to zero; thus, the orbitals on the boundaries do not enter the self-consistent field calculation. Instead, molecular mechanics type parameters are introduced that provide an accurate energy function.

Methods
QM/MM calculations were done using the software, QSite [11]. The X-ray crystal structure (1VNI) of VCPO was downloaded from the Protein Data Bank (www.rcsb.org). QSite provides a preprocessing module that selectively removes waters from outside of the active site, i.e., which are more than 5 from vanadium. Because Impact, the software component of QSite which does MM, cannot process metal bonds, the preprocessing module makes the bond order zero between the vanadium and oxygens and between V and N of His496. Manually, we adjusted the charge on vanadium to +5 and the charge on the oxygens to -1 or -2 depending on whether the oxygen is protonated. This module also optimizes the hydrogen bond network and performs a restrained minimization with OPLS-AA force field to RMSD<0.3. For the QM calculations, we used DFT with B3LYP and a basis set of LACVP using electron core potentials for the metal and 6-31G** for the rest of the atoms in the QM region.. The QM region consisted of the side chains of Lys353, His496 with a VO4 or NbO4 cofactor, Asp292, and the full residues of Arg360, Arg490, His404, Gly403, and Ser402. Arg360 is adjacent to a Pro361 residue, and QSite does not allow backbone cuts in the proline ring. Therefore following the method of the literature [3], we converted Pro361 to alanine and included this full alanine residue in the QM part of the calculation. MM constraints consisted of the protein backbone and any atoms more than 20 from the metal. QM constraints were the amino acid backbone of Arg360, Arg490, His404, Gly403, Ser402 and His496 ring as well as the hydrogens of Lys353 facing away from the cofactor.

Results and Discussion
In order to assess the viability of a niobium substituted chloroperoxidase enzyme entering the catalytic cycle, we must first ask a large number of questions such as (a) what is the resting state of the niobium substituted peroxidase? (b) what is the lowest energy protonation state of each of the proposed intermediates (c) are the proposed steps still favorable upon niobium substitution? After doing QSite calculations, Kravitz et al. [3] concluded that the resting state of VCPO is an equilibrium between 1-V and 2-V based on the energy of the optimized structures (see figure 2). In our calculations, we find that 2-V is lower in energy than 1-V by 7 kcal/mole. The optimized structure, 2-V, in comparison to 1-V and 3-V has bond lengths which are close to the bond lengths of X-ray crystal structure. In particular, for 2-V, the V-O bond lengths for equatorial oxygens are all 1.69 , and the X-Ray data reveals 1.60 for V-O10 and V-O11 and 1.61 for V-O9. One possible explanation for this difference in results from those Kravitz et al. is the following. At the time of their research in 2005, QSite could only perform calculations on 8000 atoms, and VCPO is over 9000 atoms. Therefore, they did their calculations on a truncated version of the enzyme after justifying their truncation by examining the electrostatics of the enzyme. But even Kravitz et al. [3] noticed that long range electrostatic forces play a significant role in the relative energies of various protonation states of VCPO, and hence it seems possible that their analysis of electrostatic effects was incomplete, and their truncation affected the energies that they calculated. The other QM/MM study of VCPO neglected to even do calculations on 2-V [4].
On substitution with niobium, we find a significant shift in energies of the various protonation states of the cofactor. Now, structure 3-Nb becomes lower in energy than either 1-V or 2-V. Since 2-V is lower in energy than 3-V, it is puzzling why 3-Nb should be lower in energy than 2-Nb.
We undertook a study of the natural bond orbitals (NBO) in the optimized structures. Figure 3 summarizes the results of this study regarding the bond orders between the metal and oxygens. Namely, the co-factor of VCPO, (2-V), is actually VO3. Furthermore, these results suggest the cofactor is not actually bonded to His496. Raugei et al. [4] in their QM/MM study of VCPO obtained the similar results on 1-V and 3-V by calculating Mayer's bond orders. Since the distance between V-N (2.15 ) is less than the sum of the van der Waals radii of V (2.05 ) and N (1.55 ), it is possible that there are secondary bonding interactions [12].

Conclusions
We have calculated the energies of various protonation states of VCPO and found 2-V to be the most stable. If the vanadium of VCPO were replaced with niobium, we calculated the resting state of the enzyme will be a different protonation state (i.e. 3-Nb) than in the natural enzyme. It is still thermodynamically very favourable for 3-Nb to enter a catalytic cycle analogous to the one VCPO uses to chlorinate substrate. In future work, we will examine whether the niobium substituted enzyme will improve the rate of chlorination. This work will certainly have application in industrial processes wherein chlorination plays an important role. In this conclusion, we would also like to take a step back and consider the kind of modeling that the density functional theory affords. For example, Cramer [13] in the introduction of his textbook on computational chemistry promises to present no equation without an effort to "provide an intuitive explanation for its form and the various terms within it"; however, in presenting the equations of DFT such as (5), he has to admit an "utter violation" of this promise since no physical meaning can be assigned to any term in (5). Because of this problem in DFT modeling, we should continue the search for other complementary ways of modeling matter [14].