Substrate hydroxylation in methane monooxygenase: quantitative modeling via mixed quantum mechanics/molecular mechanics techniques

Using broken-symmetry unrestricted density functional theory quantum mechanical (QM) methods in concert with mixed quantum mechanics/molecular mechanics (QM/MM) methods, the hydroxylation of methane and substituted methanes by intermediate Q in methane monooxygenase hydroxylase (MMOH) has been quantitatively modeled. This protocol allows the protein environment to be included throughout the calculations and its effects (electrostatic, van der Waals, strain) upon the reaction to be accurately evaluated. With the current results, recent kinetic data for CH3X (X = H, CH3, OH, CN, NO2) substrate hydroxylation in MMOH (Ambundo, E. A.; Friesner, R. A.; Lippard, S. J. J. Am. Chem. Soc. 2002, 124, 8770-8771) can be rationalized. Results for methane, which provide a quantitative test of the protocol, including a substantial kinetic isotope effect (KIE), are in reasonable agreement with experiment. Specific features of the interaction of each of the substrates with MMO are illuminated by the QM/MM modeling, and the resulting effects upon substrate binding are quantitatively incorporated into the calculations. The results as a whole point to the success of the QM/MM methodology and enhance our understanding of MMOH catalytic chemistry. We also identify systematic errors in the evaluation of the free energy of binding of the Michaelis complexes of the substrates, which most likely arise from inadequate sampling and/or the use of harmonic approximations to evaluate the entropy of the complex. More sophisticated sampling methods will be required to achieve greater accuracy in this aspect of the calculation.     

Electronic and spectroscopic studies of the non-heme reduced binuclear iron sites of two ribonucleotide reductase variants: comparison to reduced methane monooxygenase and contributions to O2 reactivity

Circular dichroism (CD), magnetic circular dichroism (MCD), and variable-temperature variable-field (VTVH) MCD have been used to probe the biferrous active site of two variants of ribonucleotide reductase. The aspartate to glutamate substitution (R2-D84E) at the binuclear iron site modifies the endogenous ligand set of ribonucleotide reductase to match that of the binuclear center in the hydroxylase component of methane monooxygenase (MMOH). The crystal structure of chemically reduced R2-D84E suggests that the active-site structure parallels that of MMOH. However, CD, MCD, and VTVH MCD data combined with spin-Hamiltonian analysis of reduced R2-D84E indicate a different coordination environment relative to reduced MMOH, with no mu-(1,1)(eta(1),eta(2)) carboxylate bridge. To further understand the variations in geometry of the active site, which lead to differences in reactivity, density functional theory (DFT) calculations have been carried out to identify active-site structures for R2-wt and R2-D84E consistent with these spectroscopic data. The effects of varying the ligand set, positions of bound and free waters, and additional protein constraints on the geometry and energy of the binuclear site of both R2-wt and variant R2s are also explored to identify the contributions to their structural differences and their relation to reduced MMOH.     

Energetics of oxidized and reduced methane monooxygenase active site clusters in the protein environment

Using the density functional optimized active site geometries obtained in the accompanying paper (Lovell, T.; Li, J.; Noodleman, L. Inorg. Chem. 2001, 40, 5251), a combined density functional and electrostatics approach has been applied to further address attendant uncertainties in the protonation states of the bridging ligands for MMOH(ox). The acidities (pK(a)s) associated with the bridging H(2)O ligand in Methylococcus capsulatus and corresponding energetics of each active site cluster interacting with the protein environment have been evaluated. The pK(a) calculations in combination with the results of the gas phase DFT studies allow the active site cluster in Methylosinustrichosporium to be best described as a diiron unit bridged by 2OH(-) ligands having an overall neutral net cluster charge. The presence of the exogenous acetate in M. capsulatus reveals a diiron unit bridged by 1OH(-) and 1H2O which asymmetrically shares its proton with a second-shell acetate in a very short strong AcO..H...OH hydrogen bond. For all MMOH(ox) and MMOH(red) active sites examined, significant Fe-ligand covalency is evident from the ESP atom charges, consistent with very strong ligand --> metal charge transfer from the muOH(-) and mu-carboxylato bridging ligands. The magnitude of electrostatic interaction of the individual protein residues in the active domain with the active site has been assessed via an energy decomposition scheme. Important second-shell residues are highlighted for the next level of quantum mechanics based calculations or alternatively for site-directed mutagenesis studies. Finally, from the known structural and spectroscopic evidence and the DFT studies, a possible mechanism is suggested for the conversion of MMOH(ox) into MMOH(red) that involves a combination of protein residues and solvent-derived ligands from the second coordination sphere.     

X-ray absorption spectroscopic study of the reduced hydroxylases of methane monooxygenase and toluene/o-xylene monooxygenase: differences in active site structure and effects of the coupling proteins MMOB and ToMOD

The diiron active sites of the reduced hydroxylases from methane monooxygenase (MMOH(red)) and toluene/o-xylene monooxygenase (ToMOH(red)) have been investigated by X-ray absorption spectroscopy (XAS). Results of Fe K-edge and extended X-ray absorption fine structure analysis reveal subtle differences between the hydroxylases that may be correlated to access of the active site. XAS data were also recorded for each hydroxylase in the presence of its respective coupling protein. MMOB affects the outer-shell scattering contributions in the diiron site of MMOH(red), whereas ToMOD exerts its main effect on the first-shell ligation of ToMOH(red); it also causes a slight decrease in the Fe-Fe separation. These results provide an initial step toward delineating the differences in structure and reactivity in bacterial multicomponent monooxygenase proteins.     

Density functional studies of oxidized and reduced methane monooxygenase. Optimized geometries and exchange coupling of active site clusters

The conflicting protein crystallography data for the oxidized form (MMOH(ox)) of methane monooxygenase present a dilemma regarding the identity of the solvent-derived bridging ligands within the active site: do they comprise a diiron unit bridged by 1H2O and 1OH(-) as postulated for Methylococcus capsulatus or 2OH(-) ligands as suggested for Methylosinus trichosporium? Using models derived explicitly from the M. capsulatus and M. trichosporium protein data, spin-unrestricted density functional methods have been used to study two structurally characterized forms of the hydroxylase component of methane monooxygenase. The active site geometries of the oxidized (MMOH(ox)) and two-electron-reduced (MMOH(red)) states have been geometry optimized using several quantum cluster models which take into account the antiferromagnetic (AF) and ferromagnetic (F) coupling of electron spins. Trends in cluster geometries, energetics, and Heisenberg J values have been evaluated. For the majority of models, calculated geometries are in good agreement with the X-ray analyses and appear relatively insensitive to the F or AF alignment of electron spins on adjacent Fe sites. Discrepancies between calculation and experiment appear in the orientation of the coordinated His and Glu amino acid side chains for both MMOH(ox) and MMOH(red) and also in unexpected intramolecular proton transfer in the MMOH(ox) cluster models. There is additional dispersion between (and among) calculated and experimental Fe(3+)-OH(-) distances with relevance to the correct protonation state of the solvent-derived ligands. In an accompanying paper (Lovell, T.; Li, J.; Noodleman, L. Inorg. Chem. 2001, 40, 5267), a comparison of the related energetics of the active site models examined herein is further evaluated in the full protein and solvent environment.     

QM/MM study of the absorption spectra of DsRed.M1 chromophores

We report geometries and vertical excitation energies for the red and green chromophores of the DsRed.M1 protein in the gas phase and in the solvated protein environment. Geometries are optimized using density functional theory (DFT, B3LYP functional) for the isolated chromophores and combined quantum mechanical/molecular mechanical (QM/MM) methods for the protein (B3LYP/MM). Vertical excitation energies are computed using DFT/MRCI, OM2/MRCI, and TDDFT as QM methods. In the case of the red chromophore, there is a general blue shift in the excitation energies when going from the isolated chromophore to the protein, which is caused both by structural changes and by electrostatic interactions with the environment. For the lowest ππ* transition, these two factors contribute to a similar extent to the overall DFT/MRCI shift of 0.4 eV. An enlargement of the QM region to include active‐site residues does not change the DFT/MRCI excitation energies much. The DFT/MRCI results are closest to experiment for both chromophores. OM2/MRCI and TDDFT overestimate the first vertical excitation energy by 0.3–0.5 and 0.2–0.4 eV, respectively, relative to the experimental or DFT/MRCI values. The experimental gap of 0.35 eV between the lowest ππ* excitation energies of the red (cis‐acylimine) and green (trans‐peptide) forms is well reproduced by DFT/MRCI and TDDFT (0.32 and 0.37 eV, respectively). A histogram spectrum for an equal mixture of the two forms, generated by OM2/MRCI calculations on 450 snapshots along molecular dynamics trajectories, matches the experimental spectrum quite well, with a gap of 0.23 eV and an overall blue shift of about 0.3 eV. DFT/MRCI appears as an attractive choice for calculating excitation energies in fluorescent proteins, without the shortcomings of TDDFT and computationally more affordable than CASSCF‐based approaches. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Insights into the different dioxygen activation pathways of methane and toluene monooxygenase hydroxylases

The methane and toluene monooxygenase hydroxylases (MMOH and TMOH, respectively) have almost identical active sites, yet the physical and chemical properties of their oxygenated intermediates, designated P*, H(peroxo), Q, and Q* in MMOH and ToMOH(peroxo) in a subclass of TMOH, ToMOH, are substantially different. We review and compare the structural differences in the vicinity of the active sites of these enzymes and discuss which changes could give rise to the different behavior of H(peroxo) and Q. In particular, analysis of multiple crystal structures reveals that T213 in MMOH and the analogous T201 in TMOH, located in the immediate vicinity of the active site, have different rotatory configurations. We study the rotational energy profiles of these threonine residues with the use of molecular mechanics (MM) and quantum mechanics/molecular mechanics (QM/MM) computational methods and put forward a hypothesis according to which T213 and T201 play an important role in the formation of different types of peroxodiiron(III) species in MMOH and ToMOH. The hypothesis is indirectly supported by the QM/MM calculations of the peroxodiiron(III) models of ToMOH and the theoretically computed Mo?ssbauer spectra. It also helps explain the formation of two distinct peroxodiiron(III) species in the T201S mutant of ToMOH. Additionally, a role for the ToMOD regulatory protein, which is essential for intermediate formation and protein functioning in the ToMO system, is advanced. We find that the low quadrupole splitting parameter in the Mo?ssbauer spectrum observed for a ToMOH(peroxo) intermediate can be explained by protonation of the peroxo moiety, possibly stabilized by the T201 residue. Finally, similarities between the oxygen activation mechanisms of the monooxygenases and cytochrome P450 are discussed.     

Effects of the protein environment on the structure and energetics of active sites of metalloenzymes. ONIOM study of methane monooxygenase and ribonucleotide reductase.

As the first application of our recently developed ONIOM2(QM:MM) and ONIOM3(QM:QM:MM) codes to the metalloenzymes with a large number of protein residues, two members of the non-heme protein family, methane monooxygenause and ribonucleotide reductase, have been chosen. The "active-site + four alpha-helical fragments" model was adopted which includes about 1000 atoms from 62 residues around the Fe-centered spheres. Comparison of the active-site geometries of MMOH and R2 units optimized with this model with those obtained with the "active site only" (with only 39-46 atoms) model and the X-ray results clearly demonstrates the crucial role of the active site-protein interaction in the enzymatic activities.     

Optimized quantum mechanics/molecular mechanics strategies for nitrile vibrational probes: acetonitrile and para-tolunitrile in water and tetrahydrofuran

The nitrile (Ctriple bondN) group is a powerful probe of structure and dynamics because its vibrational frequency is extraordinarily sensitive to the electrostatic and chemical characteristics of its local environment. For example, site-specific nitrile labels are useful indicators of protein structure because their infrared (IR) absorption spectra can clearly distinguish between solvent-exposed residues and residues buried in the hydrophobic core of a protein. In this work, three variants of the optimized quantum mechanics/molecular mechanics (OQM/MM) technique for computing Ctriple bondN vibrational frequencies were developed and assessed for acetonitrile in water. For the most robust variant, the transferability of the OQM/MM methodology to different solutes and solvents was evaluated by simulating the IR absorption spectra of para-tolunitrile in water and tetrahydrofuran and comparing to experiment and density functional theory (DFT) calculations. The OQM/MM frequencies compared favorably to DFT for para-tolunitrile/water, and the calculated IR absorption spectra are in qualitative agreement with experiment. This suggests that a single parametrization of the OQM/MM technique is reasonable for the calculation of nitrile line shapes when the probe is attached to different chemical moieties and when the label experiences local environments of different polarity.     

Lennard-Jones parameters for the combined QM/MM method using the B3LYP/6-31G*/AMBER potential

A combined DFT quantum mechanical and AMBER molecular mechanical potential (QM/MM) is presented for use in molecular modeling and molecular simulations of large biological systems. In our approach we evaluate Lennard-Jones parameters describing the interaction between the quantum mechanical (QM) part of a system, which is described at the B3LYP/6-31+G* level of theory, and the molecular mechanical (MM) part of the system, described by the AMBER force field. The Lennard-Jones parameters for this potential are obtained by calculating hydrogen bond energies and hydrogen bond geometries for a large set of bimolecular systems, in which one hydrogen bond monomer is described quantum mechanically and the other is treated molecular mechanically. We have investigated more than 100 different bimolecular systems, finding very good agreement between hydrogen bond energies and geometries obtained from the combined QM/MM calculations and results obtained at the QM level of theory, especially with respect to geometry. Therefore, based on the Lennard-Jones parameters obtained in our study, we anticipate that the B3LYP/6-31+G*/AMBER potential will be a precise tool to explore intermolecular interactions inside a protein environment.     

Crystal structures of the soluble methane monooxygenase hydroxylase from Methylococcus capsulatus (Bath) demonstrating geometrical variability at the dinuclear iron active site.

The oxidation of methane to methanol is performed at carboxylate-bridged dinuclear iron centers in the soluble methane monooxygenase hydroxylase (MMOH). Previous structural studies of MMOH, and the related R2 subunit of ribonucleotide reductase, have demonstrated the occurrence of carboxylate shifts involving glutamate residues that ligate the catalytic iron atoms. These shifts are thought to have important mechanistic implications. Recent kinetic and theoretical studies have also emphasized the importance of hydrogen bonding and pH effects at the active site. We report here crystal structures of MMOH from Methylococcus capsulatus (Bath) in the diiron(II), diiron(III), and mixed-valent Fe(II)Fe(III) oxidation states, and at pH values of 6.2, 7.0, and 8.5. These structures were investigated in an effort to delineate the range of possible motions at the MMOH active site and to identify hydrogen-bonding interactions that may be important in understanding catalysis by the enzyme. Our results present the first view of the diiron center in the mixed-valent state, and they indicate an increased lability for ferrous ions in the enzyme. Alternate conformations of Asn214 near the active site according to redox state and a distortion in one of the alpha-helices adjacent to the metal center in the diiron(II) state have also been identified. These changes alter the surface of the protein in the vicinity of the catalytic core and may have implications for small-molecule accessibility to the active site and for protein component interactions in the methane monooxygenase system. Collectively, these results help to explain previous spectroscopic observations and provide new insight into catalysis by the enzyme.     

DFT calculations of isomer shifts and quadrupole splitting parameters in synthetic iron-oxo complexes: applications to methane monooxygenase and ribonucleotide reductase

To predict isomer shifts and quadrupole splitting parameters of Fe atoms in the protein active sites of methane monooxygenase and ribonucleotide reductase, a correlation between experimental isomer shifts ranging 0.1-1.5 mm s(-)(1) for Fe atoms in a training set with the corresponding density functional theory (DFT) calculated electron densities at the Fe nuclei in those complexes is established. The geometries of the species in the training set, consisting of synthetic polar monomeric and dimeric iron complexes, are taken from the Cambridge structural database. A comparison of calculated M?ssbauer parameters for Fe atoms from complexes in the training set with their corresponding experimental values shows very good agreement (standard deviation of 0.11 mm/s, correlation coefficient of -0.94). However, for the Fe atoms in the active sites of the structurally characterized proteins of methane monooxygenase and ribonucleotide reductase, the calculated M?ssbauer parameters deviate more from their experimentally measured values. The high correlation that exists between calculated and observed quadrupole splitting and isomer shift parameters for the synthetic complexes leads us to conclude that the main source of the error arising for the protein active sites is due to the differing degrees of atomic-level resolution for the protein structural data, compared to the synthetic complexes in the training set. Much lower X-ray resolutions associated with the former introduce uncertainty in the accuracy of several bond lengths. This is ultimately reflected in the calculated isomer shifts and quadrupole splitting parameters of the Fe sites in the proteins. For the proteins, the closest correspondence between predicted and observed M?ssbauer isomer shifts follows the order MMOH(red), RNR(red), MMOH(ox), and RNR(ox), with average deviations from experiment of 0.17, 0.17, 0.17-0.20, and 0.32 mm/s, but this requires DFT geometry optimization of the iron-oxo dimer complexes.     

Nature of the Fe-O2 bonding in oxy-myoglobin: effect of the protein

The nature of the Fe-O2 bonding in oxy-myoglobin was probed by theoretical calculations: (a) QM/MM (hybrid quantum mechanical/molecular mechanical) calculations using DFT/MM and CASSCF/MM methods and (b) gas-phase calculations using DFT (density functional theory) and CASSCF (complete active space self-consistent field) methods. Within the protein, the O2 is hydrogen bonded by His64 and the complex feels the bulk polarity of the protein. Removal of the protein causes major changes in the complex. Thus, while CASSCF/MM and DFT/MM are similar in terms of state constitution, degree of O2 charge, and nature of the lowest triplet state, the gas-phase CASSCF(g) species is very different. Valence bond (VB) analysis of the CASSCF/MM wave function unequivocally supports the Weiss bonding mechanism. This bonding arises by electron transfer from heme-Fe(II) to O2 and the so formed species coupled then to a singlet state Fe(III)-O2(-) that possesses a dative sigma(Fe-O) bond and a weakly coupled pi(Fe-O2) bond pair. The bonding mechanism in the gas phase is similar, but now the sigma(Fe-O) bond involves higher back-donation from O2(-) to Fe(III), while the constituents of pi(Fe-O2) bond pair have greater delocalization tails. The protein thus strengthens the Fe(III)-O2(-) character of the complex and thereby affects its bonding features and the oxygen binding affinity of Mb. The VB model is generalized, showing how the protein or the axial ligand of the oxyheme complex can determine the nature of its bonding in terms of the blend of the three bonding models: Weiss, Pauling, and McClure-Goddard.     

Structure and redox properties of the protein, rubredoxin, and its ligand and metal mutants studied by electronic structure calculation

The changes in the structural and electronic properties accompanying metal ionization of the iron-containing protein, rubredoxin, and of some ligand and metal mutants, have been explored using density functional theory (DFT) calculations of the metal atom coordinated to the four immediate residues. Both isolated and embedded cluster studies have been carried out, the latter using the hybrid quantum mechanics/molecular mechanics (QM/MM) approach. The replacement of a cysteine by a serine residue has a considerable effect on both the electronic and geometric structure of the core, which can be qualitatively understood on the basis of the isolated cluster studies. The modulation of these properties caused by the protein environment is quite accurately described by the QM/MM calculations. The predicted core geometries are in good accord with both X-ray and EXAFS data, and the changes in the redox potentials are predicted, at least semiquantitatively, by considering only the core part of the protein.     

Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.     

Theoretical studies on the reactivity of molybdenum enzymes

In recent years, advances in theoretical methods and computational capabilities have made it possible to investigate reaction mechanisms in enzymes. Density functional theory (DFT) is commonly used to study reactions in model systems, while combined quantum mechanical/molecular mechanical (QM/MM) approaches allow the treatment of the complete solvated enzyme and thus provide insight into the mechanistic influence of the protein environment. This review starts with a brief overview over the available DFT and QM/MM methodology and then summarizes recent theoretical studies on biocatalysis by molybdenum-containing enzymes. It focuses on the reactions in members of the dimethylsulfoxide reductase, sulfite oxidase, and xanthine oxidase families, with special emphasis on the QM/MM studies of the latter. It concludes with a brief survey of theoretical work on some other molybdenum- and tungsten-containing enzymes.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Dynamic structures of phosphodiesterase-5 active site by combined molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical calculations

Various quantum mechanical/molecular mechanical (QM/MM) geometry optimizations starting from an x-ray crystal structure and from the snapshot structures of constrained molecular dynamics (MD) simulations have been performed to characterize two dynamically stable active site structures of phosphodiesterase-5 (PDE5) in solution. The only difference between the two PDE5 structures exists in the catalytic, second bridging ligand (BL2) which is HO- or H2O. It has been shown that, whereas BL2 (i.e. HO-) in the PDE5(BL2 = HO-) structure can really bridge the two positively charged metal ions (Zn2+ and Mg2+), BL2 (i.e. H2O) in the PDE5(BL2 = H2O) structure can only coordinate Mg2+. It has been demonstrated that the results of the QM/MM geometry optimizations are remarkably affected by the solvent water molecules, the dynamics of the protein environment, and the electronic embedding charges of the MM region in the QM part of the QMM/MM calculation. The PDE5(BL2 = H2O) geometries optimized by using the QM/MM method in different ways show strong couplings between these important factors. It is interesting to note that the PDE5(BL2 = HO-) and PDE5(BL2 = H2O) geometries determined by the QM/MM calculations neglecting these three factors are all consistent with the corresponding geometries determined by the QM/MM calculations that account for all of these three factors. These results suggest the overall effects of these three important factors on the optimized geometries can roughly cancel out. However, the QM/MM calculations that only account for some of these factors could lead to considerably different geometries. These results might be useful also in guiding future QM/MM geometry optimizations on other enzymes.