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

Resolving coupled pH titrations using alchemical free energy calculations

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured p K a . Here, we derive a formalism based on double free energy differences ( Δ Δ G ) for quantifying the individual site p K a values of coupled residues. As Δ Δ G values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue p K a s in practice. We demonstrate that our approach and a previously proposed microscopic p K a formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein p K a values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed p K a values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the p K a values of both uncoupled and coupled residues.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.

PKAD-2: New entries and expansion of functionalities of the database of experimentally measured pKa's of proteins

Almost all biological reactions are pH dependent and understanding the origin of pH dependence requires knowledge of the pKa's of ionizable groups. Here we report a new edition of PKAD, the PKAD-2, which is a database of experimentally measured pKa's of proteins, both wild type and mutant proteins. The new additions include 117 wild type and 54 mutant pKa values, resulting in total 1742 experimentally measured pKa's. The new edition of PKAD-2 includes 8 new wild type and 12 new mutant proteins, resulting in total of 220 proteins. This new edition incorporates a visual 3D image of the highlighted residue of interest within the corresponding protein or protein complex. Hydrogen bonds were identified, counted, and implemented as a search feature. Other new search features include the number of neighboring residues <4A from the heaviest atom of the side chain of a given amino acid. Here, we present PKAD-2 with the intention to continuously incorporate novel features and current data with the goal to be used as benchmark for computational methods.

Basis for Accurate Protein pKa Prediction with Machine Learning

pH regulates protein structures and the associated functions in many biological processes via protonation and deprotonation of ionizable side chains where the titration equilibria are determined by pK_a's. To accelerate pH-dependent molecular mechanism research in the life sciences or industrial protein and drug designs, fast and accurate pK_a prediction is crucial. Here we present a theoretical pK_a data set PHMD549, which was successfully applied to four distinct machine learning methods, including DeepKa, which was proposed in our previous work. To reach a valid comparison, EXP67S was selected as the test set. Encouragingly, DeepKa was improved significantly and outperforms other state-of-the-art methods, except for the constant-pH molecular dynamics, which was utilized to create PHMD549. More importantly, DeepKa reproduced experimental pK_a orders of acidic dyads in five enzyme catalytic sites. Apart from structural proteins, DeepKa was found applicable to intrinsically disordered peptides. Further, in combination with solvent exposures, it is revealed that DeepKa offers the most accurate prediction under the challenging circumstance that hydrogen bonding or salt bridge interaction is partly compensated by desolvation for a buried side chain. Finally, our benchmark data qualify PHMD549 and EXP67S as the basis for future developments of protein pK_a prediction tools driven by artificial intelligence. In addition, DeepKa built on PHMD549 has been proven an efficient protein pK_a predictor and thus can be applied immediately to, for example, pK_a database construction, protein design, drug discovery, and so on.

GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pK_a's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pK_a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pK_a's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pK_a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pK_a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

Extending the Stochastic Titration CpHMD to CHARMM36m

The impact of pH on proteins is significant but often neglected in molecular dynamics simulations. Constant-pH Molecular Dynamics (CpHMD) is the state-of-the-art methodology to deal with these effects. However, it still lacks widespread adoption by the scientific community. The stochastic titration CpHMD is one of such methods that, until now, only supported the GROMOS force field family. Here, we extend this method's implementation to include the CHARMM36m force field available in the GROMACS software package. We test this new implementation with a diverse group of proteins, namely, lysozyme, Staphylococcal nuclease, and human and E. coli thioredoxins. All proteins were conformationally stable in the simulations, even at extreme pH values. The RMSE values (pK_a prediction vs experimental) obtained were very encouraging, in particular for lysozyme and human thioredoxin. We have also identified a few residues that challenged the CpHMD simulations, highlighting scenarios where the method still needs improvement independently of the force field. The CHARMM36m all-atom implementation was more computationally efficient when compared with the GROMOS 54A7, taking advantage of a shorter nonbonded interaction cutoff and a less frequent neighboring list update. The new extension will allow the study of pH effects in many systems for which this force field is particularly suited, i.e., proteins, membrane proteins, lipid bilayers, and nucleic acids.

A Fast and Interpretable Deep Learning Approach for Accurate Electrostatics-Driven pKa Predictions in Proteins

Existing computational methods for estimating pK_a values in proteins rely on theoretical approximations and lengthy computations. In this work, we use a data set of 6 million theoretically determined pK_a shifts to train deep learning models, which are shown to rival the physics-based predictors. These neural networks managed to infer the electrostatic contributions of different chemical groups and learned the importance of solvent exposure and close interactions, including hydrogen bonds. Although trained only using theoretical data, our pKAI+ model displayed the best accuracy in a test set of ∼750 experimental values. Inference times allow speedups of more than 1000× compared to physics-based methods. By combining speed, accuracy, and a reasonable understanding of the underlying physics, our models provide a game-changing solution for fast estimations of macroscopic pK_a values from ensembles of microscopic values as well as for many downstream applications such as molecular docking and constant-pH molecular dynamics simulations.

Prediction of protein pK a with representation learning

The behavior of proteins is closely related to the protonation states of the residues. Therefore, prediction and measurement of pK a are essential to understand the basic functions of proteins. In this work, we develop a new empirical scheme for protein pK a prediction that is based on deep representation learning. It combines machine learning with atomic environment vector (AEV) and learned quantum mechanical representation from ANI-2x neural network potential (J. Chem. Theory Comput. 2020, 16, 4192). The scheme requires only the coordinate information of a protein as the input and separately estimates the pK a for all five titratable amino acid types. The accuracy of the approach was analyzed with both cross-validation and an external test set of proteins. Obtained results were compared with the widely used empirical approach PROPKA. The new empirical model provides accuracy with MAEs below 0.5 for all amino acid types. It surpasses the accuracy of PROPKA and performs significantly better than the null model. Our model is also sensitive to the local conformational changes and molecular interactions.

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.

Mapping the electrostatic potential of the nucleosome acidic patch

The nucleosome surface contains an area with negative electrostatic potential known as the acidic patch, which functions as a binding platform for various proteins to regulate chromatin biology. The dense clustering of acidic residues may impact their effective pKa and thus the electronegativity of the acidic patch, which in turn could influence nucleosome-protein interactions. We here set out to determine the pKa values of residues in and around the acidic patch in the free H2A-H2B dimer using NMR spectroscopy. We present a refined solution structure of the H2A-H2B dimer based on intermolecular distance restraints, displaying a well-defined histone-fold core. We show that the conserved histidines H2B H46 and H106 that line the acidic patch have pKa of 5.9 and 6.5, respectively, and that most acidic patch carboxyl groups have pKa values well below 5.0. For H2A D89 we find strong evidence for an elevated pKa of 5.3. Our data establish that the acidic patch is highly negatively charged at physiological pH, while protonation of H2B H106 and H2B H46 at slightly acidic pH will reduce electronegativity. These results will be valuable to understand the impact of pH changes on nucleosome-protein interactions in vitro, in silico or in vivo.

Simulation of pH-Dependent Conformational Transitions in Membrane Proteins: The CLC-ec1 Cl-/H+ Antiporter

Intracellular transport of chloride by members of the CLC transporter family involves a coupled exchange between a Cl⁻ anion and a proton (H⁺), which makes the transport function dependent on ambient pH. Transport activity peaks at pH 4.5 and stalls at neutral pH. However, a structure of the WT protein at acidic pH is not available, making it difficult to assess the global conformational rearrangements that support a pH-dependent gating mechanism. To enable modeling of the CLC-ec1 dimer at acidic pH, we have applied molecular dynamics simulations (MD) featuring a new force field modification scheme—termed an Equilibrium constant pH approach (ECpH). The ECpH method utilizes linear interpolation between the force field parameters of protonated and deprotonated states of titratable residues to achieve a representation of pH-dependence in a narrow range of physiological pH values. Simulations of the CLC-ec1 dimer at neutral and acidic pH comparing ECpH-MD to canonical MD, in which the pH-dependent protonation is represented by a binary scheme, substantiates the better agreement of the conformational changes and the final model with experimental data from NMR, cross-link and AFM studies, and reveals structural elements that support the gate-opening at pH 4.5, including the key glutamates Glu_in and Glu_ex.

Mathematical Description of the Enzymatic Activity of Proteins with Ionizable Groups Exhibiting Deviations from the Henderson-Hasselbalch Equation

The ionization equilibrium implied in the calculation of the specific activity is classically described through the Henderson-Hasselbalch equation. An extension for the description of anomalous ionization profiles using the Hill equation is presented in this communication. The proposed framework was applied to the description of the specific enzymatic activity curve as a function of pH of five enzymes presenting different ionization states in their active site. The developed equation improves the description of relative enzymatic curves that deviate from the bell curve predicted by the application of the Henderson-Hasselbalch equation, regardless of the ionization scheme related to the active site.

Titr-DMD-A Rapid, Coarse-Grained Quasi-All-Atom Constant pH Molecular Dynamics Framework

The pH-dependence of enzyme fold stability and catalytic activity is a fundamentally dynamic, structural property which is difficult to study. The challenges and expense of investigating dynamic, atomic scale behavior experimentally means that computational methods, particularly constant pH molecular dynamics (CpHMD), are well situated tools for this. However, these methods often struggle with affordable sampling of sufficiently long time scales while also obtaining accurate pKa prediction and verifying the structures they generate. We introduce Titr-DMD, an affordable CpHMD method that combines the quasi-all-atom coarse-grained discrete molecular dynamics (DMD) method for conformational sampling with Propka for pKa prediction, to circumvent these issues. The combination enables rapid sampling on limited computational resources, while simulations are still performed on the atomic scale. We benchmark the method on a set of proteins with experimentally attested pKa and on the pH triggered conformational change in a staphylococcal nuclease mutant, a rare experimental study of such behavior. Our results show Titr-DMD to be an effective and inexpensive method to study pH-coupled protein dynamics.

Efficient Quantum-Chemical Calculations of Acid Dissociation Constants from Free-Energy Relationships

The calculation of acid dissociation constants (pK_a) is an important task in computational chemistry and chemoinformatics. Theoretically and with minimal empiricism, this is possible from computed acid dissociation free energies via so-called linear free-energy relationships. In this study some modifications are introduced to the latter, providing a straightforward, broadly applicable protocol with an adjustable degree of sophistication for quantum chemistry-based calculations of pK_a in water. It targets a wide pK_a range (∼70 units) and medium-sized, flexible molecules. Herein, a focus is set on the recently published r²SCAN-3c and related efficient composite density functionals and the semiempirical GFN2-xTB method, including a newly introduced energy correction for heterolytic dissociation, both in combination with implicit solvation models. The performance is evaluated in comparison with experimental data, showing mean errors often smaller than a targeted 1 pK_a unit accuracy. Larger deviations are observed only upon inclusion of challenging highly negative (<-5) or positive (>15) pK_a values. Among all those tested, it is found that B97-3c is the best performing functional, although rather independently of the density functional theory (DFT) method used; low root-mean-square errors of 0.8-1.0 pK_a units for typical drugs are obtained. For optimal performance, it is recommended to employ DFT functional specific free-energy relationship parameters. Additionally, a significant conformational dependence of the pK_a values is revealed and quantified for some nonrigid drug molecules.

Hybrid MC/MD for protein design

Computational protein design relies on simulations of a protein structure, where selected amino acids can mutate randomly, and mutations are selected to enhance a target property, such as stability. Often, the protein backbone is held fixed and its degrees of freedom are modeled implicitly to reduce the complexity of the conformational space. We present a hybrid method where short molecular dynamics (MD) segments are used to explore conformations and alternate with Monte Carlo (MC) moves that apply mutations to side chains. The backbone is fully flexible during MD. As a test, we computed side chain acid/base constants or pK_a's in five proteins. This problem can be considered a special case of protein design, with protonation/deprotonation playing the role of mutations. The solvent was modeled as a dielectric continuum. Due to cost, in each protein we allowed just one side chain position to change its protonation state and the other position to change its type or mutate. The pK_a's were computed with a standard method that scans a range of pH values and with a new method that uses adaptive landscape flattening (ALF) to sample all protonation states in a single simulation. The hybrid method gave notably better accuracy than standard, fixed-backbone MC. ALF decreased the computational cost a factor of 13.

Continuous Constant pH Molecular Dynamics Simulations of Transmembrane Proteins

Many membrane channels, transporters, and receptors utilize a pH gradient or proton coupling to drive functionally relevant conformational transitions. Conventional molecular dynamics simulations employ fixed protonation states, thus neglecting the coupling between protonation and conformational equilibria. Here we describe the membrane-enabled hybrid-solvent continuous constant pH molecular dynamics method for capturing atomic details of proton-coupled conformational dynamics of transmembrane proteins. Example protocols from our recent application studies of proton channels and ion/substrate transporters are discussed.

Alternative proton-binding site and long-distance coupling in Escherichia coli sodium-proton antiporter NhaA

Escherichia coli NhaA is a prototypical sodium-proton antiporter responsible for maintaining cellular ion and volume homeostasis by exchanging two protons for one sodium ion; despite two decades of research, the transport mechanism of NhaA remains poorly understood. Recent crystal structure and computational studies suggested Lys300 as a second proton-binding site; however, functional measurements of several K300 mutants demonstrated electrogenic transport, thereby casting doubt on the role of Lys300. To address the controversy, we carried out state-of-the-art continuous constant pH molecular dynamics simulations of NhaA mutants K300A, K300R, K300Q/D163N, and K300Q/D163N/D133A. Simulations suggested that K300 mutants maintain the electrogenic transport by utilizing an alternative proton-binding residue Asp133. Surprisingly, while Asp133 is solely responsible for binding the second proton in K300R, Asp133 and Asp163 jointly bind the second proton in K300A, and Asp133 and Asp164 jointly bind two protons in K300Q/D163N. Intriguingly, the coupling between Asp133 and Asp163 or Asp164 is enabled through the proton-coupled hydrogen-bonding network at the flexible intersection of two disrupted helices. These data resolve the controversy and highlight the intricacy of the compensatory transport mechanism of NhaA mutants. Alternative proton-binding site and proton sharing between distant aspartates may represent important general mechanisms of proton-coupled transport in secondary active transporters.

Azo dying of α-keratin material improves microbial keratinase screening and standardization

Microbial conversion through enzymatic reactions has received a lot of attention as a cost-effective and environmentally friendly way to recover amino acids and short peptides from keratin materials. However, accurate assessment of microbial keratinase activity is not straightforward, and current available methods lack sensitivity and standardization. Here, we suggest an optimized Azokeratin assay, with substrate generated directly from azo-dyed raw keratin material. We introduced supernatant filtration in the protocol for optimal stopping of keratinase reactions instead of the widely used trichloroacetic acid (TCA), as it generated biases and impacted the sensitivity. We furthermore suggest a method for standardization of keratinase activity signals using proteinase K, a well-known keratinase, as a reference enabling reproducibility between studies. Lastly, we evaluated our developed method with several bacterial isolates through benchmarking against a commercial assay (Keratin Azure). Under different setups, the Azokeratin method was more sensitive than commonly used Keratin Azure-based assays (3-fold). We argue that this method could be applied with any type of keratin substrate, enabling more robust and sensitive results which can be used for further comparison with other studies, thus representing an important progress within the field of microbial keratin degradation.

Methods to Study the Structure and Catalytic Activity of cis-Splicing Inteins

The autocatalytic process of protein splicing is facilitated by an intein, which interrupts flanking polypeptides called exteins. The mechanism of protein splicing has been studied by overexpression in E. coli of intein fusion proteins with nonnative exteins. Inteins can be used to generate reactive α-thioesters, as well as proteins with N-terminal Cys residues, to facilitate expressed protein ligation. As such, a more detailed understanding of the function of inteins can have significant impact for biotechnology applications. Here, we provide biochemical methods to study splicing activity and NMR methods to study intein structure and the catalytic mechanism.

GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations

We present a GPU implementation of the continuous constant pH molecular dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. To test the accuracy of the tool for rapid pK_a predictions, a series of 2 ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calculated pK_a's showed a root-mean-square deviation of 0.80 and correlation coefficient of 0.83 with respect to experiment. Also, 90% of the pK_a's were converged with estimated errors below 0.1 pH units. Surprisingly, this level of accuracy is similar to our previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps^-1 (Huang, Harris, and Shen J. Chem. Inf. Model. 2018 , 58 , 1372 - 1383 ). Interestingly, for the linked titration sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the experimental macroscopic pK_a values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pK_a's is underestimated with frequent exchange attempts such as 2 ps^-1, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2 ns single-pH simulation of a 400-residue protein takes about 1 h on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, we envision that GPU-accelerated continuous CpHMD may be used in routine pK_a predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.

Major conformational changes in the structure of lysozyme obtained from a crystal with a very low solvent content

A new crystal form of lysozyme with a very low solvent content (26.35%) has been obtained in the orthorhombic space group P212121 (with unit-cell parameters a = 30.04, b = 51.68, c = 61.53 Å). The lysozyme structure obtained from these crystals does not show the typical overall fold. Instead, major conformational changes take place in some elements of the secondary structure and in the hydrophobic core of the protein. At the end of the central α-helix (α2), Glu35 is usually buried in the catalytic site and shows an abnormally high pKa value, which is key to the activity of the enzyme. The high pKa value of this glutamate residue is favoured by the hydrophobic environment, particularly by its neighbour Trp108, which is important for structural stability and saccharide binding. In this new structure, Trp108 shows a 90° rotation of its side chain, which results in the rearrangement of the hydrophobic core. Conformational changes also result in the exposure of Glu35 to the solvent, which impairs the catalytic site by increasing the distance between Glu35 and Asp52 and lowering the pKa value of the glutamate. Altogether, this new lysozyme structure reveals major conformational changes in the hydrophobic core and catalytic site that might play a role in the folding and bactericidal function of the protein.

A Hierarchical Approach to Predict Conformation-Dependent Histidine Protonation States in Stable and Flexible Proteins

Solution acidity measured by pH is an important environmental factor that affects protein structure. It influences the protonation state of protein residues, which in turn may be coupled to protein conformational changes, unfolding, and ligand binding. It remains difficult to compute and measure the p K_a of individual residues, as well as to relate them to pH-dependent protein transitions. This paper presents a hierarchical approach to compute the p K_a of individual protonatable residues, specifically histidines, coupled with underlying structural changes of a protein. A fast and efficient free energy perturbation (FEP) algorithm has also been developed utilizing a fast implementation of standard molecular dynamics (MD) algorithms. Specifically, a CUDA version of the AMBER MD engine is used in this paper. Eight histidine p K_a's are computed in a diverse set of pH stable proteins to demonstrate the proposed approach's utility and assess the predictive quality of the AMBER FF99SB force field. A reference molecule is carefully selected and tested for convergence. A hierarchical approach is used to model p K_a's of the six histidine residues of the diphtheria toxin translocation domain (DTT), which exhibits a diverse ensemble of individual conformations and pH-dependent unfolding. The hierarchical approach consists of first sampling equilibrium conformational ensembles of a protein with protonated and neutral histidine residues via long equilibrium MD simulations (Flores-Canales, J. C.; et al. bioRxiv, 2019, 572040). A clustering method is then used to identify sampled protein conformations, and p K_a's of histidines in each protein conformation are computed. Finally, an ensemble averaging formalism is developed to compute weighted average histidine p K_a's. These can be compared with an apparent experimentally measured p K_a of the DTT protein and thus allows us to propose a mechanism of pH-dependent unfolding of the DTT protein.

Side chain electrostatic interactions and pH-dependent expansion of the intrinsically disordered, highly acidic carboxyl-terminus of γ-tubulin

Intramolecular electrostatic attraction and repulsion strongly influence the conformational sampling of intrinsically disordered proteins and domains (IDPs). In order to better understand this complex relationship, we have used nuclear magnetic resonance to measure side chain pK_a values and pH-dependent translational diffusion coefficients for the unstructured and highly acidic carboxyl-terminus of γ-tubulin (γ-CT), providing insight into how the net charge of an IDP relates to overall expansion or collapse of the conformational ensemble. Many of the pK_a values in the γ-CT are shifted upward by 0.3-0.4 units and exhibit negatively cooperative ionization pH profiles, likely due to the large net negative charge that accumulates on the molecule as the pH is raised. pK_a shifts of this magnitude correspond to electrostatic interaction energies between the affected residues and the rest of the charged molecule that are each on the order of 1 kcal mol^-1 . Diffusion of the γ-CT slowed with increasing net charge, indicative of an expanding hydrodynamic radius (r_H ). The degree of expansion agreed quantitatively with what has been seen from comparisons of IDPs with different charge content, yielding the general trend that every 0.1 increase in relative charge (|Q|/res) produces a roughly 5% increase in r_H . While γ-CT pH titration data followed this trend nearly perfectly, there were substantially larger deviations for the database of different IDP sequences. This suggests that other aspects of an IDP's primary amino acid sequence beyond net charge influence the sensitivity of r_H to electrostatic interactions.

Experimental pKa Value Determination of All Ionizable Groups of a Hyperstable Protein

Electrostatic interactions significantly contribute to the stability and function of proteins. The stabilizing or destabilizing effect of local charge is reflected in the perturbation of the pK_a value of an ionizable group from the intrinsic pK_a value. Herein, the charge network of a hyperstable dimeric protein (ribbon–helix–helix (rhh) protein from plasmid pRN1 from Sulfolobus islandicus) is studied through experimental determination of the pK_a values of all ionizable groups. Transitions were monitored by multiple NMR signals per ionizable group between pH 0 and 12.5, prior to a global analysis, which accounted for the effects of neighboring residues. It is found that for several residues involved in salt bridges (four Asp and one Lys) the pK_a values are shifted in favor of the charged state. Furthermore, the pK_a values of residues C40 and Y47, both located in the hydrophobic dimer interface, are shifted beyond 13.7. The necessary energy for such a shift is about two‐thirds of the total stability of the protein, which confirms the importance of the hydrophobic core to the overall stability of the rhh protein.

Electrostatically Driven Protein Adsorption: Charge Patches versus Charge Regulation

The mechanisms of electrostatically driven adsorption of proteins on charged surfaces are studied with a new theoretical framework. The acid-base behavior, charge distribution, and electrostatic contributions to the thermodynamic properties of the proteins are modeled in the presence of a charged surface. The method is validated against experimental titration curves and apparent p K_as. The theory predicts that electrostatic interactions favor the adsorption of proteins at their isoelectric points on charged surfaces despite the fact that the protein has no net charge in solution. Two known mechanisms explain adsorption under these conditions: (i) charge regulation (the charge of the protein changes due to the presence of the surface) and (ii) charge patches (the protein orients to place charged amino acids near opposite surface charges). This work shows that both mechanisms contribute to adsorption at low ionic strengths, whereas only the charge-patch mechanism operates at high ionic strength. Interestingly, the contribution of charge regulation is insensitive to protein orientation under all conditions, which validates the use of constant-charge simulations to determine the most stable orientation of adsorbed proteins. The present study also shows that the charged surface can induce large shifts in the apparent p K_as of individual amino acids in adsorbed proteins. Our conclusions are valid for all proteins studied in this work (lysozyme, α-amylase, ribonuclease A, and β-lactoglobulin), as well as for proteins that are not isoelectric but have instead a net charge in solution of the same sign as the surface charge, i.e. the problem of protein adsorption on the "wrong side" of the isoelectric point.

Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis

Solution pH plays an important role in structure and dynamics of biomolecular systems; however, pH effects cannot be accurately accounted for in conventional molecular dynamics simulations based on fixed protonation states. Continuous constant pH molecular dynamics (CpHMD) based on the λ-dynamics framework calculates protonation states on the fly during dynamical simulation at a specified pH condition. Here we report the CPU-based implementation of the CpHMD method based on the GBNeck2 generalized Born (GB) implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. The performance of the method was tested using pH replica-exchange titration simulations of Asp, Glu and His side chains in 4 miniproteins and 7 enzymes with experimentally known p K_a's, some of which are significantly shifted from the model values. The added computational cost due to CpHMD titration ranges from 11 to 33% for the data set and scales roughly linearly as the ratio between the titrable sites and number of solute atoms. Comparison of the experimental and calculated p K_a's using 2 ns per replica sampling yielded a mean unsigned error of 0.70, a root-mean-squared error of 0.91, and a linear correlation coefficient of 0.79. Though this level of accuracy is similar to the GBSW-based CpHMD in CHARMM, in contrast to the latter, the current implementation was able to reproduce the experimental orders of the p K_a's of the coupled carboxylic dyads. We quantified the sampling errors, which revealed that prolonged simulation is needed to converge p K_a's of several titratable groups involved in salt-bridge-like interactions or deeply buried in the protein interior. Our benchmark data demonstrate that GBNeck2-CpHMD is an attractive tool for protein p K_a predictions.

Orthorhombic lysozyme crystallization at acidic pH values driven by phosphate binding

The structure of orthorhombic lysozyme has been obtained at 298 K and pH 4.5 using sodium chloride as the precipitant and in the presence of sodium phosphate at a concentration as low as 5 mM. Crystals belonging to space group P2₁2₁2₁ (unit-cell parameters a = 30, b = 56, c = 73 Å, α = β = γ = 90.00°) diffracted to a resolution higher than 1 Å, and the high quality of these crystals permitted the identification of a phosphate ion bound to Arg14 and His15. The binding of this ion produces long-range conformational changes affecting the loop containing Ser60-Asn74. The negatively charged phosphate ion shields the electrostatic repulsion of the positively charged arginine and histidine residues, resulting in higher stability of the phosphate-bound lysozyme. Additionally, a low-humidity orthorhombic variant was obtained at pH 4.5, and comparison with those previously obtained at pH 6.5 and 9.5 shows a 1.5 Å displacement of the fifth α-helix towards the active-site cavity, which might be relevant to protein function. Since lysozyme is broadly used as a model protein in studies related to protein crystallization and amyloid formation, these results indicate that the interaction of some anions must be considered when analysing experiments performed at acidic pH values.

Predicting Catalytic Proton Donors and Nucleophiles in Enzymes: How Adding Dynamics Helps Elucidate the Structure-Function Relationships

Despite the relevance of understanding structure-function relationships, robust prediction of proton donors and nucleophiles in enzyme active sites remains challenging. Here we tested three types of state-of-the-art computational methods to calculate the p K_a's of the buried and hydrogen bonded catalytic dyads in five enzymes. We asked the question what determines the p K_a order, i.e., what makes a residue proton donor vs a nucleophile. The continuous constant pH molecular dynamics simulations captured the experimental p K_a orders and revealed that the negative nucleophile is stabilized by increased hydrogen bonding and solvent exposure as compared to the proton donor. Surprisingly, this simple trend is not apparent from crystal structures and the static structure-based calculations. While the generality of the findings awaits further testing via a larger set of data, they underscore the role of dynamics in bridging enzyme structures and functions.

A Coupled Ionization-Conformational Equilibrium Is Required To Understand the Properties of Ionizable Residues in the Hydrophobic Interior of Staphylococcal Nuclease

Ionizable residues in the interior of proteins play essential roles, especially in biological energy transduction, but are relatively rare and seem incompatible with the complex and polar environment. We perform a comprehensive study of the internal ionizable residues on 21 variants of staphylococcal nuclease with internal Lys, Glu, or Asp residues. Using pH replica exchange molecular dynamics simulations, we find that, in most cases, the pK_a values of these internal ionizable residues are shifted significantly from their values in solution. Our calculated results are in excellent agreement with the experimental observations of the Garcia-Moreno group. We show that the interpretation of the experimental pK_a values requires the study of not only protonation changes but also conformational changes. The coupling between the protonation and conformational equilibria suggests a mechanism for efficient pH-sensing and regulation in proteins. This study provides new physical insights into how internal ionizable residues behave in the hydrophobic interior of proteins.

Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton

Computational free-energy correction strategies and the choice of experimental proton hydration free energy, ΔG_s^*(H⁺), are analyzed to investigate the apparent controversy in experimental thermodynamics of ionic hydration. Without corrections, the hydration free-energy (ΔG_hyd) calculations match experiments with ΔG_s^*(H⁺) = -1064 kJ/mol as reference. Using the Galvani surface potential the resulting (real) ΔG_hyd are consistent with ΔG_s^*(H⁺) = -1098 kJ/mol. When applying, in an ad hoc manner, the discrete solvent correction, ΔG_hyd matching the "consensus" ΔG_s^*(H⁺) of -1112 kJ/mol are obtained. This analysis rationalizes reports on ΔG_hyd calculations for ions using different experimental references. For neutral amino acid side chains ΔG_hyd are independent of the water model, whereas there are large differences in ΔG_hyd due to the water model for charged species, suggesting that long-range ordering of water around ions yields an important contribution to the ΔG_hyd. These differences are reduced significantly when applying consistent corrections, but to obtain the most accurate results it is recommended to use the water model belonging to the force field.

Conformational dynamics of cathepsin D and binding to a small-molecule BACE1 inhibitor

BACE1 is a major therapeutic target for prevention and treatment of Alzheimer's disease. Developing inhibitors that can selectively target BACE1 in favor of other proteases, especially cathepsin D (CatD), has presented significant challenges. Here, we investigate the conformational dynamics and protonation states of BACE1 and CatD using continuous constant pH molecular dynamics with pH replica-exchange sampling protocol. Despite similar structure, BACE1 and CatD exhibit markedly different active site dynamics. BACE1 displays pH-dependent flap dynamics that controls substrate accessibility, while the CatD flap is relatively rigid and remains open in the pH range 2.5-6. Interestingly, although each protease hydrolyzes peptide bonds, the protonation states of the catalytic dyads are different within the active pH range. The acidic and basic components of the BACE1 catalytic dyad are clear, while either aspartic acid of the CatD catalytic dyad could play the role of acid or base. Finally, we investigate binding of the inhibitor LY2811376 developed by Eli Lilly to BACE1 and CatD. Surprisingly, in the enzyme active pH range, LY2811376 forms a stronger salt bridge with the catalytic dyad in CatD than in BACE1, which might explain the retinal toxicity of the inhibitor related to off-target inhibition of CatD. This work highlights the complexity and challenge in structure-based drug design where receptor-ligand binding induces protonation state change in both the protein and the inhibitor. © 2017 Wiley Periodicals, Inc.

Intrinsically disordered inhibitor of glutamine synthetase is a functional protein with random-coil-like pKa values

The sequential action of glutamine synthetase (GS) and glutamate synthase (GOGAT) in cyanobacteria allows the incorporation of ammonium into carbon skeletons. In the cyanobacterium Synechocystis sp. PCC 6803, the activity of GS is modulated by the interaction with proteins, which include a 65-residue-long intrinsically disordered protein (IDP), the inactivating factor IF7. This interaction is regulated by the presence of charged residues in both IF7 and GS. To understand how charged amino acids can affect the binding of an IDP with its target and to provide clues on electrostatic interactions in disordered states of proteins, we measured the pK_a values of all IF7 acidic groups (Glu32, Glu36, Glu38, Asp40, Asp58, and Ser65, the backbone C-terminus) at 100 mM NaCl concentration, by using NMR spectroscopy. We also obtained solution structures of IF7 through molecular dynamics simulation, validated them on the basis of previous experiments, and used them to obtain theoretical estimates of the pK_a values. Titration values for the two Asp and three Glu residues of IF7 were similar to those reported for random-coil models, suggesting the lack of electrostatic interactions around these residues. Furthermore, our results suggest the presence of helical structure at the N-terminus of the protein and of conformational changes at acidic pH values. The overall experimental and in silico findings suggest that local interactions and conformational equilibria do not play a role in determining the electrostatic features of the acidic residues of IF7.

Determination of the pKa of the N-terminal amino group of ubiquitin by NMR

Ubiquitination regulates nearly every aspect of cellular life. It is catalysed by a cascade of three enzymes and results in the attachment of the C-terminal carboxylate of ubiquitin to a lysine side chain in the protein substrate. Chain extension occurs via addition of subsequent ubiquitin molecules to either one of the seven lysine residues of ubiquitin, or via its N-terminal α-amino group to build linear ubiquitin chains. The pK_a of lysine side chains is around 10.5 and hence E3 ligases require a mechanism to deprotonate the amino group at physiological pH to produce an effective nucleophile. In contrast, the pK_a of N-terminal α-amino groups of proteins can vary significantly, with reported values between 6.8 and 9.1, raising the possibility that linear chain synthesis may not require a general base. In this study we use NMR spectroscopy to determine the pK_a for the N-terminal α-amino group of methionine1 of ubiquitin for the first time. We show that it is 9.14, one of the highest pK_a values ever reported for this amino group, providing a rational for the observed need for a general base in the E3 ligase HOIP, which synthesizes linear ubiquitin chains.

The metalation of hen egg white lysozyme impacts protein stability as shown by ion mobility mass spectrometry, differential scanning calorimetry, and X-ray crystallography

Metalation of lysozyme with anticancer organometallics results in protein destabilisation, probably relevant in metallodrug mode of action.

All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water

Development of a pH stat to properly control solution pH in biomolecular simulations has been a long-standing goal in the community. Toward this goal recent years have witnessed the emergence of the so-called constant pH molecular dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or truncation-based electrostatic schemes. Here we report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous constant pH molecular dynamics (CpHMD) method, enabling CpHMD to be performed with a standard MD engine at a fractional added computational cost. We demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP36, BBL, HEWL, and SNase. With the sampling time of 10 ns per replica, most pK_a's are converged, yielding the average absolute and root-mean-square deviations of 0.61 and 0.77, respectively, from experiment. Linear regression of the calculated vs experimental pK_a shifts gives a correlation coefficient of 0.79, a slope of 1, and an intercept near 0. Analysis reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromolecular systems in various environments, enabling atomic insights into pH-dependent phenomena involving not only soluble proteins but also transmembrane proteins, nucleic acids, surfactants, and polysaccharides.

A Soluble, Folded Protein without Charged Amino Acid Residues

Charges are considered an integral part of protein structure and function, enhancing solubility and providing specificity in molecular interactions. We wished to investigate whether charged amino acids are indeed required for protein biogenesis and whether a protein completely free of titratable side chains can maintain solubility, stability, and function. As a model, we used a cellulose-binding domain from Cellulomonas fimi, which, among proteins of more than 100 amino acids, presently is the least charged in the Protein Data Bank, with a total of only four titratable residues. We find that the protein shows a surprising resilience toward extremes of pH, demonstrating stability and function (cellulose binding) in the pH range from 2 to 11. To ask whether the four charged residues present were required for these properties of this protein, we altered them to nontitratable ones. Remarkably, this chargeless protein is produced reasonably well in Escherichia coli, retains its stable three-dimensional structure, and is still capable of strong cellulose binding. To further deprive this protein of charges, we removed the N-terminal charge by acetylation and studied the protein at pH 2, where the C-terminus is effectively protonated. Under these conditions, the protein retains its function and proved to be both soluble and have a reversible folding-unfolding transition. To the best of our knowledge, this is the first time a soluble, functional protein with no titratable side chains has been produced.

Continuum Electrostatics Approaches to Calculating pKas and Ems in Proteins

Proteins change their charge state through protonation and redox reactions as well as through binding charged ligands. The free energy of these reactions is dominated by solvation and electrostatic energies and modulated by protein conformational relaxation in response to the ionization state changes. Although computational methods for calculating these interactions can provide very powerful tools for predicting protein charge states, they include several critical approximations of which users should be aware. This chapter discusses the strengths, weaknesses, and approximations of popular computational methods for predicting charge states and understanding the underlying electrostatic interactions. The goal of this chapter is to inform users about applications and potential caveats of these methods as well as outline directions for future theoretical and computational research.

Mimicking titration experiments with MD simulations: A protocol for the investigation of pH-dependent effects on proteins

Protein structure and function are highly dependent on the environmental pH. However, the temporal or spatial resolution of experimental approaches hampers direct observation of pH-induced conformational changes at the atomic level. Molecular dynamics (MD) simulation strategies (e.g. constant pH MD) have been developed to bridge this gap. However, one frequent problem is the sampling of unrealistic conformations, which may also lead to poor pKa predictions. To address this problem, we have developed and benchmarked the pH-titration MD (pHtMD) approach, which is inspired by wet-lab titration experiments. We give several examples how the pHtMD protocol can be applied for pKa calculation including peptide systems, Staphylococcus nuclease (SNase), and the chaperone HdeA. For HdeA, pHtMD is also capable of monitoring pH-dependent dimer dissociation in accordance with experiments. We conclude that pHtMD represents a versatile tool for pKa value calculation and simulation of pH-dependent effects in proteins.

Constant-pH MD Simulations Portray the Protonation and Structural Behavior of Four Decapeptides Designed to Coordinate Cu(2+)

The cyclic decapeptide C-Asp, containing one Asp residue and three His residues, was designed by Fragoso et al. (Chem. Eur. J. 2013, 19, 2076) to bind Cu(2+) exclusively through the side chain groups and mimic copper coordination in metalloproteins. A variant of the cyclodecapeptide where Asp is substituted by Asn (C-Asn) has also been synthesized in addition to the linear ("open") counterparts of both forms (O-Asp and O-Asn), testing the importance of cyclization and the presence of Asp in Cu(2+) coordination (Chem. Eur. J. 2013, 19, 2076; Dalton Trans. 2013, 42, 6182). All peptides formed a major species at neutral pH that was able to coordinate Cu(2+) exclusively through the neutral imidazole groups and the Asp side chain, when present, with C-Asp being the most effective. A detailed description of the protonation behavior of each histidine could help understanding the coordination species being formed in the pH range and eventually further optimizing the peptide's design. However, the standard current methods (NMR titrations) are not very suited for proximal groups titrating in the same pH range. In this work, we used the stochastic titration constant-pH molecular dynamics method to calculate the protonation curves and pKa of each titrable residue in the four decapeptides, in the absence of Cu(2+) ions. The global protonation curves obtained in our simulations are in very good agreement with the existing potentiometric titration curves. The histidines are titrating very closely, and the Asp forms abundant salt bridges with the basic residues, displaying an unusually low pKa value. In addition, we could observe that the four peptides are very unstructured in the absence of copper, and not even the cyclic forms exhibit a significant β-sheet, unlike what could be expected from the presence of β-turn inducer units in this type of scaffold.

Non-uniform sampling of NMR relaxation data

The use of non-uniform sampling of NMR spectra may give significant reductions in the data acquisition time. For quantitative experiments such as the measurement of spin relaxation rates, non-uniform sampling is however not widely used as inaccuracies in peak intensities may lead to errors in the extracted dynamic parameters. By systematic reducing the coverage of the Nyquist grid of (15)N Carr-Purcell-Meiboom-Gill (CPMG) relaxation dispersion datasets for four different proteins and performing a full data analysis of the resulting non-uniform sampled datasets, we have compared the performance of the multi-dimensional decomposition and iterative re-weighted least-squares algorithms in reconstructing spectra with accurate peak intensities. As long as a single fully sampled spectrum is included in a series of otherwise non-uniform sampled two-dimensional spectra, multi-dimensional decomposition reconstructs the non-uniform sampled spectra with high accuracy. For two of the four analyzed datasets, a coverage of only 20% results in essentially the same results as the fully sampled data. As exemplified by other data, such a low coverage is in general not enough to produce reliable results. We find that a coverage level not compromising the final results can be estimated by recording a single full two-dimensional spectrum and reducing the spectrum quality in silico.

Asymmetric protonation of EmrE

The small multidrug resistance transporter EmrE is a homodimer that uses energy provided by the proton motive force to drive the efflux of drug substrates. The pKa values of its "active-site" residues--glutamate 14 (Glu14) from each subunit--must be poised around physiological pH values to efficiently couple proton import to drug export in vivo. To assess the protonation of EmrE, pH titrations were conducted with (1)H-(15)N TROSY-HSQC nuclear magnetic resonance (NMR) spectra. Analysis of these spectra indicates that the Glu14 residues have asymmetric pKa values of 7.0 ± 0.1 and 8.2 ± 0.3 at 45°C and 6.8 ± 0.1 and 8.5 ± 0.2 at 25°C. These pKa values are substantially increased compared with typical pKa values for solvent-exposed glutamates but are within the range of published Glu14 pKa values inferred from the pH dependence of substrate binding and transport assays. The active-site mutant, E14D-EmrE, has pKa values below the physiological pH range, consistent with its impaired transport activity. The NMR spectra demonstrate that the protonation states of the active-site Glu14 residues determine both the global structure and the rate of conformational exchange between inward- and outward-facing EmrE. Thus, the pKa values of the asymmetric active-site Glu14 residues are key for proper coupling of proton import to multidrug efflux. However, the results raise new questions regarding the coupling mechanism because they show that EmrE exists in a mixture of protonation states near neutral pH and can interconvert between inward- and outward-facing forms in multiple different protonation states.