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     

Molecular modeling of green fluorescent protein: structural effects of chromophore deprotonation

Molecular dynamics (MD) simulations were carried out to study the conformational rearrangement induced by deprotonation of the fluorescent chromophore in GFP, as well as the associated changes in the hydrogen-bonding network. For both the structures with either a neutral or an anionic chromophore, it was found that the beta-barrel was stable and rigid, and the conformation of the chromophore was consistent with the available x-ray structure. The conformational change in Thr203 due to deprotonation was also found to be consistent with the three-state isomerization model. Although GFP is highly fluorescent, denatured-GFP is nonfluorescent, indicating that the environment of the protein plays an important role in its fluorescence behavior. Our MD simulations, which explore the effect of the protein shell on the conformation of the chromophore, find the flexibility of the central chromophore to be significantly restricted due to the rigid nature of the protein shell. The hydrogen-bonding between the chromophore and neighboring residues was also shown to contribute to the chromophore rigidity. In addition to the MD studies, quantum mechanics/molecular mechanics (QM/MM) ONIOM calculations were carried out to investigate the effect of the beta-barrel on the internal rotation in the chromophore. Along with providing quantitative values for torsional rotation barriers about the bridging bond in the chromophore, the ONIOM calculations also validate our MD force field parameters.     

QM/MM investigation of structure and spectroscopic properties of a vanadium-containing peroxidase

We present a comprehensive analysis of the most likely ground state configuration of the resting state of vanadium dependent chloroperoxidase (VCPO) based on quantum mechanics/molecular mechanics (QM/MM) evaluations of ground state properties, UV-vis spectra and NMR chemical shifts. Within the QM/MM framework, density functional theory (DFT) calculations are used to characterize the resting state of VCPO via time-dependent density functional theory (TD-DFT) calculations of electronic excitation energies and NMR chemical shifts. Comparison with available experimental data allows us to determine the most likely protonation state of VCPO, a state which results in a doubly protonated axial oxygen, a site largely stabilized by hydrogen bonds. We found that the bulk of the protein that is beyond the immediate layer surrounding the cofactor, has an important electrostatic effect on the absorption maximum. Through examination of frontier orbitals, we analyze the nature of two bound water molecules and the extent to which relevant residues in the active site influence the spectroscopy calculations.     

A computational study of the protein-ligand interactions in CDK2 inhibitors: using quantum mechanics/molecular mechanics interaction energy as a predictor of the biological activity

We report a combined quantum mechanics/molecular mechanics (QM/MM) study to determine the protein-ligand interaction energy between CDK2 (cyclin-dependent kinase 2) and five inhibitors with the N(2)-substituted 6-cyclohexyl-methoxy-purine scaffold. The computational results in this work show that the QM/MM interaction energy is strongly correlated to the biological activity and can be used as a predictor, at least within a family of substrates. A detailed analysis of the protein-ligand structures obtained from molecular dynamics simulations shows specific interactions within the active site that, in some cases, have not been reported before to our knowledge. The computed interaction energy gauges the strength of protein-ligand interactions. Finally, energy decomposition and multiple regression analyses were performed to check the contribution of the electrostatic and van der Waals energies to the total interaction energy and to show the capabilities of the computational model to identify new potent inhibitors.     

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.     

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.     

Vibrational analysis of a solvated green fluorescent protein chromophore

Resonance Raman (RR) spectra of green fluorescent protein (GFP) model chromophores in solution have been simulated with the CASSCF/MM methodology. Although several reports on vibrational analysis of GFP model chromophores have been recently published, the RR spectra were simulated for the first time in explicit solution with the inclusion of the counterion, as these effects are crucial for unambiguously reproducing the vibrational band assignment in the anionic form of the GFP chromophore. This strategy allows for a one-to-one correspondence of the calculated vibrational modes to the observed RR bands, concerning both the location and intensity pattern. In addition, these simulations were complemented with total energy distribution calculations to aid in the unambiguous assignment of the measured spectra. The current study helps to clarify some of the previous RR bands assignments as well as producing some new assignment for the anionic form of GFP chromophore. The explicit solvent simulations and PCM-based calculations are compared to the measured spectra, and these results demonstrate that explicit solvent simulations provide better agreement with experiment, both in terms of vibrational frequencies and intensity distribution. Figure a Correlation of explicit hydration calculations (CASSCF/6-31G*/MM) for the HBI model chromophore and experimental RR data [21]; slope = 0.982, intercept = 27.210 and regression coefficient = 0.997. b Correlation of implicit PCM calculations (CASSCF/6-31G*) for the HBI model chromophore and experimental RR data [21], slope = 1.017, intercept = −48.838 and regression coefficient = 0.984     

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.     

Multiscale methods for macromolecular simulations

In this article we review the key modeling tools available for simulating biomolecular systems. We consider recent developments and representative applications of mixed quantum mechanics/molecular mechanics (QM/MM), elastic network models (ENMs), coarse-grained molecular dynamics, and grid-based tools for calculating interactions between essentially rigid protein assemblies. We consider how the different length scales can be coupled, both in a sequential fashion (e.g. a coarse-grained or grid model using parameterization from MD simulations), and via concurrent approaches, where the calculations are performed together and together control the progression of the simulation. We suggest how the concurrent coupling approach familiar in the context of QM/MM calculations can be generalized, and describe how this has been done in the CHARMM macromolecular simulation package.     

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.     

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.     

The Role of Gln61 in HRas GTP hydrolysis: a quantum mechanics/molecular mechanics study

Activation of the water molecule involved in GTP hydrolysis within the HRas·RasGAP system is analyzed using a tailored approach based on hybrid quantum mechanics/molecular mechanics (QM/MM) simulation. A new path emerges: transfer of a proton from the attacking water molecule to a second water molecule, then a different proton is transferred from this second water molecule to the GTP. Gln(61) will stabilize the transient OH(-) and H(3)O(+) molecules thus generated. This newly proposed mechanism was generated by using, for the first time to our knowledge, the entire HRas-RasGAP protein complex in a QM/MM simulation context. It also offers a rational explanation for previous experimental results regarding the decrease of GTPase rate found in the HRas Q61A mutant and the increase exhibited by the HRas Q61E mutant.     

A quantum mechanics/molecular mechanics study of the reaction mechanism of the hepatitis C virus NS3 protease with the NS5A/5B substrate

Combined quantum mechanics and molecular mechanics (QM/MM) calculations were carried out to characterize the reaction mechanism of the NS3 protease with its preferred substrate (NS5A/5B). The main purpose of this study was to locate the barrier states and intermediates along the distinguished coordinate path (DCP) involved in this process. These structures, and in particular the one corresponding to the first barrier state and intermediate (B1 and I1), could be a starting point for the synthesis of inhibitors of this protease, which could be used to treat hepatitis C. The two first steps of the reaction mechanism were studied, i.e., the acylation step and the breaking of the peptide bond. The first step takes place through a tetracoordinated intermediate, as suggested from previous works on other Serine proteases. The importance of the different amino acid residues was also considered (perturbation study where the MM charges of each residue were set to zero independently). The residues of the oxyanion hole were confirmed as the most important for the electrostatic stabilization of the tetracoordinate intermediate. Moreover, the role of other residues, e.g., Arg-155 and Asp-79, was also explained.     

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.     

Using quantum mechanics to improve estimates of amino acid side chain rotamer energies

Amino acid side chains adopt a discrete set of favorable conformations typically referred to as rotamers. The relative energies of rotamers partially determine which side chain conformations are more often observed in protein structures and accurate estimates of these energies are important for predicting protein structure and designing new proteins. Protein modelers typically calculate side chain rotamer energies by using molecular mechanics (MM) potentials or by converting rotamer probabilities from the protein database (PDB) into relative free energies. One limitation of the knowledge‐based energies is that rotamer preferences observed in the PDB can reflect internal side chain energies as well as longer‐range interactions with the rest of the protein. Here, we test an alternative approach for calculating rotamer energies. We use three different quantum mechanics (QM) methods (second order Møller‐Plesset (MP2), density functional theory (DFT) energy calculation using the B3LYP functional, and Hartree‐Fock) to calculate the energy of amino acid rotamers in a dipeptide model system, and then use these pre‐calculated values in side chain placement simulations. Energies were calculated for over 36,000 different conformations of leucine, isoleucine, and valine dipeptides with backbone torsion angles from the helical and strand regions of the Ramachandran plot. In a subset of cases these energies differ significantly from those calculated with standard molecular mechanics potentials or those derived from PDB statistics. We find that in these cases the energies from the QM methods result in more accurate placement of amino acid side chains in structure prediction tests. Proteins 2008. © 2007 Wiley‐Liss, Inc.     

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.     

Quantum mechanics/molecular mechanics structural models of the oxygen-evolving complex of photosystem II

The annual production of 260 Gtonnes of oxygen, during the process of photosynthesis, sustains life on earth. Oxygen is produced in the thylakoid membranes of green-plant chloroplasts and the internal membranes of cyanobacteria by photocatalytic water oxidation at the oxygen-evolving complex (OEC) of photosystem II (PSII). Recent breakthroughs in X-ray crystallography and advances in quantum mechanics/molecular mechanics (QM/MM) hybrid methods have enabled the construction of chemically sensible models of the OEC of PSII. The resulting computational structural models suggest the complete ligation of the catalytic center by amino acid residues, water, hydroxide and chloride, as determined from the intrinsic electronic properties of the oxomanganese core and the perturbational influence of the surrounding protein environment. These structures are found to be consistent with available mechanistic data, and are also compatible with X-ray diffraction models and extended X-ray absorption fine structure measurements. It is therefore conjectured that these OEC models are particularly relevant for the elucidation of the catalytic mechanism of water oxidation.