Enhancing Conformation and Protonation State Sampling of Hen Egg White Lysozyme Using pH Replica Exchange Molecular Dynamics

Assessing pH-Dependent Conformational Changes in the Fusion Peptide Proximal Region of the SARS-CoV-2 Spike Glycoprotein

One of the entry mechanisms of the SARS-CoV-2 coronavirus into host cells involves endosomal acidification. It has been proposed that under acidic conditions, the fusion peptide proximal region (FPPR) of the SARS-CoV-2 spike glycoprotein acts as a pH-dependent switch, modulating immune response accessibility by influencing the positioning of the receptor binding domain (RBD). This would provide indirect coupling of RBD opening to the environmental pH. Here, we explored this possible pH-dependent conformational equilibrium of the FPPR within the SARS-CoV-2 spike glycoprotein. We analyzed hundreds of experimentally determined spike structures from the Protein Data Bank and carried out pH-replica exchange molecular dynamics to explore the extent to which the FPPR conformation depends on pH and the positioning of the RBD. A meta-analysis of experimental structures identified alternate conformations of the FPPR among structures in which this flexible regions was resolved. However, the results did not support a correlation between the FPPR conformation and either RBD position or the reported pH of the cryo-EM experiment. We calculated pKa values for titratable side chains in the FPPR region using PDB structures, but these pKa values showed large differences between alternate PDB structures that otherwise adopt the same FPPR conformation type. This hampers the comparison of pKa values in different FPPR conformations to rationalize a pH-dependent conformational change. We supplemented these PDB-based analyses with all-atom simulations and used constant-pH replica exchange molecular dynamics to estimate pKa values in the context of flexibility and explicit water. The resulting titration curves show good reproducibility between simulations, but they also suggest that the titration curves of the different FPPR conformations are the same within the error bars. In summary, we were unable to find evidence supporting the previously published hypothesis of an FPPR pH-dependent equilibrium: neither from existing experimental data nor from constant-pH MD simulations. The study underscores the complexity of the spike system and opens avenues for further exploration into the interplay between pH and SARS-CoV-2 viral entry mechanisms.

phbuilder: A Tool for Efficiently Setting up Constant pH Molecular Dynamics Simulations in GROMACS

Constant pH molecular dynamics (MD) is a powerful technique that allows the protonation state of residues to change dynamically, thereby enabling the study of pH dependence in a manner that has not been possible before. Recently, a constant pH implementation was incorporated into the GROMACS MD package. Although this implementation provides good accuracy and performance, manual modification and the preparation of simulation input files are required, which can be complicated, tedious, and prone to errors. To simplify and automate the setup process, we present phbuilder, a tool that automatically prepares constant pH MD simulations for GROMACS by modifying the input structure and topology as well as generating the necessary parameter files. phbuilder can prepare constant pH simulations from both initial structures and existing simulation systems, and it also provides functionality for performing titrations and single-site parametrizations of new titratable group types. The tool is freely available at www.gitlab.com/gromacs-constantph. We anticipate that phbuilder will make constant pH simulations easier to set up, thereby making them more accessible to the GROMACS user community.

Modeling pH-Dependent Biomolecular Photochemistry

The tuning mechanism of pH can be extremely challenging to model computationally in complex biological systems, especially with respect to the photochemical properties. This article reports a protocol aimed at modeling pH-dependent photodynamics using a combination of constant-pH molecular dynamics and semiclassical nonadiabatic molecular dynamics simulations. With retinal photoisomerization in Anabaena sensory rhodopsin (ASR) as a testbed, we show that our protocol produces pH-dependent photochemical properties, such as the isomerization quantum yield or decay rates. We decompose our results into single-titrated residue contributions, identifying some key tuning amino acids. Additionally, we assess the validity of the single protonation state picture to represent the system at a given pH and propose the most populated protein charge state as a compromise between cost and accuracy.

The role of conformational change and key glutamic acid residues in the ClC-ec1 antiporter

The triple glutamine (Q) mutant (QQQ) structure of a Cl-/H+ antiporter from Escherichia coli (ClC-ec1) displaying a novel backbone arrangement has been used to challenge the long-held notion that Cl-/H+ antiporters do not operate through large conformational motions. The QQQ mutant substitutes the glutamine residue for an external glutamate E148, an internal glutamate E203, and a third glutamate E113 that hydrogen-bonds with E203. However, it is unknown if QQQ represents a physiologically relevant state, as well as how the protonation of the wild-type glutamates relates to the global dynamics. We herein apply continuous constant-pH molecular dynamics to investigate the H+-coupled dynamics of ClC-ec1. Although any large-scale conformational rearrangement upon acidification would be due to the accumulation of excess charge within the protein, protonation of the glutamates significantly impacts mainly the local structure and dynamics. Despite the fact that the extracellular pore enlarges at acidic pHs, an occluded ClC-ec1 within the active pH range of 3.5-7.5 requires a protonated E148 to facilitate extracellular Cl- release. E203 is also involved in the intracellular H+ transfer as an H+ acceptor. The water wire connection of E148 with the intracellular solution is regulated by the charge states of the E113/E203 dyad with coupled proton titration. However, the dynamics extracted from our simulations are not QQQ-like, indicating that the QQQ mutant does not represent the behavior of the wild-type ClC-ec1. These findings reinforce the necessity of having a protonatable residue at the E203 position in ClC-ec1 and suggest that a higher level of complexity exists for the intracellular H+ transfer in Cl-/H+ antiporters.

Elucidating the Protonation State of the γ-Secretase Catalytic Dyad

γ-Secretase (GS) is an intramembrane aspartyl protease that participates in the sequential cleavage of C99 to generate different isoforms of the amyloid-β (Aβ) peptides that are associated with the development of Alzheimer's disease. Due to its importance in the proteolytic processing of C99 by GS, we performed pH replica exchange molecular dynamics (pH-REMD) simulations of GS in its apo and substrate-bound forms to sample the protonation states of the catalytic dyad. We found that the catalytic dyad is deprotonated at physiological pH in our apo form, but the presence of the substrate at the active site displaces its monoprotonated state toward physiological pH. Our results show that Asp257 acts as the general base and Asp385 as the general acid during the cleavage mechanism. We identified different amino acids such as Lys265, Arg269, and the PAL motif interacting with the catalytic dyad and promoting changes in its acid-base behavior. Finally, we also found a significant pK_a shift of Glu280 related to the internalization of TM6-CT in the GS-apo form. Our study provides critical mechanistic insight into the GS mechanism and the basis for future research on the genesis of Aβ peptides and the development of Alzheimer's disease.

Increasing the Realism of in Silico pHLIP Peptide Models with a Novel pH Gradient CpHMD Method

The pH-low insertion peptides (pHLIP) are pH-dependent membrane inserting peptides, whose function depends on the cell microenvironment acidity. Several peptide variants have been designed to improve upon the wt-sequence, particularly the state transition kinetics and the selectivity for tumor pH. The variant 3 (Var3) peptide is a 27 residue long peptide, with a key titrating residue (Asp-13) that, despite showing a modest performance in liposomes (pK^ins ∼ 5.0), excelled in tumor cell experiments. To help rationalize these results, we focused on the pH gradient in the cell membrane, which is one of the crucial properties that are not present in liposomes. We extended our CpHMD-L method and its pH replica-exchange (pHRE) implementation to include a pH gradient and mimic the pHLIP-membrane microenvironment in a cell where the internal pH is fixed (pH 7.2) and the external pH is allowed to change. We showed that, by properly modeling the pH-gradient, we can correctly predict the experimentally observed loss and gain of performance in tumor cells experiments by the wt and Var3 sequences, respectively. In sum, the pH gradient implementation allowed for more accurate and realistic pK_a estimations and was a pivotal step in bridging the in silico data and the in vivo cell experiments.

pH-dependent structural diversity of profilin allergens determines thermal stability

The family of profilin allergens is a common class of proteins found in plants, viruses and various eukaryotes including mammals. Profilins are characterized by an evolutionary conserved structural fold, which is responsible for their cross-reactive nature of Immunoglobulin E (IgE) antibodies. Despite their high overall structural similarity, they exhibit substantial differences in their biophysical properties, such as thermal and pH stability. To understand the origin of these functional differences of Amb a 8, Art v 4 and Bet v 2, we performed constant pH molecular dynamics simulation in combination with Gaussian accelerated MD simulations. Depending on the respective protonation at different pH levels, we find distinct differences in conformational flexibility, which are consistent with experimentally determined melting temperatures. These variations in flexibility are accompanied by ensemble shifts in the conformational landscape and quantified and localized by residue-wise B-factors and dihedral entropies. These findings strengthen the link between flexibility of profilin allergens and their thermal stability. Thus, our results clearly show the importance of considering protonation dependent conformational ensembles in solution to elucidate biophysical differences between these structurally similar allergens.

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.

Principles of target DNA cleavage and the role of Mg2+ in the catalysis of CRISPR-Cas9

At the core of the CRISPR-Cas9 genome-editing technology, the endonuclease Cas9 introduces site-specific breaks in DNA. However, precise mechanistic information to ameliorating Cas9 function is still missing. Here, multi-microsecond molecular dynamics, free-energy and multiscale simulations are combined with solution NMR and DNA cleavage experiments to resolve the catalytic mechanism of target DNA cleavage. We show that the conformation of an active HNH nuclease is tightly dependent on the catalytic Mg2+, unveiling its cardinal structural role. This activated Mg2+-bound HNH is consistently described through molecular simulations, solution NMR and DNA cleavage assays, revealing also that the protonation state of the catalytic H840 is strongly affected by active site mutations. Finally, ab-initio QM(DFT)/MM simulations and metadynamics establish the catalytic mechanism, showing that the catalysis is activated by H840 and completed by K866, rationalising DNA cleavage experiments. This information is critical to enhance the enzymatic function of CRISPR-Cas9 toward improved genome-editing.

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.

Generalized Born Implicit Solvent Models Do Not Reproduce Secondary Structures of De Novo Designed Glu/Lys Peptides

We test a range of standard generalized Born (GB) models and protein force fields for a set of five experimentally characterized, designed peptides comprising alternating blocks of glutamate and lysine, which have been shown to differ significantly in α-helical content. Sixty-five combinations of force fields and GB models are evaluated in >800 μs of molecular dynamics simulations. GB models generally do not reproduce the experimentally observed α-helical content, and none perform well for all five peptides. These results illustrate that these models are not usefully predictive in this context. These peptides provide a useful test set for simulation methods.

Novel US-CpHMD Protocol to Study the Protonation-Dependent Mechanism of the ATP/ADP Carrier

We have designed a protocol combining constant-pH molecular dynamics (CpHMD) simulations with an umbrella sampling (US) scheme (US-CpHMD) to study the mechanism of ADP/ATP transport (import and export) by their inner mitochondrial membrane carrier protein [ADP/ATP carrier (AAC)]. The US scheme helped overcome the limitations of sampling the slow kinetics involved in these substrates' transport, while CpHMD simulations provided an unprecedented realism by correctly capturing the associated protonation changes. The import of anionic substrates along the mitochondrial membrane has a strong energetic disadvantage due to a smaller substrate concentration and an unfavorable membrane potential. These limitations may have created an evolutionary pressure on AAC to develop specific features benefiting the import of ADP. In our work, the potential of mean force profiles showed a clear selectivity in the import of ADP compared to ATP, while in the export, no selectivity was observed. We also observed that AAC sequestered both substrates at longer distances in the import compared to the export process. Furthermore, only in the import process do we observe transient protonation of both substrates when going through the AAC cavity, which is an important advantage to counteract the unfavorable mitochondrial membrane potential. Finally, we observed a substrate-induced disruption of the matrix salt-bridge network, which can promote the conformational transition (from the C- to M-state) required to complete the import process. This work unraveled several important structural features where the complex electrostatic interactions were pivotal to interpreting the protein function and illustrated the potential of applying the US-CpHMD protocol to other transport processes involving membrane proteins.

Protein pKa Prediction by Tree-Based Machine Learning

Protonation states of ionizable protein residues modulate many essential biological processes. For correct modeling and understanding of these processes, it is crucial to accurately determine their pKa values. Here, we present four tree-based machine learning models for protein pKa prediction. The four models, Random Forest, Extra Trees, eXtreme Gradient Boosting (XGBoost), and Light Gradient Boosting Machine (LightGBM), were trained on three experimental PDB and pKa datasets, two of which included a notable portion of internal residues. We observed similar performance among the four machine learning algorithms. The best model trained on the largest dataset performs 37% better than the widely used empirical pKa prediction tool PROPKA and 15% better than the published result from the pKa prediction method DelPhiPKa. The overall root-mean-square error (RMSE) for this model is 0.69, with surface and buried RMSE values being 0.56 and 0.78, respectively, considering six residue types (Asp, Glu, His, Lys, Cys, and Tyr), and 0.63 when considering Asp, Glu, His, and Lys only. We provide pKa predictions for proteins in human proteome from the AlphaFold Protein Structure Database and observed that 1% of Asp/Glu/Lys residues have highly shifted pKa values close to the physiological pH.

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.

Intramolecular interactions play key role in stabilization of pHLIP at acidic conditions

The pH-Low Insertion Peptide (pHLIP) is a membrane-active peptide that spontaneously folds into a transmembrane α-helix upon acidification. This activity enables pHLIP to potentially act as a vector for drugs related to diseases characterized by acidosis such as cancer or heart ischemia. Presently, due to aggregation-based effects, formulations of pHLIP are only viable at near-μM concentrations. In addition, since most of pHLIP's measurable qualities involve a membrane, probing the details of pHLIP in the interstitial region is difficult. In attempts to shed light on these issues, we performed constant pH molecular dynamics simulations on pHLIP as well as P20G, a variant with increased helicity, in solution at 0 and 150 mM NaCl over a broad range of pHs. In general, the addition of ions reduced the effective pK_a of the acidic residues in pHLIP. P20G exhibits a higher helicity than pHLIP in general and is more compact than pHLIP at pH values under 4. In terms of charge effects, sodium cations localized predominantly to the C-terminus of the peptide with a high density of acidic residues. Additionally, the salt bridge between R11 and D14 is by far the most favored and particularly so with pHLIP at 150 mM NaCl. We expect that this approach will be a valuable tool to screen variants of pHLIP for favorable properties in solution, an aspect of pHLIP design that to this point has largely been neglected.

pH Effects and Cooperativity among Key Titratable Residues for Escherichia coli Glycinamide Ribonucleotide Transformylase

Human glycinamide ribonucleotide transformylase (GAR Tfase) is a regulatory enzyme in the de novo purine biosynthesis pathway that has been extensively studied as an anticancer target. To some extent, inhibition of GAR Tfase selectively targets cancer cells over normal cells and inhibits purine formation and DNA replication. In this study, we investigated E. coli GAR Tfase, which shares high sequence similarity with the human GAR Tfase, and most functional residues are conserved. Herein, we aim to predict the pH-activity curve through a computational approach. We carried out pH-replica exchange molecular dynamics (pH-REMD) simulations to investigate pH-dependent functions such as structural changes, ligand binding, and catalytic activity. To compute the pH-activity curve, we identified the catalytic residues in specific protonation states, referred to as the catalytic competent protonation states (CCPS), which maintain the structure, keep ligands bound, and facilitate catalysis. Our computed population of CCPS with respect to pH matches well with the experimental pH-activity curve. To compute the microscopic pK_a values in the catalytically active conformation, we devised a thermodynamic model that considers the coupling between protonation states of CCPS residues and conformational states. These results allow us to correctly identify the general acid and base catalysts and interpret the pH-activity curve at an atomistic level.

Charge Interactions in a Highly Charge-Depleted Protein

Electrostatic forces are important for protein folding and are favored targets of protein engineering. However, interactions between charged residues are difficult to study because of the complex network of interactions found in most proteins. We have designed a purposely simple system to investigate this problem by systematically introducing individual and pairs of charged and titratable residues in a protein otherwise free of such residues. We used constant pH molecular dynamics simulations, NMR spectroscopy, and thermodynamic double mutant cycles to probe the structure and energetics of the interaction between the charged residues. We found that the partial burial of surface charges contributes to a shift in pK_a value, causing an aspartate to titrate in the neutral pH range. Additionally, the interaction between pairs of residues was found to be highly context dependent, with some pairs having no apparent preferential interaction, while other pairs would engage in coupled titration forming a highly stabilized salt bridge. We find good agreement between experiments and simulations and use the simulations to rationalize our observations and to provide a detailed mechanistic understanding of the electrostatic interactions.

pH-Induced Local Unfolding of the Phl p 6 Pollen Allergen From cpH-MD

Susceptibility to endosomal degradation is a decisive contribution to a protein's immunogenicity. It is assumed that the processing kinetics of structured proteins are inherently linked to their probability of local unfolding. In this study, we quantify the impact of endosomal acidification on the conformational stability of the major timothy grass pollen allergen Phl p 6. We use state of the art sampling approaches in combination with constant pH MD techniques to profile pH-dependent local unfolding events in atomistic detail. Integrating our findings into the current view on type 1 allergic sensitization, we characterize local protein dynamics in the context of proteolytic degradation at neutral and acidic pH for the wild type protein and point mutants with varying proteolytic stability. We analyze extensive simulation data using Markov state models and retrieve highly reliable thermodynamic and kinetic information at varying pH levels. Thereby we capture the impact of endolysosomal acidification on the structure and dynamics of the Phl p 6 mutants. We find that upon protonation at lower pH values, the conformational flexibilities in key areas of the wild type protein, i.e., T-cell epitopes and early proteolytic cleavage sites, increase significantly. A decrease of the pH even leads to local unfolding in otherwise stable secondary structure elements, which is a prerequisite for proteolytic cleavage. This effect is even more pronounced in the destabilized mutant, while no unfolding was observed for the stabilized mutant. In summary, we report detailed structural models which rationalize the experimentally observed cleavage pattern during endosomal acidification.

On the Use of the Discrete Constant pH Molecular Dynamics to Describe the Conformational Space of Peptides

Solvent pH is an important property that defines the protonation state of the amino acids and, therefore, modulates the interactions and the conformational space of the biochemical systems. Generally, this thermodynamic variable is poorly considered in Molecular Dynamics (MD) simulations. Fortunately, this lack has been overcome by means of the Constant pH Molecular Dynamics (CPHMD) methods in the recent decades. Several studies have reported promising results from these approaches that include pH in simulations but focus on the prediction of the effective pKa of the amino acids. In this work, we want to shed some light on the CPHMD method and its implementation in the AMBER suitcase from a conformational point of view. To achieve this goal, we performed CPHMD and conventional MD (CMD) simulations of six protonatable amino acids in a blocked tripeptide structure to compare the conformational sampling and energy distributions of both methods. The results reveal strengths and weaknesses of the CPHMD method in the implementation of AMBER18 version. The change of the protonation state according to the chemical environment is presumably an improvement in the accuracy of the simulations. However, the simulations of the deprotonated forms are not consistent, which is related to an inaccurate assignment of the partial charges of the backbone atoms in the CPHMD residues. Therefore, we recommend the CPHMD methods of AMBER program but pointing out the need to compare structural properties with experimental data to bring reliability to the conformational sampling of the simulations.

pH-Dependent Conformational Changes Lead to a Highly Shifted pKa for a Buried Glutamic Acid Mutant of SNase

Ionizable residues are rarely present in the hydrophobic interior of proteins, but when they are, they play important roles in biological processes such as energy transduction and enzyme catalysis. Internal ionizable residues have anomalous experimental pK_a values with respect to their pK_a in bulk water. This work investigates the atomistic cause of the highly shifted pK_a of the internal Glu23 in the artificially mutated variant V23E of Staphylococcal Nuclease (SNase) using pH replica exchange molecular dynamics (pH-REMD) simulations. The pK_a of Glu23 obtained from our calculations is 6.55, which is elevated with respect to the glutamate pK_a of 4.40 in bulk water. The calculated value is close to the experimental pK_a of 7.10. Our simulations show that the highly shifted pK_a of Glu23 is the product of a pH-dependent conformational change, which has been observed experimentally and also seen in our simulations. We carry out an analysis of this pH-dependent conformational change in response to the protonation state change of Glu23. Using a four-state thermodynamic model, we estimate the two conformation-specific pK_a values of Glu23 and describe the coupling between the conformational and ionization equilibria.

Improved GROMOS 54A7 Charge Sets for Phosphorylated Tyr, Ser, and Thr to Deal with pH-Dependent Binding Phenomena

Phosphorylation is a ubiquitous post-translational modification in proteins, and the phosphate group is present constitutively or transiently in most biological building blocks. These phosphorylated biomolecules are involved in many high-affinity binding/unbinding events that rely predominantly on electrostatic interactions. To build accurate models of these molecules, we need an improved description of the atomic partial charges for all relevant protonation states. In this work, we showed that the commonly used protocols to derive atomic partial charges using well-solvated molecules are inadequate to model the protonation equilibria in binding events. We introduced a protocol based on PB/MC calculations with a single representative conformation (of both protonation states) and used the resulting pK_a estimations to help manually curate the atomic partial charges. The final charge set, which is fully compatible with the GROMOS 54A7 force field, proved to be very effective in modeling the protonation equilibrium in different phosphorylated peptides in the free (tetrapeptides, pentapeptides, and pY1021) and protein-complexed forms (pY1021/PLC-γ1 complex). This was particularly important in the case of the pY1021 bound to the SH2 domain of PLC-γ1, where only our curated charge set captured the correct protonation equilibrium at the neutral to slightly acidic pH range. The binding/unbinding phenomena in that pH range are biologically relevant, and to improve our models, we need to go beyond the commonly used protocols and obtain revised force field parameters for these molecules.

pH-Dependent absorption spectrum of oxyluciferin analogues in the active site of firefly luciferase

In the quest for the identification of the light emitter(s) responsible for the firefly bioluminescence, the study of oxyluciferin analogues with controlled chemical and electronic structures is of particular importance. In this article, we report the results of our experimental and computational investigation of the pH-dependent absorption spectra characterizing three analogues bound into the luciferase cavity, together with adenosine-monophosphate (AMP). While the analogue microscopic pKa values do not differ much from their reference values, it turns out that the AMP protonation state is analogue-dependent and never doubly-deprotonated. A careful analysis of the interactions evidences the main role of E344 glutamic acid, as well as the flexibility of the cavity which can accommodate any oxyluciferin analogue. The consideration of the absorption spectra suggests that the oxyluciferin enolate form has to be excluded from the list of the bioluminescence reaction products.

Exploring pH Dependent Host/Guest Binding Affinities

When the electrostatic environment surrounding binding partners changes between unbound and bound states, the net uptake or release of a proton is possible by either binding partner. This process is pH-dependent in that the free energy required to uptake or release the proton varies with pH. This pH-dependence is typically not considered in conventional free energy methods where the use of fixed protonation states is the norm. In the present paper, we apply a simple two-step approach to calculate the pH-dependent binding free energy of a model cucubit[7]uril host/guest system. By use of λ-dynamics with an enhanced sampling protocol, adaptive landscape flattening, pK_a shifts and reference binding free energies upon complexation were determined. This information enables the construction of pH-dependent binding profiles that accurately capture the pK_a shifts and reproduce binding free energies at the different pH conditions that were observed experimentally. Our calculations illustrate a general framework for computing pH-dependent binding free energies but also point to some issues in modeling the molecular charge distributions within this series of molecules with CGenFF. However, by introducing some minor charge modifications to the CGenFF force field, we saw significant improvement in accuracy of the calculated pK_a shifts.

Predicting Reactive Cysteines with Implicit-Solvent-Based Continuous Constant pH Molecular Dynamics in Amber

Cysteines existing in the deprotonated thiolate form or having a tendency to become deprotonated are important players in enzymatic and cellular redox functions and frequently exploited in covalent drug design; however, most computational studies assume cysteines as protonated. Thus, developing an efficient tool that can make accurate and reliable predictions of cysteine protonation states is timely needed. We recently implemented a generalized Born (GB) based continuous constant pH molecular dynamics (CpHMD) method in Amber for protein pKa calculations on CPUs and GPUs. Here we benchmark the performance of GB-CpHMD for predictions of cysteine pKa's and reactivities using a data set of 24 proteins with both down- and upshifted cysteine pKa's. We found that 10 ns single-pH or 4 ns replica-exchange CpHMD titrations gave root-mean-square errors of 1.2-1.3 and correlation coefficients of 0.8-0.9 with respect to experiment. The accuracy of predicting thiolates or reactive cysteines at physiological pH with single-pH titrations is 86 or 81% with a precision of 100 or 90%, respectively. This performance well surpasses the traditional structure-based methods, particularly a widely used empirical pKa tool that gives an accuracy less than 50%. We discuss simulation convergence, dependence on starting structures, common determinants of the pKa downshifts and upshifts, and the origin of the discrepancies from the structure-based calculations. Our work suggests that CpHMD titrations can be performed on a desktop computer equipped with a single GPU card to predict cysteine protonation states for a variety of applications, from understanding biological functions to covalent drug design.

Fully Atomistic Multiscale Approach for pKa Prediction

The ionization state of titratable amino acids strongly affects proteins structure and functioning in a large number of biological processes. It is therefore essential to be able to characterize the pK_a of ionizable groups inside proteins and to understand its microscopic determinants in order to gain insights into many functional properties of proteins. A big effort has been devoted to the development of theoretical approaches for the prediction of deprotonation free energies, yet the accurate theoretical/computational calculation of pK_a values is recognized as a current challenge. A methodology based on a hybrid quantum/classical approach is here proposed for the computation of deprotonation free energies. The method is applied to calculate the pK_a of formic acid, methylammonium, and methanethiol, providing results in good agreement with the corresponding experimental estimates. The pK_a is also calculated for aspartic acid and lysine as single residues in solution and for three aspartic/glutamic acids inside a well-characterized protein: hen egg white lysozyme. While for small molecules the method is able to deal with multiple protonation states of all titratable groups, this becomes computationally very expensive for proteins. The calculated pK_a values for the single amino acids (except for the zwitterionic aspartic acid) and inside the protein display a systematic shift with respect to the experimental values that suggests that the fine balance between hydrophobic and polar interactions might be not accurately reproduced by the usual classical force-fields, thus affecting the computation of deprotonation free energies. The calculated pK_a shifts inside the protein are in good agreement with the corresponding experimental ones (within 1 pK_a unit), well reproducing the pK_a changes due to the protein environment even in the case of large pK_a shifts.

Catalytic Site pKa Values of Aspartic, Cysteine, and Serine Proteases: Constant pH MD Simulations

Enzymatic function and activity of proteases is closely controlled by the pH value. The protonation states of titratable residues in the active site react to changes in the pH value, according to their pK_a, and thereby determine the functionality of the enzyme. Knowledge of the titration behavior of these residues is crucial for the development of drugs targeting the active site residues. However, experimental pK_a data are scarce, since the systems’ size and complexity make determination of these pK_a values inherently difficult. In this study, we use single pH constant pH MD simulations as a fast and robust tool to estimate the active site pK_a values of a set of aspartic, cysteine, and serine proteases. We capture characteristic pK_a shifts of the active site residues, which dictate the experimentally determined activity profiles of the respective protease family. We find clear differences of active site pK_a values within the respective families, which closely match the experimentally determined pH preferences of the respective proteases. These shifts are caused by a distinct network of electrostatic interactions characteristic for each protease family. While we find convincing agreement with experimental data for serine and aspartic proteases, we observe clear deficiencies in the description of the titration behavior of cysteines within the constant pH MD framework and highlight opportunities for improvement. Consequently, with this work, we provide a concise set of active site pK_a values of aspartic and serine proteases, which could serve as reference for future theoretical as well as experimental studies.

Exploring Coupled Redox and pH Processes with a Force-Field-Based Approach: Applications to Five Different Systems

Coupled redox and pH-driven processes are at the core of many important biological mechanisms. As the distribution of protonation and redox states in a system is associated with the pH and redox potential of the solution, having efficient computational tools that can simulate under these conditions becomes very important. Such tools have the potential to provide information that complement and drive experiments. In previous publications we have presented the implementation of the constant pH and redox potential molecular dynamics (C(pH,E)MD) method in AMBER and we have shown how multidimensional replica exchange can be used to significantly enhance the convergence efficiency of our simulations. In the current work, after an improvement in our C(pH,E)MD approach that allows a given residue to be simultaneously pH- and redox-active, we have employed our methodologies to study five different systems of interest in the literature. We present results for capped tyrosine dipeptide, two maquette systems containing one pH- and redox-active tyrosine (α₃Y and peptide A), and two proteins that contain multiple heme groups (diheme cytochrome c from Rhodobacter sphaeroides and Desulfovibrio vulgaris Hildenborough cytochrome c₃). We show that our results can provide new insights into previous theoretical and experimental findings by using a fully force-field-based and GPU-accelerated approach, which allows the simulations to be executed with high computational performance.

Effect of curcumin derivatives on hen egg white lysozyme amyloid fibrillation and their interaction study by spectroscopic methods

Two novel Boc-L-isoleucine-functionalized curcumin derivatives have been synthesized and characterized, which exhibited enhanced solubility in water compared with the natural curcumin. The solubility could reach 2.12mg/mL for the monosubstituted compound and 3.05mg/mL for the disubstituted compound, respectively. Their anti-amyloidogenic capacity on the model protein, hen egg white lysozyme (HEWL), was examined in aqueous solution. ThT fluorescence assay showed that the operation concentration was only 0.5mM when the inhibition ratio was above 70%. Meanwhile, the inhibitory capacity of monosubstituted curcumin derivative on the formation of HEWL amyloid fibrils was found to be superior to that of disubstituted derivative, suggesting that the phenolic hydroxyl group might contribute to the anti-amyloidogenic activity. Interaction study showed that both curcumin derivatives could bind with HEWL near tryptophan residues and form new ground-state complex before HEWL self-assemblies into amyloid fibrils and thus inhibits the formation of amyloid fibrils. Both of the two cucumin derivatives have displayed low cytotoxicity with HeLa cell.

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.

pH-Dependent Protonation of the Phl p 6 Pollen Allergen Studied by NMR and cpH-aMD

We use state-of-the-art NMR experiments to measure apparent pK_a values in the native protein environment and employ a cutting-edge combination of enhanced sampling and constant pH molecular dynamics (MD) simulations to rationalize strong pK_a shifts. The major timothy grass pollen allergen Phl p 6 serves as an ideal model system for both methods due to its high number of titratable residues despite its comparably small size. We present a proton transition analysis as intuitive tool to depict the captured protonation state ensemble in atomistic detail. Combining microscopic structural details from MD simulations and macroscopic ensemble averages from NMR shifts leads to a comprehensive view on pH dependencies of protonation states and tautomers. Overall, we find striking agreement between simulation-based pK_a predictions and experiment. However, our analyses suggest subtle differences in the underlying molecular origin of the observed pK_a shifts. From accelerated constant pH MD simulations, we identify immediate proximity of opposite charges, followed by vicinity of equal charges as major driving forces for pK_a shifts. NMR experiments on the other hand, suggest only a weak relation of pK_a shifts and close contacts to charged residues, while the strongest influence derives from the dipolar character of α helices. The presented study hence pinpoints opportunities for improvements concerning the theoretical description of protonation state and tautomer probabilities. However, the coherence in the resulting apparent pK_a values from simulations and experiment affirms cpH-aMD as a reliable tool to study allergen dynamics at varying pH levels.

CpHMD-Then-QM/MM Identification of the Amino Acids Responsible for the Anabaena Sensory Rhodopsin pH-Dependent Electronic Absorption Spectrum

Anabaena Sensory Rhodopsin (ASR), a microbial photoactive protein featuring the retinal chromophore in two different conformations, exhibits a pH-dependent electronic absorption spectrum. Using the recently developed CpHMD-then-QM/MM multiscale protocol applied to ASR embedded in a membrane model, the pH-induced changes in its maximum absorption wavelength have been reproduced and analyzed. While the acidic tiny red-shift is essentially correlated with the deprotonation of an aspartic acid located on the ASR extracellular side, the larger blue-shift experimentally reported at pH values larger than 5 involves a cluster of titrating residues sitting on the cytoplasmic side. The ASR pH-dependent spectrum is the consequence of the competitive stabilization of retinal ground and excited states by the protein electrostatic potential.

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.

Generalized Born Implicit Solvent Models for Biomolecules

It would often be useful in computer simulations to use an implicit description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation and can be very efficient compared to the explicit treatment of the solvent. Here, we review a particular class of so-called fast implicit solvent models, generalized Born (GB) models, which are widely used for molecular dynamics (MD) simulations of proteins and nucleic acids. These approaches model hydration effects and provide solvent-dependent forces with efficiencies comparable to molecular-mechanics calculations on the solute alone; as such, they can be incorporated into MD or other conformational searching strategies in a straightforward manner. The foundations of the GB model are reviewed, followed by examples of newer, emerging models and examples of important applications. We discuss their strengths and weaknesses, both for fidelity to the underlying continuum model and for the ability to replace explicit consideration of solvent molecules in macromolecular simulations.

Unfolding Pathways of Hen Egg-White Lysozyme in Ethanol

The aggregation of amyloid fibrils can lead to various diseases including Alzheimer's, Parkinson's disease, and transmissible spongiform encephalopathy. Amyloid fibrils can develop from a variety of proteins in the body as they misfold into a primarily β-sheet structure and aggregate. Human lysozyme has been shown to have far reaching effects in the human health-a homologous enzyme, hen egg-white lysozyme, has been shown to denature to a primarily β-sheet structure at low pH and high alcohol content solution. We have studied these systems in atomic-level detail with a combination of constant pH and microsecond long molecular dynamics simulation in explicit solvent, which cumulatively total over 10 μs of simulation time. These studies have allowed us to determine two potential unfolding pathways depending on the protonation state of a key glutamic acid residue as well as the effect of solution dynamics and pH on the unfolding process.

pH-dependent conformational dynamics of beta-secretase 1: A molecular dynamics study

Beta-secretase 1 (BACE-1) is an aspartyl protease implicated in the overproduction of β-amyloid fibrils responsible for Alzheimer disease. The process of β-amyloid genesis is known to be pH dependent, with an activity peak between solution pH of 3.5 and 5.5. We have studied the pH-dependent dynamics of BACE-1 to better understand the pH dependent mechanism. We have implemented support for graphics processor unit (GPU) accelerated constant pH molecular dynamics within the AMBER molecular dynamics software package and employed this to determine the relative population of different aspartyl dyad protonation states in the pH range of greatest β-amyloid production, followed by conventional molecular dynamics to explore the differences among the various aspartyl dyad protonation states. We observed a difference in dynamics between double-protonated, mono-protonated, and double-deprotonated states over the known pH range of higher activity. These differences include Tyr 71-aspartyl dyad proximity and active water lifetime. This work indicates that Tyr 71 stabilizes catalytic water in the aspartyl dyad active site, enabling BACE-1 activity.

Engineering Allostery into Proteins

Our ability to engineer protein structure and function has grown dramatically over recent years. Perhaps the next level in protein design is to develop proteins whose function can be regulated in response to various stimuli, including ligand binding, pH changes, and light. Endeavors toward these goals have tested and expanded on our understanding of protein function and allosteric regulation. In this chapter, we provide examples from different methods for developing new allosterically regulated proteins. These methods range from whole insertion of regulatory domains into new host proteins, to covalent attachment of photoswitches to generate light-responsive proteins, and to targeted changes to specific amino acid residues, especially to residues identified to be important for relaying allosteric information across the protein framework. Many of the examples we discuss have already found practical use in medical and biotechnology applications.

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.

Maintaining and breaking symmetry in homomeric coiled-coil assemblies

In coiled-coil (CC) protein structures α-helices wrap around one another to form rope-like assemblies. Most natural and designed CCs have two–four helices and cyclic (C_n) or dihedral (D_n) symmetry. Increasingly, CCs with five or more helices are being reported. A subset of these higher-order CCs is of interest as they have accessible central channels that can be functionalised; they are α-helical barrels. These extended cavities are surprising given the drive to maximise buried hydrophobic surfaces during protein folding and assembly in water. Here, we show that α-helical barrels can be maintained by the strategic placement of β-branched aliphatic residues lining the lumen. Otherwise, the structures collapse or adjust to give more-complex multi-helix assemblies without C_n or D_n symmetry. Nonetheless, the structural hallmark of CCs—namely, knobs-into-holes packing of side chains between helices—is maintained leading to classes of CCs hitherto unobserved in nature or accessed by design.

Higher order coiled coils with five or more helices can form α-helical barrels. Here the authors show that placing β-branched aliphatic residues along the lumen yields stable and open α-helical barrels, which is of interest for the rational design of functional proteins; whereas, the absence of β-branched side chains leads to unusual low-symmetry α-helical bundles.

Reservoir pH replica exchange

We present the reservoir pH replica exchange (R-pH-REM) method for constant pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a standard or enhanced molecular dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a constant pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from constant pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable molecules can sample two or more conformations separated by a large energy barrier, such that adequate sampling is not achieved on a time scale of standard constant pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small molecule with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calculations on this model system are in excellent agreement with a theoretical model.

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.

Membrane-Induced p Ka Shifts in wt-pHLIP and Its L16H Variant

The pH (low) insertion peptides (pHLIPs) is a family of peptides that are able to insert into a lipid bilayer at acidic pH. The molecular mechanism of pHLIPs insertion, folding, and stability in the membrane at low pH is based on multiple protonation events, which are challenging to study at the molecular level. More specifically, the relation between the experimental p K of insertion (p K^exp) of pHLIPs and the p K_a of the key residues is yet to be clarified. We carried out a computational study, complemented with new experimental data, and established the influence of (de)protonation of titrable residues on the stability of the peptide membrane-inserted state. Constant-pH molecular dynamics simulations were employed to calculate the p K_a values of these residues along the membrane normal. In the wt-pHLIP, we identified Asp14 as the key residue for the stability of the membrane-inserted state, and its p K_a value is strongly correlated with the experimental p K^exp measured in thermodynamics studies. Also, in order to narrow down the pH range at which pHLIP is stable in the membrane, we designed a new pHLIP variant, L16H, where Leu in the 16th position was replaced by a titrable His residue. Our results showed that the L16H variant undergoes two transitions. The calculated p K_a and experimentally observed p K^exp values are in good agreement. Two distinct p K^exp values delimit a pH range where the L16H peptide is stably inserted in the membrane, while, outside this range, the membrane-inserted state is destabilized and the peptide exits from the bilayer. pHLIP peptides have been successfully used to target cancer cells for the delivery of diagnostics and therapeutic agents to acidic tumors. The fine-tuning of the stability of the pHLIP inserted state and its restriction to a narrow well-defined pH range might allow the design of new peptides, able to discriminate between tissues with different extracellular pH values.

Role of Counterions in Constant-pH Molecular Dynamics Simulations of PAMAM Dendrimers

Electrostatic interactions play a pivotal role in the structure and mechanism of action of most biomolecules. There are several conceptually different methods to deal with electrostatics in molecular dynamics simulations. Ionic strength effects are usually introduced using such methodologies and can have a significant impact on the quality of the final conformation space obtained. We have previously shown that full system neutralization can lead to wrong lipidic phases in the 25% PA/PC bilayer (J. Chem. Theory Comput. 2014,10, 5483-5492). In this work, we investigate how two limit approaches to the ionic strength treatment (implicitly with GRF or using full system neutralization with either GRF or PME) can influence the conformational space of the second-generation PAMAM dendrimer. Constant-pH MD simulations were used to map PAMAM's conformational space at its full pH range (from 2.5 to 12.5). Our simulations clearly captured the coupling between protonation and conformation in PAMAM. Interestingly, the dendrimer conformational distribution was almost independent of the ionic strength treatment methods, which is in contrast to what we have observed in charged lipid bilayers. Overall, our results confirm that both GRF with implicit ionic strength and a fully neutralized system with PME are valid approaches to model charged globular systems, using the GROMOS 54A7 force field.

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.