A self-consistent, microenvironment modulated screened coulomb potential approximation to calculate pH-dependent electrostatic effects in proteins

Structural Insights of PD-1/PD-L1 Axis: An In silico Approach

BACKGROUND

Interaction of PD-1 protein (present on immune T-cell) with its ligand PD-L1 (over-expressed on cancerous cell) makes the cancerous cell survive and thrive. The association of PD-1/PD-L1 represents a classical protein-protein interaction (PPI), where receptor and ligand binding through a large flat surface. Blocking the PD-1/PDL-1 complex formation can restore the normal immune mechanism, thereby destroying cancerous cells. However, the PD-1/PDL1 interactions are only partially characterized.

OBJECTIVE

We aim to comprehend the time-dependent behavior of PD-1 upon its binding with PD-L1.

METHODS

The current work focuses on a molecular dynamics simulation (MDs) simulation study of apo and ligand bound PD-1.

RESULTS

Our simulation reveals the flexible nature of the PD-1, both in apo and bound form. Moreover, the current study also differentiates the type of strong and weak interactions which could be targeted to overcome the complex formation.

CONCLUSION

The current article could provide a valuable structural insight about the target protein (PD-1) and its ligand (PD-L1) which could open new opportunities in developing small molecule inhibitors (SMIs) targeting either PD-1 or PD-L1.

Modulation of free energy landscapes as a strategy for the design of antimicrobial peptides

Computational design of antimicrobial peptides (AMPs) is a promising area of research for developing novel agents against drug-resistant bacteria. AMPs are present naturally in many organisms, from bacteria to humans, a time-tested mechanism that makes them attractive as effective antibiotics. Depending on the environment, AMPs can exhibit α-helical or β-sheet conformations, a mix of both, or lack secondary structure; they can be linear or cyclic. Prediction of their structures is challenging but critical for rational design. Promising AMP leads can be developed using essentially two approaches: traditional modeling of the physicochemical mechanisms that determine peptide behavior in aqueous and membrane environments and knowledge-based, e.g., machine learning (ML) techniques, that exploit ever-growing AMP databases. Here, we explore the conformational landscapes of two recently ML-designed AMPs, characterize the dependence of these landscapes on the medium conditions, and identify features in peptide and membrane landscapes that mediate protein-membrane association. For both peptides, we observe greater conformational diversity in an aqueous solvent than in a less polar solvent, and one peptide is seen to alter its conformation more dramatically than the other upon the change of solvent. Our results support the view that structural rearrangement in response to environmental changes is central to the mechanism of membrane-structure disruption by linear peptides. We expect that the design of AMPs by ML will benefit from the incorporation of peptide conformational substates as quantified here with molecular simulations.

Protein pK a Prediction with Machine Learning

Protein pK _a prediction is essential for the investigation of the pH-associated relationship between protein structure and function. In this work, we introduce a deep learning-based protein pK _a predictor DeepKa, which is trained and validated with the pK _a values derived from continuous constant-pH molecular dynamics (CpHMD) simulations of 279 soluble proteins. Here, the CpHMD implemented in the Amber molecular dynamics package has been employed (Huang Y.J. Chem. Inf. Model.2018, 58, 1372-1383). Notably, to avoid discontinuities at the boundary, grid charges are proposed to represent protein electrostatics. We show that the prediction accuracy by DeepKa is close to that by CpHMD benchmarking simulations, validating DeepKa as an efficient protein pK _a predictor. In addition, the training and validation sets created in this study can be applied to the development of machine learning-based protein pK _a predictors in the future. Finally, the grid charge representation is general and applicable to other topics, such as the protein-ligand binding affinity prediction.

Automation of some macromolecular properties using a machine learning approach

In this study, we employed a newly developed method to predict macromolecular properties using a swarm artificial neural network (ANN) method as a machine learning approach. In this method, the molecular structures are represented by the feature description vectors used as training input data for a neural network. This study aims to develop an efficient approach for training an ANN using either experimental or quantum mechanics data. We aim to introduce an error model controlling the reliability of the prediction confidence interval using a bootstrapping swarm approach. We created different datasets of selected experimental or quantum mechanics results. Using this optimized ANN, we hope to predict properties and their statistical errors for new molecules. There are four datasets used in this study. That includes the dataset of 642 small organic molecules with known experimental hydration free energies, the dataset of 1475 experimental pKa values of ionizable groups in 192 proteins, the dataset of 2693 mutants in 14 proteins with given experimental values of changes in the Gibbs free energy, and a dataset of 7101 quantum mechanics heat of formation calculations. All the data are prepared and optimized using the AMBER force field in the CHARMM macromolecular computer simulation program. The bootstrapping swarm ANN code for performing the optimization and prediction is written in Python computer programming language. The descriptor vectors of the small molecules are based on the Coulomb matrix and sum over bond properties. For the macromolecular systems, they consider the chemical-physical fingerprints of the region in the vicinity of each amino acid.

Selective Stabilization of Aspartic Acid Protonation State within a Given Protein Conformation Occurs via Specific "Molecular Association"

Proteins involved in proton-/electron-transfer processes often possess "functional" aspartates/aspartic acids (Asp) with variable protonation states. The mechanism of Asp protonation-deprotonation within proteins is unclear. Two questions were asked-the possible types of determinants responsible for Asp protonation-deprotonation and the spatial arrangements of the determinants leading to selective stabilization. The questions were analyzed using nine different solvent models, which scanned the complete protein dielectric range, and four protein models, which illustrated the spatial arrangements around Asp, termed as "molecular association". The methods employed were quantum chemical calculations and constant pH simulations. The types of the determinants identified were charge-charge interaction, H bonding, dipole-π interaction, extended electronic conjugation, dielectric effect, and solvent accessibility. All solvent-exposed Asp [buried fraction (BF) less than 0.5] were aspartates, and buried Asp were either aspartic acids or aspartates, each having a different "molecular association". The exposed aspartates were stabilized via a H-bonding network with bulk water, buried aspartates via salt bridge or, minimum, two intramolecular H bonds, and buried aspartic acids via, minimum, one intramolecular H bond. An "acid-alcohol pair" (involving Ser/Thr/Tyr) was a common determinant to any "functional" buried aspartate/aspartic acid. Higher energy "molecular associations" observed within proteins compared to those within water, presumably, indicated easy molecular restructuring and alteration of the Asp protonation states during a protein-mediated proton/electron transfer.

Mass Spectrometry- and Computational Structural Biology-Based Investigation of Proteins and Peptides

Recent developments of mass spectrometry (MS) allow us to identify, estimate, and characterize proteins and protein complexes. At the same time, structural biology helps to determine the protein structure and its structure-function relationship. Together, they aid to understand the protein structure, property, function, protein-complex assembly, protein-protein interaction, and dynamics. The present chapter is organized with illustrative results to demonstrate how experimental mass spectrometry can be combined with computational structural biology for detailed studies of protein's structures. We have used tumor differentiation factor protein/peptide as ligand and Hsp70/Hsp90 as receptor protein as examples to study ligand-protein interaction. To investigate possible protein conformation, we will describe two proteins-lysozyme and myoglobin. As an application of MS-based assignment of disulfide bridges, the case of the spider venom polypeptide Phα1β will also be discussed.

Hot-spot identification on a broad class of proteins and RNA suggest unifying principles of molecular recognition

Chemically diverse fragments tend to collectively bind at localized sites on proteins, which is a cornerstone of fragment-based techniques. A central question is how general are these strategies for predicting a wide variety of molecular interactions such as small molecule-protein, protein-protein and protein-nucleic acid for both experimental and computational methods. To address this issue, we recently proposed three governing principles, (1) accurate prediction of fragment-macromolecule binding free energy, (2) accurate prediction of water-macromolecule binding free energy, and (3) locating sites on a macromolecule that have high affinity for a diversity of fragments and low affinity for water. To test the generality of these concepts we used the computational technique of Simulated Annealing of Chemical Potential to design one small fragment to break the RecA-RecA protein-protein interaction and three fragments that inhibit peptide-deformylase via water-mediated multi-body interactions. Experiments confirm the predictions that 6-hydroxydopamine potently inhibits RecA and that PDF inhibition quantitatively tracks the water-mediated binding predictions. Additionally, the principles correctly predict the essential bound waters in HIV Protease, the surprisingly extensive binding site of elastase, the pinpoint location of electron transfer in dihydrofolate reductase, the HIV TAT-TAR protein-RNA interactions, and the MDM2-MDM4 differential binding to p53. The experimental confirmations of highly non-obvious predictions combined with the precise characterization of a broad range of known phenomena lend strong support to the generality of fragment-based methods for characterizing molecular recognition.

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.

Using polarizable POSSIM force field and fuzzy-border continuum solvent model to calculate pK(a) shifts of protein residues

Our Fuzzy-Border (FB) continuum solvent model has been extended and modified to produce hydration parameters for small molecules using POlarizable Simulations Second-order Interaction Model (POSSIM) framework with an average error of 0.136 kcal/mol. It was then used to compute pKa shifts for carboxylic and basic residues of the turkey ovomucoid third domain (OMTKY3) protein. The average unsigned errors in the acid and base pKa values were 0.37 and 0.4 pH units, respectively, versus 0.58 and 0.7 pH units as calculated with a previous version of polarizable protein force field and Poisson Boltzmann continuum solvent. This POSSIM/FB result is produced with explicit refitting of the hydration parameters to the pKa values of the carboxylic and basic residues of the OMTKY3 protein; thus, the values of the acidity constants can be viewed as additional fitting target data. In addition to calculating pKa shifts for the OMTKY3 residues, we have studied aspartic acid residues of Rnase Sa. This was done without any further refitting of the parameters and agreement with the experimental pKa values is within an average unsigned error of 0.65 pH units. This result included the Asp79 residue that is buried and thus has a high experimental pKa value of 7.37 units. Thus, the presented model is capable or reproducing pKa results for residues in an environment that is significantly different from the solvated protein surface used in the fitting. Therefore, the POSSIM force field and the FB continuum solvent parameters have been demonstrated to be sufficiently robust and transferable. © 2016 Wiley Periodicals, Inc.

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.

Utility of computational structural biology in mass spectrometry

Recent developments of mass spectrometry (MS) allow us to identify, estimate, and characterize proteins and protein complexes. At the same time, structural biology helps to determine the protein structure and its structure-function relationship. Together, they aid to understand the protein structure, property, function, protein-complex assembly, protein-protein interaction and dynamics. The present chapter is organized with illustrative results to demonstrate how experimental mass spectrometry can be combined with computational structural biology for detailed studies of protein's structures. We have used tumor differentiation factor protein/peptide as ligand and Hsp70/Hsp90 as receptor protein as examples to study ligand-protein interaction. To investigate possible protein conformation we will describe two proteins, lysozyme and myoglobin.

Uncovering the determinants of a highly perturbed tyrosine pKa in the active site of ketosteroid isomerase

Within the idiosyncratic enzyme active-site environment, side chain and ligand pKa values can be profoundly perturbed relative to their values in aqueous solution. Whereas structural inspection of systems has often attributed perturbed pKa values to dominant contributions from placement near charged groups or within hydrophobic pockets, Tyr57 of a Pseudomonas putida ketosteroid isomerase (KSI) mutant, suggested to have a pKa perturbed by nearly 4 units to 6.3, is situated within a solvent-exposed active site devoid of cationic side chains, metal ions, or cofactors. Extensive comparisons among 45 variants with mutations in and around the KSI active site, along with protein semisynthesis, (13)C NMR spectroscopy, absorbance spectroscopy, and X-ray crystallography, was used to unravel the basis for this perturbed Tyr pKa. The results suggest that the origin of large energetic perturbations are more complex than suggested by visual inspection. For example, the introduction of positively charged residues near Tyr57 raises its pKa rather than lowers it; this effect, and part of the increase in the Tyr pKa from the introduction of nearby anionic groups, arises from accompanying active-site structural rearrangements. Other mutations with large effects also cause structural perturbations or appear to displace a structured water molecule that is part of a stabilizing hydrogen-bond network. Our results lead to a model in which three hydrogen bonds are donated to the stabilized ionized Tyr, with these hydrogen-bond donors, two Tyr side chains, and a water molecule positioned by other side chains and by a water-mediated hydrogen-bond network. These results support the notion that large energetic effects are often the consequence of multiple stabilizing interactions rather than a single dominant interaction. Most generally, this work provides a case study for how extensive and comprehensive comparisons via site-directed mutagenesis in a tight feedback loop with structural analysis can greatly facilitate our understanding of enzyme active-site energetics. The extensive data set provided may also be a valuable resource for those wishing to extensively test computational approaches for determining enzymatic pKa values and energetic effects.

Predicting pKa for proteins using COSMO-RS

We have used the COSMO-RS implicit solvation method to calculate the equilibrium constants, pK_a, for deprotonation of the acidic residues of the ovomucoid inhibitor protein, OMTKY3. The root mean square error for comparison with experimental data is only 0.5 pH units and the maximum error 0.8 pH units. The results show that the accuracy of pK_a prediction using COSMO-RS is as good for large biomolecules as it is for smaller inorganic and organic acids and that the method compares very well to previous pK_a predictions of the OMTKY3 protein using Quantum Mechanics/Molecular Mechanics. Our approach works well for systems of about 1000 atoms or less, which makes it useful for small proteins as well as for investigating portions of larger proteins such as active sites in enzymes.

Depth: a web server to compute depth, cavity sizes, detect potential small-molecule ligand-binding cavities and predict the pKa of ionizable residues in proteins

Residue depth accurately measures burial and parameterizes local protein environment. Depth is the distance of any atom/residue to the closest bulk water. We consider the non-bulk waters to occupy cavities, whose volumes are determined using a Voronoi procedure. Our estimation of cavity sizes is statistically superior to estimates made by CASTp and VOIDOO, and on par with McVol over a data set of 40 cavities. Our calculated cavity volumes correlated best with the experimentally determined destabilization of 34 mutants from five proteins. Some of the cavities identified are capable of binding small molecule ligands. In this study, we have enhanced our depth-based predictions of binding sites by including evolutionary information. We have demonstrated that on a database (LigASite) of ∼200 proteins, we perform on par with ConCavity and better than MetaPocket 2.0. Our predictions, while less sensitive, are more specific and precise. Finally, we use depth (and other features) to predict pK_as of GLU, ASP, LYS and HIS residues. Our results produce an average error of just <1 pH unit over 60 predictions. Our simple empirical method is statistically on par with two and superior to three other methods while inferior to only one. The DEPTH server (http://mspc.bii.a-star.edu.sg/depth/) is an ideal tool for rapid yet accurate structural analyses of protein structures.

Biomolecular electrostatics and solvation: a computational perspective

An understanding of molecular interactions is essential for insight into biological systems at the molecular scale. Among the various components of molecular interactions, electrostatics are of special importance because of their long-range nature and their influence on polar or charged molecules, including water, aqueous ions, proteins, nucleic acids, carbohydrates, and membrane lipids. In particular, robust models of electrostatic interactions are essential for understanding the solvation properties of biomolecules and the effects of solvation upon biomolecular folding, binding, enzyme catalysis, and dynamics. Electrostatics, therefore, are of central importance to understanding biomolecular structure and modeling interactions within and among biological molecules. This review discusses the solvation of biomolecules with a computational biophysics view toward describing the phenomenon. While our main focus lies on the computational aspect of the models, we provide an overview of the basic elements of biomolecular solvation (e.g. solvent structure, polarization, ion binding, and non-polar behavior) in order to provide a background to understand the different types of solvation models.

Cardiac mitochondrial matrix and respiratory complex protein phosphorylation

It has become appreciated over the last several years that protein phosphorylation within the cardiac mitochondrial matrix and respiratory complexes is extensive. Given the importance of oxidative phosphorylation and the balance of energy metabolism in the heart, the potential regulatory effect of these classical signaling events on mitochondrial function is of interest. However, the functional impact of protein phosphorylation and the kinase/phosphatase system responsible for it are relatively unknown. Exceptions include the well-characterized pyruvate dehydrogenase and branched chain α-ketoacid dehydrogenase regulatory system. The first task of this review is to update the current status of protein phosphorylation detection primarily in the matrix and evaluate evidence linking these events with enzymatic function or protein processing. To manage the scope of this effort, we have focused on the pathways involved in energy metabolism. The high sensitivity of modern methods of detecting protein phosphorylation and the low specificity of many kinases suggests that detection of protein phosphorylation sites without information on the mole fraction of phosphorylation is difficult to interpret, especially in metabolic enzymes, and is likely irrelevant to function. However, several systems including protein translocation, adenine nucleotide translocase, cytochrome c, and complex IV protein phosphorylation have been well correlated with enzymatic function along with the classical dehydrogenase systems. The second task is to review the current understanding of the kinase/phosphatase system within the matrix. Though it is clear that protein phosphorylation occurs within the matrix, based on (32)P incorporation and quantitative mass spectrometry measures, the kinase/phosphatase system responsible for this process is ill-defined. An argument is presented that remnants of the much more labile bacterial protein phosphoryl transfer system may be present in the matrix and that the evaluation of this possibility will require the application of approaches developed for bacterial cell signaling to the mitochondria.

Structural intermediates in a model of the substrate translocation path of the bacterial glutamate transporter homologue GltPh

Excitatory amino acid transporters (EAATs) are membrane proteins responsible for reuptake of glutamate from the synaptic cleft to terminate neurotransmission and help prevent neurotoxically high, extracellular glutamate concentrations. Important structural information about these proteins emerged from crystal structures of GltPh, a bacterial homologue of EAATs, in conformations facing outward and inward. These remarkably different conformations are considered to be end points of the substrate translocation path (STP), suggesting that the transport mechanism involves major conformational rearrangements that remain uncharted. To investigate possible steps in the structural transitions of the STP between the two end-point conformations, we applied a combination of computational modeling methods (motion planning, molecular dynamics simulations, and mixed elastic network models). We found that the conformational changes in the transition involve mainly the repositioning the "transport domain" and the "trimerization domain" identified previously in the crystal structures. The two domains move in opposite directions along the membrane normal, and the transport domain also tilts by ∼17° with respect to this axis. Moreover, the TM3-4 loop undergoes a flexible, "restraining bar"-like conformational change with respect to the transport domain. As a consequence of these conformational rearrangements along the transition path we calculated a significant decrease of nearly 20% in the area of the transport-to-trimerization domain interface (TTDI). Water penetrates parts of the TTDI in the modeled intermediates but very much less in the end-point conformations. We show that these characteristics of the modeled intermediate states agree with experimental results from residue-accessibility studies in individual monomers and identify specific residues that can be used to test the proposed STP. Moreover, MD simulations of complete GltPh trimers constructed from initially identical monomer intermediates suggest that asymmetry can appear in the trimer, consonant with available experimental data showing independent transport kinetics by individual monomers in the trimers.

Bluues: a program for the analysis of the electrostatic properties of proteins based on generalized Born radii

BACKGROUND

The Poisson-Boltzmann (PB) equation and its linear approximation have been widely used to describe biomolecular electrostatics. Generalized Born (GB) models offer a convenient computational approximation for the more fundamental approach based on the Poisson-Boltzmann equation, and allows estimation of pairwise contributions to electrostatic effects in the molecular context.

RESULTS

We have implemented in a single program most common analyses of the electrostatic properties of proteins. The program first computes generalized Born radii, via a surface integral and then it uses generalized Born radii (using a finite radius test particle) to perform electrostaic analyses. In particular the ouput of the program entails, depending on user's requirement: 1) the generalized Born radius of each atom; 2) the electrostatic solvation free energy; 3) the electrostatic forces on each atom (currently in a developmental stage); 4) the pH-dependent properties (total charge and pH-dependent free energy of folding in the pH range -2 to 18; 5) the pKa of all ionizable groups; 6) the electrostatic potential at the surface of the molecule; 7) the electrostatic potential in a volume surrounding the molecule;

CONCLUSIONS

Although at the expense of limited flexibility the program provides most common analyses with requirement of a single input file in PQR format. The results obtained are comparable to those obtained using state-of-the-art Poisson-Boltzmann solvers. A Linux executable with example input and output files is provided as supplementary material.

A salt-bridge controlled by ligand binding modulates the hydrolysis reaction in a GH5 endoglucanase

Cellulases, distributed in at least 15 families of glycoside hydrolases, will play a key role in biomass conversion and renewable energy challenges of the future. Cel5B from Clostridium thermocellum is a β-1,4-endoglucanase and a member of family 5 of glycoside hydrolases (GH5) and is characterized by an (α/β)(8) barrel structure. In contrast to other retaining enzymes, in which the catalytic carboxylate groups (glutamate or aspartate) are positioned ≈ 5.5 Å apart to facilitate nucleophilic attack on the anomeric carbon of the sugar substrate, these two residues in Cel5B are positioned ≈ 10 Å from each other in the unliganded wild-type structure. The structure of the enzyme solved in complex with a cleavage product (cellobiose) revealed ligand-induced conformational changes in the loop carrying Glu140 (proton donor). The reorientation of Glu140 in the complex reduces the separation of the catalytic glutamate residues to 4.3 Å. In this study, we took advantage of conventional and steered molecular dynamics (MD) simulations along with in silico and in vitro mutagenesis to investigate the ligand-induced changes of the enzyme and interactions involved in preservation of Cel5B conformations in the presence and absence of substrate. We determined that the variation in separation of catalytic glutamates in the absence and presence of substrate is due to the different protonation states of the proton donor glutamate that is largely governed by conformational changes in the β3α3 loop. In the absence of substrate, the conformation of Cel5B is preserved by an electrostatic interaction between deprotonated Glu140 and protonated His91. The ion pair is interrupted upon the binding of substrate, and the positional displacement of the β3α3 loop allows Glu140 to become oriented within the active site in a less hydrophilic microenvironment that assists in Glu140 protonation.

ACTIVE SITE PREDICTION FOR COMPARATIVE MODEL STRUCTURES WITH THEMATICS

THEMATICS (Theoretical Microscopic Titration Curves) is a simple, reliable computational predictor of the active sites of enzymes from structure. Our method, based on well-established Finite Difference Poisson–Boltzmann techniques, identifies the ionisable residues with anomalous predicted titration behavior. A cluster of two or more such perturbed residues is a very reliable predictor of the active site. The protein does not have to bear any resemblance in sequence or structure to any previously characterized protein, but the method does require the three-dimensional structure. We now present evidence that THEMATICS can also locate the active site in structures built by comparative modeling from similar structures. Results are given for a total of 21 sets of proteins, including 21 templates and 83 comparative model structures. Detailed results are presented for three sets of orthologous proteins (Triosephosphate isomerase, 6-Hydroxymethyl-7,8-dihydropterin pyrophosphokinase, and Aspartate aminotransferase) and for one set of human homologues of Aldose reductase with different functions. THEMATICS correctly locates the active site in the model structures. This suggests that the method can be applicable to a much larger set of proteins for which an experimentally determined structure is unavailable. With a few exceptions, the predicted active sites in the comparative model structures are similar to that of the corresponding template structure.

Predicting pKa values with continuous constant pH molecular dynamics

Knowledge of pK(a) values is important for understanding structure and function relationships in proteins. Over the past two decades, theoretical methods for pK(a) calculations have been mainly based on macroscopic models, in which the protein is considered as a low-dielectric cavity embedded in a high-dielectric continuum. In recent years, constant pH molecular dynamics methods have been developed based on a microscopic description of the protein. We describe here the methodology of continuous constant pH molecular dynamics (CPHMD), which has emerged as one of the most robust and accurate tools for predicting protein pK(a)s and for the study of pH-modulated conformational dynamics. We illustrate the utility of CPHMD by the calculation of pK(a)s for surface residues in ribonuclease A, buried residues in staphylococcal nuclease, and titratable groups in the intrinsically flexible protein α-lactalbumin. We will compare the CPHMD results with experimental data as well as calculations from PB-based and empirical methods. These examples demonstrate the accuracy and robustness of the CPHMD method and its ability to capture the correlation between ionization equilibria and conformational dynamics as well as the local dielectric response to structural rearrangement. Finally, we discuss future improvement of the CPHMD method.

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.

Calculation of pK(a) in proteins with the microenvironment modulated-screened coulomb potential

The MM-SCP has been applied to predict pK(a) values of titratable residues in wild type and mutants of staphylococcal nuclease (SNase). The calculations were based on crystal structures made available by the Garcia-Moreno Laboratory. In the mutants, mostly deeply buried hydrophobic residues were replaced with ionizable residues, and thus their pK(a) values could be measured and calculated using various methods. The data set used here consisted of a set of WT SNase for which His pK(a) at several ionic strengths had been measured, a set of mutants for which measured pK(a) were available and a set of 11 mutants for which the measured pK(a) were not known at the time of calculation. For this latter set, blind predictions were submitted to the protein pK(a) cooperative, 2009 workshop at Telluride, where the results of the blind predictions were discussed (the RMSD of the submitted set was 1.10 pH units). The calculations on the structures with known pK(a) indicated that in addition to weaknesses of the method, structural issues were observed that led to larger errors (>1) in pK(a) predictions. For example, different crystallographic conditions or steric clashes can lead to differences in the local environment around the titratable residue, which can produce large differences in the calculated pK(a) . To gain further insight into the reliability of the MM-SCP, pK(a) of an extended set of 54 proteins belonging to several structural classes were carried out. Here some initial results from this study are reported to help place the SNase results in the appropriate context.

CSAR scoring challenge reveals the need for new concepts in estimating protein-ligand binding affinity

The dG prediction accuracy by the Lead Finder docking software on the CSAR test set was characterized by R(2)=0.62 and rmsd=1.93 kcal/mol, and the method of preparation of the full-atom structures of the test set did not significantly affect the resulting accuracy of predictions. The primary factors determining the correlation between the predicted and experimental values were the van der Waals interactions and solvation effects. Those two factors alone accounted for R(2)=0.50. The other factors that affected the accuracy of predictions, listed in the order of decreasing importance, were the change of ligand's internal energy upon binding with protein, the electrostatic interactions, and the hydrogen bonds. It appears that those latter factors contributed to the independence of the prediction results from the method of full-atom structure preparation. Then, we turned our attention to the other factors that could potentially improve the scoring function in order to raise the accuracy of the dG prediction. It turned out that the ligand-centric factors, including Mw, cLogP, PSA, etc. or protein-centric factors, such as the functional class of protein, did not improve the prediction accuracy. Following that, we explored if the weak molecular interactions such as X-H...Ar, X-H...Hal, CO...Hal, C-H...X, stacking and π-cationic interactions (where X is N or O), that are generally of interest to the medicinal chemists despite their lack of proper molecular mechanical parametrization, could improve dG prediction. Our analysis revealed that out of these new interactions only CO...Hal is statistically significant for dG predictions using Lead FInder scoring function. Accounting for the CO...Hal interaction resulted in the reduction of the rmsd from 2.19 to 0.69 kcal/mol for the corresponding structures. The other weak interaction factors were not statistically significant and therefore irrelevant to the accuracy of dG prediction. On the basis of our findings from our participation in the CSAR scoring challenge we conclude that a significant increase of accuracy predictions necessitates breakthrough scoring approaches. We anticipate that the explicit accounting for water molecules, protein flexibility, and a more thermodynamically accurate method of dG calculation rather than single point energy calculation may lead to such breakthroughs.

Estimating binding affinities by docking/scoring methods using variable protonation states

To investigate the effects of multiple protonation states on protein-ligand recognition, we generated alternative protonation states for selected titratable groups of ligands and receptors. The selection of states was based on the predicted pK(a) of the unbound receptor and ligand and the proximity of titratable groups of the receptor to the binding site. Various ligand tautomer states were also considered. An independent docking calculation was run for each state. Several protocols were examined: using an ensemble of all generated states of ligand and receptor, using only the most probable state of the unbound ligand/receptor, and using only the state giving the most favorable docking score. The accuracies of these approaches were compared, using a set of 176 protein-ligand complexes (15 receptors) for which crystal structures and measured binding affinities are available. The best agreement with experiment was obtained when ligand poses from experimental crystal structures were used. For 9 of 15 receptors, using an ensemble of all generated protonation states of the ligand and receptor gave the best correlation between calculated and measured affinities.

Active site residue involvement in monoamine or diamine oxidation catalysed by pea seedling amine oxidase

The structures of copper amine oxidases from various sources show good similarity, suggesting similar catalytic mechanisms for all members of this enzyme family. However, the optimal substrates for each member differ, depending on the source of the enzyme and its location. The structural factors underlying substrate selectivity still remain to be discovered. With this in view, we examined the kinetic behaviour of pea seedling amine oxidase with cadaverine and hexylamine, the first bearing two, and the second only one, positively charged amino group. The dependence of K(m) and catalytic constant (k(c)) values on pH, ionic strength and temperature indicates that binding of the monoamine is driven by hydrophobic interactions. Instead, binding of the diamine is strongly facilitated by electrostatic factors, controlled by polar side-chains and two titratable residues present in the active site. The position of the docked substrate is also essential for the participation of titratable amino acid residues in the following catalytic steps. A new mechanistic model explaining the substrate-dependent kinetics of the reaction is discussed.

PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions

In this study, we have revised the rules and parameters for one of the most commonly used empirical pKa predictors, PROPKA, based on better physical description of the desolvation and dielectric response for the protein. We have introduced a new and consistent approach to interpolate the description between the previously distinct classifications into internal and surface residues, which otherwise is found to give rise to an erratic and discontinuous behavior. Since the goal of this study is to lay out the framework and validate the concept, it focuses on Asp and Glu residues where the protein pKa values and structures are assumed to be more reliable. The new and improved implementation is evaluated and discussed; it is found to agree better with experiment than the previous implementation (in parentheses): rmsd = 0.79 (0.91) for Asp and Glu, 0.75 (0.97) for Tyr, 0.65 (0.72) for Lys, and 1.00 (1.37) for His residues. The most significant advance, however, is in reducing the number of outliers and removing unreasonable sensitivity to small structural changes that arise from classifying residues as either internal or surface.

Remeasuring HEWL pK(a) values by NMR spectroscopy: methods, analysis, accuracy, and implications for theoretical pK(a) calculations

Site-specific pK(a) values measured by NMR spectroscopy provide essential information on protein electrostatics, the pH-dependence of protein structure, dynamics and function, and constitute an important benchmark for protein pK(a) calculation algorithms. Titration curves can be measured by tracking the NMR chemical shifts of several reporter nuclei versus sample pH. However, careful analysis of these curves is needed to extract residue-specific pK(a) values since pH-dependent chemical shift changes can arise from many sources, including through-bond inductive effects, through-space electric field effects, and conformational changes. We have re-measured titration curves for all carboxylates and His 15 in Hen Egg White Lysozyme (HEWL) by recording the pH-dependent chemical shifts of all backbone amide nitrogens and protons, Asp/Glu side chain protons and carboxyl carbons, and imidazole protonated carbons and protons in this protein. We extracted pK(a) values from the resulting titration curves using standard fitting methods, and compared these values to each other, and with those measured previously by ¹H NMR (Bartik et al., Biophys J 1994;66:1180–1184). This analysis gives insights into the true accuracy associated with experimentally measured pK(a) values. We find that apparent pK(a) values frequently differ by 0.5–1.0 units depending upon the nuclei monitored, and that larger differences occasionally can be observed. The variation in measured pK(a) values, which reflects the difficulty in fitting and assigning pH-dependent chemical shifts to specific ionization equilibria, has significant implications for the experimental procedures used for measuring protein pK(a) values, for the benchmarking of protein pK(a) calculation algorithms, and for the understanding of protein electrostatics in general.

Probing the structural determinants for the function of intracellular loop 2 in structurally cognate G-protein-coupled receptors

Intracellular loop 2 (IL2) in G-protein-coupled receptors (GPCRs) is functionally important, e.g., in binding to G-protein and β-arrestin. Differences in secondary structure of IL2 in the crystal structures of the very similar β(1)- and β(2)-adrenergic receptors (β(1)AR and β(2)AR, respectively), i.e., an α-helix and an L-shaped strand, respectively, emphasize the need to understand the structural basis for IL2 functionality. We studied the properties of IL2 in the context of experimental data using a Monte Carlo-based ab initio method. The procedure was validated first by verifying that the IL2 structures in β(1)AR and β(2)AR crystals were correctly reproduced, even after conformational ensemble searches at >1200 K where most secondary structure had been lost. We found that IL2 in β(1)AR and β(2)AR sampled each other's conformation but adopted different energetically preferred conformations, consistent with the crystal structures. The results indicate a persistent contextual preference for the structure of IL2, which was conserved when the IL2 sequences were interchanged between the receptors. We conclude that the protein environment, more than the IL2 sequence, regulates the IL2 structures. We extended the approach to the molecular model of 5-HT(2A)R for which no crystal structure is available and found that IL2 is predominantly helical, similar to IL2 in β(1)AR. Because the P3.57A mutation in IL2 had been shown to decrease β-arrestin binding and internalization, we predicted the effects of the mutation and found that it decreased the propensity of IL2 to form helix, identifying the helical IL2 as a component of the GPCR active form.

Cholesterol modulates the membrane effects and spatial organization of membrane-penetrating ligands for G-protein coupled receptors

The ligands of certain G-protein coupled receptors (GPCRs) are membrane soluble and reach their target from the lipid bilayer. Lipid composition and dynamics will therefore modulate the activity of these receptors, but specific roles of lipid components, including the ubiquitous cholesterol (Chol), are not clear. We have probed the organization and dynamics of such a lipid-bilayer-penetrating ligand, the endogenous ligand for the κ-opioid receptor (KOR) dynorphin A (1-17) (DynA), using molecular dynamics (MD) simulations of DynA in cholesterol-depleted and cholesterol-enriched model membranes. DynA is found to penetrate deep inside fluid dimyristoylphosphatidylcholine (DMPC) bilayers, and resides with its N-terminal helix at ∼6 Å away from the bilayer midplane, in a tilted orientation, at an ∼50° angle with respect to the membrane normal. In contrast, DynA inside DMPC/Chol membranes with 20% cholesterol (DMPC/Chol) is situated with its helical segment ∼5 A higher, i.e., closer to the lipid/water interface and in a relatively vertical orientation. The DMPC membrane shows greater thinning around the insertion and permits a stronger influx of water inside the hydrocarbon core than the DMPC/Chol membranes. Relating these results to data about key GPCR residues that have been implicated in interactions with membrane-inserting GPCR ligands, we conclude that the position of DynA in DMPC/Chol, but not in pure DMPC, correlates with generally proposed GPCR ligand entry pathways. Our predictions provide a possible mechanistic explanation as to why DynA binding to KOR, and the subsequent activation of the receptor, is facilitated in cholesterol-enriched environments. A quantitative description of DynA-induced membrane deformations is obtained with a continuum theory of membrane deformations (CTMD) that is based on hydrophobic matching. Comparison with the MD data reveals the significance of the lipid tail packing energy contribution in the DMPC/Chol mixtures in predicting equilibrium membrane shape around DynA. On this basis, specific corrections are introduced to this energy term within the CTMD framework, thereby extending the applicability of the CTMD framework to lipid raft mixtures and their interactions with GPCR proteins and their ligands.

Nucleophile activation in PD...(D/E)xK metallonucleases: an experimental and computational pK(a) study

Metallonucleases conduct metal-dependent nucleic acid hydrolysis. While metal ions serve in multiple mechanistic capacities in these enzymes, precisely how the attacking water is activated remains unclear for those lacking an obvious general base. All arguments hinge on appropriate pK(a)s for active site moieties very close to this species, and measurement of the pK(a) of a specific water molecule is difficult to access experimentally. Here we describe a computational approach for exploring the local electrostatic influences on the water-derived nucleophile in metallonucleases featuring the common PD...(D/E)xK motif. We utilized UHBD to predict the pK(a)s of active site groups, including that of a water molecule positioned to act as a nucleophile. The pK(a) of a Mg(II)-ligated water molecule hydrogen bonded to the conserved Lys70 in a Mg(II)-PvuII enzyme complex was calculated to be 6.5. The metal and the charge on the Lys group were removed in separate experiments; both resulted in the elevation of the pK(a) of this water molecule, consistent with contributions from both moieties to lowering this pK(a). This behavior is preserved among other PD...(D/E)xK metallonucleases. pK(a)s extracted from the pH dependence of the single turnover rate constant are compared to previous experimental data and the above predicted pK(a)s.

Sidechain conformational dependence of hydrogen exchange in model peptides

Peptide hydrogens that are exposed to solvent in protein X-ray structures exhibit a billion-fold range in hydroxide-catalyzed exchange rates, and these rates have previously been shown to be predictable by continuum dielectric methods to within a factor of 7, based on single protein conformations. When using a protein coil library to model the Boltzmann-weighted conformational distribution for the various N-acetyl-[X-Ala]-N-methylamides and N-acetyl-[Ala-Y]-N-methylamides, the acidity of the central amide in the individual conformers of each peptide spans nearly a million-fold range. Nevertheless, population averaging of these conformer acidities predicts the standard sidechain-dependent hydrogen exchange correction factors for nonpolar model peptides to within a factor of 30% (10(0.11)) with a correlation coefficient r=0.91. Comparison with the analogous continuum dielectric calculations for the other N-acetyl-[X-Y]-N-methylamides indicates that deviations from the isolated residue hypothesis of classical polymer theory predict appreciable errors in the exchange rates for conformationally disordered peptides when the standard sidechain-dependent hydrogen exchange rate correction factors are assumed to be independently additive. Although electronic polarizability generally dominates the dielectric shielding for the approximately 10ps lifetime of peptide ionization, evidence is presented for modest contributions from rapid intrarotamer conformational reorganization of Asn and Gln sidechains.

Titration_DB: storage and analysis of NMR-monitored protein pH titration curves

NMR-monitored pH titration experiments are routinely used to measure site-specific protein pKa values. Accurate experimental pKa values are essential in dissecting enzyme catalysis, in studying the pH-dependence of protein stability and ligand binding, in benchmarking pKa prediction algorithms, and ultimately in understanding electrostatic effects in proteins. However, due to the complex ways in which pH-dependent electrostatic and structural changes manifest themselves in NMR spectra, reported apparent pKa values are often dependent on the way that NMR pH-titration curves are analyzed. It is therefore important to retain the raw NMR spectroscopic data to allow for documentation and possible re-interpretation. We have constructed a database of primary NMR pH-titration data, which is accessible via a web interface. Here, we report statistics of the database contents and analyze the data with a global perspective to provide guidelines on best practice for fitting NMR titration curves. Titration_DB is available at http://enzyme.ucd.ie/Titration_DB. Proteins 2010. (c) 2009 Wiley-Liss, Inc.

Predicting protein pK(a) by environment similarity

A statistical method to predict protein pK(a) has been developed by using the 3D structure of a protein and a database of 434 experimental protein pK(a) values. Each pK(a) in the database is associated with a fingerprint that describes the chemical environment around an ionizable residue. A computational tool, MoKaBio, has been developed to identify automatically ionizable residues in a protein, generate fingerprints that describe the chemical environment around such residues, and predict pK(a) from the experimental pK(a) values in the database by using a similarity metric. The method, which retrieved the pK(a) of 429 of the 434 ionizable sites in the database correctly, was crossvalidated by leave-one-out and yielded root mean square error (RMSE) = 0.95, a result that is superior to that obtained by using the Null Model (RMSE 1.07) and other well-established protein pK(a) prediction tools. This novel approach is suitable to rationalize protein pK(a) by comparing the region around the ionizable site with similar regions whose ionizable site pK(a) is known. The pK(a) of residues that have a unique environment not represented in the training set cannot be predicted accurately, however, the method offers the advantage of being trainable to increase its predictive power.

A fast and accurate method for predicting pKa of residues in proteins

Predicting the pH-activities of residues in proteins is an important problem in enzyme engineering and protein design. A novel predictor called 'Pred-pK(a)' was developed based on the physicochemical properties of amino acids and protein 3D structure. The Pred-pK(a) approach considers the influence of all other residues of the protein to predict the pK(a) value of an ionizable residue. An empirical equation was formulated, in which the pK(a) value was a distance-dependent function of physicochemical parameters of 20 amino acid types, describing their electrostatic and van der Waals interaction, as well as the effects of hydrogen bonds and solvation. Two sets of coefficients, {a(alpha)} and {b(l)}, were used in the predictor: {a(alpha)} is the weight factors of 20 amino acid types and {b(l)} is the weight factors of physicochemical properties of amino acids. An iterative double least square procedure was proposed to solve the two sets of weight factors alternately and iteratively in a training set. The two coefficient sets {a(alpha)} and {b(l)} thus obtained were used to predict the pK(a) values of residues in a protein. The average predictive error is +/-0.6 pH in less than a minute in common personal computer.