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

Million-atom molecular dynamics simulations reveal the interfacial interactions and assembly of plant PSII-LHCII supercomplex

Protein-protein interface interactions dictate efficient excitation energy transfer from light-harvesting antennas to the photosystem II (PSII) core. In this work, we construct a 1.2 million atom-scale model of plant C₂S₂-type PSII-LHCII supercomplex and perform microsecond-scale molecular dynamics (MD) simulations to explore the interactions and assembly mechanisms of the sizeable PSII-LHCII supercomplex. We optimize the nonbonding interactions of the PSII-LHCII cryo-EM structure using microsecond-scale MD simulations. Binding free energy calculations with component decompositions reveal that hydrophobic interactions predominantly drive antenna-core association and the antenna-antenna interactions are relatively weak. Despite the positive electrostatic interaction energies, hydrogen bonds and salt bridges mainly provide directional or anchoring forces for interface binding. Analysis of the roles of small intrinsic subunits of PSII suggests that LHCII and CP26 first interact with small intrinsic subunits and then bind to the core proteins, whereas CP29 adopts a one-step binding process to the PSII core without the assistance of other factors. Our study provides insights into the molecular underpinnings of the self-organization and regulation of plant PSII-LHCII. It lays the framework for deciphering the general assembly principles of photosynthetic supercomplexes and possibly other macromolecular structures. The finding also has implications for repurposing photosynthetic systems to enhance photosynthesis.

SAMPL7 TrimerTrip host-guest binding poses and binding affinities from spherical-coordinates-biased simulations

Host-guest binding remains a major challenge in modern computational modelling. The newest 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge contains a new series of host-guest systems. The TrimerTrip host binds to 16 structurally diverse guests. Previously, we have successfully employed the spherical coordinates as the collective variables coupled with the enhanced sampling technique metadynamics to enhance the sampling of the binding/unbinding event, search for possible binding poses and calculate the binding affinities in all three host-guest binding cases of the 6th SAMPL challenge. In this work, we report a retrospective study on the TrimerTrip host-guest systems by employing the same protocol to investigate the TrimerTrip host in the SAMPL7 challenge. As no binding pose is provided by the SAMPL7 host, our simulations initiate from randomly selected configurations and are proceeded long enough to obtain converged free energy estimates and search for possible binding poses. The calculated binding affinities are in good agreement with the experimental reference, and the obtained binding poses serve as a nice starting point for end-point or alchemical free energy calculations. Note that as the work is performed after the close of the SAMPL7 challenge, we do not participate in the challenge and the results are not formally submitted to the SAMPL7 challenge.

Adaptive landscape flattening allows the design of both enzyme: Substrate binding and catalytic power

Designed enzymes are of fundamental and technological interest. Experimental directed evolution still has significant limitations, and computational approaches are a complementary route. A designed enzyme should satisfy multiple criteria: stability, substrate binding, transition state binding. Such multi-objective design is computationally challenging. Two recent studies used adaptive importance sampling Monte Carlo to redesign proteins for ligand binding. By first flattening the energy landscape of the apo protein, they obtained positive design for the bound state and negative design for the unbound. We have now extended the method to design an enzyme for specific transition state binding, i.e., for its catalytic power. We considered methionyl-tRNA synthetase (MetRS), which attaches methionine (Met) to its cognate tRNA, establishing codon identity. Previously, MetRS and other synthetases have been redesigned by experimental directed evolution to accept noncanonical amino acids as substrates, leading to genetic code expansion. Here, we have redesigned MetRS computationally to bind several ligands: the Met analog azidonorleucine, methionyl-adenylate (MetAMP), and the activated ligands that form the transition state for MetAMP production. Enzyme mutants known to have azidonorleucine activity were recovered by the design calculations, and 17 mutants predicted to bind MetAMP were characterized experimentally and all found to be active. Mutants predicted to have low activation free energies for MetAMP production were found to be active and the predicted reaction rates agreed well with the experimental values. We suggest the present method should become the paradigm for computational enzyme design.

Designed enzymes are of major interest. Experimental directed evolution still has significant limitations, and computational approaches are another route. Enzymes must be stable, bind substrates, and be powerful catalysts. It is challenging to design for all these properties. A method to design substrate binding was proposed recently. It used an adaptive Monte Carlo method to explore mutations of a few amino acids near the substrate. A bias energy was gradually “learned” such that, in the absence of the ligand, the simulation visited most of the possible protein mutations with comparable probabilities. Remarkably, a simulation of the protein:ligand complex, including the bias, will then preferentially sample tight-binding sequences. We generalized the method to design binding specificity. We tested it for the methionyl-tRNA synthetase enzyme, which has been engineered in order to expand the genetic code. We redesigned the enzyme to obtain variants with low activation free energies for the catalytic step. The variants proposed by the simulations were shown experimentally to be active, and the predicted activation free energies were in reasonable agreement with the experimental values. We expect the new method will become the paradigm for computational enzyme design.

Theoretical understanding of the thermodynamics and interactions in transcriptional regulator TtgR-ligand binding

The transcriptional regulator TtgR belongs to the TetR family of transcriptional repressors. It depresses the transcription of the TtgABC operon and itself and thus regulates the extrusion of noxious chemicals with efflux pumps in bacterial cells. As the ligand-binding domain of TtgR is rather flexible, it can bind with a number of structurally diverse ligands, such as antibiotics, flavonoids and aromatic solvents. In the current work, we perform equilibrium and nonequilibrium alchemical free energy simulation to predict the binding affinities of a series of ligands targeting the TtgR protein and an agreement between the theoretical prediction and the experimental result is observed. End-point methods MM/PBSA and MM/GBSA are also employed for comparison. We further study the interaction maps and contacts between the protein and the ligand and identify important interactions in the protein-ligand binding cases. The dynamics fluctuation and secondary structures are also investigated. The current work sheds light on atomic and thermodynamic understanding of the TtgR-ligand interactions.

Development of constant-pH simulation methods in implicit solvent and applications in biomolecular systems

pH is a critical parameter for biological and technological systems directly related with electrical charges. It can give rise to peculiar electrostatic phenomena, which also makes them more challenging. Due to the quantum nature of the process, involving the forming and breaking of chemical bonds, quantum methods should ideally by employed. Nevertheless, due to the very large number of ionizable sites, different macromolecular conformations, salt conditions, and all other charged species, the CPU time cost simply becomes prohibitive for computer simulations, making this a quite complex problem. Simplified methods based on Monte Carlo sampling have been devised and will be reviewed here, highlighting the updated state-of-the-art of this field, advantages, and limitations of different theoretical protocols for biomolecular systems (proteins and nucleic acids). Following a historical perspective, the discussion will be associated with the applications to protein interactions with other proteins, polyelectrolytes, and nanoparticles.

Model for screened, charge-regulated electrostatics of an eye lens protein: Bovine gammaB-crystallin

We model screened, site-specific charge regulation of the eye lens protein bovine gammaB-crystallin (γB) and study the probability distributions of its proton occupancy patterns. Using a simplified dielectric model, we solve the linearized Poisson-Boltzmann equation to calculate a 54×54 work-of-charging matrix, each entry being the modeled voltage at a given titratable site, due to an elementary charge at another site. The matrix quantifies interactions within patches of sites, including γB charge pairs. We model intrinsic pK values that would occur hypothetically in the absence of other charges, with use of experimental data on the dependence of pK values on aqueous solution conditions, the dielectric model, and literature values. We use Monte Carlo simulations to calculate a model grand-canonical partition function that incorporates both the work-of-charging and the intrinsic pK values for isolated γB molecules and we calculate the probabilities of leading proton occupancy configurations, for 4<pH<8 and Debye screening lengths from 6 to 20 Å. We select the interior dielectric value to model γB titration data. At pH 7.1 and Debye length 6.0 Å, on a given γB molecule the predicted top occupancy pattern is present nearly 20% of the time, and 90% of the time one or another of the first 100 patterns will be present. Many of these occupancy patterns differ in net charge sign as well as in surface voltage profile. We illustrate how charge pattern probabilities deviate from the multinomial distribution that would result from use of effective pK values alone and estimate the extents to which γB charge pattern distributions broaden at lower pH and narrow as ionic strength is lowered. These results suggest that for accurate modeling of orientation-dependent γB-γB interactions, consideration of numerous pairs of proton occupancy patterns will be needed.

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.

Extensive Assessment of Various Computational Methods for Aspartate's pKa Shift

A series of computational methods for pK_a shift prediction are extensively tested on a set of benchmark protein systems, aiming at identifying pitfalls and evaluating their performance on high variants. Including 19 ASP residues in 10 protein systems, the benchmark set consists of both residues with highly shifted pK_a values as well as those varying little from the reference value, with an experimental RMS free energy differences of 2.49 kcal/mol with respect to blocked amino acid, namely the RMS pK_a shift being 1.82 pK_a units. The constant pH molecular dynamics (MD), alchemical methods, PROPKA3.1, and multiconformation continuum electrostatics give RMSDs of 1.52, 2.58, 1.37, and 3.52 pK_a units, respectively, on the benchmark set. The empirical scoring method is the most accurate one with extremely low computational cost, and the pH-dependent model is also able to provide accurate results, while the accuracy of MD sampling incorporating alchemical free energy simulation is prohibited by convergence achievement and the performance of conformational search incorporating multiconformation continuum electrostatics is bad. Former research works did not define statistical uncertainty with care and yielded the questionable conclusion that alchemical methods perform well in most benchmarks. In this work the traditional alchemical methods are thoroughly tested for high variants. We also performed the first application of nonequilibrium alchemical methods to the pK_a cases.

Microscopic mechanisms that govern the titration response and pKa values of buried residues in staphylococcal nuclease mutants

To probe the microscopic mechanisms that govern the titration behavior of buried ionizable groups, microsecond explicit solvent molecular dynamics simulations are carried out for several mutants of Staphylococcal nuclease using both fixed charge and polarizable force fields. While the ionization of Asp 66, Glu 66, and Lys 125 lead to enhanced structural fluctuations and partial unfolding of adjacent α-helical regions, the ionization of Lys 25 causes local unfolding of adjacent β sheets. Using the sampled conformational ensembles, good agreement with experimental pK_a values is obtained with Poisson-Boltzmann calculations using a protein dielectric constant of 2-4 for V66D/E; slightly larger dielectric constants are needed for Lys mutants especially L25K, suggesting that structural responses beyond microseconds are involved in ionization of Lys 25. Overall, the set of unbiased simulations provides insights into the spatial and temporal scales of protein and solvent motions that dictate the diverse titration behaviors of buried protein residues. Proteins 2017; 85:268-281. © 2016 Wiley Periodicals, Inc.

Assessing the Potential Effects of Active Site Mg2+ Ions in the glmS Ribozyme-Cofactor Complex

Ribozymes employ diverse catalytic strategies in their self-cleavage mechanisms, including the use of divalent metal ions. This work explores the effects of Mg2+ ions in the active site of the glmS ribozyme-GlcN6P cofactor complex using computational methods. Deleterious and potentially beneficial effects of an active site Mg2+ ion on the self-cleavage reaction were identified. The presence of a Mg2+ ion near the scissile phosphate oxygen atoms at the cleavage site was determined to be deleterious, and thereby anticatalytic, due to electrostatic repulsion of the cofactor, disruption of key hydrogen-bonding interactions, and obstruction of nucleophilic attack. On the other hand, the presence of a Mg2+ ion at another position in the active site, the Hoogsteen face of the putative base, was found to avoid these deleterious effects and to be potentially catalytically favorable owing to the stabilization of negative charge and pKa shifting of the guanine base.

Changes in active site histidine hydrogen bonding trigger cryptochrome activation

Cryptochrome (CRY) is the principal light sensor of the insect circadian clock. Photoreduction of the Drosophila CRY (dCRY) flavin cofactor to the anionic semiquinone (ASQ) restructures a C-terminal tail helix (CTT) that otherwise inhibits interactions with targets that include the clock protein Timeless (TIM). All-atom molecular dynamics (MD) simulations indicate that flavin reduction destabilizes the CTT, which undergoes large-scale conformational changes (the CTT release) on short (25 ns) timescales. The CTT release correlates with the conformation and protonation state of conserved His378, which resides between the CTT and the flavin cofactor. Poisson-Boltzmann calculations indicate that flavin reduction substantially increases the His378 pKa Consistent with coupling between ASQ formation and His378 protonation, dCRY displays reduced photoreduction rates with increasing pH; however, His378Asn/Arg variants show no such pH dependence. Replica-exchange MD simulations also support CTT release mediated by changes in His378 hydrogen bonding and verify other responsive regions of the protein previously identified by proteolytic sensitivity assays. His378 dCRY variants show varying abilities to light-activate TIM and undergo self-degradation in cellular assays. Surprisingly, His378Arg/Lys variants do not degrade in light despite maintaining reactivity toward TIM, thereby implicating different conformational responses in these two functions. Thus, the dCRY photosensory mechanism involves flavin photoreduction coupled to protonation of His378, whose perturbed hydrogen-bonding pattern alters the CTT and surrounding regions.

Constant-pH Molecular Dynamics Study of Kyotorphin in an Explicit Bilayer

To our knowledge, we present the first constant-pH molecular dynamics study of the neuropeptide kyotorphin in the presence of an explicit lipid bilayer. The overall conformation freedom of the peptide was found to be affected by the interaction with the membrane, in accordance with previous results using different methodologies. Analysis of the interactions between the N-terminus amine group of the peptide and several lipid atoms shows that the membrane is able to stabilize both ionized and neutral forms of kyotorphin, resulting in a pKa value that is similar to the one obtained in water. This illustrates how a detailed molecular model of the membrane leads to rather different results than would be expected from simply regarding it as a low-dielectric slab.

Computation of pH-dependent binding free energies

Protein-ligand binding accompanies changes in the surrounding electrostatic environments of the two binding partners and may lead to changes in protonation upon binding. In cases where the complex formation results in a net transfer of protons, the binding process is pH-dependent. However, conventional free energy computations or molecular docking protocols typically employ fixed protonation states for the titratable groups in both binding partners set a priori, which are identical for the free and bound states. In this review, we draw attention to these important yet largely ignored binding-induced protonation changes in protein-ligand association by outlining physical origins and prevalence of the protonation changes upon binding. Following a summary of various theoretical methods for pKa prediction, we discuss the theoretical framework to examine the pH dependence of protein-ligand binding processes.

Benchmarking pKa Prediction Methods for Residues in Proteins

Methods for estimation of pKa values of residues in proteins were tested on a set of benchmark proteins with experimentally known pKa values. The benchmark set includes 80 different residues (20 each for Asp, Glu, Lys, and His), half of which consists of significantly variant cases (ΔpKa ≥ 1 pKa unit from the amino acid in solution). The method introduced by Case and co-workers [J. Am. Chem. Soc. 2004, 126, 4167-4180], referred to as the molecular dynamics/generalized-Born/thermodynamic integration (MD/GB/TI) technique, gives a root-mean-square deviation (rmsd) of 1.4 pKa units on the benchmark set. The use of explicit waters in the immediate region surrounding the residue was shown to generally reduce high errors for this method. Longer simulation time was also shown to increase the accuracy of this method. The empirical approach developed by Jensen and co-workers [Proteins 2005, 61, 704-721], PROPKA, also gives an overall rmsd of 1.4 pKa units and is more or less accurate based on residue type-the method does very well for Lys and Glu, but less so for Asp and His. Likewise, the absolute deviation is quite similar for the two methods-5.2 for PROPKA and 5.1 for MD/GB/TI. A comparison of these results with several prediction methods from the literature is presented. The error in pKa prediction is analyzed as a function of variation of the pKa from that in water and the solvent accessible surface area (SASA) of the residue. A case study of the catalytic lysine residue in 2-deoxyribose-5-phosphate aldolase (DERA) is also presented.

pKa values in proteins determined by electrostatics applied to molecular dynamics trajectories

For a benchmark set of 194 measured pKa values in 13 proteins, electrostatic energy computations are performed in which pKa values are computed by solving the Poisson-Boltzmann equation. In contrast to the previous approach of Karlsberg(+) (KB(+)) that essentially used protein crystal structures with variations in their side chain conformations, the present approach (KB2(+)MD) uses protein conformations from four molecular dynamics (MD) simulations of 10 ns each. These MD simulations are performed with different specific but fixed protonation patterns, selected to sample the conformational space for the different protonation patterns faithfully. The root-mean-square deviation between computed and measured pKa values (pKa RMSD) is shown to be reduced from 1.17 pH units using KB(+) to 0.96 pH units using KB2(+)MD. The pKa RMSD can be further reduced to 0.79 pH units, if each conformation is energy-minimized with a dielectric constant of εmin = 4 prior to calculating the electrostatic energy. The electrostatic energy expressions upon which the computations are based have been reformulated such that they do not involve terms that mix protein and solvent environment contributions and no thermodynamic cycle is needed. As a consequence, conformations of the titratable residues can be treated independently in the protein and solvent environments. In addition, the energy terms used here avoid the so-called intrinsic pKa and can therefore be interpreted without reference to arbitrary protonation states and conformations.

Role of the active site guanine in the glmS ribozyme self-cleavage mechanism: quantum mechanical/molecular mechanical free energy simulations

The glmS ribozyme catalyzes a self-cleavage reaction at the phosphodiester bond between residues A-1 and G1. This reaction is thought to occur by an acid-base mechanism involving the glucosamine-6-phosphate cofactor and G40 residue. Herein quantum mechanical/molecular mechanical free energy simulations and pKa calculations, as well as experimental measurements of the rate constant for self-cleavage, are utilized to elucidate the mechanism, particularly the role of G40. Our calculations suggest that an external base deprotonates either G40(N1) or possibly A-1(O2'), which would be followed by proton transfer from G40(N1) to A-1(O2'). After this initial deprotonation, A-1(O2') starts attacking the phosphate as a hydroxyl group, which is hydrogen-bonded to deprotonated G40, concurrent with G40(N1) moving closer to the hydroxyl group and directing the in-line attack. Proton transfer from A-1(O2') to G40 is concomitant with attack of the scissile phosphate, followed by the remainder of the cleavage reaction. A mechanism in which an external base does not participate, but rather the proton transfers from A-1(O2') to a nonbridging oxygen during nucleophilic attack, was also considered but deemed to be less likely due to its higher effective free energy barrier. The calculated rate constant for the favored mechanism is in agreement with the experimental rate constant measured at biological Mg(2+) ion concentration. According to these calculations, catalysis is optimal when G40 has an elevated pKa rather than a pKa shifted toward neutrality, although a balance among the pKa's of A-1, G40, and the nonbridging oxygen is essential. These results have general implications, as the hammerhead, hairpin, and twister ribozymes have guanines at a similar position as G40.

Microscopic basis for kinetic gating in Cytochrome c oxidase: insights from QM/MM analysis

Understanding the mechanism of vectorial proton pumping in biomolecules requires establishing the microscopic basis for the regulation of both thermodynamic and kinetic features of the relevant proton transfer steps.

Understanding the mechanism of vectorial proton pumping in biomolecules requires establishing the microscopic basis for the regulation of both thermodynamic and kinetic features of the relevant proton transfer steps. For the proton pump cytochrome c oxidase, while the regulation of thermodynamic driving force for key proton transfers has been discussed in great detail, the microscopic basis for the control of proton transfer kinetics has been poorly understood. Here we carry out extensive QM/MM free energy simulations to probe the kinetics of relevant proton transfer steps and analyze the effects of local structure and hydration level. We show that protonation of the proton loading site (PLS, taken to be a propionate of heme a₃) requires a concerted process in which a key glutamic acid (Glu286H) delivers the proton to the PLS while being reprotonated by an excess proton coming from the D-channel. The concerted nature of the mechanism is a crucial feature that enables the loading of the PLS before the cavity containing Glu286 is better hydrated to lower its pK_a to experimentally measured range; the charged rather than dipolar nature of the process also ensures a tight coupling with heme a reduction, as emphasized by Siegbahn and Blomberg. In addition, we find that rotational flexibility of the PLS allows its protonation before that of the binuclear center (the site where oxygen gets reduced to water). Together with our recent study (P. Goyal, et al., Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 18886–18891) that focused on the modulation of Glu286 pK_a, the current work suggests a mechanism that builds in a natural sequence for the protonation of the PLS prior to that of the binuclear center. This provides microscopic support to the kinetic constraints revealed by kinetic network analysis as essential elements that ensure an efficient vectorial proton transport in cytochrome c oxidase.

Electrostatic free energies in translational GTPases: Classic allostery and the rest

GTPases typically switch between an inactive, OFF conformation and an active, ON conformation when a GDP ligand is replaced by GTP. Their ON/OFF populations and activity thus depend on the stabilities of four protein complexes, two apo-protein forms, and GTP/GDP in solution. A complete characterization is usually not possible experimentally and poses major challenges for simulations. We review the most important methodological challenges and we review thermodynamic data for two GTPases involved in translation of the genetic code: archaeal Initiation Factors 2 and 5B (aIF2, aIF5B). One main challenge is the multiplicity of states and conformations, including those of GTP/GDP in solution. Another is force field accuracy, especially for interactions of GTP/GDP with co-bound divalent Mg(2+) ions. The calculation of electrostatic free energies also poses specific challenges, and requires careful protocols. For aIF2, experiments and earlier simulations showed that it is a "classic" GTPase, with distinct ON/OFF conformations that prefer to bind GTP and GDP, respectively. For aIF5B, we recently proposed a non-classic mechanism, where the ON/OFF states differ only in the protonation state of Glu81 in the nucleotide binding pocket. This model is characterized here using free energy simulations. The methodological analysis should help future studies, while the aIF2, aIF5B examples illustrate the diversity of ATPase/GTPase mechanisms. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Quantum mechanics/molecular mechanics restrained electrostatic potential fitting

We present a quantum mechanics/molecular mechanics (QM/MM) method to evaluate the partial charges of amino acid residues for use in MM potentials based on their protein environment. For each residue of interest, the nearby residues are included in the QM system while the rest of the protein is treated at the MM level of theory. After a short structural optimization, the partial charges of the central residue are fit to the electrostatic potential using the restrained electrostatic potential (RESP) method. The resulting charges and electrostatic potential account for the individual environment of the residue, although they lack the transferable nature of library partial charges. To evaluate the quality of the QM/MM RESP charges, thermodynamic integration is used to measure the pKa shift of the aspartic acid residues in three different proteins, turkey egg lysozyme, beta-cryptogein, and Thioredoxin. Compared to the AMBER ff99SB library values, the QM/MM RESP charges show better agreement between the calculated and experimental pK(a) values for almost all of the residues considered.

Changing hydration level in an internal cavity modulates the proton affinity of a key glutamate in cytochrome c oxidase

Cytochrome c oxidase contributes to the transmembrane proton gradient by removing two protons from the high-pH side of the membrane each time the binuclear center active site is reduced. One proton goes to the binuclear center, whereas the other is pumped to the low-pH periplasmic space. Glutamate 286 (Glu286) has been proposed to serve as a transiently deprotonated proton donor. Using unrestrained atomistic molecular dynamics simulations, we show that the size of and water distribution in the hydrophobic cavity that holds Glu286 is controlled by the protonation state of the propionic acid of heme a3, a group on the proton outlet pathway. Protonation of the propionate disrupts hydrogen bonding to two side chains, allowing a loop to swing open. Continuum electrostatics and atomistic free-energy perturbation calculations show that the resultant changes in hydration and electrostatic interactions lower the Glu proton affinity by at least 5 kcal/mol. These changes in the internal hydration level occur in the absence of major conformational transitions and serve to stabilize needed transient intermediates in proton transport. The trigger is not the protonation of the Glu of interest, but rather the protonation of a residue ∼10 Å away. Thus, unlike local water penetration to stabilize a new charge, this finding represents a specific role for water molecules in the protein interior, mediating proton transfers and facilitating ion transport.

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.

Progress in the prediction of pKa values in proteins

The pK(a) -cooperative aims to provide a forum for experimental and theoretical researchers interested in protein pK(a) values and protein electrostatics in general. The first round of the pK(a) -cooperative, which challenged computational labs to carry out blind predictions against pK(a) s experimentally determined in the laboratory of Bertrand Garcia-Moreno, was completed and results discussed at the Telluride meeting (July 6-10, 2009). This article serves as an introduction to the reports submitted by the blind prediction participants that will be published in a special issue of PROTEINS: Structure, Function and Bioinformatics. Here, we briefly outline existing approaches for pK(a) calculations, emphasizing methods that were used by the participants in calculating the blind pK(a) values in the first round of the cooperative. We then point out some of the difficulties encountered by the participating groups in making their blind predictions, and finally try to provide some insights for future developments aimed at improving the accuracy of pK(a) calculations.

Simulating electrostatic energies in proteins: perspectives and some recent studies of pKas, redox, and other crucial functional properties

Electrostatic energies provide what is arguably the most effective tool for structure-function correlation of biological molecules. Here, we provide an overview of the current state-of-the-art simulations of electrostatic energies in macromolecules, emphasizing the microscopic perspective but also relating it to macroscopic approaches. We comment on the convergence issue and other problems of the microscopic models and the ways of keeping the microscopic physics while moving to semi-macroscopic directions. We discuss the nature of the protein dielectric "constants" reiterating our long-standing point that the dielectric "constants" in semi-macroscopic models depend on the definition and the specific treatment. The advances and the challenges in the field are illustrated considering different functional properties including pK(a)'s, redox potentials, ion and proton channels, enzyme catalysis, ligand binding, and protein stability. We emphasize the microscopic overcharging approach for studying pK(a) 's of internal groups in proteins and give a demonstration of power of this approach. We also emphasize recent advances in coarse grained models with a physically based electrostatic treatment and provide some examples including further directions in treating voltage activated ion channels.

Measuring the successes and deficiencies of constant pH molecular dynamics: a blind prediction study

A constant pH molecular dynamics method has been used in the blind prediction of pK_a values of titratable residues in wild type and mutated structures of the Staphylococcal nuclease (SNase) protein. The predicted values have been subsequently compared to experimental values provided by the laboratory of García-Moreno. CpHMD performs well in predicting the pK_a of solvent-exposed residues. For residues in the protein interior, the CpHMD method encounters some difficulties in reaching convergence and predicting the pK_a values for residues having strong interactions with neighboring residues. These results show the need to accurately and sufficiently sample conformational space in order to obtain pK_a values consistent with experimental results.

Electrostatics in proteins and protein-ligand complexes

Accurate computational methods for predicting electrostatic energies are of major importance for our understanding of protein energetics in general for computer-aided drug design as well as for the design of novel biocatalysts and protein therapeutics. Electrostatic energies are of particular importance in such applications as virtual screening, drug design and protein-protein docking due to the high charge density of protein ligands and small-molecule drugs, and the frequent protonation state changes observed when drugs bind to their protein targets. Therefore, the development of a reliable and fast algorithm for the evaluation of electrostatic free energies, as an important contributor to the overall protein energy function, has been the focus for many scientists over the past three decades. In this review we describe the current state-of-the-art in modeling electrostatic effects in proteins and protein-ligand complexes. We focus mainly on the merits and drawbacks of the continuum methodology, and speculate on future directions in refining algorithms for calculating electrostatic energies in proteins using experimental data.

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.

Electronic continuum model for molecular dynamics simulations of biological molecules

Electronic polarizability is an important factor in molecular interactions. In the conventional force fields such as AMBER or CHARMM, however, there is inconsistency in how the effect of electronic dielectric screening of Coulombic interactions, inherent for the condensed phase media, is treated. Namely, the screening appears to be accounted for via effective charges only for neutral moieties, whereas the charged residues are treated as if they were in vacuum. As a result, the electrostatic interactions between ionized groups are exaggerated in molecular simulations by the factor of about 2. The discussed here MDEC (Molecular Dynamics in Electronic Continuum) model provides a theoretical framework for modification of the standard non-polarizable force fields to make them consistent with the idea of uniform electronic screening of partial atomic charges. The present theory states that the charges of ionized groups and ions should be scaled; i.e. reduced by a factor of about 0.7. In several examples, including the interaction between Na (+) ions, which is of interest for ion-channel simulations, and the dynamics of an important salt-bride in Cytochrome c Oxidase, we compared the standard non-polarizable MD simulations with MDEC simulations, and demonstrated that MDEC charge scaling procedure results in more accurate interactions. The inclusion of electronic screening for charged moieties is shown to result in significant changes in protein dynamics and can give rise to new qualitative results compared with the traditional non-polarizable force fields simulations.

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.

Protein-ion binding process on finite macromolecular concentration. A Poisson-Boltzmann and Monte Carlo study

Electrostatic interactions are one of the key driving forces for protein-ligands complexation. Different levels for the theoretical modeling of such processes are available on the literature. Most of the studies on the Molecular Biology field are performed within numerical solutions of the Poisson-Boltzmann Equation and the dielectric continuum models framework. In such dielectric continuum models, there are two pivotal questions: (a) how the protein dielectric medium should be modeled, and (b) what protocol should be used when solving this effective Hamiltonian. By means of Monte Carlo (MC) and Poisson-Boltzmann (PB) calculations, we define the applicability of the PB approach with linear and nonlinear responses for macromolecular electrostatic interactions in electrolyte solution, revealing some physical mechanisms and limitations behind it especially due the raise of both macromolecular charge and concentration out of the strong coupling regime. A discrepancy between PB and MC for binding constant shifts is shown and explained in terms of the manner PB approximates the excess chemical potentials of the ligand, and not as a consequence of the nonlinear thermal treatment and/or explicit ion-ion interactions as it could be argued. Our findings also show that the nonlinear PB predictions with a low dielectric response well reproduce the pK shifts calculations carried out with an uniform dielectric model. This confirms and completes previous results obtained by both MC and linear PB calculations.

Electronic continuum model for molecular dynamics simulations

A simple model for accounting for electronic polarization in molecular dynamics (MD) simulations is discussed. In this model, called molecular dynamics electronic continuum (MDEC), the electronic polarization is treated explicitly in terms of the electronic continuum (EC) approximation, while the nuclear dynamics is described with a fixed-charge force field. In such a force-field all atomic charges are scaled to reflect the screening effect by the electronic continuum. The MDEC model is rather similar but not equivalent to the standard nonpolarizable force-fields; the differences are discussed. Of our particular interest is the calculation of the electrostatic part of solvation energy using standard nonpolarizable MD simulations. In a low-dielectric environment, such as protein, the standard MD approach produces qualitatively wrong results. The difficulty is in mistreatment of the electronic polarizability. We show how the results can be much improved using the MDEC approach. We also show how the dielectric constant of the medium obtained in a MD simulation with nonpolarizable force-field is related to the static (total) dielectric constant, which includes both the nuclear and electronic relaxation effects. Using the MDEC model, we discuss recent calculations of dielectric constants of alcohols and alkanes, and show that the MDEC results are comparable with those obtained with the polarizable Drude oscillator model. The applicability of the method to calculations of dielectric properties of proteins is discussed.

Binding of tetracyclines to elongation factor Tu, the Tet repressor, and the ribosome: a molecular dynamics simulation study

Tetracycline (Tc) is a broad-spectrum antibiotic that kills bacteria by interrupting protein biosynthesis. It is thought that the bacteriostatic action of Tc is associated with its binding to the acceptor site (or A site) in the bacterial ribosome, interfering with the attachment of aminoacyl-tRNA. Recently, however, the crystal structure of a complex between Tc and trypsin-modified elongation factor Tu (tm-EF-Tu) was determined, raising the question of whether Tc binding to EF-Tu has a role in its inhibition of protein synthesis. We address this question using computer simulations. As controls, we first compute relative ribosome binding free energies for seven Tc variants for which experimental data are available, obtaining good agreement. We then consider the binding of Tc to both the trypsin-modified and unmodified EF-Tu-GDP complexes. We show that the direct contribution of EF-Tu to the binding free energy is negligible; rather, the binding can be solely attributed to interactions of Tc with a bridging Mg(2+) ion and the GDP phosphate groups. The effects of trypsin modification are modest. Further, our calculations show that EF-Tu does not exhibit any binding preference for Tc over the nonantibiotic, 4-dedimethyl-Tc, and EF-Tu does not bind the Tc analogue tigecycline, which is a potent antibiotic. In contrast, both the ribosome and the Tet Repressor protein (involved in Tc resistance) do show a binding preference for Tc over 4-dedimethyl-Tc, and the ribosome prefers to bind tigecycline over Tc. Overall, our results provide insights into the binding properties of tetracyclines and support the idea that EF-Tu is not their primary target.

Computer simulations reveal a novel nucleotide-type binding orientation for ellipticine-based anticancer c-kit kinase inhibitors

Receptor tyrosine kinase (RTK) enzymes regulate cell signaling pathways and so are an important target for cancer chemotherapy. Current inhibitors of c-kit, a key RTK stem cell factor receptor, are inactive against the most common mutated variant Asp816Val, associated with highly malignant cancers. Recent combined experimental/simulation work has highlighted the utility of the ellipticine pharmacore in inhibiting mutant c-kit, and the present simulation study applies a combination of high-level simulation tools to probe further the binding of ellipticine-based derivatives to c-kit. We find a large preference for protonation of bound ellipticine, which stabilizes the negative protein residues that coordinated ADP.Mg (2+) in the native complex. The resulting ellipticine inhibitor binding mode resembles the native nucleotide complex and serves to explain some existing experimental data on binding specificities, indicating that functionalization at the C4/C5 sites of ellipticine derivatives may be important for the design of novel nucleotide analogues that inhibit mutant c-kit.

Probing electrostatic interactions and ligand binding in aspartyl-tRNA synthetase through site-directed mutagenesis and computer simulations

Faithful genetic code translation requires that each aminoacyl-tRNA synthetase recognise its cognate amino acid ligand specifically. Aspartyl-tRNA synthetase (AspRS) distinguishes between its negatively-charged Asp substrate and two competitors, neutral Asn and di-negative succinate, using a complex network of electrostatic interactions. Here, we used molecular dynamics simulations and site-directed mutagenesis experiments to probe these interactions further. We attempt to decrease the Asp/Asn binding free energy difference via single, double and triple mutations that reduce the net positive charge in the active site of Escherichia coli AspRS. Earlier, Glutamine 199 was changed to a negatively-charged glutamate, giving a computed reduction in Asp affinity in good agreement with experiment. Here, Lysine 198 was changed to a neutral leucine; then, Lys198 and Gln199 were mutated simultaneously. Both mutants are predicted to have reduced Asp binding and improved Asn binding, but the changes are insufficient to overcome the initial, high specificity of the native enzyme, which retains a preference for Asp. Probing the aminoacyl-adenylation reaction through pyrophosphate exchange experiments, we found no detectable activity for the mutant enzymes, indicating weaker Asp binding and/or poorer transition state stabilization. The simulations show that the mutations' effect is partly offset by proton uptake by a nearby histidine. Therefore, we performed additional simulations where the nearby Histidines 448 and 449 were mutated to neutral or negative residues: (Lys198Leu, His448Gln, His449Gln), and (Lys198Leu, His448Glu, His449Gln). This led to unexpected conformational changes and loss of active site preorganization, suggesting that the AspRS active site has a limited structural tolerance for electrostatic modifications. The data give insights into the complex electrostatic network in the AspRS active site and illustrate the difficulty in engineering charged-to-neutral changes of the preferred ligand.

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.

Role of electrostatic interactions in the assembly of empty spherical viral capsids

We examine the role of electrostatic interactions in the assembly of empty spherical viral capsids. The charges on the protein subunits that make the viral capsid mutually interact and are expected to yield electrostatic repulsion acting against the assembly of capsids. Thus, attractive protein-protein interactions of nonelectrostatic origin must act to enable the capsid formation. We investigate whether the interplay of repulsive electrostatic and attractive interactions between the protein subunits can result in the formation of spherical viral capsids of a preferred radius. For this to be the case, we find that the attractive interactions must depend on the angle between the neighboring protein subunits (i.e., on the mean curvature of the viral capsid) so that a particular angle(s) is (are) preferred energywise. Our results for the electrostatic contributions to energetics of viral capsids nicely correlate with recent experimental determinations of the energetics of protein-protein contacts in the hepatitis B virus [P. Ceres A. Zlotnick, Biochemistry 41, 11525 (2002)].

Ammonium Scanning in an Enzyme Active Site

D-amino acids are largely excluded from protein synthesis, yet they are of great interest in biotechnology. Aspartyl-tRNA synthetase (AspRS) can misacylate tRNA(Asp) with D-aspartate instead of its usual substrate, L-Asp. We investigate how the preference for L-Asp arises, using molecular dynamics simulations. Asp presents a special problem, having pseudosymmetry broken only by its ammonium group, and AspRS must protect not only against D-Asp, but against an "inverted" orientation where the two substrate carboxylates are swapped. We compare L-Asp and D-Asp, in both orientations, and succinate, where the ammonium group is removed and the ligand has an additional negative charge. All possible ammonium positions on the ligand are thus scanned, providing information on electrostatic interactions. As controls, we simulate a Q199E mutation, obtaining a reduction in binding free energy in agreement with experiment, and we simulate TyrRS, which can misacylate tRNA(Tyr) with D-Tyr. For both TyrRS and AspRS, we obtain a moderate binding free energy difference DeltaDeltaG between the L- and D-amino acids, in agreement with their known ability to misacylate their tRNAs. In contrast, we predict that AspRS is strongly protected against inverted L-Asp binding. For succinate, kinetic measurements reveal a DeltaDeltaG of over 5 kcal/mol, favoring L-Asp. The simulations show how chiral discriminations arises from the structures, with two AspRS conformations acting in different ways and proton uptake by nearby histidines playing a role. A complex network of charges protects AspRS against most binding errors, making the engineering of its specificity a difficult challenge.

Attraction of iodide ions by the free water surface, revealed by simulations with a polarizable force field based on Drude oscillators

Recent theoretical and experimental studies have shown that polarizable anions, such as iodide and bromide, preferentially accumulate close to the surface of electrolyte solutions. This finding is in sharp contrast to the previously prevailing idea that salts are dielectrically excluded from the free water surface and opens up new avenues for research in specific salt effects. In this work, we have verified the ability of a recently introduced polarizable water model, SWM4-DP, to reproduce this behavior, by simulations of a NaI/water slab, corresponding to a 1.2 M solution. The water and ion polarizabilities are modeled by classical Drude oscillator particles. As revealed by the simulations, a double layer is formed close to the free water surface, with the iodide ions located closer to the interface and the sodium ions at a neighboring, interior layer. Near the surface, all solution species acquire an induced dipole moment, that is perpendicular to the surface and points toward the exterior. The double charge layer causes ordering of water at a subsurface region. Simulations with a simpler system of a single iodide ion in a water slab show that the surface position is stabilized by induced charge interactions; in contrast, the charge-dipole interactions between the iodide permanent charge and the water permanent dipole moment favor the bulk position. Thus, the polarizabilities of ion and water are essential for explaining the increased preference of iodide for the air-water interface, in accordance with other studies.

Monte Carlo and modified Tanford-Kirkwood results for macromolecular electrostatics calculations

The understanding of electrostatic interactions is an essential aspect of the complex correlation between structure and function of biological macromolecules. It is also important in protein engineering and design. Theoretical studies of such interactions are predominantly done within the framework of Debye-Hückel theory. A classical example is the Tanford-Kirkwood (TK) model. Besides other limitations, this model assumes an infinitesimally small macromolecule concentration. By comparison to Monte Carlo (MC) simulations, it is shown that TK predictions for the shifts in ion binding constants upon addition of salt become less reliable even at moderately macromolecular concentrations. A simple modification based on colloidal literature is suggested to the TK scheme. The modified TK models suggested here satisfactorily predict MC and experimental shifts in the calcium binding constant as a function of protein concentration for the calbindin D(9k) mutant and calmodulin.

Protonation patterns in tetracycline:tet repressor recognition: simulations and experiments

Resistance to the antibiotic tetracycline (Tc) is regulated by its binding as a Tc:Mg2+ complex to the Tet Repressor protein (TetR). Tc:TetR recognition is a complex problem, with the protein and ligand each having several possible conformations and protonation states, which are difficult to elucidate by experiment alone. We used a combination of free-energy simulations and crystallographic analysis to investigate the electrostatic interactions between protein and ligand and the possible role of induced fit in Tc binding. Tc in solution was described quantum mechanically, while Tc:TetR interactions were described by a recent, high-quality molecular-mechanics model. The orientations of the amide and imidazole groups were determined experimentally by a careful analysis of Debye-Waller factors in alternate crystallographic models. The agreement with experiment for these orientations suggested that the simulations and their more detailed, thermodynamic predictions were reliable. We found that the ligand prefers an extended, zwitterionic state both in solution and in complexation with the protein. Tc is thus preorganized for binding, while the protein combines lock-and-key behavior for regions close to the ligand's amide, enolate, and ammonium groups, with an induced fit for regions close to the Mg2+ ion. These insights and the modeling techniques employed should be of interest for engineering improved TetR ligands and improved TetR proteins for gene regulation, as well as for drug design.

Computational sidechain placement and protein mutagenesis with implicit solvent models

Structure prediction and computational protein design should benefit from accurate solvent models. We have applied implicit solvent models to two problems that are central to this area. First, we performed sidechain placement for 29 proteins, using a solvent model that combines a screened Coulomb term with an Accessible Surface Area term (CASA model). With optimized parameters, the prediction quality is comparable with earlier work that omitted electrostatics and solvation altogether. Second, we computed the stability changes associated with point mutations involving ionized sidechains. For over 1000 mutations, including many fully or partly buried positions, we compared CASA and two generalized Born models (GB) with a more accurate model, which solves the Poisson equation of continuum electrostatics numerically. CASA predicts the correct sign and order of magnitude of the stability change for 81% of the mutations, compared to 97% with the best GB. We also considered 140 mutations for which experimental data are available. Comparing to experiment requires additional assumptions about the unfolded protein structure, protein relaxation in response to the mutations, and contributions from the hydrophobic effect. With a simple, commonly-used unfolded state model, the mean unsigned error is 2.1 kcal/mol with both CASA and the best GB. Overall, the electrostatic model is not important for sidechain placement; CASA and GB are equivalent for surface mutations, while GB is far superior for fully or partly buried positions. Thus, for problems like protein design that involve all these aspects, the most recent GB models represent an important step forward. Along with the recent discovery of efficient, pairwise implementations of GB, this will open new possibilities for the computational engineering of proteins.

Recognition of ribonuclease A by 3'-5'-pyrophosphate-linked dinucleotide inhibitors: a molecular dynamics/continuum electrostatics analysis

The proteins of the pancreatic ribonuclease A (RNase A) family catalyze the cleavage of the RNA polymer chain. The development of RNase inhibitors is of significant interest, as some of these compounds may have a therapeutic effect in pathological conditions associated with these proteins. The most potent low molecular weight inhibitor of RNase reported to date is the compound 5'-phospho-2'-deoxyuridine-3-pyrophosphate (P-->5)-adenosine-3-phosphate (pdUppA-3'-p). The 3',5'-pyrophosphate group of this compound increases its affinity and introduces structural features which seem to be unique in pyrophosphate-containing ligands bound to RNase A, such as the adoption of a syn conformation by the adenosine base at RNase subsite B(2) and the placement of the 5'-beta-phosphate of the adenylate (instead of the alpha-phosphate) at subsite P(1) where the phosphodiester bond cleavage occurs. In this work, we study by multi-ns molecular dynamics simulations the structural properties of RNase A complexes with the ligand pdUppA-3'-p and the related weaker inhibitor dUppA, which lacks the 3' and 5' terminal phosphate groups of pdUppA-3'-p. The simulations show that the adenylate 5'-beta-phosphate binding position and the adenosine syn orientation constitute robust structural features in both complexes, stabilized by persistent interactions with specific active-site residues of subsites P(1) and B(2). The simulation structures are used in conjunction with a continuum-electrostatics (Poisson-Boltzmann) model, to evaluate the relative binding affinity of the two complexes. The computed relative affinity of pdUppA-3'-p varies between -7.9 kcal/mol and -2.8 kcal/mol for a range of protein/ligand dielectric constants (epsilon(p)) 2-20, in good agreement with the experimental value (-3.6 kcal/mol); the agreement becomes exact with epsilon(p) = 8. The success of the continuum-electrostatics model suggests that the differences in affinity of the two ligands originate mainly from electrostatic interactions. A residue decomposition of the electrostatic free energies shows that the terminal phosphate groups of pdUppA-3'-p make increased interactions with residues Lys(7) and Lys(66) of the more remote sites P(2) and P(0), and His(119) of site P(1).