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

Constant pH Simulation with FMM Electrostatics in GROMACS. (A) Design and Applications

The structural dynamics of biological macromolecules, such as proteins, DNA/RNA, or complexes thereof, are strongly influenced by protonation changes of their typically many titratable groups, which explains their sensitivity to pH changes. Conversely, conformational and environmental changes of the biomolecule affect the protonation state of these groups. With few exceptions, conventional force field-based molecular dynamics (MD) simulations neither account for these effects nor do they allow for coupling to a pH buffer. Here, we present design decisions and applications of a rigorous Hamiltonian interpolation λ-dynamics constant pH method in GROMACS, which rests on GPU-accelerated Fast Multipole Method (FMM) electrostatics. Our implementation supports both CHARMM36m and Amber99sb*-ILDN force fields and is largely automated to enable seamless switching from regular MD to constant pH MD, involving minimal changes to the input files. Here, the first of two companion papers describes the underlying constant pH protocol and sample applications to several prototypical benchmark systems such as cardiotoxin V, lysozyme, and staphylococcal nuclease. Enhanced convergence is achieved through a new dynamic barrier height optimization method, and high pKa accuracy is demonstrated. We use Functional Mode Analysis (FMA) and Mutual Information (MI) to explore the complex intra- and intermolecular couplings between the protonation states of titratable groups as well as those between protonation states and conformational dynamics. We identify striking conformation-dependent pKa variations and unexpected inter-residue couplings. Conformation-protonation coupling is identified as a primary cause of the slow protonation convergence notorious to constant pH simulations involving multiple titratable groups, suggesting enhanced sampling methods to accelerate convergence.

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.

Benchmarking In Silico Tools for Cysteine pKa Prediction

Accurate estimation of the pK_a's of cysteine residues in proteins could inform targeted approaches in hit discovery. The pK_a of a targetable cysteine residue in a disease-related protein is an important physiochemical parameter in covalent drug discovery, as it influences the fraction of nucleophilic thiolate amenable to chemical protein modification. Traditional structure-based in silico tools are limited in their predictive accuracy of cysteine pK_a's relative to other titratable residues. Additionally, there are limited comprehensive benchmark assessments for cysteine pK_a predictive tools. This raises the need for extensive assessment and evaluation of methods for cysteine pK_a prediction. Here, we report the performance of several computational pK_a methods, including single-structure and ensemble-based approaches, on a diverse test set of experimental cysteine pK_a's retrieved from the PKAD database. The dataset consisted of 16 wildtype and 10 mutant proteins with experimentally measured cysteine pK_a values. Our results highlight that these methods are varied in their overall predictive accuracies. Among the test set of wildtype proteins evaluated, the best method (MOE) yielded a mean absolute error of 2.3 pK units, highlighting the need for improvement of existing pK_a methods for accurate cysteine pK_a estimation. Given the limited accuracy of these methods, further development is needed before these approaches can be routinely employed to drive design decisions in early drug discovery efforts.

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.