Calculation of chromophore excited state energy shifts in response to molecular dynamics of pigment-protein complexes

The absorption and energy transfer properties of photosynthetic pigments are strongly influenced by their local environment or “site.” Local electrostatic fields vary in time with protein and chromophore molecular movement and thus transiently influence the excited state transition properties of individual chromophores. Site-specific information is experimentally inaccessible in many light-harvesting pigment–proteins due to multiple chromophores with overlapping spectra. Full quantum mechanical calculations of each chromophores excited state properties are too computationally demanding to efficiently calculate the changing excitation energies along a molecular dynamics trajectory in a pigment–protein complex. A simplified calculation of electrostatic interactions with each chromophores ground to excited state transition, the so-called charge density coupling (CDC) for site energy, CDC, has previously been developed to address this problem. We compared CDC to more rigorous quantum chemical calculations to determine its accuracy in computing excited state energy shifts and their fluctuations within a molecular dynamics simulation of the bacteriochlorophyll containing light-harvesting Fenna–Mathews–Olson (FMO) protein. In most cases CDC calculations differed from quantum mechanical (QM) calculations in predicting both excited state energy and its fluctuations. The discrepancies arose from the inability of CDC to account for the differing effects of charge on ground and excited state electron orbitals. Results of our study show that QM calculations are indispensible for site energy computations and the quantification of contributions from different parts of the system to the overall site energy shift. We suggest an extension of QM/MM methodology of site energy shift calculations capable of accounting for long-range electrostatic potential contributions from the whole system, including solvent and ions.     

Chromophore Structure of Cyanobacterial Phytochrome Cph1 in the Pr State: Reconciling Structural and Spectroscopic Data by QM/MM Calculations

A quantum mechanics (QM)/molecular mechanics (MM) hybrid method was applied to the Pr state of the cyanobacterial phytochrome Cph1 to calculate the Raman spectra of the bound PCB cofactor. Two QM/MM models were derived from the atomic coordinates of the crystal structure. The models differed in the protonation site of His²⁶⁰ in the chromophore-binding pocket such that either the δ-nitrogen (M-HSD) or the ɛ-nitrogen (M-HSE) carried a hydrogen. The optimized structures of the two models display small differences specifically in the orientation of His²⁶⁰ with respect to the PCB cofactor and the hydrogen bond network at the cofactor-binding site. For both models, the calculated Raman spectra of the cofactor reveal a good overall agreement with the experimental resonance Raman (RR) spectra obtained from Cph1 in the crystalline state and in solution, including Cph1 adducts with isotopically labeled PCB. However, a distinctly better reproduction of important details in the experimental spectra is provided by the M-HSD model, which therefore may represent an improved structure of the cofactor site. Thus, QM/MM calculations of chromoproteins may allow for refining crystal structure models in the chromophore-binding pocket guided by the comparison with experimental RR spectra. Analysis of the calculated and experimental spectra also allowed us to identify and assign the modes that sensitively respond to chromophore-protein interactions. The most pronounced effect was noted for the stretching mode of the methine bridge A-B adjacent to the covalent attachment site of PCB. Due a distinct narrowing of the A-B methine bridge bond angle, this mode undergoes a large frequency upshift as compared with the spectrum obtained by QM calculations for the chromophore in vacuo. This protein-induced distortion of the PCB geometry is the main origin of a previous erroneous interpretation of the RR spectra based on QM calculations of the isolated cofactor.Abbreviations: Agp1, phytochrome from Agrobacterium tumefaciens; α-CPC, α-subunit of C-phycocyanin; BV, biliverdin IXα; B3LYP, three-parameter exchange functional according to Becke, Lee, Yang, and Parr; DFT, density functional theory; DrBphP, phytochrome from Deinococcus radiodurans; GAF, domain found in cGMP-specific phosphodiesterases; MM, molecular mechanics; MD, molecular dynamics; N-H ip, N-H in-plane bending; PCB, phycocyanobilin; PED, potential energy distribution; phyA, plant phytochrome; Pr, Pfr, red- and far-red absorbing parent states of phytochrome; PΦB, phytochromobilin; QM, quantum mechanics; RMSD, root mean-square deviation; RR, resonance Raman     

The reaction mechanism of paraoxon hydrolysis by phosphotriesterase from combined QM/MM simulations

Molecular dynamics simulations employing combined quantum mechanical and molecular mechanical (QM/MM) potentials have been carried out to investigate the reaction mechanism of the hydrolysis of paraoxon by phosphotriesterase (PTE). We used a dual-level QM/MM approach that synthesizes accurate results from high-level electronic structure calculations with computational efficiency of semiempirical QM/MM potentials for free energy simulations. In particular, the intrinsic (gas-phase) energies of the active site in the QM region are determined by using density functional theory (B3LYP) and second-order M?ller-Plesset perturbation theory (MP2) and the molecular dynamics free energy simulations are performed by using the mixed AM1:CHARMM potential. The simulation results suggest a revised mechanism for the phosphotriester hydrolysis mechanism by PTE. The reaction free energy profile is mirrored by structural motions of the binuclear metal center in the active site. The two zinc ions occupy a compact conformation with an average zinc-zinc distance of 3.5 +/- 0.1 A in the Michaelis complex, whereas it is elongated to 5.3 +/- 0.3 A at the transition state and product state. The substrate is loosely bound to the more exposed zinc ion (Znbeta2+) at an average distance of 3.8 A +/- 0.3 A. The P=O bond of the substrate paraoxon is activated by adopting a tight coordination to the Znbeta2+, releasing the coordinate to the bridging hydroxide ion and increasing its nucleophilicity. It was also found that a water molecule enters into the binding pocket of the loosely bound binuclear center, originally occupied by the nucleophilic hydroxide ion. We suggest that the proton of this water molecule is taken up by His254 at low pH or released to the solvent at high pH, resulting in a hydroxide ion that pulls the Znbeta2+ ion closer to form the compact configuration and restores the resting state of the enzyme.     

Density functional theory calculations on entire proteins for free energies of binding: Application to a model polar binding site

In drug optimization calculations, the molecular mechanics Poisson‐Boltzmann surface area (MM‐PBSA) method can be used to compute free energies of binding of ligands to proteins. The method involves the evaluation of the energy of configurations in an implicit solvent model. One source of errors is the force field used, which can potentially lead to large errors due to the restrictions in accuracy imposed by its empirical nature. To assess the effect of the force field on the calculation of binding energies, in this article we use large‐scale density functional theory (DFT) calculations as an alternative method to evaluate the energies of the configurations in a “QM‐PBSA” approach. Our DFT calculations are performed with a near‐complete basis set and a minimal parameter implicit solvent model, within the self‐consistent calculation, using the ONETEP program on protein–ligand complexes containing more than 2600 atoms. We apply this approach to the T4‐lysozyme double mutant L99A/M102Q protein, which is a well‐studied model of a polar binding site, using a set of eight small aromatic ligands. We observe that there is very good correlation between the MM and QM binding energies in vacuum but less so in the solvent. The relative binding free energies from DFT are more accurate than the ones from the MM calculations, and give markedly better agreement with experiment for six of the eight ligands. Furthermore, in contrast to MM‐PBSA, QM‐PBSA is able to correctly predict a nonbinder. Proteins 2014; 82:3335–3346. © 2014 Wiley Periodicals, Inc.     

Combined quantum and molecular mechanics calculations on metalloproteins

The combination of quantum mechanics and molecular mechanics (QM/MM) methods is one of the most promising approaches to study the structure, function and properties of proteins. The number of QM/MM applications on metalloproteins is steadily increasing, especially studies with density functional methods on redox-active metal centres. Recent developments include new parameterised methods to treat covalent bonds between the quantum and classical systems, methods to obtain free energy from QM/MM results, and the combination of quantum chemistry and protein crystallography.     

QM/MM theoretical study of the pentacoordinate Mn(III) and resting states of manganese-reconstituted cytochrome P450<Subscript>cam</Subscript>

Quantum mechanical/molecular mechanical (QM/MM) theoretical calculations were performed for the pentacoordinate Mn(III) and water-bound resting states of the Mn-reconstituted mutant of cytochrome P450_cam (Mn-P450_cam) in order to obtain insights into their characters, especially, their spin state ordering. The QM/MM study was carried out by use of the B3LYP and BLYP density functional theory (DFT) methods coupled to the CHARMM force field. Although the relative energies of possible spin states for the Mn-P450_cam species varied depending on the functional, this dependence was less significant compared with previous calculations on the corresponding intermediates of wild-type P450_cam. The results suggested that both Mn-P450_cam intermediates have quintet ground states. Additional time-dependent DFT (TDDFT) calculations were carried out for the quintet states of these species using the B3LYP and BP86 functionals with the electrostatic environmental effect included. The TDDFT results enabled us to assign the origins of the peaks observed in optical absorption spectra (Makris et al. in J. Inorg. Biochem. 100:507–518, 2006). Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Determinants of reactivity and selectivity in soluble epoxide hydrolase from quantum mechanics/molecular mechanics modeling

Soluble epoxide hydrolase (sEH) is an enzyme involved in drug metabolism that catalyzes the hydrolysis of epoxides to form their corresponding diols. sEH has a broad substrate range and shows high regio- and enantioselectivity for nucleophilic ring opening by Asp333. Epoxide hydrolases therefore have potential synthetic applications. We have used combined quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics (MD) simulations (at the AM1/CHARMM22 level) and high-level ab initio (SCS-MP2) QM/MM calculations to analyze the reactions, and determinants of selectivity, for two substrates: trans-stilbene oxide (t-SO) and trans-diphenylpropene oxide (t-DPPO). The calculated free energy barriers from the QM/MM (AM1/CHARMM22) umbrella sampling MD simulations show a lower barrier for phenyl attack in t-DPPO, compared with that for benzylic attack, in agreement with experiment. Activation barriers in agreement with experimental rate constants are obtained only with the highest level of QM theory (SCS-MP2) used. Our results show that the selectivity of the ring-opening reaction is influenced by several factors, including proximity to the nucleophile, electronic stabilization of the transition state, and hydrogen bonding to two active site tyrosine residues. The protonation state of His523 during nucleophilic attack has also been investigated, and our results show that the protonated form is most consistent with experimental findings. The work presented here illustrates how determinants of selectivity can be identified from QM/MM simulations. These insights may also provide useful information for the design of novel catalysts for use in the synthesis of enantiopure compounds.     

How does the protein environment optimize the thermodynamics of thiol sulfenylation? Insights from model systems to QM/MM calculations on human 2-Cys peroxiredoxin

Protein thiol/sulfenic acid oxidation potentials provide a tool to select specific oxidation agents, but are experimentally difficult to obtain. Here, insights into the thiol sulfenylation thermodynamics are obtained from model calculations on small systems and from a quantum mechanics/molecular mechanics (QM/MM) analysis on human 2-Cys peroxiredoxin thioredoxin peroxidase B (Tpx-B). To study thiol sulfenylation in Tpx-B, our recently developed computational method to determine reduction potentials relatively compared to a reference system and based on reaction energies reduction potential from electronic energies is updated. Tpx-B forms a sulfenic acid (R-SO^?) on one of its active site cysteines during reactive oxygen scavenging. The observed effect of the conserved active site residues is consistent with the observed hydrogen bond interactions in the QM/MM optimized Tpx-B structures and with free energy calculations on small model systems. The ligand effect could be linked to the complexation energies of ligand L with CH₃S^? and CH₃SO^?. Compared to QM only calculations on Tpx-B’s active site, the QM/MM calculations give an improved understanding of sulfenylation thermodynamics by showing that other residues from the protein environment other than the active site residues can play an important role.     

QM/MM model study on properties and structure of some antibiotics in gas phase: Comparison of energy and NMR chemical shift

The combination of Quantum Mechanics (QM) and Molecular Mechanics (MM) methods has become an alternative tool for many applications for which pure QM and MM are not suitable. The QM/MM method has been used for different types of problems, for example: structural biology, surface phenomena, and liquid phase. In this paper, we have used these methods for antibiotics and then we compare results. The calculations were done by the full ab initio method (HF/3-21G) and the (HF/STO-3G) and QM/MM (ONIOM) method with HF (3-21G)/AM1/UFF and HF (STO-3G)/AM1/UFF. We found the geometry that has obtained by the QM/MM method to be very accurate, and we can use this rapid method in place of time consuming ab initio methods for large molecules. Comparison of energy values in the QM/MM and QM methods is given. In the present work, we compare chemical shifts and conclude that the QM/MM method is a perturbed full QM method. The work has been done on penicillin, streptomycin, benzyl penicillin, neomycin, kanamycin, gentamicin, and amoxicillin.     

Recent developments in QM/MM methods towards open-boundary multi-scale simulations

Combined quantum mechanics/molecular mechanics (QM/MM) methods have been widely used in multi-scale modelling and simulations of physical, chemical and biological processes in complex environments. In this review, we provide an overview of the recently developed QM/MM algorithms, with emphasis on our works, towards the ultimate goal of establishing an open boundary between the QM and MM subsystems. The open boundary is characterised by on-the-fly exchanges of partial charges and atoms between the QM and MM subsystems, allowing us to focus on the small QM subsystem of primary interest in dynamics simulations. An open-boundary scheme has the promise to the utilisations of small QM subsystems, high-levels of QM theory and long simulation times, which can potentially lead to new insights.     

Low-barrier hydrogen bond hypothesis in the catalytic triad residue of serine proteases: correlation between structural rearrangement and chemical shifts in the acylation process

To elucidate the catalytic advantage of the low-barrier hydrogen bond (LBHB), we analyze the hydrogen bonding network of the catalytic triad (His57-Asp102-Ser195) of serine protease trypsin, one of the best examples of the LBHB reaction mechanism. Especially, we focus on the correlation between the change of the chemical shifts and the structural rearrangement of the active site in the acylation process. To clarify LBHB, we evaluate the two complementary properties. First, we calculate the NMR chemical shifts of the imidazole ring of His57 by the gauge-including atomic orbital (GIAO) approach within the ab initio QM/MM framework. Second, the free energy profile of the proton transfer from His57 to Asp102 in the tetrahedral intermediate is obtained by ab initio QM/MM calculations combined with molecular dynamics free energy perturbation (MD-FEP) simulations. The present analyses reveal that the calculated shifts reasonably reproduce the observed values for (1)H chemical shift of H(epsilon)(1) and H(delta)(1) in His57. The (15)N and (13)C chemical shifts are also consistent with the experiments. It is also shown that the proton between His57 and Asp102 is localized at the His57 side. This largely downfield chemical shift is originated from the strong electrostatic interaction, not a covalent-like bonding character between His57 and Asp102. Also, it is proved that a slight downfield character of H(epsilon)(1) is originated from a electrostatic interaction between His57 and the backbone carbonyl group of Val213 and Ser214. These downfield chemical shifts are observed only when the tetrahedral intermediate is formed in the acylation process.     

Application of a SCC-DFTB QM/MM approach to the investigation of the catalytic mechanism of fatty acid amide hydrolase

Self-consistent charge density functional tight binding (SCC-DFTB) is a promising method for hybrid quantum mechanics/molecular mechanics (QM/MM) simulations of enzyme-catalyzed reactions. The acylation reaction of fatty acid amide hydrolase (FAAH), a promising drug target, was investigated by applying a SCC-DFTB/CHARMM27 scheme. Calculated potential energy barriers resulted in reasonable agreement with experiments for oleamide (OA) and oleoylmethyl ester (OME) substrates, outperforming previous calculations performed at the PM3/CHARMM22 level. Furthermore, the experimental preference of FAAH in hydrolyzing OA faster than OME was adequately reproduced by calculations. All these findings indicate that the SCC-DFTB/CHARMM27 approach can be successfully applied to mechanistic investigations of FAAH-catalyzed reactions.     

Binding free energy calculation with QM/MM hybrid methods for Abl-Kinase inhibitor

We report a Quantum mechanics/Molecular Mechanics–Poisson-Boltzmann/ Surface Area (QM/MM-PB/SA) method to calculate the binding free energy of c-Abl human tyrosine kinase by combining the QM and MM principles where the ligand is treated quantum mechanically and the rest of the receptor by classical molecular mechanics. To study the role of entropy and the flexibility of the protein ligand complex in a solvated environment, molecular dynamics calculations are performed using a hybrid QM/MM approach. This work shows that the results of the QM/MM approach are strongly correlated with the binding affinity. The QM/MM interaction energy in our reported study confirms the importance of electronic and polarization contributions, which are often neglected in classical MM-PB/SA calculations. Moreover, a comparison of semi-empirical methods like DFTB-SCC, PM3, MNDO, MNDO-PDDG, and PDDG-PM3 is also performed. The results of the study show that the implementation of a DFTB-SCC semi-empirical Hamiltonian that is derived from DFT gives better results than other methods. We have performed such studies using the AMBER molecular dynamic package for the first time. The calculated binding free energy is also in agreement with the experimentally determined binding affinity for c-Abl tyrosine kinase complex with Imatinib.     

QM/MM model study on properties and structure of some antibiotics in gas phase: Comparison of energy and NMR chemical shift

The combination of Quantum Mechanics (QM) and Molecular Mechanics (MM) methods has become an alternative tool for many applications for which pure QM and MM are not suitable. The QM/MM method has been used for different types of problems, for example: structural biology, surface phenomena, and liquid phase. In this paper, we have used these methods for antibiotics and then we compared results. The calculations were done by the full ab initio method (HF/3-21G) and the (HF/STO-3G) and QM/MM (ONIOM) method with HF (3-21G)/AM1/UFF and HF (STO-3G)/AM1/UFF. We found the geometry obtained by the QM/MM method to be very accurate, and we can use this rapid method in place of time consuming ab initio methods for large molecules. Comparison of energy values in the QM/MM and QM methods is given. In the present work, we compare chemical shifts and conclude that the QM/MM method is a perturbed full QM method. The work has been done on penicillin, streptomycin, benzyl penicillin, neomycin, kanamycin, gentamicin, and amoxicillin.     

Validation of semiempirical methods for modeling of corrinoid systems

Several semiempirical methods (MNDO-d, PM3tm, PM3-d, PM5, PM6, and AM1-d) have been tested against experimental data and density functional theory (DFT) results in search for the best methods that can be used for quantum-mechanical-molecular mechanics (QM/MM) modeling of corrinoid systems of vitamin B(12) co-factor. It has been found that the PM6 parametrization in its present form gives results closest to hybrid DFT calculations that are most widely used thus far. In comparison with pure DFT and experimental data the best agreement is obtained for PM3tm parametrization, while PM6 yields slightly worse results. AM1-d yields bad geometry of the corrin moiety. The worst performance was observed for MNDO-d, which has severe problem with position and orientation of the alpha-ligands.     

Modeling reaction routes from rhodopsin to bathorhodopsin

The quantum mechanical–molecular mechanical (QM/MM) theory was applied to calculate accurate structural parameters, vibrational and optical spectra of bathorhodopsin (BATHO), one of the primary photoproducts of the functional cycle of the visual pigment rhodopsin (RHO), and to characterize reaction routes from RHO to BATHO. The recently resolved crystal structure of BATHO (PDBID: 2G87) served as an initial source of coordinates of heavy atoms. Protein structures in the ground electronic state and vibrational frequencies were determined by using the density functional theory in the PBE0/cc‐pVDZ approximation for the QM part and the AMBER force field parameters in the MM part. Calculated and assigned vibrational spectra of both model protein systems, BATHO and RHO, cover three main regions referring to the hydrogen‐out‐of‐plan (HOOP) motion, the C?C ethylenic stretches, and the C? C single‐bond stretches. The S₀–S₁ electronic excitation energies of the QM part, including the chromophore group in the field of the protein matrix, were estimated by using the advanced quantum chemistry methods. The computed structural parameters as well as the spectral bands match perfectly the experimental findings. A structure of the transition state on the S₀ potential energy surface for the ground electronic state rearrangement from RHO to BATHO was located proving a possible route of the thermal protein activation to the primary photoproduct. Proteins 2010. © 2009 Wiley‐Liss, Inc.     

Comparison of a QM/MM force field and molecular mechanics force fields in simulations of alanine and glycine "dipeptides" (Ace-Ala-Nme and Ace-Gly-Nme) in water in relation to the problem of modeling the unfolded peptide backbone in solution

We compare the conformational distributions of Ace-Ala-Nme and Ace-Gly-Nme sampled in long simulations with several molecular mechanics (MM) force fields and with a fast combined quantum mechanics/molecular mechanics (QM/MM) force field, in which the solute's intramolecular energy and forces are calculated with the self-consistent charge density functional tight binding method (SCCDFTB), and the solvent is represented by either one of the well-known SPC and TIP3P models. All MM force fields give two main states for Ace-Ala-Nme, beta and alpha separated by free energy barriers, but the ratio in which these are sampled varies by a factor of 30, from a high in favor of beta of 6 to a low of 1/5. The frequency of transitions between states is particularly low with the amber and charmm force fields, for which the distributions are noticeably narrower, and the energy barriers between states higher. The lower of the two barriers lies between alpha and beta at values of psi near 0 for all MM simulations except for charmm22. The results of the QM/MM simulations vary less with the choice of MM force field; the ratio beta/alpha varies between 1.5 and 2.2, the easy pass lies at psi near 0, and transitions between states are more frequent than for amber and charmm, but less frequent than for cedar. For Ace-Gly-Nme, all force fields locate a diffuse stable region around phi = pi and psi = pi, whereas the amber force field gives two additional densely sampled states near phi = +/-100 degrees and psi = 0, which are also found with the QM/MM force field. For both solutes, the distribution from the QM/MM simulation shows greater similarity with the distribution in high-resolution protein structures than is the case for any of the MM simulations.