Reliable and Accurate Prediction of Single-Residue pKa Values through Free Energy Perturbation Calculations

Improving pKa Predictions with Reparameterized Force Fields and Free Energy Calculations

Given the growing interest in designing targeted covalent inhibitors, methods for rapidly and accurately probing pKas─and, by extension, the reactivities─of target cysteines are highly desirable. Complementary to cysteine, histidine is similarly relevant due to its frequent presence in protein active sites and its unique ability to exist in two tautomeric states. Here, we demonstrate that nonequilibrium free energy calculations can accurately determine the pKa values of both residues, often outperforming conventional predictors. Importantly, we find that (1) increasing the van der Waals radius of cysteine's sulfur atom, (2) modifying the backbone charges of histidine, and (3) introducing effective polarization by downscaling the side chain partial charges of both residues can all significantly improve pKa prediction accuracy. Using the modified CHARMM36m force field on the full dataset reduces the prediction error from 2.12 ± 0.27 pK to 1.28 ± 0.15 pK and increases the correlation with experiment from 0.25 ± 0.09 to 0.58 ± 0.08. Similarly, using the modified Amber14SB force field decreases the error from 3.21 ± 0.29 pK to 1.69 ± 0.23 pK and improves the correlation from 0.15 ± 0.10 to 0.36 ± 0.10.

Methods for Theoretical Treatment of Local Fields in Proteins and Enzymes

Electric fields generated by protein scaffolds are crucial in enzymatic catalysis. This review surveys theoretical approaches for detecting, analyzing, and comparing electric fields, electrostatic potentials, and their effects on the charge density within enzyme active sites. Pioneering methods like the empirical valence bond approach rely on evaluating ionic and covalent resonance forms influenced by the field. Strategies employing polarizable force fields also facilitate field detection. The vibrational Stark effect connects computational simulations to experimental Stark spectroscopy, enabling direct comparisons. We highlight how protein dynamics induce fluctuations in local fields, influencing enzyme activity. Recent techniques assess electric fields throughout the active site volume rather than only at specific bonds, and machine learning helps relate these global fields to reactivity. Quantum theory of atoms in molecules captures the entire electron density landscape, providing a chemically intuitive perspective on field-driven catalysis. Overall, these methodologies show protein-generated fields are highly dynamic and heterogeneous, and understanding both aspects is critical for elucidating enzyme mechanisms. This holistic view empowers rational enzyme engineering by tuning electric fields, promising new avenues in drug design, biocatalysis, and industrial applications. Future directions include incorporating electric fields as explicit design targets to enhance catalytic performance and biochemical functionalities.

Using AlphaFold2 to Predict the Conformations of Side Chains in Folded Proteins

AlphaFold has revolutionized protein structure prediction by accurately creating 3D structures from just the amino acid sequence. However, even with extensive research validating its overall accuracy, a key question remains: Can AlphaFold predict the conformation of individual amino acid residue side chains within a folded protein? This is important for the field of molecular modeling, particularly when predicting the effects of mutations on protein stability and ligand binding. AlphaFold generates a set of atomic coordinates not just for the mutated side chain but also for potential rearrangements across the entire protein structure. In this study we investigate the ability of ColabFold, an online implementation of AlphaFold2 (AF2), to predict the conformations of residue side chains in folded proteins. We find that over a set of 10 benchmark proteins, the side chain conformation prediction error of ColabFold is on average ~ 14 % forχ 1 dihedral angles, and increases to ~ 48 % forχ 3 dihedral angles. The prediction error is smaller for non-polar side chains and is somewhat improved using structural templates. ColabFold demonstrates a bias towards the most prevalent rotamer states in the protein data bank (PDB), potentially limiting its ability to capture rare side chain conformations effectively. As an application of AlphaFold to explore the structural consequences of strongly cooperative mutations on side chain rearrangements, we employ a Potts sequence-based statistical energy model to perform large scale mutational scans of two proteins ABL1 and PIM1 kinase, searching for the most strongly cooperative mutational pairs, and then use ColabFold to predict the structural signatures of this cooperativity on the interacting side chains. Our results demonstrate that integration of the sequence-based Potts model with AlphaFold into a single pipeline provides a new tool that can be used to explore the fundamental relationship between protein mutations, cooperative changes in structure, and fitness.

Comparative Assessment of Water Models in Protein-Glycan Interaction: Insights from Alchemical Free Energy Calculations and Molecular Dynamics Simulations

Accurate computational simulations of protein-glycan dynamics are crucial for a comprehensive understanding of critical biological mechanisms, including host-pathogen interactions, immune system defenses, and intercellular communication. The accuracy of these simulations, including molecular dynamics (MD) simulation and alchemical free energy calculations, critically relies on the appropriate parameters, including the water model, because of the extensive hydrogen bonding with glycan hydroxyl groups. However, a systematic evaluation of water models' accuracy in simulating protein-glycan interaction at the molecular level is still lacking. In this study, we used full atomistic MD simulations and alchemical absolute binding free energy (ABFE) calculations to investigate the performance of five distinct water models in six protein-glycan complex systems. We evaluated water models' impact on structural dynamics and binding affinity through over 5.8 μs of simulation time per system. Our results reveal that most protein-glycan complexes are stable in the overall structural dynamics regardless of the water model used, while some show obvious fluctuations with specific water models. More importantly, we discover that the stability of the binding motif's conformation is dependent on the water model chosen when its residues form weak hydrogen bonds with the glycan. The water model also influences the conformational stability of the glycan in its bound state according to density functional theory (DFT) calculations. Using alchemical ABFE calculations, we find that the OPC water model exhibits exceptional consistency with experimental binding affinity data, whereas commonly used models such as TIP3P are less accurate. The findings demonstrate how different water models affect protein-glycan interactions and the accuracy of binding affinity calculations, which is crucial in developing therapeutic strategies targeting these interactions.

Enhanced and Efficient Predictions of Dynamic Ionization through Constant-pH Adiabatic Free Energy Dynamics

Dynamic or structurally induced ionization is a critical aspect of many physical, chemical, and biological processes. Molecular dynamics (MD) based simulation approaches, specifically constant pH MD methods, have been developed to simulate ionization states of molecules or proteins under experimentally or physiologically relevant conditions. While such approaches are now widely utilized to predict ionization sites of macromolecules or to study physical or biological phenomena, they are often computationally expensive and require long simulation times to converge. In this article, using the principles of adiabatic free energy dynamics, we introduce an efficient technique for performing constant pH MD simulations within the framework of the adiabatic free energy dynamics (AFED) approach. We call the new approach pH-AFED. We show that pH-AFED provides highly accurate predictions of protein residue pKa values, with a MUE of 0.5 pKa units when coupled with driven adiabatic free energy dynamics (d-AFED), while reducing the required simulation times by more than an order of magnitude. In addition, pH-AFED can be easily integrated into most constant pH MD codes or implementations and flexibly adapted to work in conjunction with enhanced sampling algorithms that target collective variables. We demonstrate that our approaches, with both pH-AFED standalone as well as pH-AFED combined with collective variable based enhanced sampling, provide promising predictive accuracy, with a MUE of 0.6 and 0.5 pKa units respectively, on a diverse range of proteins and enzymes, ranging up to 186 residues and 21 titratable sites. Lastly, we demonstrate how this approach can be utilized to understand the in vivo performance engineered antibodies for immunotherapy.

Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery

Computational methods constitute efficient strategies for screening and optimizing potential drug molecules. A critical factor in this process is the binding affinity between candidate molecules and targets, quantified as binding free energy. Among various estimation methods, alchemical transformation methods stand out for their theoretical rigor. Despite challenges in force field accuracy and sampling efficiency, advancements in algorithms, software, and hardware have increased the application of free energy perturbation (FEP) calculations in the pharmaceutical industry. Here, we review the practical applications of FEP in drug discovery projects since 2018, covering both ligand-centric and residue-centric transformations. We show that relative binding free energy calculations have steadily achieved chemical accuracy in real-world applications. In addition, we discuss alternative physics-based simulation methods and the incorporation of deep learning into free energy calculations.

Robust Prediction of Relative Binding Energies for Protein-Protein Complex Mutations Using Free Energy Perturbation Calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how Protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

A Comprehensive Study of Factors Affecting the Prediction of the pKa Shift of Asp26 in Thioredoxin Protein

The stable protonation state of ionizable amino acids in a protein can be predicted by computing the pKa shift of that residue within the protein environment. Thermodynamic Integration (TI) is an ideal molecular dynamics-based approach for predicting the pKa shift of ionizable protein residues. Here, we probe TI-based simulation protocols for their ability to accurately predict the pKa shift of Asp26 in thioredoxin. While implicit solvent models can predict the pKa shift accurately, explicit solvent models result in substantial errors. To understand the underlying reason for this surprising discrepancy, we investigate the role of various factors such as solvent models, conformational sampling, background charges, and polarization.

Mining Medicinally Relevant Bioreduction Substrates Inspired by Ligand-Based Drug Design

Exploring the scope of biocatalytic transformations in the absence of enzyme structures without extensive experimentation is a challenging task. To expand the limited substrate capacity of carrot-mediated bioreduction and hunt for new medicinally relevant ketones with minimum cost of labor and time, we deployed a practical method inspired by ligand-based drug design. Through analyzing collected literature data and building pharmacophore and reactivity prediction models, we screened a self-built virtual library of >8000 ketones bearing the most frequently used N,O,S-heterocycles and functional groups in drug discovery. Representative examples were validated, expanding the bioreduction substrate scope. The public availability of our models alongside the straightforward screening workflow makes it time-, labor-, and cost-saving to evaluate unknown bioreduction substrates for medicinal chemistry applications, especially for a large set of structurally differentiated ketones. Our studies also showcase the novelty of utilizing medicinal chemistry principles to solve a general biocatalysis problem.

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.

Robust prediction of relative binding energies for protein-protein complex mutations using free energy perturbation calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Titratable residues that drive RND efflux: Insights from molecular simulations

The resistance-nodulation-division efflux machinery confers antimicrobial resistance to Gram-negative bacteria by actively pumping antibiotics out of the cell. The protein complex is powered by proton motive force; however, the proton transfer mechanism itself and indeed even its stoichiometry is still unclear. Here we review computational studies from the last decade that focus on elucidating the number of protons transferred per conformational cycle of the pump. Given the difficulties in studying proton movement using even state-of-the-art structural biology methods, the contributions from computational studies have been invaluable from a mechanistic perspective.

Potts Hamiltonian Models and Molecular Dynamics Free Energy Simulations for Predicting the Impact of Mutations on Protein Kinase Stability

Single-point mutations in kinase proteins can affect their stability and fitness, and computational analysis of these effects can provide insights into the relationships among protein sequence, structure, and function for this enzyme family. To assess the impact of mutations on protein stability, we used a sequence-based Potts Hamiltonian model trained on a kinase family multiple-sequence alignment (MSA) to calculate the statistical energy (fitness) effects of mutations and compared these against relative folding free energies (ΔΔGs) calculated from all-atom molecular dynamics free energy perturbation (FEP) simulations in explicit solvent. The fitness effects of mutations in the Potts model (ΔEs) showed good agreement with experimental thermostability data (Pearson r = 0.68), similar to the correlation we observed with ΔΔGs predicted from structure-based relative FEP simulations. Recognizing the possible advantages of using Potts models to rapidly estimate protein stability effects of kinase mutations seen in cancer genomics data, we used the Potts statistical energy model to estimate the stability effects of 65 conservative and nonconservative mutations across three distinct kinases (Wee1, Abl1, and Cdc7) with somatic mutations reported in the Genomic Data Commons (GDC) database. The ΔEs of these mutations calculated from the Potts model are consistent with the corresponding ΔΔGs from FEP simulations (Pearson ratio of 0.72). The agreement between these methods suggests that the Potts model may be used as a sequence-based tool for high-throughput screening of mutational effects as part of a computational pipeline for predicting the stability effects of mutations. We also demonstrate how the scalability of the fitness-based Potts model calculations permits analyses that are not easily accessed using FEP simulations. To this end, we employed site-saturation mutagenesis in the Potts model in order to investigate the relative stability effects of mutations seen in different cancer evolutionary scenarios. We used this approach to analyze the effects of drug pressure in Abl kinase by contrasting the relative fitness penalties of somatic mutations seen in miscellaneous cancer types with those calculated for mutations associated with cancer drug resistance. We observed that, in contrast to somatic mutations of Abl seen in various tumors that appear to have evolved neutrally, cancer mutations that evolved under drug pressure in Abl-targeted therapies tend to preserve enzyme stability.

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.

Direct measurements of biomolecular electrostatics through experiments

Biomolecular electrostatics has been a subject of computational investigations based on 3D structures. This situation is changing because emerging experimental tools allow us to quantitatively investigate biomolecular electrostatics without any use of structure information. Now, electrostatic potentials around biomolecules can directly be measured for many residues simultaneously by nuclear magnetic resonance (NMR) spectroscopy. This NMR method can be used to study electrostatic aspects of various processes, including macromolecular association and liquid-liquid phase separation. Applications to structurally flexible biomolecules such as intrinsically disordered proteins are particularly useful. The new tools also facilitate examination of theoretical models and methods for biomolecular electrostatics.