Scalable Constant pH Molecular Dynamics in GROMACS

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.

Simulating the Skin Permeation Process of Ionizable Molecules

It is commonly assumed that ionizable molecules, such as drugs, permeate through the skin barrier in their neutral form. By using molecular dynamics simulations of the charged and neutral states separately, we can study the dynamic protonation behavior during the permeation process. We have studied three weak acids and three weak bases and conclude that the acids are ionized to a larger extent than the bases, when passing through the headgroup region of the lipid barrier structure, at pH values close to their pKa. It can also be observed that even if these dynamic protonation simulations are informative, in the cases studied herein they are not necessary for the calculation of permeability coefficients. It is sufficient to base the calculations only on the neutral form, as is commonly done.

Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures

Classical molecular dynamics (MD) simulations provide unmatched spatial and time resolution of protein structure and function. However, the accuracy of MD simulations often depends on the quality of force field parameters and the time scale of sampling. Another limitation of conventional MD simulations is that the protonation states of titratable amino acid residues remain fixed during simulations, even though protonation state changes coupled to conformational dynamics are central to protein function. Due to the uncertainty in selecting protonation states, classical MD simulations are sometimes performed with all amino acids modeled in their standard charged states at pH 7. Here, we performed and analyzed classical MD simulations on high-resolution cryo-EM structures of two large membrane proteins that transfer protons by catalyzing protonation/deprotonation reactions. In simulations performed with titratable amino acids modeled in their standard protonation (charged) states, the structure diverges far from its starting conformation. In comparison, MD simulations performed with predetermined protonation states of amino acid residues reproduce the structural conformation, protein hydration, and protein-water and protein-protein interactions of the structure much better. The results support the notion that it is crucial to perform basic protonation state calculations, especially on structures where protonation changes play an important functional role, prior to the launch of any conventional MD simulations. Furthermore, the combined approach of fast protonation state prediction and MD simulations can provide valuable information about the charge states of amino acids in the cryo-EM sample. Even though accurate prediction of protonation states in proteinaceous environments currently remains a challenge, we introduce an approach of combining pKa prediction with cryo-EM density map analysis that helps in improving not only the protonation state predictions but also the atomic modeling of density data.

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.

On the nature of screening in charge-regulated macroion solutions

We present a derivation of the screening length for a solution containing a charge-regulated macroion, e.g. protein, with its counterions. We show that it can be obtained directly from the second derivatives of the total free energy by taking recourse to the "uncertainty relation" of the Legendre transform, which connects the Hessians or the local curvatures of the free energy as a function of density and its Legendre transform, i.e., osmotic pressure, as a function of chemical potentials. Based on the Fowler-Guggenheim-Frumkin model of charge regulation, we then analyze the "screening resonance" and the "overscreening" of the screening properties of the charge-regulated macroion solution.

Efficient Assessment of 'Instantaneous pK' Values from Molecular Dynamics Simulations

We present a molecular simulation approach to studying the role of local and momentary molecular environment for potential acid-base reactions. For this, we combine thermodynamic considerations on the pK of ionic species with rapid sampling of energy changes related to (de)protonation. Using dispersed carbonate ions in water as a reference, our approach aims at the fast assessment of the momentary protonation energy, and thus the 'instantaneous pK', of calcium-carbonate ion aggregates. The latter include transient complexes that are elusive to long sampling runs. This motivated the elaboration of approximate, yet particularly fast assessable sampling strategies. Along this line, we were able to characterize instantaneous pK values at a statistical accuracy of 0.4 pK units within sampling runs of only 10 ps duration, whereas statistical errors reduce to 0.1 pK units in 75 ps sampling runs, respectively. This readily enabled the required time resolution for the characterization of [Cax (CO3 )y ]2(x-y) aggregates with x=1,2 and y=1,2,3, respectively. In turn, the analysis of the pH-dependent nature of calcite-water interfaces and dynamically ordered liquid-like oxyanion polymers (dollop) domains is outlined at 10 ps resolution.

Computational approaches to lipid-based nucleic acid delivery systems

Nucleic acid-based therapies have shown enormous effectiveness as vaccines against the recent COVID19 pandemics and hold great promises in the fight of a broad spectrum of diseases ranging from viral infections to cancer up to genetically transmitted pathologies. Due to their highly degradable polyanionic nature, nucleic acids need to be packed in sophisticate delivery vehicles which compact them up, protect them from early degradation and help delivery them to the right tissue/cells. Lipid-based nanoparticles (LNP) represent, at present, the main solution for nucleic acid delivery. They are made of a mixture of lipids whose key ingredient is an ionizable cationic lipid. Indeed, the interactions between the polyanionic nucleic acids and the ionizable cationic lipids, and their pH-dependent regulation in the life cycle of the nanoparticle, from production to cargo delivery, mostly determine the effectiveness of the therapeutic approach. Notwithstanding the large improvements in the delivery efficiency of LNPs in the last two decades, it is estimated that only a small fraction of the cargo is actually delivered, stimulating further research for the design of more effective LNP formulations. A rationally driven design would profit from the knowledge of the precise molecular structure of these materials, which is however still either missing or characterized by poor spatial resolution. Computational approaches have often been used as a molecular microscope either to enrich the available experimental data and provide a molecular-level picture of the LNPs or even simulate specific processes involving the formation and/or the molecular mechanisms of action of the LNP. Here, I review the recent literature in the field.

Electrostatic interactions between charge regulated spherical macroions

We study the interaction between two charge regulating spherical macroions with dielectric interior and dissociable surface groups immersed in a monovalent electrolyte solution. The charge dissociation is modelled via the Frumkin-Fowler-Guggenheim isotherm, which allows for multiple adsorption equilibrium states. The interactions are derived from the solutions of the mean-field Poisson-Boltzmann type theory with charge regulation boundary conditions. For a range of conditions we find symmetry breaking transitions from symmetric to asymmetric charge distribution exhibiting annealed charge patchiness, which results in like-charge attraction even in a univalent electrolyte-thus fundamentally modifying the nature of electrostatic interactions in charge-stabilized colloidal suspensions.

ATP-Bound State of the Uncoupling Protein 1 (UCP1) from Molecular Simulations

The uncoupling protein 1 (UCP1) dissipates the transmembrane (TM) proton gradient in the inner mitochondrial membrane (IMM) by leaking protons across the membrane and producing heat in the process. Such a nonshivering production of heat in the brown adipose tissue can combat obesity-related diseases. UCP1-associated proton leak is activated by free fatty acids and inhibited by purine nucleotides. The mechanism of proton leak and the binding sites of the activators (fatty acids) remain unknown, while the binding site of the inhibitors (nucleotides) was described recently. Using molecular dynamics simulations, we generated a conformational ensemble of UCP1. Using metadynamics-based free energy calculations, we obtained the most likely ATP-bound conformation of UCP1. Our conformational ensemble provides a molecular basis for a breadth of prior biochemical data available for UCP1. Based on the simulations, we make the following testable predictions about the mechanisms of activation of proton leak and proton leak inhibition by ATP: (1) R277 plays the dual role of stabilizing ATP at the binding site for inhibition and acting as a proton surrogate for D28 in the absence of a proton during proton transport, (2) the binding of ATP to UCP1 is mediated by residues R84, R92, R183, and S88, (3) R92 shuttles ATP from the E191-R92 gate in the intermembrane space to the nucleotide binding site and serves to increase ATP affinity, (4) ATP can inhibit proton leak by controlling the ionization states of matrix facing lysine residues such as K269 and K56, and (5) fatty acids can bind to UCP1 from the IMM either via the cavity between TM1 and TM2 or between TM5 and TM6. Our simulations set the platform for future investigations into the proton transport and inhibition mechanisms of UCP1.

Pragmatic Coarse-Graining of Proteins: Models and Applications

The molecular details involved in the folding, dynamics, organization, and interaction of proteins with other molecules are often difficult to assess by experimental techniques. Consequently, computational models play an ever-increasing role in the field. However, biological processes involving large-scale protein assemblies or long time scale dynamics are still computationally expensive to study in atomistic detail. For these applications, employing coarse-grained (CG) modeling approaches has become a key strategy. In this Review, we provide an overview of what we call pragmatic CG protein models, which are strategies combining, at least in part, a physics-based implementation and a top-down experimental approach to their parametrization. In particular, we focus on CG models in which most protein residues are represented by at least two beads, allowing these models to retain some degree of chemical specificity. A description of the main modern pragmatic protein CG models is provided, including a review of the most recent applications and an outlook on future perspectives in the field.

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.

Fighting Celiac Disease: Improvement of pH Stability of Cathepsin L In Vitro by Computational Design

Roughly 1% of the global population is susceptible to celiac disease (CD)-inheritable autoimmune inflammation of the small intestine caused by intolerance to gliadin proteins present in wheat, rye, and barley grains, and called gluten in wheat. Classical treatment is a life-long gluten-free diet, which is constraining and costly. An alternative approach is based upon the development and oral reception of effective peptidases that degrade in the stomach immunogenic proline- and glutamine-rich gliadin peptides, which are the cause of the severe reaction in the intestine. In previous research, we have established that the major digestive peptidase of an insect Tribolium castaneum-cathepsin L-hydrolyzes immunogenic prolamins after Gln residues but is unstable in the extremely acidic environment (pH 2-4) of the human stomach and cannot be used as a digestive aid. In this work, using molecular dynamics simulations, we discover the probable cause of the pH instability of cathepsin L-loss of the catalytically competent rotameric state of one of the active site residues, His 275. To "fix" the correct orientation of this residue, we designed a V277A mutant variant, which extends the range of stability of the peptidase in the acidic environment while retaining most of its activity. We suggest this protein as a lead glutenase for the development of oral medical preparation that fights CD and gluten intolerance in susceptible people.

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.

Basis for Accurate Protein pKa Prediction with Machine Learning

pH regulates protein structures and the associated functions in many biological processes via protonation and deprotonation of ionizable side chains where the titration equilibria are determined by pK_a's. To accelerate pH-dependent molecular mechanism research in the life sciences or industrial protein and drug designs, fast and accurate pK_a prediction is crucial. Here we present a theoretical pK_a data set PHMD549, which was successfully applied to four distinct machine learning methods, including DeepKa, which was proposed in our previous work. To reach a valid comparison, EXP67S was selected as the test set. Encouragingly, DeepKa was improved significantly and outperforms other state-of-the-art methods, except for the constant-pH molecular dynamics, which was utilized to create PHMD549. More importantly, DeepKa reproduced experimental pK_a orders of acidic dyads in five enzyme catalytic sites. Apart from structural proteins, DeepKa was found applicable to intrinsically disordered peptides. Further, in combination with solvent exposures, it is revealed that DeepKa offers the most accurate prediction under the challenging circumstance that hydrogen bonding or salt bridge interaction is partly compensated by desolvation for a buried side chain. Finally, our benchmark data qualify PHMD549 and EXP67S as the basis for future developments of protein pK_a prediction tools driven by artificial intelligence. In addition, DeepKa built on PHMD549 has been proven an efficient protein pK_a predictor and thus can be applied immediately to, for example, pK_a database construction, protein design, drug discovery, and so on.

Exploring the Potential Hormonal Effects of Tire Polymers (TPs) on Different Species Based on a Theoretical Computational Approach

Tire polymers (TPs) are the most prevalent type of microplastics and are of great concern due to their potential environmental risks. This study aims to determine the toxicity of TPs with the help of molecular-dynamics simulations of their interactions with receptors and to highlight the differences in the toxicity characteristics of TPs in different environmental media (marine environment, freshwater environment, soil environment). For this purpose, five TPs—natural rubber, styrene–butadiene rubber (SBR), butadiene rubber, nitrile–butadiene rubber, and isobutylene–isoprene rubber—were analyzed. Molecular-dynamics calculations were conducted on their binding energies to neurotoxic, developmental, and reproductive receptors of various organisms to characterize the toxic effects of the five TPs. The organisms included freshwater species (freshwater nematodes, snails, shrimp, and freshwater fish), marine species (marine nematodes, mussels, crab, and marine fish), and soil species (soil nematodes, springtails, earthworms, and spiders). A multilevel empowerment method was used to determine the bio-toxicity of the TPs in various environmental media. A coupled-normalization method–principal-component analysis–factor-analysis weighting method—was used to calculate the weights of the TP toxicity (first level) categories. The results revealed that the TPs were the most biologically neurotoxic to three environmental media (20.79% and 10.57% higher compared with developmental and reproductive toxicity, respectively). Regarding the effects of TPs on organisms in various environmental media (second level), using a subjective empowerment approach, a gradual increase in toxicity was observed with increasing trophic levels due to the enrichment of TPs and the feeding behavior of organisms. TPs had the greatest influence in the freshwater-environment organisms according to the subjective empowerment approach employed to weight the three environmental media (third level). Therefore, using the minimum-value method coupled with the feature-aggregation method, the interval-deflation method coupled with the entropy-weighting method, and the standard-deviation normalization method, the three toxicity characteristics of SBR in three environmental media and four organisms were determined. SBR was found to have the greatest impact on the overall toxicity of the freshwater environment (12.38% and 9.33% higher than the marine and soil environments, respectively). The greatest contribution to neurotoxicity (26.01% and 15.95% higher than developmental and reproductive toxicity, respectively) and the greatest impact on snails and shrimp among organisms in the freshwater environment were observed. The causes of the heterogeneity of SBR’s toxicity were elucidated using amino-acid-residue analysis. SBR primarily interacted with toxic receptors through van der Waals, hydrophobic, π-π, and π-sigma interactions, and the more stable the binding, the more toxic the effect. The toxicity characteristics of TMPs to various organisms in different environments identified in this paper provide a theoretical basis for subsequent studies on the prevention and control of TMPs in the environment.

A refined pH-dependent coarse-grained model for peptide structure prediction in aqueous solution

Introduction: Peptides carry out diverse biological functions and the knowledge of the conformational ensemble of polypeptides in various experimental conditions is important for biological applications. All fast dedicated softwares perform well in aqueous solution at neutral pH. Methods: In this study, we go one step beyond by combining the Debye-Hückel formalism for charged-charged amino acid interactions and a coarse-grained potential of the amino acids to treat pH and salt variations. Results: Using the PEP-FOLD framework, we show that our approach performs as well as the machine-leaning AlphaFold2 and TrRosetta methods for 15 well-structured sequences, but shows significant improvement in structure prediction of six poly-charged amino acids and two sequences that have no homologous in the Protein Data Bank, expanding the range of possibilities for the understanding of peptide biological roles and the design of candidate therapeutic peptides.

Molecular determinants of acrylamide neurotoxicity through covalent docking

Acrylamide (ACR) is formed during food processing by Maillard reaction between sugars and proteins at high temperatures. It is also used in many industries, from water waste treatment to manufacture of paper, fabrics, dyes and cosmetics. Unfortunately, cumulative exposure to acrylamide, either from diet or at the workplace, may result in neurotoxicity. Such adverse effects arise from covalent adducts formed between acrylamide and cysteine residues of several neuronal proteins via a Michael addition reaction. The molecular determinants of acrylamide reactivity and its impact on protein function are not completely understood. Here we have compiled a list of acrylamide protein targets reported so far in the literature in connection with neurotoxicity and performed a systematic covalent docking study. Our results indicate that acrylamide binding to cysteine is favored in the presence of nearby positively charged amino acids, such as lysines and arginines. For proteins with more than one reactive Cys, docking scores were able to discriminate between the primary ACR modification site and secondary sites modified only at high ACR concentrations. Therefore, docking scores emerge as a potential filter to predict Cys reactivity against acrylamide. Inspection of the ACR-protein complex structures provides insights into the putative functional consequences of ACR modification, especially for non-enzyme proteins. Based on our study, covalent docking is a promising computational tool to predict other potential protein targets mediating acrylamide neurotoxicity.

Assessing transition rates as functions of environmental variables

We present a method to estimate the transition rates of molecular systems under different environmental conditions that cause the formation or the breaking of bonds and require the sampling of the Grand Canonical Ensemble. For this purpose, we model the molecular system in terms of probable "scenarios," governed by different potential energy functions, which are separately sampled by classical MD simulations. Reweighting the canonical distribution of each scenario according to specific environmental variables, we estimate the grand canonical distribution, then use the Square Root Approximation method to discretize the Fokker-Planck operator into a rate matrix and the robust Perron Cluster Cluster Analysis method to coarse-grain the kinetic model. This permits efficiently estimating the transition rates of conformational states as functions of environmental variables, for example, the local pH at a cell membrane. In this work, we formalize the theoretical framework of the procedure, and we present a numerical experiment comparing the results with those provided by a constant-pH method based on non-equilibrium Molecular Dynamics Monte Carlo simulations. The method is relevant for the development of new drug design strategies that take into account how the cellular environment influences biochemical processes.

Best Practices in Constant pH MD Simulations: Accuracy and Sampling

Various approaches have been proposed to include the effect of pH in molecular dynamics (MD) simulations. Among these, the λ-dynamics approach proposed by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414−2423] can be performed with little computational overhead and hfor each typeence be used to routinely perform MD simulations at microsecond time scales, as shown in the accompanying paper [Aho, N. et al. J. Chem. Theory Comput.2022, DOI: 10.1021/acs.jctc.2c00516]. At such time scales, however, the accuracy of the molecular mechanics force field and the parametrization becomes critical. Here, we address these issues and provide the community with guidelines on how to set up and perform long time scale constant pH MD simulations. We found that barriers associated with the torsions of side chains in the CHARMM36m force field are too high for reaching convergence in constant pH MD simulations on microsecond time scales. To avoid the high computational cost of extending the sampling, we propose small modifications to the force field to selectively reduce the torsional barriers. We demonstrate that with such modifications we obtain converged distributions of both protonation and torsional degrees of freedom and hence consistent pK_a estimates, while the sampling of the overall configurational space accessible to proteins is unaffected as compared to normal MD simulations. We also show that the results of constant pH MD depend on the accuracy of the correction potentials. While these potentials are typically obtained by fitting a low-order polynomial to calculated free energy profiles, we find that higher order fits are essential to provide accurate and consistent results. By resolving problems in accuracy and sampling, the work described in this and the accompanying paper paves the way to the widespread application of constant pH MD beyond pK_a prediction.