constant-pH Simulations Research Articles

Force Field Limitations of All-Atom Continuous Constant pH Molecular Dynamics.

All-atom constant pH molecular dynamics simulations offer a powerful tool for understanding pH-mediated and proton-coupled biological processes. As the protonation equilibria of protein side chains are shifted by electrostatic interactions and desolvation energies, pKa values calculated from the constant pH simulations may be sensitive to the underlying protein force field and water model. Here we investigated the force field dependence of the all-atom particle mesh Ewald (PME) continuous constant pH (PME-CpHMD) simulations of a mini-protein BBL. The replica-exchange titration simulations based on the Amber ff19sb and ff14sb force fields with the respective water models showed significantly overestimated pKa downshifts for a buried histidine (His166) and for two glutamic acids (Glu141 and Glu161) that are involved in salt-bridge interactions. These errors (due to undersolvation of neutral histidines and overstabilization of salt bridges) are consistent with the previously reported pKa's based on the CHARMM c22/CMAP force field, albeit in larger magnitudes. The pKa calculations also demonstrated that ff19sb with OPC water is significantly more accurate than ff14sb with TIP3P water, and the salt-bridge related pKa downshifts can be partially alleviated by the atom-pair specific Lennard-Jones corrections (NBFIX). Together, these data suggest that the accuracies of the protonation equilibria of proteins from constant pH simulations can significantly benefit from improvements of force fields.

PH-Dependent Assembly and Stability of Toll-Like Receptor 3/dsRNA Signaling Complex: Insights from Constant pH Molecular Dynamics and Metadynamics Simulations.

The pH-dependent assembly of Toll-like receptors (TLRs), which triggers a threshold-like response, is a key principle in immune signaling. While crystallography has revealed the intricate structure of these assembly complexes, the mechanisms underlying their pH dependency remain unclear. Herein, constant pH simulations and metadynamics are employed to investigate the pH-dependent assembly and stability of the TLR3/dsRNA signaling complex. The findings demonstrate that system pH regulates complex assembly and stability by modulating the protonation and charge states of histidines. Histidines in TLR3 act as pH-dependent, positively charged binding sites that capture negatively charged dsRNA. Additionally, these histidines form a [H682⁺]-[E626⁻] dipole, facilitating the assembly of two TLR3 molecules into an antisymmetric dimer through dipole-dipole interactions. Surprisingly, TLR3 can shift the pKa values of key histidines from their model pKa of 6.5, increasing protonation likelihood and enhancing ligand binding. Notably, the aromatic residue Phe84, located within the dsRNA binding site [His39⁺-His60⁺-Phe84-His108⁺], alters the pKa of His60 through cation-π interactions with its protonated state. This study offers new insights into the molecular mechanisms underlying pH-dependent immune signaling via higher-order assemblies and suggests potential applications for histidine in self-assembling biomaterials.

Force Field Limitations of All-Atom Continuous Constant pH Molecular Dynamics.

All-atom constant pH molecular dynamics simulations offer a powerful tool for understanding pH-mediated and proton-coupled biological processes. As the protonation equilibria of protein sidechains are shifted by electrostatic interactions and desolvation energies, pK a values calculated from the constant pH simulations may be sensitive to the underlying protein force field and water model. Here we investigated the force field dependence of the all-atom particle mesh Ewald (PME) continuous constant pH (PME-CpHMD) simulations of a mini-protein BBL. The replica-exchange titration simulations based on the Amber ff19SB and ff14SB force fields with the respective water models showed significantly overestimated pK a downshifts for a buried histidine (His166) and for two glutamic acids (Glu141 and Glu161) that are involved in salt-bridge interactions. These errors (due to undersolvation of neutral histidines and overstabilization of salt bridges) are consistent with the previously reported pK a's based on the CHARMM c22/CMAP force field, albeit in larger magnitudes. The pK a calculations also demonstrated that ff19SB with OPC water is significantly more accurate than ff14SB with TIP3P water, and the salt-bridge related pK a downshifts can be partially alleviated by the atom-pair specific Lennard-Jones corrections (NBFIX). Together, these data suggest that the accuracies of the protonation equilibria of proteins from constant pH simulations can significantly benefit from improvements of force fields.

Effect of cationic surfactants on titration behavior of isotactic and atactic poly(methacrylic acid)

By forming strong complexes with oppositely charged polyelectrolytes, cationic surfactants significantly affect their solution properties. We investigated the effect of bound surfactant on the titration behavior of two isomer forms of poly(methacrylic acid), PMA, atactic and isotactic PMA, aPMA and iPMA, respectively, which are known for the cooperative change in chain conformation in water. The bound surfactant micelles increase the initial degree of ionization, α, of carboxyl groups on PMA, shift the conformational transition to higher α and make it narrower, and exclude the complexed chains from the solution. The titration curves were analyzed in the framework of the Henderson-Hasselbalch theory and the Gibbs free energy change of the conformational transition was calculated. The results were discussed in the context of a complex influence of surfactant on the PMA chain conformation in solution. It is proposed that the PMA-surfactant interaction changes from a so-called hydrophobic mode, which is responsible for the solubilization of the unionized and water-insoluble iPMA at low α, to a predominantly electrostatic one at high α, which leads to the precipitation of both aPMA and iPMA from solution. Constant pH simulations using simple coarse-grained models for the PMA polymers and individual surfactants were performed to rationalize the titration behavior of PMAs in the absence and presence of surfactant. By inducing an attractive interaction between PMA monomers, an excellent agreement between the experimental and numerical titration curves was obtained. Such agreement becomes poorer upon the addition of surfactant; however, the main features of the PMA titration curves are preserved, which supports the role of hydrophobic interactions in PMA-surfactant association at low pH values and that of electrostatics at higher pH.

Constant-pH Simulations with the Polarizable Atomic Multipole AMOEBA Force Field.

Accurately predicting protein behavior across diverse pH environments remains a significant challenge in biomolecular simulations. Existing constant-pH molecular dynamics (CpHMD) algorithms are limited to fixed-charge force fields, hindering their application to biomolecular systems described by permanent atomic multipoles or induced dipoles. This work overcomes these limitations by introducing the first polarizable CpHMD algorithm in the context of the Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) force field. Additionally, our implementation in the open-source Force Field X (FFX) software has the unique ability to handle titration state changes for crystalline systems including flexible support for all 230 space groups. The evaluation of constant-pH molecular dynamics (CpHMD) with the AMOEBA force field was performed on 11 crystalline peptide systems that span the titrating amino acids (Asp, Glu, His, Lys, and Cys). Titration states were correctly predicted for 15 out of the 16 amino acids present in the 11 systems, including for the coordination of Zn2+ by cysteines. The lone exception was for a HIS-ALA peptide where CpHMD predicted both neutral histidine tautomers to be equally populated, whereas the experimental model did not consider multiple conformers and diffraction data are unavailable for rerefinement. This work demonstrates the promise polarizable CpHMD simulations for pKa predictions, the study of biochemical mechanisms such as the catalytic triad of proteases, and for improved protein-ligand binding affinity accuracy in the context of pharmaceutical lead optimization.

Constant-pH Simulations of a Coarse-Grained Model of Polyfunctional Weak Charged Biopolymers

A coarse-grained model of linear polyfunctional weak charged biopolymers was implemented, formed of different proportions of acid-base groups resembling the composition of humic substances. These substances are mainly present in dissolved organic matter in natural water. The influence of electrostatic interactions computing methods, factors concerning the structure of the chain, different functional groups, and the ionic strength on polyelectrolytes were studied. Langevin dynamics with constant pH simulations were performed using the ESPResSO package and the Python-based Molecule Builder for ESPResSo (pyMBE) library. The coverage was fitted to a polyfunctional Frumkin isotherm, with a mean-field interaction between charged beads. The composition of the chain affects the charge while ionic strength affects both the charge and the radius of gyration. Additionally, the parameters intrinsic to the polyelectrolyte model were well reproduced by fitting the polyfunctional Frumkin isotherm. In contrast, the non-intrinsic parameters depended on the ionic strength. The method developed and applied to a polyfunctional polypeptide model, that resembles a humic acid, will be very useful for characterizing biopolymers with several acid-base functional groups, where their structure, the composition of the different functional groups, and the determination of the main intrinsic proton binding constants and their proportion are not exactly known.

On the Validity of Constant pH Simulations.

Constant pH (cpH) simulations are now a standard tool for investigating charge regulation in coarse-grained models of polyelectrolytes and colloidal systems. Originally developed for studying solutions with implicit ions, extending this method to systems with explicit ions or solvents presents several challenges. Ensuring proper charge neutrality within the simulation cell requires performing titration moves in sync with the insertion or deletion of ions, a crucial aspect often overlooked in the literature. Contrary to the prevailing views, cpH simulations are inherently grand-canonical, meaning that the controlled pH is that of the reservoir. The presence of the Donnan potential between the implicit reservoir and the simulation cell introduces significant differences between titration curves calculated for open and closed systems; the pH of an isolated (closed) system is different from the pH of the reservoir for the same protonation state of the polyelectrolyte. To underscore this point, in this paper, we will compare the titration curves calculated using the usual cpH algorithm with those from the exact canonical simulation algorithm. In the latter case, titration moves adhere to the correct detailed balance condition, and pH is calculated using the recently introduced surface Widom insertion algorithm. Our findings reveal a very significant difference between the titration isotherms obtained using the standard cpH algorithm and the canonical titration algorithm, emphasizing the importance of using the correct simulation approach when studying charge regulation of polyelectrolytes, proteins, and colloidal particles.

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.

A new method for reactive constant pH simulations.

We present a simulation method that allows us to calculate the titration curves for systems undergoing protonation/deprotonation reactions-such as charged colloidal suspensions with acidic/basic surface groups, polyelectrolytes, polyampholytes, and proteins. The new approach allows us to simultaneously obtain titration curves both for systems in contact with salt and acid reservoir (semi-grand canonical ensemble) and for isolated suspensions (canonical ensemble). To treat the electrostatic interactions, we present a new method based on Ewald summation-which accounts for the existence of both Bethe and Donnan potentials within the simulation cell. We show that the Donnan potential dramatically affects the pH of a suspension. Counterintuitively, we find that in suspensions with a large volume fraction of nanoparticles and low ionic strength, the number of deprotonated groups can be 100% larger in an isolated system, compared to a system connected to a reservoir by a semi-permeable membrane-both systems being at exactly the same pH.

Histidine Protonation and Conformational Switching in Diphtheria Toxin Translocation Domain.

Protonation of key histidine residues has been long implicated in the acid-mediated cellular action of the diphtheria toxin translocation (T-) domain, responsible for the delivery of the catalytic domain into the cell. Here, we use a combination of computational (constant-pH Molecular Dynamics simulations) and experimental (NMR, circular dichroism, and fluorescence spectroscopy along with the X-ray crystallography) approaches to characterize the initial stages of conformational change happening in solution in the wild-type T-domain and in the H223Q/H257Q double mutant. This replacement suppresses the acid-induced transition, resulting in the retention of a more stable protein structure in solutions at pH 5.5 and, consequently, in reduced membrane-disrupting activity. Here, for the first time, we report the pKa values of the histidine residues of the T-domain, measured by NMR-monitored pH titrations. Most peaks in the histidine side chain spectral region are titrated with pKas ranging from 6.2 to 6.8. However, the two most up-field peaks display little change down to pH 6, which is a limiting pH for this protein in solution at concentrations required for NMR. These peaks are absent in the double mutant, suggesting they belong to H223 and H257. The constant-pH simulations indicate that for the T-domain in solution, the pKa values for histidine residues range from 3.0 to 6.5, with those most difficult to protonate being H251 and H257. Taken together, our experimental and computational data demonstrate that previously suggested cooperative protonation of all six histidines in the T-domain does not occur.

A pH-dependent cluster of charges in a conserved cryptic pocket on flaviviral envelopes.

Flaviviruses are enveloped viruses which include human pathogens that are predominantly transmitted by mosquitoes and ticks. Some, such as dengue virus, exhibit the phenomenon of antibody-dependent enhancement (ADE) of disease, making vaccine-based routes of fighting infections problematic. The pH-dependent conformational change of the envelope (E) protein required for fusion between the viral and endosomal membranes is an attractive point of inhibition by antivirals as it has the potential to diminish the effects of ADE. We examined six flaviviruses by employing large-scale molecular dynamics (MD) simulations of raft systems that represent a substantial portion of the flaviviral envelope. We utilised a benzene-mapping approach that led to a discovery of shared hotspots and conserved cryptic sites. A cryptic pocket previously shown to bind a detergent molecule exhibited strain-specific characteristics. An alternative conserved cryptic site at the E protein domain interfaces showed a consistent dynamic behaviour across flaviviruses and contained a conserved cluster of ionisable residues. Constant-pH simulations revealed cluster and domain-interface disruption under low pH conditions. Based on this, we propose a cluster-dependent mechanism that addresses inconsistencies in the histidine-switch hypothesis and highlights the role of cluster protonation in orchestrating the domain dissociation pivotal for the formation of the fusogenic trimer.

Inferring the transport stoichiometry of the proton/zinc transporter YiiP from microscale thermophoresis and simulations.

The transport stoichiometry of the Zn2+/H+ antiporter YiiP remains an open question. We used experimental microscale thermophoresis (MST) zinc binding data in combination with a microscopic thermodynamic model and constant pH simulation in order to infer the likely microstates of H+ and Zn2+ binding. The pH dependence of Zn2+ binding in different binding sites was measured with MST. For these sites, we constructed microscopic, thermodynamically consistent models based on the pKas of the binding residues and zinc binding affinities. Such a model relates the macroscopically measured MST binding curves to microscopic protonation states and zinc binding in a concentration and pH-dependent manner. We inferred microscopic model parameters from the experimental data with a Monte Carlo approach. We compared the MST-inferred microscopic parameters to results from atomistic constant pH MD simulations in the absence of Zn2+. For zinc binding site B on a solvent-exposed loop between transmembrane helices 2 and 3, MST and simulation data broadly agree with each other. For the transport site A in the transmembrane domain, the CpHMD simulations indicate that the binding site His residue is protonated before any Asp. Using the CpHMD results as a starting point for the MST inference approach, we find that the microstate distribution that best fits the experimental data indicates competitive binding between either one Zn2+ ion or two or three protons bound, with the one-proton state suppressed, possibly due to high cooperativity between the Asp and His residues. The microscopic model indicates a transport stoichiometry of 1:2, or possibly 1:3 at lower pH.

Residues involved in proton-coupled transport of antimalarial drugs by PfCRT.

A major setback in the fight against malaria has been the emergence of P. falciparum parasites resistant to chloroquine (CQ), the primary antimalarial drug. Neutral CQ enters the parasite's digestive vacuole (DV) and obtains two protons preventing its passive efflux through the membrane. The trapped CQ2+ kills the parasite by disrupting the detoxification process of heme released from proteolyzed host hemoglobin. CQ-resistant parasites show mutations in a DV membrane transporter, P. falciparum CQ resistance transporter (PfCRT), that enables active, proton-coupled efflux of CQ2+. While biochemical measurements established the coupling of PfCRT's transport function to the DV membrane proton gradient, the proton-coupled transport mechanism has remained elusive. Particularly it is still unclear what the proton binding site(s) is(are). Recently, the cryo-EM structure of PfCRT from a CQ-resistant variant was determined with a bound Fab fragment occupying its binding pocket. Using us-long all-atom molecular dynamics (MD) simulations we examined the stability of the transporter in the DV-open conformation in the absence of bound Fab fragments, discovering surprising protein stability upon removal of Fab. Next, utilizing constant pH (cpHMD) simulations at 7 different pH values, and considering Galvani potential offset corrections, we determined the pKa of the four titratable residues (H97, D137, D326, and D329) lining the binding pocket. These simulations allowed unambiguous determination of the primary site responding to DV's acidic environment. Finally, with protonation states compatible with the acidic pH of the DV for the aforementioned residues, we studied binding of CQ2+ to the transporter via a series of driven MD simulations. We then repeated the cpHMD simulations and determined the proton binding site in the presence of the bound CQ2+. Our results shed light on PfCRT's transport mechanism and open new possibilities to design new effective drugs against malaria.

Canonical titration simulations.

We present a Monte Carlo approach for performing titration simulations in the canonical ensemble. The standard constant pH (cpH) simulation methods are intrinsically grand canonical, allowing us to study the protonation state of molecules only as a function of pH in the reservoir. Due to the Donnan potential between a system and an (implicit) reservoir of a semi-grand canonical simulation, the pH of the reservoir can be significantly different from that of an isolated system, for an identical protonation state. The new titration method avoids this difficulty by using the canonical reactive Monte Carlo algorithm to calculate the protonation state of macromolecules as a function of the total number of protons present inside the simulation cell. The pH of an equilibrated system is then calculated using a new surface insertion Widom algorithm, which bypasses the difficulties associated with the bulk Widom particle insertion for intermediate and high pH values. To properly treat the long range Coulomb force, we use the Ewald summation method, showing the importance of the Bethe potential for calculating the pH of canonical systems.

Galvani Offset Potential and Constant-pH Simulations of Membrane Proteins.

A central problem in computational biophysics is the treatment of titratable residues in molecular dynamics simulations of large biological macromolecular systems. Conventional simulation methods ascribe a fixed ionization state to titratable residues in accordance with their pKa and the pH of the system, assuming that an effective average model will be able to capture the predominant behavior of the system. While this assumption may be justifiable in many cases, it is certainly limited, and it is important to design alternative methodologies allowing a more realistic treatment. Constant-pH simulation methods provide powerful approaches to handle titratable residues more realistically by allowing the ionization state to vary statistically during the simulation. Extending the molecular mechanical (MM) potential energy function to a family of potential functions accounting for different ionization states, constant-pH simulations are designed to sample all accessible configurations and ionization states, properly weighted according to their Boltzmann factor. Because protonation and deprotonation events correspond to a change in the total charge, difficulties arise when the long-range Coulomb interaction is treated on the basis of an idealized infinite simulation model and periodic boundary conditions with particle-mesh Ewald lattice sums. Charging free-energy calculations performed under these conditions in aqueous solution depend on the Galvani potential of the bulk water phase. This has important implications for the equilibrium and nonequilibrium constant-pH simulation methods grounded in the relative free-energy difference corresponding to the protonated and unprotonated residues. Here, the effect of the Galvani potential is clarified, and a simple practical solution is introduced to address this issue in constant-pH simulations of the acid-sensing ion channel (ASIC).

Self-association features of NS1 proteins from different flaviviruses

Flaviviruses comprise a large group of arboviral species that are distributed in several countries of the tropics, neotropics, and some temperate zones. Since they can produce neurological pathologies or vascular damage, there has been intense research seeking better diagnosis and treatments for their infections in the last decades. The flavivirus NS1 protein is a relevant clinical target because it is involved in viral replication, immune evasion, and virulence. Being a key factor in endothelial and tissue-specific modulation, NS1 has been largely studied to understand the molecular mechanisms exploited by the virus to reprogram host cells. A central part of the viral maturation processes is the NS1 oligomerization because many stages rely on these protein-protein assemblies. In the present study, the self-associations of NS1 proteins from Zika, Dengue, and West Nile viruses are examined through constant-pH coarse-grained biophysical simulations. Free energies of interactions were estimated for different oligomeric states and pH conditions. Our results show that these proteins can form both dimers and tetramers under conditions near physiological pH even without the presence of lipids. Moreover, pH plays an important role mainly controlling the regimes where van der Waals interactions govern their association. Finally, despite the similarity at the sequence level, we found that each flavivirus has a well-characteristic protein-protein interaction profile. These specific features can provide new hints for the development of binders both for better diagnostic tools and the formulation of new therapeutic drugs.

Theoretical prediction of pH-dependent electronic spectra in aqueous solution: A combinational application of QM/MM calculations and constant-pH simulations with configuration-selection scheme

We propose a self-contained theoretical methodology to predict the pH effects on the molecular electronic spectra in aqueous solution, combining our recently developed the configuration-selection (CS) constant-pH (CpH) simulation (CS-CpH) method with quantum mechanics/molecular mechanics (QM/MM) method. By computationally sampling many protonation state transitions by CS-CpH simulations and estimating the electronic energies by QM/MM calculations, we can obtain the theoretical description of the pH dependent electronic spectra. As an example, we applied the present methodology to para-nitrophenol (pNP) in aqueous solution, and demonstrated that the method could successfully reproduce the electronic transition energies and the isosbestic point.

BLaDE: A Basic Lambda Dynamics Engine for GPU-Accelerated Molecular Dynamics Free Energy Calculations.

There is an accelerating interest in practical applications of alchemical free energy methods to problems in protein design, constant pH simulations, and especially computer-aided drug design. In the present paper, we describe a basic lambda dynamics engine (BLaDE) that enables alchemical free energy simulations, including multisite λ dynamics (MSλD) simulations, on graphical processor units (GPUs). We find that BLaDE is 5 to 8 times faster than the current GPU implementation of MSλD-based free energy calculations in CHARMM. We also demonstrate that BLaDE running standard molecular dynamics attains a performance competitive with and sometimes exceeding that of the highly optimized OpenMM GPU code. BLaDE is available as a standalone program and through an API in CHARMM.

Inducing pH control over the critical micelle concentration of zwitterionic surfactants via polyacids adsorption: Effect of chain length and structure

HypothesisThe critical concentration above which micelles form from zwitterionic surfactant solutions and their thermodynamic stability is affected by the interaction with weak Brønsted polyacid chains (An) via the formation of charged hydrogen bonds between the latter and anionic moieties. ExperimentsThe interaction between zwitterionic micelles and polyacids capable of forming hydrogen bonds, and its dependence on the environmental pH and polymer structure, has been studied with constant–pH simulations and a restricted primitive model for all electrolytes. FindingsAt low pH, the formation of polyacid/micelle complexes is witnessed independently of the polymer size or structure, so that the concentration above which micelles form is substantially decreased compared to polyacid–free cases. Upon rising pH, polymer desorption takes place within a narrow range of pH values, its location markedly depending on the size and structure of polyacids, and on the relative disposition between headgroup charged moieties. Thus, the desorption onset for long linear polyacids (A60) interacting with sulphobetainic headgroups is roughly two pH units higher than for six decameric chains (6A10) adsorbed onto micelles bearing phosphorylcholinic headgroups. This effect, together with the preferential desorption of chain ends at intermediate pH, may be exploited for drug delivery purposes or building advanced metamaterials.

The Transport Cycle of a Sodium/Proton Antiporter

Sodium/proton antiporters, also known as sodium/proton exchangers, are essential secondary active transporters that are involved in pH homeostasis and regulation of sodium concentration and cell volume. They belong to the monovalent cation/proton antiporter families (CPA1 and CPA2) which all appear to share a common fold, spanning prokaryotic to mammalian family members. They obey general principles of secondary active transporters, namely they switch between at least two conformations in accordance with the alternating access model, the major states (outward facing (OF) and inward facing (OF)) are rooted in structurally breaking an inverted repeat symmetry, and they transduce energy from a driving ion gradient to vectorial substrate transport through conformational cycles. However, these broad principles constrain but do not immediately provide insights into the molecular mechanism of sodium/proton antiport. We have been using molecular dynamics simulations (long equilibrium MD, free energy calculations, enhanced sampling for rare events, constant pH simulations) in combination with experimental techniques such as X-ray crystallography and functional measurements to investigate key steps in the transport cycle of the prokaryotic CPA2 transporter TtNapA. A competitive binding mechanism between protons and sodium emerges from an analysis of the thermodynamics of the OF and IF states and naturally explains the experimentally observed activity peak in a narrow pH range. Direct computational sampling of the conformational transitions shows that NapA can undergo fast (sub-microsecond) conformational changes in either the proton or sodium-loaded states and reveals a distinct meta-stable intermediate state. Taken together, the information from all the major steps in the complete transport cycle can then form the basis for multi-scale kinetic models of the physiological function of the antiporter.