Calculation of higher protonation states and of a new resting state for vanadium chloroperoxidase using QM/MM, with an Atom-in-Molecules analysis.

Earlier QM/MM studies of the resting state of vanadium chloroperoxidase (VCPO) focused on the diprotonated states of the vanadate cofactor. Herein, we report a new extensive QM/MM study that includes the tri- and quadprotonated states of VCPO at neutral pH. We identify certain di- and triprotonated states as being candidates for the resting state based on a comparison of relative energies. The quadprotonated states as well as some of the triprotonated states are ruled out as the resting state. An Atoms-in-Molecules (AIM) analysis of the complex hydrogen bonding around the vanadate cofactor helps to explain the relative energies of the protonation states considered herein, and it also indicates new hydrogen bonding which has not been recognized previously. A Natural Bond Orbital (NBO) study is presented to give a better understanding of the electronic structure of the vanadate co-factor.

ABSTRACT. Earlier QM/MM studies of the resting state of vanadium chloroperoxidase (VCPO) focused on the diprotonated states of the vanadate cofactor. Herein, we report a new extensive QM/MM study that includes the tri-and quadprotonated states of VCPO at neutral pH.
We identify certain di-and triprotonated states as being candidates for the resting state based on a comparison of relative energies. The quadprotonated states as well as some of the triprotonated states are ruled out as the resting state. An Atoms-in-Molecules (AIM) analysis of the complex hydrogen bonding around the vanadate cofactor helps to explain the relative energies of the protonation states considered herein, and it also indicates new hydrogen bonding which has not been recognized previously. A Natural Bond Orbital (NBO) study is presented to give a better understanding of the electronic structure of the vanadate co-factor.

Introduction:
Vanadium haloperoxidases are found in many algae, lichens, fungi and bacteria, and they also have practical applications in the biotechnology industry. [1,2] In Nature, they have a crucial role in biosynthetic pathways that lead to antibiotics as well as the emission of greenhouse gases (i.e., chloroform and bromoform) into the atmosphere. [3] Within this class of enzymes, we have focused on vanadium chloroperoxidase (VCPO), which catalyzes the two-electron oxidation of halide to hypohalous acid (Scheme 1); the hypohalous acid diffuses out of the active site and reacts in a nonspecific way with electrophilic substrates. In the fungus, Curvularia inequalis, the synthesized hypohalous acid breaks down the cell membrane of plants which the fungus is attacking. [3] In order to begin the study of the catalytic mechanism [4] of this enzyme, one first has to decide what is its resting state. In many cases, the investigation of the mechanism of catalysis of enzymes requires a determination of the resting state. [5] Scheme 1. VCPO catalyzes the two-electron oxidation of halide.
In the case of VCPO, the results from X-ray diffraction studies give only partial information about the resting state. Nevertheless, it is clear that the active site of the enzyme consists of a vanadate in a trigonal bipyramidal configuration with His496 covalently bonded to one of the axial positions of the vanadate, and the vanadate is hydrogen bonded to a sphere of amino acid residues surrounding it, namely, Lys353, Ser402, Gly403, His404, Arg490, and Arg360 ( Figure   1). [6] Furthermore, Asp 292 forms a salt bridge to Arg490. The resolution of the X-ray crystal structure ( 0.24 Å) is insufficient to reveal the protonation state of the vanadate, and computational studies have been carried out to determine the resting state of VCPO. [7][8][9][10][11] .
Spectroscopic and biochemical studies have also established that the oxidation state of vanadium remains +5 throughout the catalytic cycle [12][13][14]; therefore, it has been thought that the protonation state of the oxo groups coordinating vanadium may be responsible for redox tuning of the cofactor. [15] Recently, Gupta et al. performed 51 V magic angle spinning (MAS) NMR spectroscopy and determined certain NMR parameters (i.e., the reduced anisotropy, the asymmetry parameter, and the electric field gradient (EFG) tensor parameters) which are influenced by the protonation state of vanadate in the resting state. [15] Then, they constructed models of the active site with the vanadate in various protonation states, and calculated these NMR parameters by density functional theory using the B3LYP hybrid functional and 6-311G(d,p) basis set while constraining the position of all atoms except for the hydrogens and vanadate. They argued that the protonation state is dependent on the pH in the following way: diprotonated 3 (pH 8.3), triprotonated 10 or 13 (pH 7.3) and quadprotonated 16 (pH 6.3) (Figure 2). At pH 7.3, the observed NMR parameters were consistent with either triprotonated 10 or 13.  [7] Another advantage of QM/MM modeling of VCPO is that many heavy atoms do not need to be constrained in order to obtain convergence as they do in a QM model.
We describe the results of a QM/MM study starting with input structures, 1-23, in the mono, di-, tri-and quadprotonated states at pH 7. This QM/MM study is conducted with two QM regions consisting of 8 residues and 13 residues. Thus far, no QM/MM study of VCPO with 13 residues in the QM region has been reported. This study provides structural information about VCPO as well as the stability of the tri and quadprotonated states. It also provides the relative energies of the protonation states. All of this information combined together allows us to indicate what is likely to be the resting state of VCPO.
Next, we apply Bader's Atoms-in-Molecules (AIM) approach [16] to elucidate the complex hydrogen bonding to the vanadate cofactor. The AIM approach has been increasingly used to determine noncovalent interactions. [17][18][19][20][21][22] In the literature on VCPO, almost all knowledge of the hydrogen bonding is based on distance considerations, but distance considerations can be misleading. The AIM approach involves a topological analysis of the probability density function. Atoms that are interacting (bonding or secondary interactions) have a single line of locally maximum electron density, and this line is referred to as a bond path. A point on the bond path between two interacting atoms with the lowest value of the electron density is called the bond critical point, and we report herein the values of the electron density, ρ, at the bond critical points in units of electrons/bohr 3 . The bond critical point along a covalent bond, will typically have ρ= 0.2 electrons/bohr 3 , but for a hydrogen bond the ρ value at the bond critical point will be less than 0.1 electrons/bohr 3 . Knowing the precise hydrogen bonding in VCPO will aid the study of the mechanism, which has so far been limited to the study of the cofactor in isolation or at most by taking into consideration the hydrogen bonding between Lys353 and vanadate. [4,23] In addition to this AIM study, we report a natural bond orbital (NBO) study of the tri-and quadprotonated states of vanadate to gain a better understanding of their electronic structure.   . QM Region I with atom numbering. Hydrogen bonding which is present in most of the optimized structures is indicated by dotted lines. The hydrogen bonding between H6 or H7 and the crystallographic water molecules varies between structures and is not indicated in this Figure; this hydrogen bonding is discussed in the text and can be found from the Tables S3-S10 in the Supplementary Material.

Methods:
QM/MM calculations were carried out with QSite. [24,25] The topological analysis of the electron density, involved in the AIM approach, is built into the QSite code and is invoked with iplotnoncov=1. Natural Bond Orbital calculations were done with NBO 6.0 which is implemented in Jaguar 10.1. [25,26] The X-Ray crystal structure, 1IDQ, was preprocessed with Schrodinger Protein Preparation Wizard.
[27,28] The Schrodinger Protein Preparation Wizard may flip the orientation of Gln, His and Asn residues in the enzyme structure to optimize the hydrogen-bonding network because the orientation of these residues is unclear from the X-ray crystal structure. The titratable residues (i.e., histidines, aspartates, glutamates) outside of the QM region were protonated assuming pH 7.0 because the PROPKA method [29] for calculating the pKa values of amino acid residues failed for this structure. [30] The Protein Preparation Wizard detected that the following residues were incomplete in the X-ray crystal structure: Glu83, Glu119, Gln120, Pro121, Asn122, Pro123, Asn124, Pro125, Asn128, His185. Although they were on the surface of the enzyme and unlikely to affect the calculations, they were repaired. Furthermore, the hydroxyl hydrogen (H28) of Ser402 was rotated to point to the equatorial O3 of the vanadate. (See Figure 3 for atom numbering.) The resulting structure was then relaxed through a restrained minimization to RMSD 0.3 Å on hydrogens to relieve any steric clashes. (Restrained minimization on the heavy atoms was not done because it was found to result in significant and irreproducible changes in the positions of atoms.) This structure was then manually altered to generate input structures, 1-23, which are grouped into four series: mono-(1-2), di-(3-8), tri-(9-14) and quad- (15)(16)(17)(18)(19)(20)(21)(22)(23) protonated states. (See Figure 2.) Two QM regions were used. The first QM region, which is shown in Figure 3, consisted of His496 (side chain), Lys353 (side chain), Arg360 (full), Arg490 (full), Ser402 (full), Gly403 (full), HIE404 (full) [31], Asp292 (side chain), H2O1911, H2O1778 and vanadate (VO4) in various protonation states. We will refer to this first QM region as "QM Region I" (Figure 3).
The second QM region, "QM Region II" (Figure 4), consisted of these 8 residues and vanadate as well as all of the atoms in the following residues: Ala399-Pro398-Phe397-Pro396-Pro395. It is also included two additional crystallographic water molecules: H2O1832, H2O1628. The series of three water molecules-namely, H2O1911, H2O1778, and H2O1832-form a hydrogen bonded bridge between the axial oxygen (O4) of vanadate and Asp292 (i.e., O13). QM Region I consisted of 149-153 atoms depending on the protonation state of vanadate. QM Region II consisted of 226-230 atoms.
QM optimizations in this QM/MM study were carried out with density functional theory with the B3LYP functional and the lacvp* basis set. [32] In the case of calculations concerning 3-8, the following atoms were constrained for both QM Region I and QM Region II: the backbone atoms of Ser402, Gly403, His404, the methylene atoms of Arg360 and Arg 490, H50 and H51 of the guanidine of Arg360, as well as H52 of HIE404. For the other input structures (1-2, 9-23), all of these atoms as well as the side chain atoms of Asp292 were constrained. These constraints were used to speed up convergence of the geometry optimization and follow the literature. [7] As stated above, the entire Arg360 residue, which is next to Pro361, was included in the QM region, but the only parameterized cut for proline in QSite is along C-Cα; hence, the entire Pro361 would have had to be included. As a result, Pro361 was simplified to alanine following the example of Kravitz et al., and this alanine was constrained during the geometry optimizations. [7] In the case of QM Region II, all of the atoms in the sequence Ala399-Pro395 were constrained except for the phenyl atoms (i.e., CH2-Ph) of Phe397. The MM region consisted of the remainder of the enzyme; the protein backbone was frozen during the MM optimizations. Furthermore, any MM atoms beyond 20 Å from vanadium were frozen. MM optimization was done with the OPLS_2005 force field. [33] Solvent effects were modeled by neutralizing charged surface residues. [34] The solvent accessible surface area (SASA) of each residue was determined by the method of Lee and Richards [35] which is implemented in Bioluminate, a software of Schrodinger Corporation, and Bioluminate was used to determine the SASA of each residue. [36] Any charged surface residue with a SASA of greater than 1 Å 2 was neutralized.

Justification for the protonation state of amino acids in QM region
The justification for choosing the protonation states of the various amino acid residues in the QM region of the input structures for the QM/MM calculations was the following. In a study of the kinetics of VCPO, van Schijndel et al. observed that at pH<4 H2O2 will not bind to VCPO and chloride inhibits VCPO. From this observation, they argued that a histidine within the active site must be monoprotonated at pH>4. [37,38] Zhang and Gascon argued from this observation that His404 is becoming diprotonated at low pH, and diprotonation of His404 leads to the enzyme being inhibited by halide at low pH. [8] As a result, in their QM/MM study of VCPO, Zhang and Gascon kept His404 in the monoprotonated state, HIE, as shown in Figure 3   In all of the optimized quadprotonated structures (15-O to 23-O), one finds that a proton has transferred to His404 (i.e., to N26). One can argue that these structures may not enter the catalytic cycle since, as discussed earlier, diprotonation of His404 results in inhibition of the enzyme by substrate. [8,37,38] In the other cases (15, 16, 18, 21, and 22), an axial proton (H6 or H7) transfers from axial oxygen, O4, to N26 (e.g., 16-O, Figure 6). Whenever the equatorial oxygen, O3, is protonated (with H10) as in the cases of 17, 19, 20, and 23, the optimized structures show that H10 has transferred to N26 (eg., 17-O, Figure 7).

QM/MM Optimizations with QM Region II
In order to confirm that the proton transfers are not an artifact of the size of the QM region being used, we carried out a set of QM/MM calculations on input structures 1-23 with the much larger 13 residue, QM Region II ( Figure 4). QM Region II included three additional residues (i.e., Pro 395, Phe397, Ala399) modeled in the QM study of Gupta et al. [15] Two more proline residues (Pro396, Pro398) had to be included because QSite allows parameterized cuts of proline only along the C-Cα bond. The resulting optimized structures have many similarities with the optimized structures found using the 8 residue QM Region I. Table S2 in

Bond lengths and bond angles of optimized structures
The calculated lengths of the V-O bonds and the V-N bonds in the optimized structures with QM Region I are reported in Table 1. (The calculated bond lengths for the optimized structures found with QM Region II are given in Table S1 of the Supplementary Material; these bond lengths do not significantly differ from the results given in Table 1  It is not possible to compare the calculated bond lengths (Table 1) with the X-ray crystallographic data found in the PDB because these data have a large error of 0.24 Å. [6] However, the Extended X-ray Absorption Fine Structure (EXAFS) regions obtained at pH 6.0 at 77 K from vanadium K-edge X-ray Absorption spectroscopy led to more precise bond lengths ( 0.02 Å). Namely, two V-O bonds are 1.69 Å, one V-O bond is 1.54 Å, and one V-O bond is 1.93 Å. [13] Based on the calculated bond lengths, one can conclude that this short bond length of 1.54 Å must correspond to an unprotonated O4. Furthermore, the longer bond length of 1.93 Å corresponds to either mono-or diprotonated O4. The bond lengths of 1.69 Å in this spectroscopic data correspond to an unprotonated or monoprotonated oxygen. Table 1. Calculated bond lengths (in Å) and bond angles (O4-V-N5) for optimized structures (QM Region I) where 1-23 are the input structure as well as bond lengths and bond angle for the X-ray crystal structure data (1IDQ) of VCPO.

NBO study of Vanadate
To understand the electronic structure of the vanadate protonation states, we undertook a natural bond orbital (NBO) study of an isolated triprotonated 24 and quadprotonated 25 ( Figure 9). First, the NBOs of vanadate for 24 and 25 are as follows: Although there is no NBO between V and O3 for 24, there is electron density between V and O3 from delocalization as found in the second order perturbative mixing. Namely there is mixing between the lone pair of O3 with the antibonding orbital between V-O1 and with an extremely large E2 stabilization energy of 171 kcal/mol (Table 2) Figure 10 in which the overlap is illustrated for the case of V-O1 bond orbital as the donor and V-O2 antibonding orbital as the acceptor. It is also interesting that although there is no NBO between axial O4 and V, there is significant delocalization shown in the perturbative mixing. Figure 11 shows this mixing between the donor orbital (i.e., the lone pair on O4) and the acceptor orbital (i.e., the antibonding orbital V-N5) which has significant density in the area of V-O4 bond with stabilization energy (E2) of 60.70 kcal/mol. Table 2. NBO result for 24: V-O bond orbitals involved in perturbative mixing between donor and acceptor bond orbital with stabilization energy, E2. There is also significant mixing of the lone pair on O3 with an antibonding orbital between V1-O1.

Atoms-in-Molecules Study
The Atoms-in-Molecules (AIM) approach was applied to clearly determine the complex Some other structurally significant hydrogen bonds become evident from the AIM study.
Namely, the carbonyl of Pro401 is strongly hydrogen bonded (ρ=0.04178) to one of the amino hydrogens of Lys353 (i.e., O18-H32). We expect that this hydrogen bond will stabilize the conformation of Lys353 so that Lys353 is strongly hydrogen bonded to the equatorial oxygen (O2) of vanadate; a QM study of the mechanism has shown that this hydrogen bond between Lys353 and vanadate plays a significant role. [4,23] Using QM Region II, we find that the other amino hydrogen of Lys353 is strong hydrogen bonded to the carbonyl of Pro398 (i.e., ρ(O58-H33)=0.033 on average). There are also two significant dihydrogen bonds involving Arg490, namely between H43 and H44 and between H35 and H45; these dihydrogen bonds should restrict the conformation of Arg490. Furthermore, with QM Region II, Pro396 has two hydrogen bonds to Arg360 (ρ(O60-H50)=0.013 and ρ(O60-H66)=0.026 on average); thus, the conformation of Arg360 is also constrained.

Relative Energies of the Optimized Structures with QM Region I
The relative energies of the QM/MM optimized structures (with QM Region I) in each of the series of mono-, di-, tri-and quadprotonated states of VCPO are shown in Table 4. In the diprotonated series (3)(4)(5)(6)(7)(8), optimized structure 5-O ( Figure 14) with protonation (H6) of axial O4 and protonation (H10) of equatorial O3 is the lowest in energy. An earlier QM/MM study [7] obtained different results for these relative energies, but it should be noted that this earlier study did not include two crystallographic water molecules in the QM region. Another earlier study used a much smaller QM region and obtained different relative energies. [8] The present results show that these water molecules are involved in strong hydrogen bonds to hydrogens protonating axial O4. In 5-O, H10 is strongly hydrogen bonded to His404 (i.e., ρ(H10-N26)= 0.03733) and H6 is hydrogen bonded to a water molecule with ρ=0.06714 electrons/bohr 3 . In comparing      kcal/mol although O2 is unprotonated and O1 is diprotonated. In the case of 18-O, there are also steric clashes between H8 and H23 and a hydrogen bond forms between H8 and N40 (ρ(H8-N40)=0.03398) of Arg490. This steric clash becomes evident with space filling models ( Figure  16A). With any model, one can notice how H23 has been pushed out of the guanidine plane ( Figure 16B).

Relative Energies with QM Region II.
When the much larger QM Region II was used in the geometry optimization of 1-23, we found that the lowest and highest energy structures remain the same in the di-and quadprotonated series. Table 5 compares the relative energies found with QM Region I and II with input structures, 1-8 and 15-23. One significant difference is that with the larger QM Region II we find

What is the resting state?
In the above discussion, it has been argued that the quadprotonated forms of VCPO are unlikely  Across the various series, we find that the state in which the rear O3 is protonated is often the lowest energy structure because of strong hydrogen bonding to His404 and because of minimal steric interactions which are present in forms where O1 or O2 is protonated. When O1 or O2 is protonated, steric interactions with Lys353 and Arg360/490 destabilize these protonation states based on an examination of space filling models.
We have also reported the results of an AIM study to ascertain the complex hydrogen   Table S3. ρ values in units of e/bohr 3 for monoprotonated (1-2) and diprotonated (3)(4)(5)(6)(7)(8) forms [43] calculated using QM Region I. Atom numbering is given in Figure 3 n/a n/a 0.00000 0.03891 n/a 0.04854 n/a n/a O16-H8 n/a n/a 0.00000 n/a n/a n/a 0.00978 0.02776 O15-H6 0.05746 n/a n/a 0.06006 0.06714 n/a n/a n/a O18-H32 0 n/a n/a 0.00000 0.00000 n/a n/a n/a N26-H10 n/a 0.04017 n/a n/a 0.03733 n/a 0.05815 n/a H45-H35 0.01255 0.01417 0.01553 0.01524 0.01434 0.01249 0.01578 0.01276 O17-H34 0.0217 0.02023 0.02077 0.01934 0.02010 0.01849 0.01670 0.01988 N5-H22 0 0 0.00000 0.00000 0.00000 0.00000 0.01274 0.01042 O15-H7 n/a n/a 0.08131 0.00000 n/a 0.00000 n/a n/a H21-O2 0 0 0.00000 0.00000 0.00000 0.02063 0.00000 0.00000 a n/a means "not applicable": one or both of the atoms were not present. "NC" means ρ value was not calculated, and it is expected to be greater than 0.1. A value of ρ=0 indicates that there is no density between these atoms, and hence there is no bond between these atoms. b H22 transfers from N19 to O1; H22-O2 distance is 1.000 Å and ρ(N19-H22)=0.04895 c H22 transfers from N19 to O1; H22-O2 distance is 1.002 Å and ρ(N19-H22)=0.04922 Table S4. ρ values in units of e/bohr 3 for triprotonated forms, 9-14 a calculated using QM Region I 0.00000 0.00000 0.00000 0.00095 0.00000 0.01157 a n/a means "not applicable": one or both of the atoms were not present. NC means ρ value was "not calculated" and is expected to be greater than 0. O4-H10 0.02554 n/a 0.00000 0.00000 a n/a means "not applicable": one or both of the atoms were not present. NC means ρ value was not calculated and is expected to be greater than 0.1. A value of ρ=0 indicates that there is no density between these atoms, and hence there is no bond between these atoms. O16-H8 n/a n/a 0.00000 n/a n/a n/a 0.05718 0.04715 O15-H6 0 n/a n/a 0.04845 0 n/a n/a n/a H43-H44 N26-H6 0.03228 n/a n/a 0.00000 0.00000 n/a n/a n/a N26-H7 n/a n/a 0.04203 0.04982 n/a 0.06149 n/a n/a N26-H10 n/a 0.03610 n/a n/a 0. 0.00000 n/a n/a 0.00000 0.06497 N/A n/a n/a H75-O16 0 0 0 0 0.01264 0 0 0 O15-H11 n/a n/a n/a n/a n/a 0.02525 n/a n/a O4-O16

H9-N40
n/a n/a n/a n/a O4-H73 There are no conflicts of interest to declare.