A Poisson−Boltzmann Study of Charge Insertion in an Enzyme Active Site: The Effect of Dielectric Relaxation

Deciphering the selectivity of inhibitor MKC9989 towards residue K907 in IRE1α; a multiscale in silico approach

The selectivity of the ligand MKC9989, as inhibitor of the Inositol-Requiring Enzyme 1α (IRE1α) transmembrane kinase/ribonuclease protein, towards the residue K907 in the context of Schiff base formation, has been investigated by employing an array of in silico techniques including Multi-Conformation Continuum Electrostatics (MCCE) simulations, Quantum Mechanics/Molecular Mechanics (QM/MM) calculations, covalent docking, and Molecular Dynamics (MD) simulations. According to the MCCE results, K907 displays the lowest pK_a value among all 23 lysine residues in IRE1α. The MMCE simulations also indicate a critical interaction between K907 and D885 within the hydrophobic pocket which increases significantly at low protein dielectric constants. The QM/MM calculations reveal a spontaneous proton transfer from K907 to D885, consistent with the low pK_a value of K907. A Potential Energy Surface (PES) scan confirms the lack of energy barrier and transition state associated with this proton transfer reaction. Covalent docking and MD simulations verify that the protein pocket containing K907 can effectively stabilize the inhibitor by strong π–π and hydrogen bonding interactions. In addition, Radial Distribution Function (RDF) analysis shows that the imine group formed in the chemical reaction between MKC9989 and K907 is inaccessible to water molecules and thus the probability of imine hydrolysis is almost zero. The results of the current study explain the high selectivity of the MKC9989 inhibitor towards the K907 residue of IRE1α.

The high selectivity of inhibitor MKC9989 towards Lys907 of IRE1α is explained by the unique pK_a properties of the lysine.

Comparing pairwise-additive and many-body generalized Born models for acid/base calculations and protein design

Generalized Born (GB) solvent models are common in acid/base calculations and protein design. With GB, the interaction between a pair of solute atoms depends on the shape of the protein/solvent boundary and, therefore, the positions of all solute atoms, so that GB is a many-body potential. For compute-intensive applications, the model is often simplified further, by introducing a mean, native-like protein/solvent boundary, which removes the many-body property. We investigate a method for both acid/base calculations and protein design that uses Monte Carlo simulations in which side chains can explore rotamers, bind/release protons, or mutate. The fluctuating protein/solvent dielectric boundary is treated in a way that is numerically exact (within the GB framework), in contrast to a mean boundary. Its originality is that it captures the many-body character while retaining the residue-pairwise complexity given by a fixed boundary. The method is implemented in the Proteus protein design software. It yields a slight but systematic improvement for acid/base constants in nine proteins and a significant improvement for the computational design of three PDZ domains. It eliminates a source of model uncertainty, which will facilitate the analysis of other model limitations. © 2017 Wiley Periodicals, Inc.

Dielectric Allostery of Protein: Response of Myosin to ATP Binding

Protein uses allostery to execute biological function. The physical mechanism underlying the allostery has long been studied, with the focus on the mechanical response by ligand binding. Here, we highlight the electrostatic response, presenting an idea of "dielectric allostery". We conducted molecular dynamics simulations of myosin, a motor protein with allostery, and analyzed the response to ATP binding which is a crucial step in force-generating function, forcing myosin to unbind from the actin filament. We found that the net negative charge of ATP causes a large-scale, anisotropic dielectric response in myosin, altering the electrostatic potential in the distant actin-binding region and accordingly retracting a positively charged actin-binding loop. A large-scale rearrangement of electrostatic bond network was found to occur upon ATP binding. Since proteins are dielectric and ligands are charged/polar in general, the dielectric allostery might underlie a wide spectrum of functions by proteins.

The MOD-QM/MM Method: Applications to Studies of Photosystem II and DNA G-Quadruplexes

Quantum mechanics/molecular mechanics (QM/MM) hybrid methods are currently the most powerful computational tools for studies of structure/function relations and catalytic sites embedded in macrobiomolecules (eg, proteins and nucleic acids). QM/MM methodologies are highly efficient since they implement quantum chemistry methods for modeling only the portion of the system involving bond-breaking/forming processes (QM layer), as influenced by the surrounding molecular environment described in terms of molecular mechanics force fields (MM layer). Some of the limitations of QM/MM methods when polarization effects are not explicitly considered include the approximate treatment of electrostatic interactions between QM and MM layers. Here, we review recent advances in the development of computational protocols that allow for rigorous modeling of electrostatic interactions in biomacromolecules and structural refinement, beyond the common limitations of QM/MM hybrid methods. We focus on photosystem II (PSII) with emphasis on the description of the oxygen-evolving complex (OEC) and its high-resolution extended X-ray absorption fine structure spectra (EXAFS) in conjunction with Monte Carlo structural refinement. Furthermore, we review QM/MM structural refinement studies of DNA G4 quadruplexes with embedded monovalent cations and direct comparisons to NMR data.

An Overview of Electrostatic Free Energy Computations for Solutions and Proteins

Free energy simulations for electrostatic and charging processes in complex molecular systems encounter specific difficulties owing to the long-range, 1/r Coulomb interaction. To calculate the solvation free energy of a simple ion, it is essential to take into account the polarization of nearby solvent but also the electrostatic potential drop across the liquid-gas boundary, however distant. The latter does not exist in a simulation model based on periodic boundary conditions because there is no physical boundary to the system. An important consequence is that the reference value of the electrostatic potential is not an ion in a vacuum. Also, in an infinite system, the electrostatic potential felt by a perturbing charge is conditionally convergent and dependent on the choice of computational conventions. Furthermore, with Ewald lattice summation and tinfoil conducting boundary conditions, the charges experience a spurious shift in the potential that depends on the details of the simulation system such as the volume fraction occupied by the solvent. All these issues can be handled with established computational protocols, as reviewed here and illustrated for several small ions and three solvated proteins.

Understanding photosynthetic light-harvesting: a bottom up theoretical approach

We discuss a bottom up approach for modeling photosynthetic light-harvesting. Methods are reviewed for a full structure-based parameterization of the Hamiltonian of pigment-protein complexes (PPCs). These parameters comprise (i) the local transition energies of the pigments in their binding sites in the protein, the site energies; (ii) the couplings between optical transitions of the pigments, the excitonic couplings; and (iii) the spectral density characterizing the dynamic modulation of pigment transition energies and excitonic couplings by protein vibrations. Starting with quantum mechanics perturbation theory, we provide a microscopic foundation for the standard PPC Hamiltonian and relate the expressions obtained for its matrix elements to quantities that can be calculated with classical molecular mechanics/electrostatics approaches including the whole PPC in atomic detail and using charge and transition densities obtained with quantum chemical calculations on the isolated building blocks of the PPC. In the second part of this perspective, the Hamiltonian is utilized to describe the quantum dynamics of excitons. Situations are discussed that differ in the relative strength of excitonic and exciton-vibrational coupling. The predictive power of the approaches is demonstrated in application to different PPCs, and challenges for future work are outlined.

Constant pH molecular dynamics (CpHMD) and molecular docking studies of CquiOBP1 pH-induced ligand releasing mechanism

The odorant binding protein of Culex quinquefasciatus (CquiOBP1), expressed on the insect antenna, is crucial for the investigation of trapping baited with oviposition semi-chemicals and controlling mosquito populations. The acidic titratable residues pKa prediction and the ligand binding poses investigation in two systems (pH 7 and pH 5) are studied by constant pH molecular dynamics (CpHMD) and molecular docking methods. Research results reveal that the change of the protonation states would disrupt some important H-bonds, such as Asp 66-Asp 70, Glu 105-Asn 102, etc. The cleavage of these H-bonds leads to the movement of the relative position of hydrophobic tunnel, N- and C- termini loops and pH-sensing triad (His23-Tyr54-Val125) in acid solution. Ligand MOP has lower affinity and shows different binding poses to protein CquiOBP1 at pH 5. This ligand may be released from another tunnel between helices α3 and α4 in acidic environment. However, it would bind to the protein with high affinity in neutral environment. This work could provide more penetrating understanding of the pH-induced ligand-releasing mechanism.

Nucleotide recognition by the initiation factor aIF5B: free energy simulations of a neoclassical GTPase

The GTPase aIF5B is a universally conserved initiation factor that assists ribosome assembly. Crystal structures of its nucleotide complexes, X-ray(GTP) and X-ray(GDP), are similar in the nucleotide vicinity, but differ in the orientation of a distant domain IV. This has led to two, contradictory, mechanistic models. One postulates that X-ray(GTP) and X-ray(GDP) are, respectively, the active, "ON" and the inactive, "OFF" states; the other postulates that both structures are OFF, whereas the ON state is still uncharacterized. We study GTP/GDP binding using molecular dynamics and a continuum electrostatic free energy method. We predict that X-ray(GTP) has a ≈ 3 kcal/mol preference to bind GDP, apparently contradicting its assignment as ON. However, the preference arises mainly from a single, nearby residue from the switch 2 motif: Glu81, which becomes protonated upon GTP binding, with a free energy cost of about 4 kcal/mol. We then propose a different model, where Glu81 protonation/deprotonation defines the ON/OFF states. With this model, the X-ray(GTP):GTP complex, with its protonated Glu81, is ON, whereas X-ray(GTP):GDP is OFF. The model postulates that distant conformational changes such as domain IV rotation are "uncoupled" from GTP/GDP exchange and do not affect the relative GTP/GDP binding affinities. We analyze the model using a general thermodynamic framework for GTPases. It yields rather precise predictions for the nucleotide specificities of each state, and the state specificities of each nucleotide, which are roughly comparable to the homologues IF2 and aIF2, despite the lack of any conformational switching in the model.

QM/MM Prediction of the Stark Shift in the Active Site of a Protein

Recent developments in the biophysical characterization of proteins have provided a means of directly measuring electrostatic fields by introducing a probe molecule to the system of interest and interpreting photon absorption in the context of the Stark effect. To fully account for this effect, the development of accurate atomistic models is of paramount importance. However, suitable computational protocols for evaluating Stark shifts in proteins are yet to be established. In this work, we present a comprehensive computational method to predict the change in absorption frequency of a probe functional group as a direct result of a perturbation in its surrounding electrostatic field created by a protein environment, i.e., the Stark shift. We apply the method to human aldose reductase, a key protein enzyme that catalyzes the reduction of monosaccharides. We develop a protocol based on a combination of molecular dynamics and moving-domain QM/MM methods, which achieves quantitative agreement with experiment. We outline the difficulties in predicting localized electrostatic field changes within a protein environment, and by extension the Stark shift, due to a protein site mutation. Furthermore, the combined use of Stark effect spectroscopy and computational modeling is used to predict the protonation state of ionizable residues in the vicinity of the electrostatic probe.

Constant pH molecular dynamics (CpHMD) and mutation studies: insights into AaegOBP1 pH-induced ligand releasing mechanism

AaegOBP1, isolated from the male and female antenna of yellow fever mosquitoes, may serve as crucial molecular targets for the development of mosquitoes' attractants and for the control of mosquito populations. Nowadays crystal structures of AaegOBP1 in the neutral environment have been obtained, whereas little research is focused on the conformational change of AaegOBP1 in the acid solution. In our study, the conformational change and the ligand bound poses in different solution pH were investigated using constant pH molecular dynamics (CpHMD) as well as mutation studies. Results demonstrate that the protein changes dramatically in low pH solution and that the pH-sensing triad (Arg23-Tyr54-Ile125) plays an indispensable role in the structural stability and ligand binding. In addition, we have proved that the residue Arg23 is the most important one of the pH-sensing triad. This work could provide more penetrating understanding of the pH-induced ligand-releasing mechanism.

Effects of two solvent conditions on the free energy landscape of the BBL peripheral subunit binding domain

BBL is a small independently folding domain with two main parallel helices. The experiment of C(α) secondary shifts has shown that changing the pH from ~7 to ~5 results in the reduced helicity at the C-terminus of helix 2. Combining constant pH molecular dynamics with replica exchange, we sampled the protein conformation space and protonation states extensively under a neutral pH condition and an acidic condition. Our results reveal that the solvent conditions influence the free energy landscape. Under the neutral pH condition, the denatured state and the native state are well separated. The condition of the acidic pH reshapes the free energy surface, leading to a broadly populated denatured-state basin and a low free energy barrier between the denatured state and the native state. The acidic pH shifts the equilibrium between the denatured state and the native state in favor of the denatured state. Caution must be used to interpret experimental data under the acidic condition because the contribution of the denatured state is significant. Our simulation results are supported by the fact that the calculated chemical shifts are in good agreement with the experiment data.

An extrapolation method for computing protein solvation energies based on density fragmentation of a graphical surface tessellation

Modeling chemical events inside proteins often require the incorporation of solvent effects via continuum polarizable models. One of these approaches is based on the assumption that the interface between solute and solvent acts as a conductor. Image charges are added on the molecular surface to satisfy the appropriate conductor boundary conditions in the presence of solute charges. As in the case of other polarizable continuum models that are based on surface tessellation, the simplest implementation of this approach is often limited to several hundred atoms due to a matrix inversion, which scales as the cube of the number or tesserae. For larger systems, approaches that use iterative matrix solvers coupled to fast summation methods must be used. In the present work, we develop a self-consistent approach to obtain conductor-like screening charges suitable for applications in proteins. The approach is based on a density fragmentation of a graphical surface tessellation. This method, although approximate, provides a straightforward scheme of parallelization, which can in principle be added to existing linear scaling implementations of conductor-like models. We implement this method in conjunction with a fixed charge model for the protein, as well as with a moving domain QM/MM description of the protein. In the latter case, the overall result leads to a charge distribution within the protein determined by self-polarization and polarization due to solvent.

Virtual screening using molecular simulations

Effective virtual screening relies on our ability to make accurate prediction of protein-ligand binding, which remains a great challenge. In this work, utilizing the molecular-mechanics Poisson-Boltzmann (or Generalized Born) surface area approach, we have evaluated the binding affinity of a set of 156 ligands to seven families of proteins, trypsin β, thrombin α, cyclin-dependent kinase (CDK), cAMP-dependent kinase (PKA), urokinase-type plasminogen activator, β-glucosidase A, and coagulation factor Xa. The effect of protein dielectric constant in the implicit-solvent model on the binding free energy calculation is shown to be important. The statistical correlations between the binding energy calculated from the implicit-solvent approach and experimental free energy are in the range of 0.56-0.79 across all the families. This performance is better than that of typical docking programs especially given that the latter is directly trained using known binding data whereas the molecular mechanics is based on general physical parameters. Estimation of entropic contribution remains the barrier to accurate free energy calculation. We show that the traditional rigid rotor harmonic oscillator approximation is unable to improve the binding free energy prediction. Inclusion of conformational restriction seems to be promising but requires further investigation. On the other hand, our preliminary study suggests that implicit-solvent based alchemical perturbation, which offers explicit sampling of configuration entropy, can be a viable approach to significantly improve the prediction of binding free energy. Overall, the molecular mechanics approach has the potential for medium to high-throughput computational drug discovery.

Alchemical free energy simulations for biological complexes: powerful but temperamental

Free energy simulations compare multiple ligand:receptor complexes by "alchemically" transforming one into another, yielding binding free energy differences. Since their introduction in the 1980s, many technical and theoretical obstacles were surmounted, and the method ("MDFE," since molecular dynamics are often used) has matured into a powerful tool. We describe its current status, its effectiveness, and the challenges it faces. MDFE has provided chemical accuracy for many systems but remains expensive, with significant human overhead costs. The bottlenecks have shifted, partly due to increased computer power. To study diverse sets of ligands, force field availability and accuracy can be a major difficulty. Another difficulty is the frequent need to consider multiple states, related to sidechain protonation or buried waters, for example. Sophisticated, automated methods to sample these states are maturing, such as constant pH simulations. Meanwhile, combinations of MDFE and simpler approaches, like continuum dielectric models, can be very effective. As illustrations, we show how, with careful force field parameterization, MDFE accurately predicts binding specificities between complex tetracycline ligands and their targets. We describe substrate binding to the aspartyl-tRNA synthetase enzyme, where many distinct electrostatic states play a role, and a histidine and a Mg(2+) ion act as coupled switches that help enforce a strict preference for the aspartate substrate, relative to several analogs. Overall, MDFE has achieved a predictive status, where novel ligands can be studied and molecular recognition elucidated in depth. It should play an increasing role in the analysis of complex cellular processes and biomolecular engineering.

Constant pH replica exchange molecular dynamics in biomolecules using a discrete protonation model

A constant pH replica exchange molecular dynamics (REMD) method is proposed and implemented to improve coupled protonation and conformational state sampling. By mixing conformational sampling at constant pH (with discrete protonation states) with a temperature ladder, this method avoids conformational trapping. Our method was tested and applied to seven different biological systems. The constant pH REMD not only predicted pKa correctly for small, model compounds but also converged faster than constant pH molecular dynamics (MD). We further tested our constant pH REMD on a heptapeptide from ovomucoid third domain (OMTKY3). Although constant pH REMD and MD produced very close pKa values, the constant pH REMD showed its advantage in the efficiency of conformational and protonation state samplings.

Redox potential of the Rieske iron-sulfur protein quantum-chemical and electrostatic study

Quantum-chemical study of structures, energies, and effective partial charge distribution for several models of the Rieske protein redox center is performed in terms of the B3LYP density functional method in combination with the broken symmetry approach using three different atomic basis sets. The structure of the redox complex optimized in vacuum differs markedly from that inside the protein. This means that the protein matrix imposes some stress on the active site resulting in distortion of its structure. The redox potentials calculated for the real active site structure are in a substantially better agreement with the experiment than those calculated for the idealized structure. This shows an important role of the active site distortion in tuning its redox potential. The reference absolute electrode potential of the standard hydrogen electrode is used that accounts for the correction caused by the water surface potential. Electrostatic calculations are performed in the framework of the polarizable solute model. Two dielectric permittivities of the protein are employed: the optical permittivity for calculation of the intraprotein electric field, and the static permittivity for calculation of the dielectric response energy. Only this approach results in a reasonable agreement of the calculated and experimental redox potentials.

A single molecule as a dielectric medium

For three molecules with weak or negligible charge overlap, we prove that the three-body interaction energy obtained from quantum perturbation theory (to leading order) fits a dielectric model with a nonlocal electronic screening function. The electronic charge cloud of each molecule acts as a dielectric medium for the interaction of the remaining two with the nonlocal dielectric function epsilon(r,r') obtained by O. S. Jenkins and K. L. C. Hunt [J. Chem. Phys. 119, 8250 (2003)], by considering the charge redistribution induced in a single molecule by an external perturbation. The dielectric function depends parametrically on the coordinates of the nuclei, within the Born-Oppenheimer approximation. We also prove that the force on each nucleus in molecule A depends on intramolecular dielectric screening within A. The potential from the charge distribution of B, screened by C acting as a dielectric medium, is further screened linearly within A; and similarly, with the roles of B and C reversed. In addition, the potential due to the unperturbed charge distribution of B and the potential due to the unperturbed charge distribution of C, acting simultaneously, are screened nonlinearly within A. The results show that nonlocal dielectric theory holds on the molecular level, provided that the overlap of the electronic charge distributions is weak.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Dielectric relaxation of cytochrome c oxidase: Comparison of the microscopic and continuum models

We have studied a charge-insertion process that models the deprotonation of a histidine side chain in the active site of cytochrome c oxidase (CcO) using both the continuum electrostatic calculations and the microscopic simulations. The group of interest is a ligand to Cu(B) center of CcO, which has been previously suggested to play the role of the proton pumping element in the enzyme; the group is located near a large internal water cavity in the protein. Using the nonpolarizable Amber-99 force field in molecular dynamics (MD) simulations, we have calculated the nuclear part of the reaction-field energy of charging of the His group and combined it with the electronic part, which we estimated in terms of the electronic continuum (EC) model, to obtain the total reaction-field energy of charging. The total free energy obtained in this MDEC approach was then compared with that calculated using pure continuum electrostatic model with variable dielectric parameters. The dielectric constant for the "dry" protein and that of the internal water cavity of CcO were determined as those parameters that provide best agreement between the continuum and microscopic MDEC model. The nuclear (MD) polarization alone (without electronic part) of a dry protein was found to correspond to an unphysically low dielectric constant of only about 1.3, whereas the inclusion of electronic polarizability increases the protein dielectric constant to 2.6-2.8. A detailed analysis is presented as to how the protein structure should be selected for the continuum calculations, as well as which probe and atomic radii should be used for cavity definition. The dielectric constant of the internal water cavity was found to be 80 or even higher using "standard" parameters of water probe radius, 1.4 A, and protein atomic radii from the MD force field for cavity description; such high values are ascribed to the fact that the standard procedure produces unphysically small cavities. Using x-ray data for internal water in CcO, we have explored optimization of the parameters and the algorithm of cavity description. For Amber radii, the optimal probe size was found to be 1.25 A; the dielectric of water cavity in this case is in the range of 10-16. The most satisfactory cavity description, however, was achieved with ProtOr atomic radii, while keeping the probe radius to be standard 1.4 A. In this case, the value of cavity dielectric constant was found to be in the range of 3-6. The obtained results are discussed in the context of recent calculations and experimental measurements of dielectric properties of proteins.

Accurate Molecular Polarizabilities Based on Continuum Electrostatics

A novel approach for representing the intramolecular polarizability as a continuum dielectric is introduced to account for molecular electronic polarization. It is shown, using a finite-difference solution to the Poisson equation, that the Electronic Polarization from Internal Continuum (EPIC) model yields accurate gas-phase molecular polarizability tensors for a test set of 98 challenging molecules composed of heteroaromatics, alkanes and diatomics. The electronic polarization originates from a high intramolecular dielectric that produces polarizabilities consistent with B3LYP/aug-cc-pVTZ and experimental values when surrounded by vacuum dielectric. In contrast to other approaches to model electronic polarization, this simple model avoids the polarizability catastrophe and accurately calculates molecular anisotropy with the use of very few fitted parameters and without resorting to auxiliary sites or anisotropic atomic centers. On average, the unsigned error in the average polarizability and anisotropy compared to B3LYP are 2% and 5%, respectively. The correlation between the polarizability components from B3LYP and this approach lead to a R2 of 0.990 and a slope of 0.999. Even the F2 anisotropy, shown to be a difficult case for existing polarizability models, can be reproduced within 2% error. In addition to providing new parameters for a rapid method directly applicable to the calculation of polarizabilities, this work extends the widely used Poisson equation to areas where accurate molecular polarizabilities matter.

Semi-continuum electrostatic calculations of redox potentials in photosystem I

The midpoint redox potentials (E(m)) of all cofactors in photosystem I from Synechococcus elongatus as well as of the iron-sulfur (Fe(4)S(4)) clusters in two soluble ferredoxins from Azotobacter vinelandii and Clostridium acidiurici were calculated within the framework of a semi-continuum dielectric approach. The widely used treatment of proteins as uniform media with single dielectric permittivity is oversimplified, particularly, because permanent charges are considered both as a source for intraprotein electric field and as a part of dielectric polarizability. Our approach overcomes this inconsistency by using two dielectric constants: optical epsilon(o)=2.5 for permanent charges pre-existing in crystal structure, and static epsilon(s) for newly formed charges. We also take into account a substantial dielectric heterogeneity of photosystem I revealed by photoelectric measurements and a liquid junction potential correction for E(m) values of relevant redox cofactors measured in aprotic solvents. We show that calculations based on a single permittivity have the discrepancy with experimental data larger than 0.7 V, whereas E(m) values calculated within our approach fall in the range of experimental estimates. The electrostatic analysis combined with quantum chemistry calculations shows that (i) the energy decrease upon chlorophyll dimerization is essential for the downhill mode of primary charge separation between the special pair P(700) and the primary acceptor A(0); (ii) the primary donor is apparently P(700) but not a pair of accessory chlorophylls; (iii) the electron transfer from the A branch quinone Q(A) to the iron-sulfur cluster F(X) is most probably downhill, whereas that from the B branch quinone Q(B) to F(X) is essentially downhill.

Dielectric relaxation in proteins: the computational perspective

In photoexcitation and electron transfer, a new dipole or charge is introduced, and the structure is adjusted. This adjustment represents dielectric relaxation, which is the focus of this review. We concentrate on a few selected topics. We discuss linear response theory, as a unifying framework and a tool to describe non-equilibrium states. We review recent, molecular dynamics simulation studies that illustrate the calculation of dynamic and thermodynamic properties, such as Stokes shifts or reorganization free energies. We then turn to the macroscopic, continuum electrostatic view. We recall the physical definition of a dielectric constant and revisit the decomposition of the free energy into a reorganization and a static term. We review some illustrative continuum studies and discuss some difficulties that can arise with the continuum approach. In conclusion, we consider recent developments that will increase the accuracy and broaden the scope of all these methods.

Structure-based pKa calculations using continuum electrostatics methods

Electrostatic free energy is useful for correlating structure with function in proteins in which ionizable groups play essential functional roles. To this end, the pK(a) values of ionizable groups must be known and their molecular determinants must be understood. Structure-based calculations of electrostatic energies and pK(a) values are necessary for this purpose. This unit describes protocols for pK(a) calculations with continuum electrostatics methods based on the numerical solution of the linearized Poisson-Boltzmann equation by the method of finite differences. Critical discussion of key parameters, approximations, and shortcomings of these methods is included. Two protocols are described for calculations with methods modified empirically to maximize agreement between measured and calculated pK(a) values. Applied judiciously, these methods can contribute useful and novel insight into properties of surface ionizable groups in proteins.

Proton pathways in a [NiFe]-hydrogenase: A theoretical study

We present here a theoretical study to investigate possible proton pathways in the [NiFe]-hydrogenase from Desulfovibrio gigas. The approach used in this study consists of a combination of Poisson-Boltzmann and Monte Carlo simulations together with a distance-based network analysis to find possible groups involved in the proton transfer. Results obtained at different pH values show a reasonable number of proton active residues distributed by the protein interior and surface, with a concentration around the metal centres. The electrostatic interactions in this protein are strong, as shown by the unusual shape of the titration curves of several sites. Some residue pairs show strongly correlated protonations, indicating the sharing and probably exchange of a proton between them. The conjugation of the PB and MC simulations with the distance-based analysis allows a detailed characterization of the possible proton pathways. We discuss previous suggestions and propose a new complete pathway for the proton transfer between the active site and the surface. This pathway is mainly composed of histidines and glutamic acid residues.

pH modulation of transport properties of alamethicin oligomers inserted in zwitterionic-based artificial lipid membranes

Electric features of biological membranes are major determinants of the function and physiological manifestation of membrane-penetrating peptides, and such features are prone to be modulated by the properties of the surrounding aqueous medium. In this work, we demonstrate that pH plays crucial roles in modulating electric characteristics of zwitterionic-based artificial lipid membranes. The effect of pH on electrical properties of such membranes was probed by evaluating the transport properties of embedded alamethicin oligomers over a wide range of pH values (i.e., 0.65, 2.08, 2.94, 7 and 10.1). Our data strongly support the paradigm of a pH-dependent variation of the surface and membrane dipole potential which, in conjunction with possible lateral pressure effects within the lipid membrane, lead to a non-monotonic modulation of the electrical conductance of alamethicin oligomers. As expected, pH modulation of transport properties through the alamethicin oligomer is more visible for narrower pores (that is, the 1st conductive state) with slightly better cation selectivity as compared to larger oligomers.

Spatial distribution of dielectric shielding in the interior of Pyrococcus furiosus rubredoxin as sampled in the subnanosecond timeframe by hydrogen exchange

Experimental pK values of ionizable sidechains provide the most direct test for models representing dielectric shielding within the interior of a protein. However, only the strongly shifted pK values are particularly useful for discriminating among models. NMR titration studies have usually found only one or two such shifted pK values in each protein, so that the fitting of the experimental data to a uniform internal dielectric (epsilon(int)) model is not well constrained. The observed variation among proteins for such epsilon(int) estimates may reflect nonuniformity of dielectric shielding within each protein interior or qualitative differences between individual proteins. The differential amide kinetic acidities for a series of metal-substituted rubredoxins are shown to be consistent with Poisson-Boltzmann predictions of dielectric shielding that is relatively uniform for all of the amides that are sensitive to the metal charge, a region which corresponds to roughly 1/3 of the internal volume. The effective epsilon(int) values near 6 that are found in this study are significantly lower than many such estimates derived from sidechain pK measurements. The differing timeframes in which dielectric relaxation can respond to the highly transient peptide anion as compared to the longer lived states of the charged sidechains offers an explanation for the lower apparent dielectric constant deduced from these measurements.

The pH-dependent conformational states of kyotorphin: a constant-pH molecular dynamics study

An extensive conformational study of the analgesic dipeptide kyotorphin (L-Tyr-L-Arg) at different pH values was performed using a constant-pH molecular dynamics method. This dipeptide showed a remarkable pH-dependent conformational variety. The protonation of the N-terminal amine was identified as a key element in the transition between the more extended and the more packed conformational states, as monitored by the dihedral angle defined by the atoms 1Cbeta-1Calpha-2Calpha-2Cbeta. The principal-component analysis of kyotorphin identified two major conformational populations (the extended trans and the packed cis) together with conformations that occur exclusively at extreme pH values. Other, less stable conformations were also identified, which help us to understand the transitions between the predominant populations. The fitting of kyotorphin's conformational space to the structure of morphine resulted in a set of conformers that were able to fulfill most of the constraints for the mu-receptor. These results suggest that there may be strong similarities between the kyotorphin receptor and the structural family of opioid receptors.

Free-energy simulations and experiments reveal long-range electrostatic interactions and substrate-assisted specificity in an aminoacyl-tRNA synthetase

Specific recognition of their cognate amino acid substrates by the aminoacyl-tRNA synthetase enzymes is essential for the correct translation of the genetic code. For aspartyl-tRNA synthetase (AspRS), electrostatic interactions are expected to play an important role, since its three substrates (aspartate, ATP, tRNA) are all electrically charged. We used molecular-dynamics free-energy simulations and experiments to compare the binding of the substrate Asp and its electrically neutral analogue Asn to AspRS. The preference for Asp is found to be very strong, with good agreement between simulations and experiment. The simulations reveal long-range interactions that electrostatically couple the amino acid ligand, ATP, and its associated Mg2+ cations, a histidine side chain (His448) next to the amino acid ligand and a flexible loop that closes over the active site in response to amino acid binding. Closing this loop brings a negatively charged glutamate into the active site; this causes His448 to recruit a labile proton, which interacts favorably with Asp and accounts for most of the Asp/Asn discrimination. Cobinding of the second substrate, ATP, increases specificity for Asp further and makes the system robust towards removal of His448, which is mutated to a neutral amino acid in many organisms. Thus, AspRS specificity is assisted by a labile proton and a cosubstrate, and ATP acts as a mobile discriminator for specific Asp binding to AspRS. In asparaginyl-tRNA synthetase, a close homologue of AspRS, a few binding-pocket differences modify the charge balance so that asparagine binding predominates.

Free energy simulations: use of reverse cumulative averaging to determine the equilibrated region and the time required for convergence

A method is proposed for improving the accuracy and efficiency of free energy simulations. The essential idea is that the convergence of the relevant measure (e.g., the free energy derivative in thermodynamic integration) is monitored in the reverse direction starting from the last frame of the trajectory, instead of the usual approach, which begins with the first frame and goes in the forward direction. This simple change in the use of the simulation data makes it straightforward to eliminate the contamination of the averages by contributions from the equilibrating region. A statistical criterion is introduced for distinguishing the equilibrated (production) region from the equilibrating region. The proposed method, called reverse cumulative averaging, is illustrated by its application to the well-studied case of the alchemical free energy simulation of ethane to methanol.

Probing molecular docking in a charged model binding site

A model binding site was used to investigate charge-charge interactions in molecular docking. This simple site, a small (180A(3)) engineered cavity in cyctochrome c peroxidase (CCP), is negatively charged and completely buried from solvent, allowing us to explore the balance between electrostatic energy and ligand desolvation energy in a system where many of the common approximations in docking do not apply. A database with about 5300 molecules was docked into this cavity. Retrospective testing with known ligands and decoys showed that overall the balance between electrostatic interaction and desolvation energy was captured. More interesting were prospective docking scre"ens that looked for novel ligands, especially those that might reveal problems with the docking and energy methods. Based on screens of the 5300 compound database, both high-scoring and low-scoring molecules were acquired and tested for binding. Out of 16 new, high-scoring compounds tested, 15 were observed to bind. All of these were small heterocyclic cations. Binding constants were measured for a few of these, they ranged between 20microM and 60microM. Crystal structures were determined for ten of these ligands in complex with the protein. The observed ligand geometry corresponded closely to that predicted by docking. Several low-scoring alkyl amino cations were also tested and found to bind. The low docking score of these molecules owed to the relatively high charge density of the charged amino group and the corresponding high desolvation penalty. When the complex structures of those ligands were determined, a bound water molecule was observed interacting with the amino group and a backbone carbonyl group of the cavity. This water molecule mitigates the desolvation penalty and improves the interaction energy relative to that of the "naked" site used in the docking screen. Finally, six low-scoring neutral molecules were also tested, with a view to looking for false negative predictions. Whereas most of these did not bind, two did (phenol and 3-fluorocatechol). Crystal structures for these two ligands in complex with the cavity site suggest reasons for their binding. That these neutral molecules do, in fact bind, contradicts previous results in this site and, along with the alkyl amines, provides instructive false negatives that help identify weaknesses in our scoring functions. Several improvements of these are considered.

Constant-pH Molecular Dynamics with Ionic Strength Effects: Protonation−Conformation Coupling in Decalysine

A new implementation of the stochastic titration method for constant-pH molecular dynamics is presented, which introduces ionic strength effects in the simulations. In addition, the new implementation uses a faster molecular dynamics algorithm and an improved treatment of protonation events and of their effect on force field parameters. This new methodology is applied to a decalysine peptide, yielding very good quantitative agreement with experiments, both in terms of titration and helix-coil transition. The results show a significant dependence on ionic strength, illustrating the importance of including this parameter in constant-pH molecular dynamics simulations. Overall, the method seems to properly capture the protonation-conformation coupling and its dependence on ionic strength.

The role of backbone stability near Ala44 in the high reduction potential class of rubredoxins

Rubredoxins may be separated into high and low reduction potential classes, with reduction potentials differing by approximately 50 mV. Our previous work showed that a local shift in the polar backbone due to an A(44) versus V(44) side-chain size causes this reduction potential difference. However, this work also indicated that in the low potential Clostridium pasteurianum (Cp) rubredoxin, a V(44) --> A(44) mutation causes larger local backbone flexibility, because the V(44) side-chain present in the wild-type (wt) is no longer present to interlock with neighboring residues to stabilize the subsequent G(45). Since Pyrococcus furiosus (Pf) and other high potential rubredoxins generally have a P(45), it was presumed that a G(45) --> P(45) mutation might stabilize a V(44) --> A(44) mutation in Cp rubredoxin. Here crystal structure analysis, energy minimization, and molecular dynamics (MD) were performed for wt V(44)G(45), single mutant A(44)G(45) and double mutant A(44)P(45) Cp, and for wt A(44)P(45) Pf rubredoxins. The local structural, dynamical, and electrostatic properties of Cp gradually approach wt Pf in the order wt Cp to single to double mutant because of greater sequence similarity, as expected. The double mutant A(44)P(45) Cp exhibits increased backbone stability near residue 44 and thus enhances the probability that the backbone dipoles point toward the redox site, which favors an increase in the electrostatic contribution to the reduction potential. It appears that the electrostatic potential of residue 44 and the solvent accessibility to the redox are both determinants for the reduction potentials of homologous rubredoxins. Overall, these results indicate that an A(44) in a rubredoxin may require a P(45) for backbone stability whereas a V(44) can accommodate a G(45), since the valine side-chain can interlock with its neighbors.

On the Use of Different Dielectric Constants for Computing Individual and Pairwise Terms in Poisson−Boltzmann Studies of Protein Ionization Equilibrium

Poisson-Boltzmann (PB) models are a fast and common tool for studying electrostatic processes in proteins, particularly their ionization equilibrium (protonation and/or reduction), often yielding quite good results when compared with more detailed models. Yet, they are conceptually very simple and necessarily approximate, their empirical character being most evident when it comes to the choice of the dielectric constant assigned to the protein region. The present study analyzes several factors affecting the ability of PB-based methods to model protein ionization equilibrium. We give particular attention to a suggestion made by Warshel and co-workers (e.g., Sham et al. J. Phys. Chem. B 1997, 101, 4458) of using different protein dielectric constants for computing the individual (site) and the pairwise (site-site) terms of the ionization free energies. Our prediction of pK(a) values for several proteins indicates that no advantage is obtained by such a procedure, even for sites that are buried and/or display large pK(a) shifts relative to the solution values. In particular, the present methodology gives the best predictions using a dielectric constant around 20, for shifted/buried and nonshifted/exposed sites alike. The similarities and differences between the PB model and Warshel's PDLD/S model are discussed, as well as the reasons behind their apparently discrepant results. The present PB model is shown to predict also good reduction potentials in redox proteins.

Proton binding to proteins: a free-energy component analysis using a dielectric continuum model

Proton binding plays a critical role in protein structure and function. We report pK(a) calculations for three aspartates in two proteins, using a linear response approach, as well as a "standard" Poisson-Boltzmann approach. Averaging over conformations from the two endpoints of the proton-binding reaction, the protein's atomic degrees of freedom are explicitly modeled. Treating macroscopically the protein's electronic polarizability and the solvent, a meaningful model is obtained, without adjustable parameters. It reproduces qualitatively the electrostatic potentials, proton-binding free energies, Marcus reorganization free energies, and pK(a) shifts from explicit solvent molecular dynamics simulations, and the pK(a) shifts from experiment. For thioredoxin Asp-26, which has a large pK(a) upshift, we correctly capture the balance between unfavorable carboxylate desolvation and favorable interactions with a nearby lysine; similarly for RNase A Asp-14, which has a large pK(a) downshift. For the unshifted thioredoxin Asp-20, desolvation by the protein cavity is overestimated by 2.9 pK(a) units; several effects could explain this. "Standard" Poisson-Boltzmann methods sidestep this problem by using a large, ad hoc protein dielectric; but protein charge-charge interactions are then incorrectly downscaled, giving an unbalanced description of the reaction and a large error for the shifted pK(a) values of Asp-26 and Asp-14.

Accurate pKa Determination for a Heterogeneous Group of Organic Molecules

Single-molecule studies that allow to compute pKa values, proton affinities (gas-phase acidity/basicity) and the electrostatic energy of solvation have been performed for a heterogeneous set of 26 organic compounds. Quantum mechanical density functional theory (DFT) using the Becke-half&half and B3LYP functionals on optimized molecular geometries have been carried out to investigate the energetics of gas-phase protonation. The electrostatic contribution to the solvation energies of protonated and deprotonated compounds were calculated by solving the Poisson equation using atomic charges generated by fitting the electrostatic potential derived from the molecular wave functions in vacuum. The combination of gas-phase and electrostatic solvation energies by means of the thermodynamic cycle enabled us to compute pKa values for the 26 compounds, which cover six distinct chemical groups (carboxylic acids, benzoic acids, phenols, imides, pyridines and imidazoles). The computational procedure for determining pKa values is accurate and transferable with a root-mean-square deviation of 0.53 and 0.57 pKa units and a maximum error of 1.0 pKa and 1.3 pKa units for Becke-half&half and B3LYP DFT functionals, respectively.

Escherichia coli glutaminyl-tRNA synthetase is electrostatically optimized for binding of its cognate substrates

Natural evolution has resulted in protein molecules displaying a wide range of binding properties that include extremes of affinity and specificity. A detailed understanding of the principles underlying protein structure-function relationships, particularly with respect to binding properties, would greatly enhance molecular engineering and ligand design studies. Here, we have analyzed the interactions of an aminoacyl-tRNA synthetase for which strong evolutionary pressure has enforced high specificity for substrate binding and catalysis. Electrostatic interactions have been identified as one efficient mechanism for enhancing binding specificity; as such, the effects of charged and polar groups were the focus of this study. The binding of glutaminyl-tRNA synthetase from Escherichia coli to several ligands, including the natural substrates, was analyzed. The electrostatic complementarity of the enzyme to its ligands was assessed using measures derived from affinity optimization theory. The results were independent of the details of the calculational parameters, including the value used for the protein dielectric constant. Glutamine and ATP, two of the natural ligands, were found to be extremely complementary to their binding sites, particularly in regions seen to make electrostatic interactions in the structure. These data suggest that the optimization of electrostatic interactions has played an important role in guiding the evolution of this enzyme. The results also show that the enzyme is able to effectively select for high affinity and specificity for the same chemical moieties both in the context of smaller substrates, and in that of a larger reactive intermediate. The regions of greatest non-complementarity between the enzyme and ligands are the portions of the ligand that make few polar contacts with the binding site, as well as the sites of chemical reaction, where overly strong electrostatic binding interactions with the substrate could hinder catalysis. The results also suggest that the negative charge on the phosphorus center of glutaminyl-adenylate plays an important role in the tight binding of this intermediate, and thus that adenylate analogs that preserve the negative charge in this region may bind substantially tighter than analogs where this group is replaced with a neutral group, such as the sulfamoyl family, which can make similar hydrogen bonds but is uncharged.

A quantum mechanical study on phosphotyrosyl peptide binding to the SH2 domain of p56lck tyrosine kinase with insights into the biochemistry of intracellular signal transduction events

A study on the interaction between a phosphotyrosyl peptide with the SH2 domain of Lck kinase has been undertaken with the aid of semiempirical linear-scaling quantum mechanical methods. The structure of this complex has been solved at atomic resolution and, hence, it represents the ideal candidate for studying the charge deformation effects induced by the phosphopeptide on the binding site. Substantial changes in the charge of amino acid residues located in the binding pocket of the protein are observed upon ligand binding. More specifically, our quantum chemical calculations indicate that H-bonds involving charged side-chains are subject to consistent charge deformation effects whereas those forming salt bridges are unaffected by ligand binding. Furthermore, ligand binding has the effect of changing both the magnitude and direction of the protein's macrodipole, which rotates approximately 150 degrees with respect that of the unliganded protein. This suggests that a change in the polarization state of the protein might acts as a switch during the transmission of intracellular signals. The binding energy calculated with the aid of the COSMO solvation model corresponds to about -200 kcal/mol, most of which is attributed to the interaction of the phosphotyrosine head with the amino acid chains located in the binding site of the SH2 domain.

Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models

Ionizable residues play important roles in protein structure and activity, and proton binding is a valuable reporter of electrostatic interactions in these systems. We use molecular dynamics free energy simulations (MDFE) to compute proton pKa shifts, relative to a model compound in solution, for three aspartate side chains in two proteins. Simulations with explicit solvent and with an implicit, dielectric continuum solvent are reported. The implicit solvent simulations use the generalized Born (GB) model, which provides an approximate, analytical solution to Poisson's equation. With explicit solvent, the direction of the pKa shifts is correct in all three cases with one force field (AMBER) and in two out of three cases with another (CHARMM). For two aspartates, the dielectric response to ionization is found to be linear, even though the separate protein and solvent responses can be nonlinear. For thioredoxin Asp26, nonlinearity arises from the presence of two substates that correspond to the two possible orientations of the protonated carboxylate. For this side chain, which is partly buried and has a large pKa upshift, very long simulations are needed to correctly sample several slow degrees of freedom that reorganize in response to the ionization. Thus, nearby Lys57 rotates to form a salt bridge and becomes buried, while three waters intercalate along the opposite edge of Asp26. Such strong and anisotropic reorganization is very difficult to predict with Poisson-Boltzmann methods that only consider electrostatic interactions and employ a single protein structure. In contrast, MDFE with a GB dielectric continuum solvent, used for the first time for pKa calculations, can describe protein reorganization accurately and gives encouraging agreement with experiment and with the explicit solvent simulations.