|Home | About | Journals | Submit | Contact Us | Français|
Quantum mechanics/molecular mechanics (QM/MM) hybrid methods are currently the most powerful computational tools for studies of structure/function relations and structural refinement of macrobiomolecules (e.g., proteins and nucleic acids). These methods are highly efficient since they implement quantum chemistry techniques for modeling only the small part of the system (QM layer) that undergoes chemical modifications, charge transfer, etc., under the influence of the surrounding environment. The rest of the system (MM layer) is described in terms of molecular mechanics force fields, assuming that its influence on the QM layer can be roughly decomposed in terms of electrostatic interactions and steric hindrance. Common limitations of QM/MM methods include inaccuracies in the MM force fields when polarization effects are not explicitly considered, and the approximate treatment of electrostatic interactions at the boundaries between QM and MM layers. This paper reviews recent advances in the development of computational protocols that allow for rigorous modeling of electrostatic interactions in extended systems beyond the common limitations of QM/MM hybrid methods. We focus on the Moving-Domain QM/MM (MoD-QM/MM) methodology that partitions the system into many molecular domains and obtains the electrostatic and structural properties of the whole system from an iterative self-consistent treatment of the constituent molecular fragments. We illustrate the MoD-QM/MM method as applied to the description of photosystem II as well as in conjunction with the application of spectroscopically constrained QM/MM optimization methods, based on high-resolution spectroscopic data (extended X-ray absorption fine structure spectra, and exchange coupling constants).
In recent years, the development of computer hardware has allowed computational chemistry methods to establish themselves as powerful tools for studies of structure/function relations and structural refinement of macrobiomolecules (e.g., proteins and nucleic acids). In particular, quantum mechanics/molecular mechanics (QM/MM) hybrid methods have become the most popular approaches for modeling prosthetic groups that undergo significant structural or electronic modifications when embedded in macromolecular structures (Maseras and Morokuma 1995;Vreven and Morokuma 2000a, b;Warshel and Levitt 1976). The main aspect, common to all QM/MM methods, is the partitioning of the system into two fragments: a small part (QM layer) that includes the prosthetic group as described by quantum chemistry methods and the rest of the system (MM layer) that is modeled in terms of molecular mechanics force fields, since its influence can be roughly decomposed in terms of steric hindrance and electrostatic interactions. This basic idea is at the heart of the Moving-Domain QM/MM (MoD-QM/MM) methodology and the applications studies reviewed in this paper (Gascon et al. 2006a;Menikarachchi and Gascon 2008;Newcomer et al. 2009).
The various different versions of QM/MM methods differ in the way the energy is evaluated (e.g., additive and subtractive schemes are available); in the way the boundary between QM and MM layers is treated (e.g., link-atoms, localized bond orbitals and the pseudopotential approach); on the specific level of quantum chemistry and MM force field; and on how the electrostatic interactions between QM and MM layers are described, including mechanical embedding, electronic embedding, and polarized embedding schemes.
QM/MM methodologies require a starting structure, which is often obtained from X-ray diffraction models. However, the available X-ray diffraction structures might have only moderate resolution, or suffer from radiation damage due to the high doses of X-ray necessary for data collection. In addition, the potential energy landscape of complex molecular systems is often so rugged that the obtained minimum energy structures are typically strongly correlated to the initial guess configuration. Another major problem with traditional QM/MM modeling is the limited accuracy of the electrostatic potentials provided by molecular mechanics force fields. These limitations often introduce inaccuracies in the description of the interactions between the QM and MM layers in QM/MM hybrid models.
A structural refinement scheme has been developed to circumvent the risk of obtaining minimum energy structures that are inconsistent with high-resolution spectroscopic data (Sproviero et al. 2008a). The method adjusts the overall configuration to minimize the total QM/MM energy subject to the constraint of a minimum χ2 deviation between the calculated and experimental properties of the system as described by an independent set of experimental data. This methodology has been successfully applied in conjunction with high-resolution polarized extended X-ray absorption fine structure (P-EXAFS) spectra to refine structural models of the Mn4Ca active site of PSII (Sproviero et al. 2008a).
The MoD-QM/MM methodology has been developed to provide an accurate description of the electrostatic interactions between MM and QM layers in QM/MM hybrid models. The MoD-QM/MM approach partitions the system into many molecular domains and obtains the electrostatic and structural properties of the whole system from an iterative self-consistent treatment of the constituent molecular fragments (Gascon et al. 2006a;Menikarachchi and Gascon 2008;Newcomer et al. 2009). This involves a polarization scheme rendering a description of prosthetic groups at the polarized embedding QM/MM level. Self-consistent optimized structures and electrostatic-Potential (ESP) atomic charges are obtained, modeling each fragment as a QM layer polarized by the distribution of atomic charges in the surrounding molecular fragments. The resulting computational task scales linearly with the size of the system, bypassing the enormous demands of memory and computational resources that would otherwise be required by ‘brute-force’ full QM calculations of the complete system. Computations for several benchmark systems that allow for full QM calculations have shown that the MoD-QM/MM method produces ab initio quality electrostatic potentials and optimized structures in quantitative agreement with full QM results (Gascon et al. 2006a;Menikarachchi and Gascon 2008;Newcomer et al. 2009).
The MoD-QM/MM approach is not limited to any specific level of quantum chemistry and has been applied in conjunction with Density Functional Theory (DFT) in studies of photosystem II and other macrobiological systems (Gascon et al. 2006a;Gascon et al. 2006b;Menikarachchi and Gascon 2008;Sproviero et al. 2006b, 2008d). The main advantage of DFT is its efficiency for the description of electron correlation at a computational cost comparable to single-determinant calculations, and the efficient description of metal centers and inorganic cores of metalloproteins in various possible spin-electronic states. MoD-QM/MM calculations have modeled charge redistribution, polarization effects and changes in the ligation scheme of metal centers induced by oxidation state transitions, bond breaking and forming processes associated with chemical reactivity and proton transfer. In addition, the MoD-QM/MM method has provided a reliable description of electronic and magnetic properties of prosthetic groups, embedded in complex molecular environments, and has allowed for the application of spectroscopically constrained QM/MM optimization methods based on simulations of X-ray absorption spectra and direct comparisons with high-resolution spectroscopic data (e.g. FT-IR, NMR, EXAFS, and EPR) (Sproviero et al. 2008c, d).
The work reviewed in this paper illustrates the MoD-QM/MM methodology as applied to the description of macromolecules of biological interest, including photosystem II and in conjunction with spectroscopically constrained QM/MM structural refinement schemes based on simulations of X-ray absorption spectra and direct comparisons with experimental data. The application studies show that the MoD-QM/MM methodology is well suited to study equilibrium and transition states of biological systems beyond the limitations of model structures obtained by X-ray diffraction or NMR spectroscopy, including reaction intermediates generated by water exchange, proton transfer and the rearrangement of ligands during the catalytic cycle of water oxidation at the oxygen-evolving complex (OEC) of photosystem II (PSII) (Sproviero et al. 2006b;Sproviero et al. 2007;Sproviero et al. 2008b, d). The paper is organized as follows. First, we review the QM/MM methodology, including a discussion of its specific implementation to studies of high-valent metal complexes in metalloproteins. Then, we review the MoD-QM/MM method and a Spectroscopically Constrained (SC) QM/MM refinement scheme based on simulations of X-ray absorption spectra and direct comparisons with experimental measurements.
The MoD-QM/MM methodology (Gascon et al. 2006a;Menikarachchi and Gascon 2008;Newcomer et al. 2009;Sproviero et al. 2008e) has been applied in conjunction with the two-layer ONIOM Electronic-Embedding (EE) link-hydrogen atom method (Dapprich et al. 1999), as implemented in Gaussian 03 (Frisch et al. 2004). The ONIOM approach divides the system into two partitions (Fig. 1), including a reduced system X (the QM layer) and the rest of the system (region Y). The total energy E is obtained from three independent calculations:
where EMM,X+Y is the energy of the complete system computed at the molecular mechanics level of theory, while EQM,X and EMM,X correspond to the energy of the reduced-system X computed at the QM and MM levels of theory, respectively. Electrostatic interactions between regions X and Y are included in the calculation of both EQM,X and EMM,X at the quantum mechanical and molecular mechanics levels, respectively. Thus, the electrostatic interactions computed at the MM level in EMM,X and EMM,full cancel and the resulting DFT QM/MM evaluation of the total energy involves a quantum mechanical description of the polarization of the reduced system, due to the electrostatic influence of the surrounding protein environment.
Figure 2 illustrates the MoD-QM/MM computational protocol where the system is partitioned into n molecular domains by using a simple space-domain decomposition scheme (Gascon et al. 2006a;Menikarachchi and Gascon 2008;Newcomer et al. 2009). ElectroStatic-Potential (ESP) atomic-charges of the constituent molecular domains are sequentially computed by using QM/MM hybrid methods (e.g., the two-layer ONIOM-EE scheme) considering the updated distribution of atomic charges in the molecular domains previously modeled as QM layers in the cycle. The whole cycle of QM/MM calculations is iterated several times until obtaining a self-consistent point-charge model of the electrostatic potential that explicitly includes mutual polarization interactions between the constituent molecular domains. The underlying computational effort usually scales linearly with the size of the system since the iterative scheme typically converges within a few iteration-cycles (i.e., 4 or 5 cycles). The method, thus, bypasses the enormous demands of memory and computational resources that would be required by a ‘brute-force’ quantum chemistry calculation of the complete system.
The MoD-QM/MM method typically corrects the atomic charges prescribed by popular molecular mechanics force fields (Cornell et al. 1995;Dykstra 1993) by only 10–20%. However, summing the corrections from multiple domains over the entire QM/MM interface often introduces significant corrections (10–15 kcal/mol) to the total interaction between the QM and MM layers. These energy corrections are often critical since they are comparable to energy differences associated with different protonation states, hydrogen-bonding structures, or the energy splitting between high-spin and low-spin states of the metal cluster.
The accuracy of the resulting molecular electrostatic potential comes at the expense of transferability since the atomic charges obtained for one configuration are, in principle, nontransferable to other configurations of the system. Therefore, while accurate, the computed electrostatic potential is most useful for applications where the system does not involve large conformational changes, such as intermediate structures of the OEC of PSII where there are no significant temperature-dependent variations in structure, protonation state, or charge localization, as recently reported by studies of X-ray absorption at 20 K and room-temperature (Haumann et al. 2005).
The resulting generalized computational protocol requires the space-domain fragmentation of the system into partially overlapping molecular domains and an iterative cycle of QM/MM calculations as depicted in Fig. 2. Each step in the cycle considers a different molecular domain as a QM layer and updates its distribution of ESP atomic charges, after relaxing its geometry by energy minimization. The rest of the system is considered at the MM level, with an updated distribution of atomic charges for all molecular domains previously considered as QM layers in the iteration cycle
The geometry optimization is performed with standard methods (Fan and Ziegler 1991;Versluis and Ziegler 1988), using the Newton-Raphson method and updating the second derivative (Hessian matrix) with the Broyden-Fletcher-Goldfarb-Shanno strategy (Schlegel 1987). Starting with the first molecular domain R1, its geometry is relaxed as influenced by the surrounding environment. The ESP charges of the QM layer in its minimum energy configuration are computed and updated for subsequent steps in the cycle (see Fig. 2). The relaxed configuration and ESP atomic-charges of the next domain (R2) are computed analogously, considering the updated charges and configuration of the previously considered molecular domain. The procedure is subsequently applied to all fragments and the entire computational cycle is iterated several times until reaching self-consistency in the description of the geometry and distribution of atomic charges as influenced by mutual polarization effects. The capabilities of the resulting protocol have been recently demonstrated at the ONIOM-EE level, as applied to the description of benchmark models systems that allowed for direct comparisons with full QM calculations as well as in the study of the structural characterization of the DNA Oxytricha nova guanine quadruplex (G4) (Newcomer et al. 2009).
The typical energy landscape of macrobiological systems is very rugged, with a manifold of configurations of comparable energy (e.g., within the ~5 kcal/mol energy range of error of DFT calculations). Therefore, it is rather difficult to determine the most likely configuration of the system solely from the energetic analysis of possible configurations. In order to address these limitations, a spectroscopically constrained MoD-QM/MM structural refinement scheme has been introduced. The method allows one to assess the likelihood of structures of similar energies in terms of simulations of high-resolution spectroscopy (e.g., EXAFS spectra) and direct comparisons with experimental data (Sproviero et al. 2008c). These MoD-QM/MM methods have been successfully applied to the development of structural models of the OEC of PSII (Sproviero et al. 2008a) and other systems (Soderhjelm and Ryde 2006) and proved particularly valuable to refine configuration of prosthetic groups with multiple high-valent metal centers beyond the current limitations of DFT-QM/MM methods.
The refinement algorithm starts with a reference configuration of the system (e.g., the MoD-QM/MM minimum energy structure) and iteratively adjusts its configuration to minimize the mean square deviation between the simulated and experimental EXAFS spectra, subject to the constraint of minimum displacements relative to the reference configuration (see Fig. 3). Upon convergence, the resulting structures are often consistent with the reference MoD-QM/MM minimum energy structures (with energies and root-mean-square deviations within the range of error of DFT QM/MM minimization methods) and predict high-resolution spectroscopic data in quantitative agreement with experiments. The refinement scheme allows one to narrow the range of possible structures through direct comparisons with high-resolution spectroscopic data. As mentioned earlier, this method is particularly valuable for metalloproteins for which the most accurate and reliable experimental information on the structure of the prosthetic group metal clusters comes from high-resolution spectroscopy, e.g., EXAFS (Yano et al. 2006).
The MoD-QM/MM refinement protocol has been applied in conjunction with simulations of EXAFS spectra for the development of models of the OEC of PSII (Sproviero et al. 2008c). The EXAFS spectra were computed as a function of the nuclear coordinates and compared to the experimental spectra. The nuclear coordinates of the Mn cluster and coordinated ligands were adjusted to minimize the mean-square deviation , where Aicalc and Aiexp correspond to the calculated and experimental EXAFS amplitudes for a grid of electron kinetic energies Ki with i = 1, 2, …, N. The χ2 minimization is performed subject to the constraints of minimum displacements of nuclear positions relative to their coordinates in the reference MoD-QM/MM configuration. The minimization algorithm implements a multidimensional variable metric method through the Broyden-Fletcher-Goldfarb-Shanno (BFGS) approach (Broyden 1970;Fletcher 1970;Goldfarb 1970;Shanno 1970), using the Cartesian coordinates xi of the nuclei as input variables and an external simulation procedure that returns the corresponding spectrum as computed by the program FEFF8 (version 8.2) (Ankudinov et al. 2002;Ankudinov et al. 1998;Bouldin et al. 2001).
Ab initio simulations of EXAFS spectra involve solving the multiscattering problem associated with the photodetached electrons emitted upon X-ray absorption by the Mn centers. The quantum mechanical interference of outgoing photoelectrons with the backscattering waves from atoms surrounding the Mn ions, gives rise to oscillations of EXAFS intensities as a function of the kinetic energy (or momentum) of the photoelectron. The Fourier transform of these oscillations provide information on the Mn-Mn distances, the coordination of Mn ions, and changes in the Mn coordination sphere as determined by changes in the oxidation state of the metal centers. Calculations are typically carried out according to the Real Space Green function approach as implemented in the program FEFF8 (version 8.2) (Ankudinov et al. 2002;Ankudinov et al. 1998;Bouldin et al. 2001), using the theory of the oscillatory structure due to multiple-scattering originally proposed by Kronig (Kronig 1931;Kronig and Penney 1931) and worked out in detail by Sayers (Sayers et al. 1971), Stern (Stern 1974) Lee and Pendry (Lee and Pendry 1975), and by Ashley and Doniach (Ashley and Doniach 1975).
One of the advantages of the structural refinement scheme based on EXAFS simulations and comparisons with experimental data is that EXAFS experiments require doses of X-ray radiation that are substantially lower than those required by X-ray diffraction (XRD) measurements (Yano et al. 2005). In addition, EXAFS provides metal-to-metal and metal-to-ligand distances with high accuracy (~0.02 Å) and a resolution of ~0.1 Å. For polarized EXAFS experiments on one-dimensionally oriented membranes or three-dimensionally oriented single crystals, the EXAFS amplitude is orientation dependent and proportional to 2 (θ), where θ is the angle between the electric field vector of the polarized X-ray beam and the absorber back-scatterer vector. Therefore, the structural refinement can provide additional geometric information about the metal site in metalloprotein crystals (Flank et al. 1986;Scott et al. 1982;Yano et al. 2005).
The SC-QM/MM analysis of MoD-QM/MM structural models can be based on calculations of magnetic properties and direct comparisons with experimental data. Biomimetic oxomanganese complexes (Cooper and Calvin 1977;Cooper et al. 1978;Limburg et al. 1999;Manchanda et al. 1994;Ruettinger et al. 1999;Sarneski et al. 1990;Stebler et al. 1986;Thorp and Brudvig 1991) and models of the Mn4Ca complex of PSII extracted from complete MoD-QM/MM models (Sproviero et al. 2006b) have been recently analyzed (Pantazis et al. 2009;Sproviero et al. 2006a). The exchange coupling constants Jij between centers i, j were computed assuming that the centers interact according to the Heisenberg Hamiltonian,
where Si is the spin corresponding to metal center i. In the weak-coupling limit, the spin states are essentially localized on individual metal centers and the exchange coupling constants are determined by splittings between energy levels as follows,
where and are the magnetic energy contributions from the pair of metal centers i and j (with Si < Sj) in the broken-symmetry (BS) and high-spin (HS) configurations, respectively. The HS configuration corresponds to equally oriented spins while in the BS state the α and β densities are allowed to localize on different atomic centers, providing a multiconfigurational character to the spin state (Noodleman 1981;Noodleman and Case 1992;Noodleman and Davidson 1986;Noodleman et al. 1995).
The physical picture supporting this calculation is that in the weak-coupling limit, atomic d-orbitals remain singly occupied and the spin states are essentially localized on individual metal centers. Under those conditions, the ground state of a multinuclear complex involves anti-parallel coupling between spins localized on different centers. However, metal-metal interactions induce electronic delocalization over multiple centers, forming bonding and antibonding orbitals by either direct overlap of atomic orbitals or superexchange. In the BS state, α and β electronic densities can be localized on different centers. These relaxed symmetry requirements allow for substantial interactions mediated by d-orbitals of the metal centers and the p-orbitals of the ligands.
To perform ab initio calculations of J, it is convenient to define spin configurations that differ only in the orientation of the spin of one center with respect to the HS configuration. The energy difference associated with flipping a spin of center k with respect to the HS configuration is,
where k = 1, …, N. Considering that there are N + 1 relevant spin configurations, including one with the highest spin and N with a spin flipped in only one center, the coupling constants can be readily obtained since the number of equations is equal or greater than the number of unknown Js.
MoD-QM/MM models of the OEC of PSII have been based on the 1S5L X-ray crystal structure of PSII from the cyanobacterium Thermosynechococcus elongatus 1S5L (Ferreira et al. 2004) (see Fig. 4). The models explicitly considered 1987 atoms of PSII, including the proposed Mn3CaO4Mn unit and all amino acid residues with α-carbons within 15 Å from any atom in the OEC metal ion cluster. The coordination spheres of the Mn centers were completed by hydration, assuming a minimum displacement of the proteinaceous ligands from their crystallographic positions. The usual coordination of 5 or 6 ligands, for high-valent Mn ions with oxidation states III and IV, was also satisfied. In addition, the coordination of calcium (typically with 6–8 ligands) was satisfied by coordination of water molecules and negative counterions.
The MoD-QM/MM molecular structures were hydrated by including water molecules that could bind to polar amino acid residues, or existing water molecules in the model. The hydration process involved inserting the model system in a large box containing a thermal distribution of water molecules. After adding to the model the water molecules that did not interfer sterically with the protein residues, the structure of hydrogen bonds was optimized by energy minimization. The hydration and relaxation processes were repeated several times until the number of water molecules in the system converged, typically with the addition of 85–100 water molecules. A few of the additional water molecules (up to six molecules) attached to metal centers in the cuboidal Mn3CaO4Mn cluster, including 2 waters molecules that were considered as potential substrate water molecules responsible for O-O bond formation in the S4 to S0 transition.
MoD-QM/MM computations of the OEC of PSII could only be efficiently performed, according to the ONIOM-EE level of theory, after obtaining high-quality initial-guess states of the QM layer. These correspond to spin-electronic states of the cluster of Mn ions based on ligand field theory as implemented in Jaguar 5.5 (Jaguar 5.5; Schroedinger). The electronic states of the Mn cluster involve antiferromagnetic couplings between manganese centers and are defined by broken-symmetry (BS) states, in which the α and β densities are allowed to localize on different atomic fragments (Noodleman 1981;Noodleman and Case 1992;Noodleman and Davidson 1986;Noodleman et al. 1995). Each fragment is simply a group of atoms that would be bonded together even if all bonds between transition metal atoms and non-transition metal atoms were broken. Typically, the system is broken into ligand fragments and transition metal fragments. The BS approach, in conjunction with DFT, is the most accurate approach available for these types of high-valent metal clusters. BS-DFT also facilitates the ligand field analysis, including metal d→d, charge transfer (ligand→metal, metal→ligand), and intervalence charge transfer (metal→metal or ligand→ligand) transitions.
From a theoretical point of view, changes in the spin state will affect the spatial distribution and energetics of the molecular orbitals. As different electronic distributions will give different relaxed geometries, having a proper description of the spin state of the system allows a correct description of the corresponding relaxed geometry (Lundberg and Siegbahn 2004).
Once the initial state was properly defined, the combined approach exploits important capabilities of ONIOM as implemented in Gaussian 03 (Frisch et al. 2004), including both the link-hydrogen atom scheme for efficient and flexible definitions of QM layers and the possibility of modeling open-shell systems by performing unrestricted DFT (e.g. UB3LYP) calculations in a variety of spin-electronic states.
The QM layer of MoD-QM/MM calculations (region X) included the Mn3CaO4-Mn complex and the directly ligated proteinaceous carboxylate groups of D1-D189, CP43-E354, D1-A344, D1- E333, D1-D170, D1-D342. In addition, the imidazole ring of D1- H332 as well as bound water molecules, hydroxide and chloride ions were included in the QM layer. The rest of the system defined the MM layer (region Y). The boundary between the QM and MM layers was defined by cutting covalent bonds in the corresponding amino acid residues (i.e., D1-D189, CP43-E354, D1-A344, D1- E333, D1-D170, D1-D342, and D1-H332) and completing the covalency of frontier atoms according to the standard link-hydrogen atom scheme (see Fig. 1).
The efficiency of MoD-QM/MM calculations has been optimized by using a combination of basis sets for the QM layer (Sproviero et al. 2006c). Such a combination of basis sets has been validated through extensive benchmark calculations on high-valent manganese complexes (Sproviero et al. 2006a). The molecular structure beyond the QM layer is described by the Amber MM force field (Cornell et al. 1995).
Fully relaxed MoD-QM/MM model structures were obtained by energy minimization of the QM layer, surrounded by the MM layer of amino acid residues with α-carbons within 15–20 Å from any atom in the OEC ion cluster. The buffer shell is subjected to harmonic constraints in order to preserve the natural shape of the system even in the absence of the membrane. Upon convergence, the resulting optimized structures are analyzed and assessed not only in terms of the total energy of the system but also in terms of the comparison of structural, electronic and mechanistic features that should be consistent with experimental data. Such a comprehensive analysis is crucial for these complex biological systems often with many possible configurations with similar energies (e.g., within the ~5 kcal/mol range of error of DFT QM/MM calculations).
The MoD-QM/MM minimum energy structure of the OEC of PSII in the S1 state was used as the initial reference configuration for the EXAFS structural refinement algorithm based on P-EXAFS data from Yano et al. (Yano et al. 2008). The algorithm searched for a configuration with the smallest possible displacement relative to the initial structure with the constraint of a minimum χ2 deviation between the calculated and experimental EXAFS data. Upon convergence of the structural refinement procedure, the isotropic and polarized EXAFS spectra of the resulting refined (R)-QM/MM model were found to be in quantitative agreement with high-resolution EXAFS measurements (Haumann et al. 2005;Yano et al. 2005) (Figs. 5 and and6)6) as well as with the description of the protein environment suggested by XRD models (Ferreira et al. 2004;Loll et al. 2005) (Fig. 7).
The predicted metal-metal and metal-ligand distances are consistent with EXAFS studies of frozen solutions of PSII that have provided accurate distances for Mn-Mn, Mn-Ca, and Mn/Ca-ligand vectors in the Mn4Ca cluster of PSII (Dau et al. 2001;Haumann et al. 2005;Robblee et al. 2001;Yano et al. 2005). In addition, the R-QM/MM model provides insights on structural details currently unresolved by XRD spectroscopy, including the spatial orientation of the cluster, which has been inherit from the experimental models (Yano et al. 2006), as well as the approximate placement of the models within the protein environment, which is achieved in this study (Sproviero et al. 2008a). However, it is important to emphasize that the R-QM/MM model does not rule out other possible structures since it is not the result of an exhaustive refinement of all possible structures of the system but rather a local minimum solution relative to the reference MoD-QM/MM structure.
The oxomanganese inorganic core of the R-QM/MM model involves a cuboidal core Mn3CaO4 with a dangler manganese ligated to a corner μ4-oxide ion. Dangler Mn models were originally proposed (Britt et al. 2004;Britt et al. 2000) within a set of 11 possible structural motifs (Derose et al. 1994), suggested by solution EXAFS studies, from which the dimer-of-dimers model was highly favored (Britt et al. 2004;Britt et al. 2000;Yachandra 2002). Similar dangler Mn models were also preferred by 55Mn electronuclear double resonance (ENDOR) studies that strongly disfavored the dimer-of-dimers motif over models with a trinuclear Mn core and a fourth Mn set off from the core by a longer Mn-Mn internuclear distance (Britt et al. 2000). However, recent polarized-EXAFS studies of single crystals have disfavored dangler models (Yano et al. 2006) and have suggested instead four alternative models (models I, II, IIa and III (Yano et al. 2006)) for which the simulated EXAFS spectra (in the absence of ligands) match the experimental data. These empirical models, however, remain uncertain since structural differences between some of them (e.g., models I and II) are large. Furthermore, placing any of the four models into the XRD structures of PSII results in unsatisfactory metal-ligand distances, coordination numbers and geometries. The underlying inconsistencies might be due to the coordinate error in the X-ray crystal models, or due to the intrinsic limitations of polarized-EXAFS models built by neglecting the contributions of electron scattering paths from the ligands (i.e., assuming that the EXAFS amplitudes result solely from scattering paths associated with metal centers and oxo-bridges in the inorganic core). These issues should be addressed in the development of empirical models since calculations of polarized-EXAFS spectra suggest that the contributions of backscattering from the ligands could introduce amplitude corrections of the order of 10–15 % (Sproviero et al. 2008c).
Benchmark calculations of exchange coupling constants involved the analysis of of di-, tri- and tetra-nuclear Mn complexes (Sproviero et al. 2006a), previously characterized both chemically and spectroscopically (Cooper and Calvin 1977;Cooper et al. 1978) (see Table 1). These include the di-μ-oxo bridged dimers 1: [MnIVMnIII(μ-O)2-(H2O)2(terpy)2]3+ (terpy = 2,2’:6,2”-terpyridine) and 2: [MnIIIMnIV(μ-O)2(phen)2]3+ (phen = 1,10-phenanthroline), and the trimer 3: [Mn3O4(bpy)4(OH2)2]4+ (bpy = 2,2’-bipyridine), shown in Fig. 8. These synthetic complexes have been previously characterized by electrochemical and high-resolution EPR and X-ray diffraction methods (Cooper and Calvin 1977;Cooper et al. 1978).
Table 1 compares exchange-coupling constants for synthetic oxomanganese complexes 1–3 to readily available experimental data (Cooper and Calvin 1977;Cooper et al. 1978;Sarneski et al. 1990;Stebler et al. 1986). Exchange coupling constants were obtained at the DFT/B3LYP level and show very good agreement with experimental data. These results serve to validate the analyzed model structure and its electronic description, as applied to the description of magnetic properties of oxomanganese complexes with pre-defined spin-electronic states. The reported agreement between calculated and experimental exchange coupling constants was consistent with anti-ferromagnetic couplings of high-valent metal centers in high-spin states. The analysis of the electronic structures show that the most important exchange interactions between Mn centers involved symmetric super-exchange pathways Jxz/xz and Jyz/zz mediated by the overlap of orbitals dxz and dyz of Mn with out-of-the-plane px orbitals of the μ-oxo bridges (Sproviero et al. 2006a). These results are particularly relevant to the description of active sites of metalloproteins with common prosthetic groups and provide insight on the origin of electronic delocalization over multiple centers (Pantazis et al. 2009).
The MoD-QM/MM methodology has been recently applied to develop structural models of transition-state structures the OEC of PSII generated by water exchange, in an effort to provide insights on the nature of water exchange rates, and rigorous validations of the structural models by direct comparisons of water-exchange rates with time-resolved mass spectrometry (MS) measurements (Sproviero et al. 2008e).
Time-resolved MS studies have determined different water exchange rates (kex) of two substrate water molecules of the OEC with bulk 18O-labelled water. The more slowly exchanging water (Wslow) was found to be associated with Ca2+, implying that the fast-exchanging water (Wfast) would be bound to a manganese ion (Hendry and Wydrzynski 2002;Hillier and Wydrzynski 2008;Messinger et al. 1995). These findings were rather surprising since high-valent manganese ions (e.g. Mn3+ or Mn4+) were expected to bind water more strongly than Ca2+. In addition, it was observed that the exchange rate of Wslow (kex slow) increased by two orders of magnitude upon S1 to S2 oxidation, with kex(S1) = 0.02 s−1 and kex(S2) = 2.0 s−1 (Hendry and Wydrzynski 2003;Hillier and Wydrzynski 2004). These rates correspond approximately to activation barriers of 20 and 17 kcal/mol, respectively. This faster exchange of Wslow in the S2 state was also intriguing since it implied that the oxidation of a manganese center should indirectly affect the exchange rate of a calcium-bound water molecule.
The QM/MM structural models allowed for a quantitative analysis of specific metal-water interactions and potential-energy profiles associated along minimum energy paths (MEPs) for water detachment and exchange (Sproviero et al. 2008e). The MEPs were found by energy minimization with respect to nuclear and electronic coordinates, while progressively detaching substrate water molecules from their corresponding coordination metal centers. Since elongation of the metal-oxygen bond is the primary step in water exchange, and presumably rate determining in this case, the resulting structural rearrangements provided insights on the detailed mechanisms of water exchange.
The MoD-QM/MM models are consistent with coordination of substrate water molecules as terminal ligands, in agreement with earlier proposals (Haumann and Junge 1999;Hoganson and Babcock 1997;McEvoy and Brudvig 2004;Messinger 2004;Schlodder and Witt 1999;Sproviero et al. 2006b;Sproviero et al. 2007), although in contrast to other suggested models that involve coordination as oxo-bridges between Mn ions (Brudvig and Crabtree 1986;Messinger 2004;Nugent et al. 2001;Pecoraro et al. 1994;Robblee et al. 2001;Yachandra et al. 1996). The energy barriers for water exchange in the MoD-QM/MM models were estimated to be 19.3, 16.6, 8.4 and 7.9 kcal/mol for exchange of water ligands from Ca2+(S1), Ca2+(S2), Mn(4)3+(S2) and Mn(4)3+(S1), respectively. Figs. 9–11 show the initial and TS configurations along the MEPs for water exchange from Ca2+ and Mn(4), as described by the DFT-QM/MM structural models of the OEC, including a detailed description of coordination bond lengths. In the initial states (Fig. 9), substrate water molecules are shown attached to the metal centers forming hydrogen bonds with the exchanging water molecules. In the TS configurations (Figs. 10 and and11),11), the coordination bonds of the substrate water molecules are stretched and the substituting water molecules are displaced relative to their original configurations. Further displacement beyond the transition state leads to coordination of the exchanging water molecule to the metal center. The underlying process involves converted exchange with dissociative character (i.e. without formation of any intermediate state of lower or higher coordination number).
The comparison of the coordination bond lengths of water ligands in TS configurations indicates that the coordination bonds are more stretched when water is exchange from Ca2+ in the S2 state, suggesting that the exchange mechanism gains dissociative character upon S1/S2 oxidation of the OEC. The electronic origin of these differences has been traced to the redistribution of charge induced by charge-transfer interactions with D1-A344 leading to a net reduction of the atomic charge of Ca2+. The electrostatic interactions with negatively charged ligands are strenghen upon S1 to S2 oxidation since the metal complex becomes positively charged. Most of the electronic density acquired by Ca2+ upon S1/S2 oxidation (Δq = −0.21) is transferred from the ligated carboxylate group of D1-A344. These results suggest that oxidation state transitions of the oxomanganese complex can indirectly regulate substrate water binding to the metal cluster by modulating charge-transfer interactions.
Water exchange in the coordination sphere of the dangling manganese is also predicted to have dissociative character in both the S1 and the S2 states, probably due to the reduced charge of the metal center in the presence of charge transfer interactions between Mn and oxo-bridges. The calculated energy barriers for water exchange from Mn(4) (7.9 and 8.4 kcal/mol in the S1 and S2 states, respectively) are comparable to the barriers in hydrated Mn complexes, such as [(H2O)2(OH)2 MnIV(μ-O)2MnIV(H2O)2(OH)2] and [MnIV(H2O)2(OH)4], where the exchange of terminal water ligands requires 8.6–9.6 kcal/mol (Lundberg et al. 2003;Tsutsui et al. 1999). In contrast, the DFT-QM/MM energy barriers for water exchange from Ca2+ are much higher (19.3 and 16.6 kcal/mol for the OEC in the S1 and S2 states, respectively). These larger barriers stabilizing substrate water molecules bound to the catalytic metal center for much longer times (ms) than typically bound to similar complexes in solutions.
The analysis of MoD-QM/MM models suggests that the higher energy barriers are mainly determined by the nature of the interactions in the OEC-binding pocket and the incomplete solvation of dissociated water ligands. These calculations, thus, provided theoretical support for the rather slow exchange rates and for the surprising experimental finding that the slow-exchanging substrate water is associated with calcium, while the fast-exchanging substrate water is coordinated with manganese. Furthermore, the QM/MM models provided a rationale for the opposite effect of the S1/S2 transition on the two water-exchange rate constants. It was concluded that the mechanism and energetics of water binding to metal centers in the OEC is not simply correlated to the formal oxidation states of the ions but rather to their corresponding partial charges, as modulated by charge-transfer interactions with coordinated ligands.
A wide range of computational methods has been developed during the last few decades to model electrostatic interactions in complex molecular systems, including the influence of polarization effects on structural refinement and the analysis of structure/function relations. However, rigorous and practical approaches to obtain ab initio quality electrostatic potentials have yet to be established. Several techniques have emerged within the field of linear scaling methods, based on the fragmentation of extended system into small molecular domains and the subsequent calculation of the properties of the entire system from accurate quantum chemistry calculations of the constituent fragments. The recently developed MoD-QM/MM hybrid method belongs to this general class of linear scaling approaches, and is based on a simple space domain decomposition scheme and the iterative self-consistent treatment of the constituent molecular fragments according to state-of-the-art QM/MM methods. The resulting methodology complements previous approaches (Dapprich et al. 1999) and extends the capabilities of current QM/MM hybrid methods to the description of electrostatic interactions that explicitly include mutual polarization effects in both the QM and MM layers.
MoD-QM/MM structural models have been used as reference configurations in conjunction with the development and application of SC-QM/MM structural refinement schemes. The SC-QM/MM refinement methodology allowed for the assessment of MoD-QM/MM molecular structures of similar energy (e.g., within ~5 kcal mol−1 from the minimum energy configuration) in terms of direct comparisons with high-resolution spectroscopic data. The method has been developed and implemented in terms of simulations of EXAFS spectra and the minimization of χ2-deviations between simulated and experimental spectra by adjustment of nuclear coordinates, subject to the constraint of minimum displacements relative to the reference MOD-QMM minimum energy configuration. The same methodology could be implemented, more generally, with other types of spectroscopic data (e.g., FTIR, NMR) whenever the experimental data are available so long as the spectroscopic properties can be readily computed. Such SC-QM/MM refinement methods are particularly valuable for the analysis of complex systems with rugged potential energy landscapes, where many configurations have similar energy, including metalloproteins for which the most accurate and reliable information on the structure of prosthetic metal clusters comes from high-resolution spectroscopy.
The reviewed applications addressed fundamental questions that for many years have remained largely elusive to rigorous first-principles examinations, showing that QM/MM hybrid methods where polarization effects are explicitly considered are powerful tools of modern computational chemistry when applied in conjunction with high-resolution spectroscopy. It is, therefore, expected that these QM/MM techniques will continue to make many more important contributions in the characterization of structure and reactivity of biological macromolecules.
V.S.B. acknowledges a generous allocation of DOE supercomputer time from the National Energy Research Scientific Computing Center (NERSC) and financial support from the Grants NSF ECCS # 0404191, NSF CHE # 0345984, DOE DE-FG02-07ER15909, NIH 2R01-GM043278 and the US-Israel BSF R08164. G.W.B. acknowledges support from NIH GM32715.
Eduardo M. Sproviero, Yale University, Department of Chemistry, P. O. Box 208107, New Haven Connecticut 06520-8107 USA.
Michael B. Newcomer, Yale University, Department of Chemistry, P. O. Box 208107, New Haven Connecticut 06520-8107 USA.
José A. Gascón, Yale University, Department of Chemistry, P. O. Box 208107, New Haven Connecticut 06520-8107 USA.
Enrique R. Batista, Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545.
Gary W. Brudvig, Yale University, Department of Chemistry, P. O. Box 208107, New Haven Connecticut 06520-8107 USA.
Victor S. Batista, Yale University, Department of Chemistry, P. O. Box 208107, New Haven Connecticut 06520-8107 USA.