Approach to Study pH-Dependent Protein Association Using Constant-pH Molecular Dynamics: Application to the Dimerization of β-Lactoglobulin

Multistate Method to Efficiently Account for Tautomerism and Protonation in Alchemical Free-Energy Calculations

The majority of drug-like molecules contain at least one ionizable group, and many common drug scaffolds are subject to tautomeric equilibria. Thus, these compounds are found in a mixture of protonation and/or tautomeric states at physiological pH. Intrinsically, standard classical molecular dynamics (MD) simulations cannot describe such equilibria between states, which negatively impacts the prediction of key molecular properties in silico. Following the formalism described by de Oliveira and co-workers (J. Chem. Theory Comput. 2019, 15, 424-435) to consider the influence of all states on the binding process based on alchemical free-energy calculations, we demonstrate in this work that the multistate method replica-exchange enveloping distribution sampling (RE-EDS) is well suited to describe molecules with multiple protonation and/or tautomeric states in a single simulation. We apply our methodology to a series of eight inhibitors of factor Xa with two protonation states and a series of eight inhibitors of glycogen synthase kinase 3β (GSK3β) with two tautomeric states. In particular, we show that given a sufficient phase-space overlap between the states, RE-EDS is computationally more efficient than standard pairwise free-energy methods.

Anions and Cations Affect Amino Acid Dissociation Equilibria via Distinct Mechanisms

Salts reduce the pKa of weak acids by a mechanism sensitive to ion identity and concentration via charge screening of the deprotonated state. In this study, we utilize constant pH molecular dynamics simulations to understand the molecular mechanism behind the salt-dependent dissociation of aspartic acid (Asp). We calculate the pKa of Asp in the presence of a monovalent salt and investigate Hofmeister ion effects by systematically varying the ionic radii. We observe that increasing the anion size leads to a monotonic decrease in Asp pKa. Conversely, the cation size affects the pKa nonmonotonically, interpretable in the context of the law of matching water affinity. The net effect of salt on Asp acidity is governed by an interplay of solvation and competing ion interactions. The proposed mechanism is rather general and can be applicable to several problems in Hofmeister ion chemistry, such as pH effects on protein stability and soft matter interfaces.

Tools for analyzing protonation states and for tracing proton transfer pathways with examples from the Rb. sphaeroides photosynthetic reaction centers

Protons participate in many reactions. In proteins, protons need paths to move in and out of buried active sites. The vectorial movement of protons coupled to electron transfer reactions establishes the transmembrane electrochemical gradient used for many reactions, including ATP synthesis. Protons move through hydrogen bonded chains of waters and hydroxy side chains via the Grotthuss mechanism and by proton binding and release from acidic and basic residues. MCCE analysis shows that proteins exist in a large number of protonation states. Knowledge of the equilibrium ensemble can provide a rational basis for setting protonation states in simulations that fix them, such as molecular dynamics (MD). The proton path into the Q_B site in the bacterial reaction centers (RCs) of Rb. sphaeroides is analyzed by MD to provide an example of the benefits of using protonation states found by the MCCE program. A tangled web of side chains and waters link the cytoplasm to Q_B. MCCE analysis of snapshots from multiple trajectories shows that changing the input protonation state of a residue in MD biases the trajectory shifting the proton affinity of that residue. However, the proton affinity of some residues is more sensitive to the input structure. The proton transfer networks derived from different trajectories are quite robust. There are some changes in connectivity that are largely restricted to the specific residues whose protonation state is changed. Trajectories with Q_B^•- are compared with earlier results obtained with Q_B [Wei et. al Photosynthesis Research volume 152, pages153-165 (2022)] showing only modest changes. While introducing new methods the study highlights the difficulty of establishing the connections between protein conformation.

Constant pH molecular dynamics simulations: Current status and recent applications

Many important protein functions are carried out through proton-coupled conformational dynamics. Thus, the ability to accurately model protonation states dynamically has wide-ranging implications. Over the past two decades, two main types of constant pH methods (discrete and continuous) have been developed to enable proton-coupled molecular dynamics (MD) simulations. In this short review, we discuss the current status of the development and highlight recent applications that have advanced our understanding of protein structure-function relationships. We conclude the review by outlining the remaining challenges in the method development and projecting important areas for future applications.

Computational Study of the pH-Dependent Ionic Environment around β-Lactoglobulin

Ions are involved in multiple biological processes and may exist bound to biomolecules or may be associated with their surface. Although the presence of ions in nucleic acids has traditionally gained more interest, ion-protein interactions, often with a marked dependency on pH, are beginning to gather attention. Here we present a detailed analysis on the binding and distribution of ions around β-lactoglobulin using a constant-pH MD (CpHMD) method, at a pH range 3-8, and compare it with the more traditional Poisson-Boltzmann (PB) model and the existing experimental data. Most analyses used ion concentration maps built around the protein, obtained from either the CpHMD simulations or PB calculations. The requirements of approximate charge neutrality and ionic strength equal to bulk, imposed on the MD box, imply that the absolute value of the ion excess should be half the protein charge, which is in agreement with experimental observation on other proteins ( Proc. Natl. Acad. Sci. U.S.A. 2021, 118, e2015879118) and lends support to this protocol. In addition, the protein total charge (including territorially bound ions) estimated with MD is in excellent agreement with electrophoretic measurements. Overall, the CpHMD simulations show good agreement with the nonlinear form of the PB (NLPB) model but not with its linear form, which involves a theoretical inconsistency in the calculation of the concentration maps. In several analyses, the observed pH-dependent trends for the counterions and co-ions are those generally expected, and the ion concentration maps correctly converge to the bulk ionic strength as one moves away from the protein. Despite the overall similarity, the CpHMD and NLPB approaches show some discrepancies when analyzed in more detail, which may be related to an apparent overestimation of counterion excess and underestimation of co-ion exclusion by the NLPB model, particularly at short distances from the protein.

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.

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.

FixBox: A General Algorithm to Fix Molecular Systems in Periodic Boxes

Periodic boundary conditions (PBCs) are a standard feature of molecular simulations, and their mathematical and computational aspects are well-understood and relatively straightforward. However, they can in practice be a nuisance when simulating heterogeneous systems, especially when different types of molecules change their relative positions during the simulation. Although the translation required to fix a broken molecular complex of interest can in most cases be easily inferred by visual inspection, it typically depends on the type of system, its configuration, and the box geometry, making automated procedures problematic. We present here a general algorithm, named FixBox, that can fix a molecular complex of interest from a minimal set of definitions of its assembling parts and intended arrangement in the simulation box. It uses a unified triclinic framework for the box geometric periodicity, does not require a full molecular topology, and is applicable to various types of systems and configurations, making it possible to fully and easily automate the fixing of a broken molecular complex. The performance of the algorithm is illustrated with problematic configurations of various types of simulated systems. The presented formal framework can generally be useful for algorithms that need to perform geometrical transformations on systems with PBCs.