Calculation of protonation patterns in proteins with structural relaxation and molecular ensembles – application to the photosynthetic reaction center

The conventional method to determine protonation patterns of proteins was extended by explicit consideration of structural relaxation. The inclusion of structural relaxation was achieved by alternating energy minimization with the calculation of protonation pattern in an iterative manner until consistency of minimized structure and protonation pattern was reached. We applied this method to the bacterial photosynthetic reaction center (bRC) of Rps. viridis and could show that the relaxation procedure accounts for the nuclear polarization and therefore allows one to lower the dielectric constant for the protein from the typically chosen value of ɛ _p = 4 to a value of ɛ _p = 2 without fundamentally changing the results. Owing to the lower dielectric shielding at ɛ _p = 2, the charges of the titratable groups interact more strongly, which leads to sampling problems during Monte Carlo titration. We solved this problem by introducing triple moves in addition to the conventional single and double moves. We also present a new method that considers ensembles of protein conformations for the calculation of protonation patterns. Our method was successfully applied to calculate the redox potential differences of the quinones in the bRC using the relaxed structures for the different redox states of the quinones. Received: 30 October 1997 / Revised version: 2 March 1998 / Accepted: 7 March 1998     

Calculating pH-dependent free energy of proteins by using Monte Carlo protonation probabilities of ionizable residues

Protein folding, stability, and function are usually influenced by pH. And free energy plays a fundamental role in analysis of such pH-dependent properties. Electrostatics-based theoretical framework using dielectric solvent continuum model and solving Poisson-Boltzmann equation numerically has been shown to be very successful in understanding the pH-dependent properties. However, in this approach the exact computation of pH-dependent free energy becomes impractical for proteins possessing more than several tens of ionizable sites (e.g.>30), because exact evaluation of the partition function requires a summation over a vast number of possible protonation microstates. Here we present a method which computes the free energy using the average energy and the protonation probabilities of ionizable sites obtained by the well-established Monte Carlo sampling procedure. The key feature is to calculate the entropy by using the protonation probabilities. We used this method to examine a well-studied protein (lysozyme) and produced results which agree very well with the exact calculations. Applications to the optimum pH of maximal stability of proteins and protein–DNA interactions have also resulted in good agreement with experimental data. These examples recommend our method for application to the elucidation of the pH-dependent properties of proteins.     

Electrostatic interactions play an essential role in DNA repair and cold-adaptation of Uracil DNA glycosylase

Life has adapted to most environments on earth, including low and high temperature niches. The increased catalytic efficiency and thermoliability observed for enzymes from organisms living in constantly cold regions when compared to their mesophilic and thermophilic cousins are poorly understood at the molecular level. Uracil DNA glycosylase (UNG) from cod (cUNG) catalyzes removal of uracil from DNA with an increased k_cat and reduced K_m relative to its warm-active human (hUNG) counterpart. Specific issues related to DNA repair and substrate binding/recognition (K_m) are here investigated by continuum electrostatics calculations, MD simulations and free energy calculations. Continuum electrostatic calculations reveal that cUNG has surface potentials that are more complementary to the DNA potential at and around the catalytic site when compared to hUNG, indicating improved substrate binding. Comparative MD simulations combined with free energy calculations using the molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) method show that large opposing energies are involved when forming the enzyme-substrate complexes. Furthermore, the binding free energies obtained reveal that the Michaelis-Menten complex is more stable for cUNG, primarily due to enhanced electrostatic properties, suggesting that energetic fine-tuning of electrostatics can be utilized for enzymatic temperature adaptation. Energy decomposition pinpoints the residual determinants responsible for this adaptation. Figure Electrostatic isosurfaces of cod uracil DNA glycosylase in complex with double stranded DNA     

Computational analysis of photosynthetic systems

The use of various computational techniques for the study of photosynthetic systems is described ranging from genome analysis to density functional simulations of the oxygen evolving complex of PSII. The use of simulations for analyzing protein structures can aid in clarifying ambiguous and incomplete experimental results to identifying underlying rules to create efficient light-initiated charge separation at high efficiency.     

Acetylcholinesterase: mechanisms of covalent inhibition of H447I mutant determined by computational analyses

The reaction mechanisms of two inhibitor TFK(+) and TFK(0) binding to H447I mutant mouse acetylcholinesterase (mAChE) have been investigated by using a combined ab initio quantum mechanical/molecular mechanical (QM/MM) approach and classical molecular dynamics (MD) simulations. TFK(+) binding to the H447I mutant may proceed with a different reaction mechanism from the wild-type. A water molecule takes over the role of His447 and participates in the bond breaking and forming as a "charge relayer". Unlike in the wild-type mAChE case, Glu334, a conserved residue from the catalytic triad, acts as a catalytic base in the reaction. The calculated energy barrier for this reaction is about 8kcal/mol. These predictions await experimental verification. In the case of the neutral ligand TFK(0), however, multiple MD simulations on the TFK(0)/H447I complex reveal that none of the water molecules can be retained in the active site as a "catalytic" water. Taken together our computational studies confirm that TFK(0) is almost inactive in the H447I mutant, and also provide detailed mechanistic insights into the experimental observations.     

Energetics of charge-charge interactions in proteins

Electrostatic interactions between pairs of atoms in proteins are calculated with a model based on the linearized Poisson-Boltzmann equation. The equation is solved accurately by a method that takes into account the detailed shape of the protein. This paper presents applications to several systems. Experimental data for the interaction of ionized residues with an active site histidine in subtilisin BPN' allow the model to be tested, using various assumptions for the electrical properties of the protein and solvent. The electrostatic stabilization of the active site thiolate of rhodanese is analyzed, with attention to the influence of alpha-helices. Finally, relationships between electrostatic potential and charge-charge distance are reported for large and small globular proteins. The above results are compared with those of simpler electrostatic models, including Coulomb's law with both a distance-dependent dielectric constant (epsilon = R) and a fixed dielectric constant (epsilon = 2), and Tanford-Kirkwood theory. The primary conclusions are as follows: 1) The Poisson-Boltzmann model agrees with the subtilisin data over a range of ionic strengths; 2) two alpha-helices generate a large potential in the active site of rhodanese; 3) epsilon = R overestimates weak electrostatic interactions but yields relatively good results for strong ones; 4) Tanford-Kirkwood theory is a useful approximation to detailed solutions of the linearized Poisson-Boltzmann equation in globular proteins; and 5) the modified Tanford-Kirkwood theory over-screens the measured electrostatic interactions in subtilisin.     

Dielectric relaxation in hydroxyethyl cellulose

This paper presents the results of the dielectric analysis of hydroxyethyl cellulose (HEC) in the temperature range 100–350 K and in the frequency range 20 Hz–1 MHz. The results show a distinct broad relaxation process in the temperature range 150–250 K with activation energy of about 33.4 MJ/kmol. The strength of the relaxation, ε_s−ε∞, increases only slightly with temperature and the apparent increase of the height of the relaxation maximum is attributed to an increase of the co-operativity parameter . This increase of is interpreted as decrease of the co-operativity of dipole motions with increasing temperature. The origin of the relaxation is ascribed to the reorientation of double ethylene oxide groups or part of these.     

A comparison of the electromechanical properties of structurally diverse proteins by molecular dynamics simulation

Proteins are subjected to electric fields both within the cell and during routine biochemical analysis. We have used atomistic molecular dynamics simulations to study conformational changes within three structurally diverse proteins subjected to high electric fields. At electric fields in excess of .5?V/nm, major structural changes were observed in all three proteins due to charge redistribution within the biomolecule. However, the electromechanical resilience was found to be highly dependent on the protein secondary structure, with α-helices showing a particularly high susceptibility to deformation by the applied electric field.     

Focusing of electric fields in the active site of Cu-Zn superoxide dismutase: effects of ionic strength and amino-acid modification

In this paper we report the implementation of a finite-difference algorithm which solves the linearized Poisson-Boltzmann equation for molecules of arbitrary shape and charge distribution and which includes the screening effects of electrolytes. The microcoding of the algorithm on an ST-100 array processor allows us to obtain electrostatic potential maps in and around a protein, including the effects of ionic strength, in about 30 minutes. We have applied the algorithm to a dimer of the protein Cu-Zn superoxide dismutase (SOD) and compared our results to those obtained from uniform dielectric models based on coulombic potentials. We find that both the shape of the protein-solvent boundary and the ionic strength of the solvent have a profound effect on the potentials in the solvent. For the case of SOD, the cluster of positive charge at the bottom of the active site channel produces a strongly enhanced positive potential due to the focusing of field lines in the channel-a result that cannot be obtained with any uniform dielectric model. The remainder of the protein is surrounded by a weak negative potential. The electrostatic potential of the enzyme seems designed to provide a large cross-sectional area for productive collisions. Based on the ionic strength dependence of the size of the positive potential region emanating from the active site and the repulsive negative potential barrier surrounding the protein, we are able to suggest an explanation for the ionic strength dependence of the activity of the native and chemically modified forms of the enzyme.     

Molecular dynamics simulations of membrane proteins

Membrane proteins control the traffic across cell membranes and thereby play an essential role in cell function from transport of various solutes to immune response via molecular recognition. Because it is very difficult to determine the structures of membrane proteins experimentally, computational methods have been increasingly used to study their structure and function. Here we focus on two classes of membrane proteins—ion channels and transporters—which are responsible for the generation of action potentials in nerves, muscles, and other excitable cells. We describe how computational methods have been used to construct models for these proteins and to study the transport mechanism. The main computational tool is the molecular dynamics (MD) simulation, which can be used for everything from refinement of protein structures to free energy calculations of transport processes. We illustrate with specific examples from gramicidin and potassium channels and aspartate transporters how the function of these membrane proteins can be investigated using MD simulations.     

Electrostatics of the FeS clusters in respiratory complex I

Respiratory complex I couples the transfer of electrons from NADH to ubiquinone and the translocation of protons across the mitochondrial membrane. A detailed understanding of the midpoint reduction potentials (E_m) of each redox center and the factors which influence those potentials are critical in the elucidation of the mechanism of electron transfer in this enzyme. We present accurate electrostatic interaction energies for the iron-sulfur (FeS) clusters of complex I to facilitate the development of models and the interpretation of experiments in connection to electron transfer (ET) in this enzyme. To calculate redox titration curves for the FeS clusters it is necessary to include interactions between clusters, which in turn can be used to refine E_m values and validate spectroscopic assignments of each cluster. Calculated titration curves for clusters N4, N5, and N6a are discussed. Furthermore, we present some initial findings on the electrostatics of the redox centers of complex I under the influence of externally applied membrane potentials. A means of determining the location of the FeS cofactors within the holo-complex based on electrostatic arguments is proposed. A simple electrostatic model of the protein/membrane system is examined to illustrate the viability of our hypothesis.     

Automated computational framework for the analysis of electrostatic similarities of proteins

Charge plays an important role in protein-protein interactions. In the case of excessively charged proteins, their electrostatic potentials contribute to the processes of recognition and binding with other proteins or ligands. We present an automated computational framework for determining the contribution of each charged amino acid to the electrostatic properties of proteins, at atomic resolution level. This framework involves computational alanine scans, calculation of Poisson-Boltzmann electrostatic potentials, calculation of electrostatic similarity distances (ESDs), hierarchical clustering analysis of ESDs, calculation of solvation free energies of association, and visualization of the spatial distributions of electrostatic potentials. The framework is useful to classify families of mutants with similar electrostatic properties and to compare them with the parent proteins in the complex. The alanine scan mutants introduce perturbations in the local electrostatic properties of the proteins and aim in delineating the contribution of each mutated amino acid in the spatial distribution of electrostatic potential, and in biological function when electrostatics is a dominant contributing factor in protein-protein interactions. The framework can be used to design new proteins with tailored electrostatic properties, such as immune system regulators, inhibitors, and vaccines, and in guiding experimental studies. We present an example for the interaction of the immune system protein C3d (the d-fragment of complement protein C3) with its receptor CR2, and we discuss our data in view of a binding site controversy.     

Conformational selection through electrostatics: Free energy simulations of GTP and GDP binding to archaeal initiation factor 2

Archaeal Initiation Factor 2 is a GTPase involved in protein biosynthesis. In its GTP-bound, "ON" conformation, it binds an initiator tRNA and carries it to the ribosome. In its GDP-bound, "OFF" conformation, it dissociates from tRNA. To understand the specific binding of GTP and GDP and their dependence on the conformational state, molecular dynamics free energy simulations were performed. The ON state specificity was predicted to be weak, with a GTP/GDP binding free energy difference of -1 kcal/mol, favoring GTP. The OFF state specificity is larger, 4 kcal/mol, favoring GDP. The overall effects result from a competition among many interactions in several complexes. To interpret them, we use a simpler, dielectric continuum model. Several effects are robust with respect to the model details. Both nucleotides have a net negative charge, so that removing them from solvent into the binding pocket carries a desolvation penalty, which is large for the ON state, and strongly disfavors GTP binding compared to GDP. Short-range interactions between the additional GTP phosphate group and ionized sidechains in the binding pocket offset most, but not all of the desolvation penalty; more distant groups also contribute significantly, and the switch 1 loop only slightly. The desolvation penalty is lower for the more open, wetter OFF state, and the GTP/GDP difference much smaller. Short-range interactions in the binding pocket and with more distant groups again make a significant contribution. Overall, the simulations help explain how conformational selection is achieved with a single phosphate group.     

Spatial optimization of electrostatic interactions between the ionized groups in globular proteins

A model approach is suggested to estimate the degree of spatial optimization of the electrostatic interactions in protein molecules. The method is tested on a set of 44 globular proteins, representative of the available crystallographic data. The theoretical model is based on macroscopic computation of the contribution of charge–charge interactions to the electrostatic term of the free energy for the native proteins and for a big number of virtual structures with randomly distributed on protein surface charge consetellations (generated by a Monte-Carlo technique). The statistical probability of occurrence of random structures with electrostatic energies lower than the energy of the native protein is suggested as a criterion for spatial optimization of the electrostatic interactions. The results support the hypothesis that the folding process optimizes the stabilizing effect of electrostatic interactions, but to very different degree for different proteins. A parallel analysis of ion pairs shows that the optimization of the electrostatic term in globular proteins has increasingly gone in the direction of rejecting the repulsive short contacts between charges of equal sign than of creating of more salt bridges (in comparison with the statistically expected number of shortrange ion pairs in the simulated random structures). It is observed that the decrease in the spatial optimization of the electrostatic interactions is usually compensated for by an appearance of disulfide bridges in the covalent structure of the examined proteins. © 1994 Wiley-Liss, Inc.     

A single amino acid substitution in the active site of Escherichia coli aspartate transcarbamoylase prevents the allosteric transition

Modeling of the tetrahedral intermediate within the active site of Escherichia coli aspartate transcarbamoylase revealed a specific interaction with the side-chain of Gln137, an interaction not previously observed in the structure of the X-ray enzyme in the presence of N-phosphonacetyl-L-aspartate (PALA). Previous site-specific mutagenesis experiments showed that when Gln137 was replaced by alanine, the resulting mutant enzyme (Q137A) exhibited approximately 50-fold less activity than the wild-type enzyme, exhibited no homotropic cooperativity, and the binding of both carbamoyl phosphate and aspartate were extremely compromised. To elucidate the structural alterations in the mutant enzyme that might lead to such pronounced changes in kinetic and binding properties, the Q137A enzyme was studied by time-resolved, small-angle X-ray scattering and its structure was determined in the presence of PALA to 2.7 angstroms resolution. Time-resolved, small-angle X-ray scattering established that the natural substrates, carbamoyl phosphate and L-aspartate, do not induce in the Q137A enzyme the same conformational changes as observed for the wild-type enzyme, although the scattering pattern of the Q137A and wild-type enzymes in the presence of PALA were identical. The overall structure of the Q137A enzyme is similar to that of the R-state structure of wild-type enzyme with PALA bound. However, there are differences in the manner by which the Q137A enzyme coordinates PALA, especially in the side-chain positions of Arg105 and His134. The replacement of Gln137 by Ala also has a dramatic effect on the electrostatics of the active site. These data taken together suggest that the side-chain of Gln137 in the wild-type enzyme is required for the binding of carbamoyl phosphate in the proper orientation so as to induce conformational changes required for the creation of the high-affinity aspartate-binding site. The inability of carbamoyl phosphate to create the high-affinity binding site in the Q137A enzyme results in an enzyme locked in the low-activity low-affinity T state. These results emphasize the absolute requirement of the binding of carbamoyl phosphate for the creation of the high-affinity aspartate-binding site and for inducing the homotropic cooperativity in aspartate transcarbamoylase.     

Thermodynamic and electrostatic properties of ternary Oxytricha nova TEBP-DNA complex

Telomeres constitute the nucleoprotein ends of eukaryotic chromosomes which are essential for their proper function. Telomere end binding protein (TEBP) from Oxytricha nova was among the first telomeric proteins, which were well characterized biologically. TEBP consists of two protein subunits (alpha, beta) and forms a ternary complex with single stranded telomeric DNA containing tandem repeats TTTTGGGG. This work presents the characterization of the thermodynamic and electrostatic properties of this complex by computational chemistry methods (continuum Poisson-Boltzmann and solvent accessible surface calculations). Our calculations give a new insight into molecular properties of studied system. Based on the thermodynamic analysis we provide a rationale for the experimental observation that alpha and ssDNA forms a binary complex and the beta subunit joins alpha:ssDNA complex only after the latter is formed. Calculations of distribution of the molecular electrostatic potential for protein subunits alone and for all possible binary complexes revealed the important role of the "guiding funnel" potential generated by alpha:ssDNA complex. This potential may help the beta subunit to dock to the already formed alpha:DNA intermediate in highly steric and electrostatic favorable manner. Our pK(a) calculations of TEBP are able to explain the experimental mobility shifts of the complex in electrophoretic non-denaturating gels.     

Interplay between the electrostatic membrane potential and conformational changes in membrane proteins

Transmembrane electrostatic membrane potential is a major energy source of the cell. Importantly, it determines the structure as well as function of charge‐carrying membrane proteins. Here, we discuss the relationship between membrane potential and membrane proteins, in particular whether the conformation of these proteins is integrally connected to the membrane potential. Together, these concepts provide a framework for rationalizing the types of conformational changes that have been observed in membrane proteins and for better understanding the electrostatic effects of the membrane potential on both reversible as well as unidirectional dynamic processes of membrane proteins.     

A Comparison of Seven Fast but Approximate Methods to Compute the Free Energy of Deprotonation for Amino Acids in Aqueous Solution

Abstract

In recent years, a variety of methods based on statistical mechanics have been successfully applied to calculate free energy differences of chemical reactions from molecular simulation. The accuracy and computational efficiency vary strongly between these methods. Seven approximate but fast methods to calculate free energy differences are compared in terms of accuracy and efficiency with the accurate but expensive thermodynamic integration method as reference, using 28 protonation and deprotonation reactions of aspartic acid in aqueous solution as test cases. At least two simulations are required to obtain an accurate free energy difference between two states of the system. Both, the averaged one-step perturbation method and the linear response method yield the most accurate results, while the latter method shows the fastest convergence.     

Kinase conformations: a computational study of the effect of ligand binding.

Protein function is often controlled by ligand-induced conformational transitions. Yet, in spite of the increasing number of three-dimensional crystal structures of proteins in different conformations, not much is known about the driving forces of these transitions. As an initial step toward exploring the conformational and energetic landscape of protein kinases by computational methods, intramolecular energies and hydration free energies were calculated for different conformations of the catalytic domain of cAMP-dependent protein kinase (cAPK) with a continuum (Poisson) model for the electrostatics. Three protein kinase crystal structures for ternary complexes of cAPK with the peptide inhibitor PKI(5-24) and ATP or AMP-PNP were modeled into idealized intermediate and open conformations. Concordant with experimental observation, we find that the binding of PKI(5-24) is more effective in stabilizing the closed and intermediate forms of cAPK than ATP. PKI(5-24) seems to drive the final closure of the active site cleft from intermediate to closed state because ATP does not distinguish between these two states. Binding of PKI(5-24) and ATP is energetically additive.     

Insights into the local residual entropy of proteins provided by NMR relaxation.

A simple model is used to illustrate the relationship between the dynamics measured by NMR relaxation methods and the local residual entropy of proteins. The expected local dynamic behavior of well-packed extended amino acid side chains are described by employing a one-dimensional vibrator that encapsulates both the spatial and temporal character of the motion. This model is then related to entropy and to the generalized order parameter of the popular "model-free" treatment often used in the analysis of NMR relaxation data. Simulations indicate that order parameters observed for the methyl symmetry axes in, for example, human ubiquitin correspond to significant local entropies. These observations have obvious significance for the issue of the physical basis of protein structure, dynamics, and stability.