Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent

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.

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.

Biomolecular dynamics in the 21st century

The relevance of motions in biological macromolecules has been clear since the early structural analyses of proteins by X-ray crystallography. Computer simulations have been applied to provide a deeper understanding of the dynamics of biological macromolecules since 1976, and are now a standard tool in many labs working on the structure and function of biomolecules. In this mini-review we highlight some areas of current interest and active development for simulations, in particular all-atom molecular dynamics simulations.

Electromagnetic Field Stimulation Therapy for Alzheimer's Disease

Alzheimer's disease (AD) is the most common neurodegenerative dementia worldwide. AD is a multifactorial disease that causes a progressive decline in memory and function precipitated by toxic beta-amyloid (Aβ) proteins, a key player in AD pathology. In 2022, 6.5 million Americans lived with AD, costing the nation $321billion. The standard of care for AD treatment includes acetylcholinesterase inhibitors (AchEIs), NMDA receptor antagonists, and monoclonal antibodies (mAbs). However, these methods are either: 1) ineffective in improving cognition, 2) unable to change disease progression, 3) limited in the number of therapeutic targets, 4) prone to cause severe side effects (brain swelling, microhemorrhages with mAb, and bradycardia and syncope with AchEIs), 5) unable to effectively cross the blood-brain barrier, and 6) lack of understanding of the aging process on the disease. mAbs are available to lower Aβ, but the difficulties of reducing the levels of the toxic Aβ proteins in the brain without triggering brain swelling or microhemorrhages associated with mAbs make the risk-benefit profile of mAbs unclear. A novel multitarget, effective, and safe non-invasive approach utilizing Repeated Electromagnetic Field Stimulation (REMFS) lowers Aβ levels in human neurons and memory areas, prevents neuronal death, stops disease progression, and improves memory without causing brain edema or bleeds in AD mice. This REMFS treatment has not been developed for humans because current EMF devices have poor penetration depth and inhomogeneous E-field distribution in the brain. Here, we discussed the biology of these effects in neurons and the design of optimal devices to treat AD.

The Electronic Structure of Genome Editors from the First Principles

Genome editing based on the CRISPR-Cas9 system has paved new avenues for medicine, pharmaceutics, biotechnology, and beyond. This article reports the role of first-principles (ab-initio) molecular dynamics (MD) in the CRISPR-Cas9 revolution, achieving a profound understanding of the enzymatic function and offering valuable insights for enzyme engineering. We introduce the methodologies and explain the use of ab-initio MD simulations to characterize the two-metal dependent mechanism of DNA cleavage in the RuvC domain of the Cas9 enzyme, and how a second catalytic domain, HNH, cleaves the target DNA with the aid of a single metal ion. A detailed description of how ab-initio MD is combined with free-energy methods - i.e., thermodynamic integration and metadynamics - to break and form chemical bonds is given, explaining the use of these methods to determine the chemical landscape and establish the catalytic mechanism in CRISPR-Cas9. The critical role of classical methods is also discussed, explaining theory and application of constant pH MD simulations, used to accurately predict the catalytic residues' protonation states. Overall, first-principles methods are shown to unravel the electronic structure of the Cas9 enzyme, providing valuable insights that can serve for the design of genome editing tools with improved catalytic efficiency or controllable activity.

GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pK_a's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pK_a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pK_a's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pK_a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pK_a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

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.

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.

Multidimensional insights into the repeated electromagnetic field stimulation and biosystems interaction in aging and age-related diseases

We provide a multidimensional sequence of events that describe the electromagnetic field (EMF) stimulation and biological system interaction. We describe this process from the quantum to the molecular, cellular, and organismal levels. We hypothesized that the sequence of events of these interactions starts with the oscillatory effect of the repeated electromagnetic stimulation (REMFS). These oscillations affect the interfacial water of an RNA causing changes at the quantum and molecular levels that release protons by quantum tunneling. Then protonation of RNA produces conformational changes that allow it to bind and activate Heat Shock Transcription Factor 1 (HSF1). Activated HSF1 binds to the DNA expressing chaperones that help regulate autophagy and degradation of abnormal proteins. This action helps to prevent and treat diseases such as Alzheimer’s and Parkinson’s disease (PD) by increasing clearance of pathologic proteins. This framework is based on multiple mathematical models, computer simulations, biophysical experiments, and cellular and animal studies. Results of the literature review and our research point towards the capacity of REMFS to manipulate various networks altered in aging (Reale et al. PloS one 9, e104973, 2014), including delay of cellular senescence (Perez et al. 2008, Exp Gerontol 43, 307-316) and reduction in levels of amyloid-β peptides (Aβ) (Perez et al. 2021, Sci Rep 11, 621). Results of these experiments using REMFS at low frequencies can be applied to the treatment of patients with age-related diseases. The use of EMF as a non-invasive therapeutic modality for Alzheimer’s disease, specifically, holds promise. It is also necessary to consider the complicated and interconnected genetic and epigenetic effects of the REMFS-biological system’s interaction while avoiding any possible adverse effects.

Multidimensional redox potential/pKa coupling in multicopper oxidases from molecular dynamics: implications for the proton transfer mechanism

Bilirubin oxidases (BOD) are metalloenzymes that catalyze the conversion of O₂ and bilirubin to biliverdin and water in the metabolism of chlorophyll and porphyrin. In this work we have used the CpHMD method to analyze the effects of the different oxidation states on the BOD trinuclear cluster (TNC). Our results demonstrate that there is a link between the different oxidation states of copper ions and the protonation capacity of nearby titratable residues. Each configuration affects pK_a differently, creating proton gradients within the enzyme that act in an extremely orderly manner. This order is closely linked to the catalytic mechanism and leads us to the conclusion of the entry of the O₂ molecule and its reduction in water molecules is associated with the probability of the release of protons from nearby acid groups. With this information, we deduce that under the initial reaction conditions the acidic side chains of nearby residues can be protonated; this allows the enzyme to reduce the activation energy of the reaction by coupling the proton transfer to oxidation state changes in the metallic center.

Free-Energy Landscape of a pH-Modulated G·C Base Pair Transition from Watson-Crick to Hoogsteen State in Duplex DNA

In double-helical DNAs, the most stable Watson-Crick (WC) base pair (bp) can be in thermal equilibrium with much less abundant Hoogsteen (HG) bp by the spontaneous rotation of the glycosidic angle in purine bases. Previous experimental studies showed that in the case of a G·C bp, the population of the transient HG is enhanced as a protonated form (HG⁺) through the protonation of the cytosine base under weakly acidic conditions. Hence, pH is a key factor that can modulate this transition event from the WC to HG⁺ bp. In this study, to computationally probe the overall free-energy landscapes of this pH-modulated G·C HG breathing, a comprehensive classical molecular dynamics (MD) simulation protocol is proposed using an enhanced sampling MD in conjunction with the standard thermodynamic integration method. From this MD protocol proposed, the free-energy surfaces of the G·C bp transition from the WC to HG bp were constructed successfully at any pH range, producing pH-dependent free-energy quantities in close agreement with previously reported experimental results. The simulation protocol is expected to provide valuable atomistic insight into the DNA bp transition events coupled with protonation or tautomeric shift in a target bp.

A Constant-pH Hybrid Monte Carlo Method with a Configuration-Selection Scheme Using the Zero Energy Difference Condition: Elucidation of Molecular Diffusivity Correlated with a pH-Dependent Solvation Shell

We have proposed a new constant-pH (CpH) hybrid Monte Carlo (MC) method with a configuration-selection (CS) scheme, called the CS-CpH method, to obtain pH-dependent physical properties within a framework of atomistic molecular simulation. The CS-CpH method consists of carrying out a short equilibrium molecular dynamics (MD) and a searching MD coupled with thermostats and barostats to generate physically plausible configurations with changed protonation states (PSs) that are subsequently accepted or rejected according to the Metropolis MC procedure. As an example, we have applied it to glutamic acid in aqueous solution and have demonstrated that it can work to generate reasonably the pH-dependent microscopic configuration ensemble compatible with the experimental pK_a value and also to show interestingly the molecular diffusivity correlated with pH-dependent solvation shell. In conclusion, we believe that the present CS-CpH method becomes a quite useful tool to study the microscopic origin of various pH-dependent phenomena, interpreting them in the atomistic chemical processes.

Biomolecular QM/MM Simulations: What Are Some of the "Burning Issues"?

QM/MM simulations have become an indispensable tool in many chemical and biochemical investigations. Considering the tremendous degree of success, including recognition by a 2013 Nobel Prize in Chemistry, are there still "burning challenges" in QM/MM methods, especially for biomolecular systems? In this short Perspective, we discuss several issues that we believe greatly impact the robustness and quantitative applicability of QM/MM simulations to many, if not all, biomolecules. We highlight these issues with observations and relevant advances from recent studies in our group and others in the field. Despite such limited scope, we hope the discussions are of general interest and will stimulate additional developments that help push the field forward in meaningful directions.

Continuous Constant pH Molecular Dynamics Simulations of Transmembrane Proteins

Many membrane channels, transporters, and receptors utilize a pH gradient or proton coupling to drive functionally relevant conformational transitions. Conventional molecular dynamics simulations employ fixed protonation states, thus neglecting the coupling between protonation and conformational equilibria. Here we describe the membrane-enabled hybrid-solvent continuous constant pH molecular dynamics method for capturing atomic details of proton-coupled conformational dynamics of transmembrane proteins. Example protocols from our recent application studies of proton channels and ion/substrate transporters are discussed.

Titration of Adenine in a GA Mismatch with Grand Canonical Simulations

Computational prediction of the pKa of ionizable groups remains a central challenge in biomolecular modeling. Although all-atom fixed-charge force fields could be accurate to describe the interaction network within the biomolecules, proper sampling techniques are required to obtain the thermodynamic information in the (de)protonation event. Sufficient sampling requires an ensemble of structures from simulations with proper treatments of the acid-base equilibria, and the grand canonical simulation technique could be used to model the growth/annihilation of hydrogen atoms by merging Hamiltonians of different protonation states into one simulation ensemble. The electrostatic feature of nucleotide systems is especially difficult to model, and the situation becomes more challenging when the ionizable site is highly perturbed. Although there are many successful predictions obtained from the grand canonical constant pH simulations, few reports focus on highly perturbed nucleotide systems with unconventional base-pair features. In this work, with the discrete constant pH method, we investigate the titration thermodynamics of an adenine in the catalytic triad in a 35-nucleotide single-stranded RNA hairpin, featuring an unconventional GA mismatch and a substantially shifted pKa value. Validation tests are performed with two system setups, both of which provide pKa predictions in good agreement with the experimental value. A single-configuration-based technique is used to calculate the pKa for comparison. The current success indicates the predictive power of the current nucleotide modeling framework.

Computing pKa Shifts Using Traditional Molecular Dynamics: Example of Acid-Base Indicator Dyes in Organized Solutions

A compound's acidity constant (K_a) in a given medium determines its protonation state and, thus, its behavior and physicochemical properties. Therefore, it is among the key characteristics considered during the design of new compounds for the needs of advanced technology, medicine, and biological research, a notable example being pH sensors. The computational prediction of K_a for weak acids and bases in homogeneous solvents is presently rather well developed. However, it is not the case for more complex media, such as microheterogeneous solutions. The constant-pH molecular dynamics (MD) method is a notable contribution to the solution of the problem, but it is not commonly used. Here, we develop an approach for predicting K_a changes of weak small-molecule acids upon transfer from water to colloid solutions by means of traditional classical molecular dynamics. The approach is based on free energy (ΔG) computations and requires limited experiment data input during calibration. It was successfully tested on a series of pH-sensitive acid-base indicator dyes in micellar solutions of surfactants. The difficulty of finite-size effects affecting ΔG computation between states with different total charges is taken into account by evaluating relevant corrections; their impact on the results is discussed, and it is found non-negligible (0.1-0.4 pK_a units). A marked bias is found in the ΔG values of acid deprotonation, as computed from MD, which is apparently caused by force-field issues. It is hypothesized to affect the constant-pH MD and reaction ensemble MD methods as well. Consequently, for these methods, a preliminary calibration is suggested.

Exploring pH Dependent Host/Guest Binding Affinities

When the electrostatic environment surrounding binding partners changes between unbound and bound states, the net uptake or release of a proton is possible by either binding partner. This process is pH-dependent in that the free energy required to uptake or release the proton varies with pH. This pH-dependence is typically not considered in conventional free energy methods where the use of fixed protonation states is the norm. In the present paper, we apply a simple two-step approach to calculate the pH-dependent binding free energy of a model cucubit[7]uril host/guest system. By use of λ-dynamics with an enhanced sampling protocol, adaptive landscape flattening, pK_a shifts and reference binding free energies upon complexation were determined. This information enables the construction of pH-dependent binding profiles that accurately capture the pK_a shifts and reproduce binding free energies at the different pH conditions that were observed experimentally. Our calculations illustrate a general framework for computing pH-dependent binding free energies but also point to some issues in modeling the molecular charge distributions within this series of molecules with CGenFF. However, by introducing some minor charge modifications to the CGenFF force field, we saw significant improvement in accuracy of the calculated pK_a shifts.

Electrostatics, Hydrogen Bonding, and Molecular Structure at Polycation and Peptide:Lipid Membrane Interfaces

Polycation and peptide-modified surfaces represent opportunities for developing potentially novel biocidal materials in a growing effort to combat bacterial resistance to traditional bactericides. It is well-known that the positive charge of these compounds is crucial to their function in biofouling prevention and as antimicrobials; however, methods for quantifying the number of positive charges on surface-bound polycations and peptides are necessary to predict, control, and optimize the design and therefore the utility of these compounds. This Spotlight on Applications reports on such an approach that combines second harmonic generation (SHG) spectroscopy, quartz crystal microbalance with dissipation monitoring (QCM-D), and atomistic simulations to obtain mechanistic insight into polycation-membrane interactions using supported lipid bilayers (SLBs) as our model system. We find that at high surface coverage, the large polycations we surveyed feature a considerably smaller percentage of ionization when compared to the smaller polycations and peptides. At these high charge densities, we suspect a pK_a shift of the charged groups to lower charge-charge repulsion as well as the formation of a looplike conformation such that less monomeric units form contact-ion pairs with the bilayer. Our sum frequency generation (SFG) spectroscopy results complement our understanding of the polycation-membrane interaction. At a high density of the polycation poly(allylamine hydrochloride) (PAH), second-order spectral line shapes are consistent with the expulsion of interfacial water molecules possibly due to contact-ion pair formation between PAH and the lipid bilayer. This finding could be essential for understanding the underlying first steps of cell lysis and penetration by polycations and should be explored further.

Probing the Accuracy of Explicit Solvent Constant pH Molecular Dynamics Simulations for Peptides

Protonation states of titratable amino acids play a key role in many biomolecular processes. Knowledge of protonatable residue charges at a given pH is essential for a correct understanding of protein catalysis, inter- and intramolecular interactions, substrate binding, and protein dynamics for instance. However, acquiring experimental values for individual amino acid protonation states of complex systems is not straightforward; therefore, several in silico approaches have been developed to tackle this issue. In this work, we assess the accuracy of our previously developed constant pH MD approach by comparing our theoretically obtained pK_a values for titratable residues with experimental values from an equivalent NMR study. We selected a set of four pentapeptides, of adequately small size to ensure comprehensive sampling, but concurrently, due to their charge composition, posing a challenge for protonation state calculation. The comparison of the pK_a values shows good agreement of the experimental and the theoretical approach with a largest difference of 0.25 pK_a units. Further, the corresponding titration curves are in fair agreement, although the shift of the Hill coefficient from a value of 1 was not always reproduced in simulations. The phase space overlap in Cartesian space between trajectories generated in constant pH and standard MD simulations is fair and suggests that our constant pH MD approach reasonably well preserves the dynamics of the system, allowing dynamic protonation MD simulations without introducing structural artifacts.

M2 amphipathic helices facilitate pH-dependent conformational transition in influenza A virus

The matrix-2 (M2) protein from influenza A virus is a tetrameric, integral transmembrane (TM) protein that plays a vital role in viral replication by proton flux into the virus. The His37 tetrad is a pH sensor in the center of the M2 TM helix that activates the channel in response to the low endosomal pH. M2 consists of different regions that are believed to be involved in membrane targeting, packaging, nucleocapsid binding, and proton transport. Although M2 has been the target of many experimental and theoretical studies that have led to significant insights into its structure and function under differing conditions, the main mechanism of proton transport, its conformational dynamics, and the role of the amphipathic helices (AHs) on proton conductance remain elusive. To this end, we have applied explicit solvent constant pH molecular dynamics using the multisite λ-dynamics approach (CpHMD^MSλD) to investigate the buried ionizable residues comprehensively and to elucidate their effect on the conformational transition. Our model recapitulates the pH-dependent conformational transition of M2 from closed to open state when the AH domain is included in the M2 construct, revealing the role of the amphipathic helices on this transition and shedding light on the proton-transport mechanism. This work demonstrates the importance of including the amphipathic helices in future experimental and theoretical studies of ion channels. Finally, our work shows that explicit solvent CpHMD^MSλD provides a realistic pH-dependent model for membrane proteins.

GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations

We present a GPU implementation of the continuous constant pH molecular dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. To test the accuracy of the tool for rapid pK_a predictions, a series of 2 ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calculated pK_a's showed a root-mean-square deviation of 0.80 and correlation coefficient of 0.83 with respect to experiment. Also, 90% of the pK_a's were converged with estimated errors below 0.1 pH units. Surprisingly, this level of accuracy is similar to our previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps^-1 (Huang, Harris, and Shen J. Chem. Inf. Model. 2018 , 58 , 1372 - 1383 ). Interestingly, for the linked titration sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the experimental macroscopic pK_a values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pK_a's is underestimated with frequent exchange attempts such as 2 ps^-1, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2 ns single-pH simulation of a 400-residue protein takes about 1 h on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, we envision that GPU-accelerated continuous CpHMD may be used in routine pK_a predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.

pH-Dependent Protonation of the Phl p 6 Pollen Allergen Studied by NMR and cpH-aMD

We use state-of-the-art NMR experiments to measure apparent pK_a values in the native protein environment and employ a cutting-edge combination of enhanced sampling and constant pH molecular dynamics (MD) simulations to rationalize strong pK_a shifts. The major timothy grass pollen allergen Phl p 6 serves as an ideal model system for both methods due to its high number of titratable residues despite its comparably small size. We present a proton transition analysis as intuitive tool to depict the captured protonation state ensemble in atomistic detail. Combining microscopic structural details from MD simulations and macroscopic ensemble averages from NMR shifts leads to a comprehensive view on pH dependencies of protonation states and tautomers. Overall, we find striking agreement between simulation-based pK_a predictions and experiment. However, our analyses suggest subtle differences in the underlying molecular origin of the observed pK_a shifts. From accelerated constant pH MD simulations, we identify immediate proximity of opposite charges, followed by vicinity of equal charges as major driving forces for pK_a shifts. NMR experiments on the other hand, suggest only a weak relation of pK_a shifts and close contacts to charged residues, while the strongest influence derives from the dipolar character of α helices. The presented study hence pinpoints opportunities for improvements concerning the theoretical description of protonation state and tautomer probabilities. However, the coherence in the resulting apparent pK_a values from simulations and experiment affirms cpH-aMD as a reliable tool to study allergen dynamics at varying pH levels.

Modeling and Simulation of Oligonucleotide Hybrids: Outlining a Strategy

Molecular dynamics simulations with a state-of-the-art force field provide an atomistic detailed description of the structural and thermodynamic features of biomolecules. Effects of chemical modifications and of the environment such as sequence, solvent, and ionic strength can explicitly be taken into account. Molecular simulation techniques can also provide insight in change in binding affinity, in protonation (pK_a shift) and tautomeric propensity due to changes in the environment or in the molecular system. The quality and reliability of a simulation depend on the quality of the force field and on the reproducibility of the data, and validation depends on the availability of suitable experimental data. Here, we describe the workflow to investigate oligonucleotide hybrids using molecular simulation including hardware and software information.

Atomistic simulations indicate the functional loop-to-coiled-coil transition in influenza hemagglutinin is not downhill

Influenza hemagglutinin (HA) mediates viral entry into host cells through a large-scale conformational rearrangement at low pH that leads to fusion of the viral and endosomal membranes. Crystallographic and biochemical data suggest that a loop-to-coiled-coil transition of the B-loop region of HA is important for driving this structural rearrangement. However, the microscopic picture for this proposed "spring-loaded" movement is missing. In this study, we focus on understanding the transition of the B loop and perform a set of all-atom molecular dynamics simulations of the full B-loop trimeric structure with the CHARMM36 force field. The free-energy profile constructed from our simulations describes a B loop that stably folds half of the postfusion coiled coil in tens of microseconds, but the full coiled coil is unfavorable. A buried hydrophilic residue, Thr59, is implicated in destabilizing the coiled coil. Interestingly, this conserved threonine is the only residue in the B loop that strictly differentiates between the group 1 and 2 HA molecules. Microsecond-scale constant temperature simulations revealed that kinetic traps in the structural switch of the B loop can be caused by nonnative, intramonomer, or intermonomer β-sheets. The addition of the A helix stabilized the postfusion state of the B loop, but introduced the possibility for further β-sheet structures. Overall, our results do not support a description of the B loop in group 2 HAs as a stiff spring, but, rather, it allows for more structural heterogeneity in the placement of the fusion peptides during the fusion process.

Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis

Solution pH plays an important role in structure and dynamics of biomolecular systems; however, pH effects cannot be accurately accounted for in conventional molecular dynamics simulations based on fixed protonation states. Continuous constant pH molecular dynamics (CpHMD) based on the λ-dynamics framework calculates protonation states on the fly during dynamical simulation at a specified pH condition. Here we report the CPU-based implementation of the CpHMD method based on the GBNeck2 generalized Born (GB) implicit-solvent model in the pmemd engine of the Amber molecular dynamics package. The performance of the method was tested using pH replica-exchange titration simulations of Asp, Glu and His side chains in 4 miniproteins and 7 enzymes with experimentally known p K_a's, some of which are significantly shifted from the model values. The added computational cost due to CpHMD titration ranges from 11 to 33% for the data set and scales roughly linearly as the ratio between the titrable sites and number of solute atoms. Comparison of the experimental and calculated p K_a's using 2 ns per replica sampling yielded a mean unsigned error of 0.70, a root-mean-squared error of 0.91, and a linear correlation coefficient of 0.79. Though this level of accuracy is similar to the GBSW-based CpHMD in CHARMM, in contrast to the latter, the current implementation was able to reproduce the experimental orders of the p K_a's of the coupled carboxylic dyads. We quantified the sampling errors, which revealed that prolonged simulation is needed to converge p K_a's of several titratable groups involved in salt-bridge-like interactions or deeply buried in the protein interior. Our benchmark data demonstrate that GBNeck2-CpHMD is an attractive tool for protein p K_a predictions.

Modeling p K Shift in DNA Triplexes Containing Locked Nucleic Acids

The protonation states for nucleic acid bases are difficult to assess experimentally. In the context of DNA triplex, the protonation state of cytidine in the third strand is particularly important, because it needs to be protonated in order to form Hoogsteen hydrogen bonds. A sugar modification, locked nucleic acid (LNA), is widely used in triplex forming oligonucleotides to target sites in the human genome. In this study, the parameters for LNA are developed in line with the CHARMM nucleic acid force field and validated toward the available structural experimental data. In conjunction, two computational methods were used to calculate the protonation state of the third strand cytidine in various DNA triplex environments: λ-dynamics and multiple pH regime. Both approaches predict p K of this cytidine shifted above physiological pH when cytidine is in the third strand in a triplex environment. Both methods show an upshift due to cytidine methylation, and a small downshift when the sugar configuration is locked. The predicted p K values for cytidine in DNA triplex environment can inform the design of better-binding oligonucleotides.

Structural Articulation of Biochemical Reactions Using Restrained Geometries and Topology Switching

A strategy named "restrained geometries and topology switching" (RGATS) is presented to obtain detailed trajectories for complex biochemical reactions using molecular mechanics (MM) methods. It enables prediction of realistic dynamical pathways for chemical reactions, especially for accurately characterizing the structural adjustments of highly complex environments to any proximal biochemical reaction. It can be used to generate reactive conformations, model stepwise or concerted reactions in complex environments, and probe the influence of changes in the environment. Its ability to take reactively nonoptimal conformations and generate favorable starting conformations for a biochemical reaction is illustrated for a proton transfer between two model compounds. Its ability to study concerted reactions in explicit solvent is illustrated using proton transfers between an ammonium ion and two conserved histidines in an ammonia transporter channel embedded in a lipid membrane. Its ability to characterize the changes induced by subtle differences in the active site environment is illustrated using nucleotide addition by a DNA polymerase in the presence of two versus three Mg²⁺ ions. RGATS can be employed within any MM program and requires no additional software implementation. This allows the full assortment of computational methods implemented in all available MM programs to be used to tackle virtually any question about biochemical reactions that is answerable without using a quantum mechanical (QM) model. It can also be applied to generate reasonable starting structures for more detailed and expensive QM or QM/MM methods. In particular, this strategy enables rapid prediction of reactant, intermediary, or product state structures in any macromolecular context, with the only requirement being that the structure in any one of these states is either known or can be accurately modeled.

Computational Study of Uracil Tautomeric Forms in the Ribosome: The Case of Uracil and 5-Oxyacetic Acid Uracil in the First Anticodon Position of tRNA

Tautomerism is important in many biomolecular interactions, not least in RNA biology. Crystallographic studies show the possible presence of minor tautomer forms of transfer-RNA (tRNA) anticodon bases in the ribosome. The hydrogen positions are not resolved in the X-ray studies, and we have used ab initio calculations and molecular dynamics simulations to understand if and how the minor enol form of uracil (U), or the modified uracil 5-oxyacetic acid (cmo⁵U), can be accommodated in the tRNA-messenger-RNA interactions in the ribosome decoding center. Ab initio calculations on isolated bases show that the modification affects the keto-enol equilibrium of the uracil base only slightly; the keto form is dominant (>99.99%) in both U and cmo⁵U. Other factors such as interactions with the surrounding nucleotides or ions would be required to shift the equilibrium toward the enol tautomer. Classical molecular simulations show a better agreement with the X-ray structures for the enol form, but free energy calculations indicate that the most stable form is the keto. In the ribosome, the enol tautomers of U and cmo⁵U pair with a guanine forming two hydrogen bonds, which do not involve the enol group. The oxyacetic acid modification has a minor effect on the keto-enol equilibrium.

Sampling the protonation states: the pH-dependent UV absorption spectrum of a polypeptide dyad

When a chromophore interacts with several titratable molecular sites, the modeling of its photophysical properties requires to take into account all their probable protonation states.

CHARMM-GUI 10 years for biomolecular modeling and simulation

CHARMM-GUI, http://www.charmm-gui.org, is a web-based graphical user interface that prepares complex biomolecular systems for molecular simulations. CHARMM-GUI creates input files for a number of programs including CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM. Since its original development in 2006, CHARMM-GUI has been widely adopted for various purposes and now contains a number of different modules designed to set up a broad range of simulations: (1) PDB Reader & Manipulator, Glycan Reader, and Ligand Reader & Modeler for reading and modifying molecules; (2) Quick MD Simulator, Membrane Builder, Nanodisc Builder, HMMM Builder, Monolayer Builder, Micelle Builder, and Hex Phase Builder for building all-atom simulation systems in various environments; (3) PACE CG Builder and Martini Maker for building coarse-grained simulation systems; (4) DEER Facilitator and MDFF/xMDFF Utilizer for experimentally guided simulations; (5) Implicit Solvent Modeler, PBEQ-Solver, and GCMC/BD Ion Simulator for implicit solvent related calculations; (6) Ligand Binder for ligand solvation and binding free energy simulations; and (7) Drude Prepper for preparation of simulations with the CHARMM Drude polarizable force field. Recently, new modules have been integrated into CHARMM-GUI, such as Glycolipid Modeler for generation of various glycolipid structures, and LPS Modeler for generation of lipopolysaccharide structures from various Gram-negative bacteria. These new features together with existing modules are expected to facilitate advanced molecular modeling and simulation thereby leading to an improved understanding of the structure and dynamics of complex biomolecular systems. Here, we briefly review these capabilities and discuss potential future directions in the CHARMM-GUI development project. © 2016 Wiley Periodicals, Inc.

Role of Pseudoisocytidine Tautomerization in Triplex-Forming Oligonucleotides: In Silico and in Vitro Studies

Pseudoisocytidine (^ΨC) is a synthetic cytidine analogue that can target DNA duplex to form parallel triplex at neutral pH. Pseudoisocytidine has mainly two tautomers, of which only one is favorable for triplex formation. In this study, we investigated the effect of sequence on ^ΨC tautomerization using λ-dynamics simulation, which takes into account transitions between states. We also performed in vitro binding experiments with sequences containing ^ΨC and furthermore characterized the structure of the formed triplex using molecular dynamics simulation. We found that the neighboring methylated or protonated cytidine promotes the formation of the favorable tautomer, whereas the neighboring thymine or locked nucleic acid has a poor effect, and consecutive ^ΨC has a negative influence. The deleterious effect of consecutive ^ΨC in a triplex formation was confirmed using in vitro binding experiments. Our findings contribute to improving the design of ^ΨC-containing triplex-forming oligonucleotides directed to target G-rich DNA sequences.

Gibbs Sampler-Based λ-Dynamics and Rao-Blackwell Estimator for Alchemical Free Energy Calculation

λ-dynamics is a generalized ensemble method for alchemical free energy calculations. In traditional λ-dynamics, the alchemical switch variable λ is treated as a continuous variable ranging from 0 to 1 and an empirical estimator is utilized to approximate the free energy. In the present article, we describe an alternative formulation of λ-dynamics that utilizes the Gibbs sampler framework, which we call Gibbs sampler-based λ-dynamics (GSLD). GSLD, like traditional λ-dynamics, can be readily extended to calculate free energy differences between multiple ligands in one simulation. We also introduce a new free energy estimator, the Rao-Blackwell estimator (RBE), for use in conjunction with GSLD. Compared with the current empirical estimator, the advantage of RBE is that RBE is an unbiased estimator and its variance is usually smaller than the current empirical estimator. We also show that the multistate Bennett acceptance ratio equation or the unbinned weighted histogram analysis method equation can be derived using the RBE. We illustrate the use and performance of this new free energy computational framework by application to a simple harmonic system as well as relevant calculations of small molecule relative free energies of solvation and binding to a protein receptor. Our findings demonstrate consistent and improved performance compared with conventional alchemical free energy methods.

Accurate Three States Model for Amino Acids with Two Chemically Coupled Titrating Sites in Explicit Solvent Atomistic Constant pH Simulations and pK(a) Calculations

Correct protonation of titratable groups in biomolecules is crucial for their accurate description by molecular dynamics simulations. In the context of constant pH simulations, an additional protonation degree of freedom is introduced for each titratable site, allowing the protonation state to change dynamically with changing structure or electrostatics. Here, we extend previous approaches for an accurate description of chemically coupled titrating sites. A second reaction coordinate is used to switch between two tautomeric states of an amino acid with chemically coupled titratable sites, such as aspartate (Asp), glutamate (Glu), and histidine (His). To this aim, we test a scheme involving three protonation states. To facilitate charge neutrality as required for periodic boundary conditions and Particle Mesh Ewald (PME) electrostatics, titration of each respective amino acid is coupled to a "water" molecule that is charged in the opposite direction. Additionally, a force field modification for Amber99sb is introduced and tested for the description of carboxyl group protonation. Our three states model is tested by titration simulations of Asp, Glu, and His, yielding a good agreement, reproducing the correct geometry of the groups in their different protonation forms. We further show that the ion concentration change due to the neutralizing "water" molecules does not significantly affect the protonation free energies of the titratable groups, suggesting that the three states model provides a good description of biomolecular dynamics at constant pH.

Probing the Structure of the Escherichia coli Periplasmic Proteins HdeA and YmgD by Molecular Dynamics Simulations

HdeA and YmgD are structurally homologous proteins in the periplasm of Escherichia coli. HdeA has been shown to represent an acid-activated chaperone, whereas the function of YmgD has not yet been characterized. We performed pH-titrating molecular dynamics simulations (pHtMD) to investigate the structural changes of both proteins and to assess whether YmgD may also exhibit an unfolding behavior similar to that of HdeA. The unfolding pathway of HdeA includes partially unfolded dimer structures, which represent a prerequisite for subsequent dissociation. In contrast to the coupled unfolding and dissociation of HdeA, YmgD displays dissociation of the folded subunits, and the subunits do not undergo significant unfolding even at low pH values. The differences in subunit stability between HdeA and YmgD may be explained by the structural features of helix D, which represents the starting point of unfolding in HdeA. In summary, the present study suggests that YmgD either is not an acid-activated chaperone or, at least, does not require unfolding for activation.

All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water

Development of a pH stat to properly control solution pH in biomolecular simulations has been a long-standing goal in the community. Toward this goal recent years have witnessed the emergence of the so-called constant pH molecular dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or truncation-based electrostatic schemes. Here we report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous constant pH molecular dynamics (CpHMD) method, enabling CpHMD to be performed with a standard MD engine at a fractional added computational cost. We demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP36, BBL, HEWL, and SNase. With the sampling time of 10 ns per replica, most pK_a's are converged, yielding the average absolute and root-mean-square deviations of 0.61 and 0.77, respectively, from experiment. Linear regression of the calculated vs experimental pK_a shifts gives a correlation coefficient of 0.79, a slope of 1, and an intercept near 0. Analysis reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromolecular systems in various environments, enabling atomic insights into pH-dependent phenomena involving not only soluble proteins but also transmembrane proteins, nucleic acids, surfactants, and polysaccharides.

Specific Binding of Cholesterol to C99 Domain of Amyloid Precursor Protein Depends Critically on Charge State of Protein

Recent NMR chemical shift measurements of the 99 residue C-terminal fragment of amyloid precursor protein (APP-C99) in the presence of cholesterol provide evidence of binary complex formation between C99 and cholesterol in membrane mimetic environments. It has also been observed that the production of Aβ protein is enhanced under conditions of high cholesterol concentration. In this study, we investigated the impact of the charge state of C99 on the structure and stability of the C99-cholesterol complex. We observed that the binding of C99 to cholesterol depends critically on the charge state of Glu 693 (E22) and Asp 694 (D23). Evaluation of the pKa values of the Asp and Glu side chains suggests that these residues may be predominantly neutral in existing experimental observations of a stable C99-cholesterol complex at lower pH (characteristic of the endosomal environment), while binding is destabilized near neutral pH (characteristic of the cytoplasm). These observations suggest that specific binding of cholesterol to C99 is a sensitive function of the pH encountered in vivo, with key E22 and D23 residues serving as a "pH switch" controlling C99-cholesterol binding.

Constant-pH Hybrid Nonequilibrium Molecular Dynamics-Monte Carlo Simulation Method

A computational method is developed to carry out explicit solvent simulations of complex molecular systems under conditions of constant pH. In constant-pH simulations, preidentified ionizable sites are allowed to spontaneously protonate and deprotonate as a function of time in response to the environment and the imposed pH. The method, based on a hybrid scheme originally proposed by H. A. Stern (J. Chem. Phys. 2007, 126, 164112), consists of carrying out short nonequilibrium molecular dynamics (neMD) switching trajectories to generate physically plausible configurations with changed protonation states that are subsequently accepted or rejected according to a Metropolis Monte Carlo (MC) criterion. To ensure microscopic detailed balance arising from such nonequilibrium switches, the atomic momenta are altered according to the symmetric two-ends momentum reversal prescription. To achieve higher efficiency, the original neMD-MC scheme is separated into two steps, reducing the need for generating a large number of unproductive and costly nonequilibrium trajectories. In the first step, the protonation state of a site is randomly attributed via a Metropolis MC process on the basis of an intrinsic pKa; an attempted nonequilibrium switch is generated only if this change in protonation state is accepted. This hybrid two-step inherent pKa neMD-MC simulation method is tested with single amino acids in solution (Asp, Glu, and His) and then applied to turkey ovomucoid third domain and hen egg-white lysozyme. Because of the simple linear increase in the computational cost relative to the number of titratable sites, the present method is naturally able to treat extremely large systems.

Diverse mixtures of 2,4-dihydroxy tautomers and O4 protonated conformers of uridine and 2'-deoxyuridine coexist in the gas phase

The gas-phase conformations of protonated uridine, [Urd+H](+), and its modified form, protonated 2'-deoxyuridine, [dUrd+H](+), generated by electrospray ionization are investigated using infrared multiple photon dissociation (IRMPD) action spectroscopy techniques. IRMPD action spectra of [Urd+H](+) and [dUrd+H](+) are measured over the IR fingerprint and hydrogen-stretching regions. [Urd+H](+) and [dUrd+H](+) exhibit very similar IRMPD spectral profiles. However, the IRMPD yields of [Urd+H](+) exceed those of [dUrd+H](+) in both the IR fingerprint and hydrogen-stretching regions. The measured spectra are compared to the linear IR spectra predicted for the stable low-energy structures of these species computed at the B3LYP/6-311+G(d,p) level of theory to determine the tautomeric conformations populated by electrospray ionization. Both B3LYP and MP2 methods find O4 and O2 protonated canonical as well as 2,4-dihydroxy tautomers among the stable low-energy structures of [Urd+H](+) and [dUrd+H](+). Comparison between the measured IRMPD and calculated linear IR spectra suggests that these species exist in their ring-closed forms and that both 2,4-dihydroxy tautomers as well as O4 protonated canonical conformers coexist in the population generated by electrospray ionization for both [Urd+H](+) and [dUrd+H](+). The 2'-deoxy modification of [dUrd+H](+) reduces the variety of 2,4-dihydroxy tautomers populated in the experiments vs. those of [Urd+H](+).

Coarse-grained modeling of RNA 3D structure

Functional RNA molecules depend on three-dimensional (3D) structures to carry out their tasks within the cell. Understanding how these molecules interact to carry out their biological roles requires a detailed knowledge of RNA 3D structure and dynamics as well as thermodynamics, which strongly governs the folding of RNA and RNA-RNA interactions as well as a host of other interactions within the cellular environment. Experimental determination of these properties is difficult, and various computational methods have been developed to model the folding of RNA 3D structures and their interactions with other molecules. However, computational methods also have their limitations, especially when the biological effects demand computation of the dynamics beyond a few hundred nanoseconds. For the researcher confronted with such challenges, a more amenable approach is to resort to coarse-grained modeling to reduce the number of data points and computational demand to a more tractable size, while sacrificing as little critical information as possible. This review presents an introduction to the topic of coarse-grained modeling of RNA 3D structures and dynamics, covering both high- and low-resolution strategies. We discuss how physics-based approaches compare with knowledge based methods that rely on databases of information. In the course of this review, we discuss important aspects in the reasoning process behind building different models and the goals and pitfalls that can result.

The Role of Electrostatic Interactions in Folding of β-Proteins

Atomic-level molecular dynamic simulations are capable of fully folding structurally diverse proteins; however, they are limited in their ability to accurately represent electrostatic interactions. Here we have experimentally tested the role of charged residues on stability and folding kinetics of one of the most widely simulated β-proteins, the WW domain. The folding of wild type Pin1 WW domain, which has two positively charged residues in the first turn, was compared to the fast folding mutant FiP35 Pin1, which introduces a negative charge into the first turn. A combination of FTIR spectroscopy and laser-induced temperature-jump coupled with infrared spectroscopy was used to probe changes in the amide I region. The relaxation dynamics of the peptide backbone, β-sheets and β-turns, and negatively charged aspartic acid side chain of FiP35 were measured independently by probing the corresponding bands assigned in the amide I region. Folding is initiated in the turns and the β-sheets form last. While the global folding mechanism is in good agreement with simulation predictions, we observe changes in the protonation state of aspartic acid during folding that have not been captured by simulation methods. The protonation state of aspartic acid is coupled to protein folding; the apparent pKa of aspartic acid in the folded protein is 6.4. The dynamics of the aspartic acid follow the dynamics of the intermediate phase, supporting assignment of this phase to formation of the first hairpin. These results demonstrate the importance of electrostatic interactions in turn stability and formation of extended β-sheet structures.

Computation of pH-dependent binding free energies

Protein-ligand binding accompanies changes in the surrounding electrostatic environments of the two binding partners and may lead to changes in protonation upon binding. In cases where the complex formation results in a net transfer of protons, the binding process is pH-dependent. However, conventional free energy computations or molecular docking protocols typically employ fixed protonation states for the titratable groups in both binding partners set a priori, which are identical for the free and bound states. In this review, we draw attention to these important yet largely ignored binding-induced protonation changes in protein-ligand association by outlining physical origins and prevalence of the protonation changes upon binding. Following a summary of various theoretical methods for pKa prediction, we discuss the theoretical framework to examine the pH dependence of protein-ligand binding processes.

Simulations of outer membrane channels and their permeability

Channels in the outer membrane of Gram-negative bacteria provide essential pathways for the controlled and unidirectional transport of ions, nutrients and metabolites into the cell. At the same time the outer membrane serves as a physical barrier for the penetration of noxious substances such as antibiotics into the bacteria. Most antibiotics have to pass through these membrane channels to either reach cytoplasmic bound targets or to further cross the hydrophobic inner membrane. Considering the pharmaceutical significance of antibiotics, understanding the functional role and mechanism of these channels is of fundamental importance in developing strategies to design new drugs with enhanced permeation abilities. Due to the biological complexity of membrane channels and experimental limitations, computer simulations have proven to be a powerful tool to investigate the structure, dynamics and interactions of membrane channels. Considerable progress has been made in computer simulations of membrane channels during the last decade. The goal of this review is to provide an overview of the computational techniques and their roles in modeling the transport across outer membrane channels. A special emphasis is put on all-atom molecular dynamics simulations employed to better understand the transport of molecules. Moreover, recent molecular simulations of ion, substrate and antibiotics translocation through membrane pores are briefly summarized. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.

Role of Deprotonation Free Energies in pKa Prediction and Molecule Ranking

A computationally efficient classical molecular simulation technique is derived for ranking the pKa values of a set of chemically similar congeneric molecules in an implicit solvent model of water. This uses the deprotonation free energy of the titratable group in the gas and aqueous phases obtained by thermodynamic integration (TI). For a series of alcohols and acids a strong linear correlation is demonstrated between the experimental pKa and the deprotonation free energy difference in the gas and liquid phases. These calculations also show that classical TI is more efficient than slow-growth TI in calculating deprotonation free energies for the series of molecules considered herein.

A Virtual Mixture Approach to the Study of Multistate Equilibrium: Application to Constant pH Simulation in Explicit Water

Chemical and thermodynamic equilibrium of multiple states is a fundamental phenomenon in biology systems and has been the focus of many experimental and computational studies. This work presents a simulation method to directly study the equilibrium of multiple states. This method constructs a virtual mixture of multiple states (VMMS) to sample the conformational space of all chemical states simultaneously. The VMMS system consists of multiple subsystems, one for each state. The subsystem contains a solute and a solvent environment. The solute molecules in all subsystems share the same conformation but have their own solvent environments. Transition between states is implicated by the change of their molar fractions. Simulation of a VMMS system allows efficient calculation of relative free energies of all states, which in turn determine their equilibrium molar fractions. For systems with a large number of state transition sites, an implicit site approximation is introduced to minimize the cost of simulation. A direct application of the VMMS method is for constant pH simulation to study protonation equilibrium. Applying the VMMS method to a heptapeptide of 3 ionizable residues, we calculated the pKas of those residues both with all explicit states and with implicit sites and obtained consistent results. For mouse epidermal growth factor of 9 ionizable groups, our VMMS simulations with implicit sites produced pKas of all 9 ionizable groups and the results agree qualitatively with NMR measurement. This example demonstrates the VMMS method can be applied to systems of a large number of ionizable groups and the computational cost scales linearly with the number of ionizable groups. For one of the most challenging systems in constant pH calculation, SNase Δ+PHS/V66K, our VMMS simulation shows that it is the state-dependent water penetration that causes the large deviation in lysine66’s pKa.

Computer simulation plays an important role to understand molecular systems and has been applied to problems of increasing complexity. Multistate equilibrium is a fundamental concept behind the structure and function of biological systems. Due to the limit in computing resources and lack of good alternative methods, computer simulation has been conducted for systems in a single state, sampling from one state to another to infer equilibrium properties. This sequential approach has been successful in many cases such as protonation equilibrium with implicit solvation model. However, state transition is difficult when explicit solvent is used for more accurate solvation description. Many efforts have been dedicated to overcome this difficulty. Analogous to real multistate systems, we proposed a virtual mixture of multiple states (VMMS) to directly simulate the equilibrium. State transitions are replaced by changes in state molar fractions. Mimicking a test tube environment, all states are simulated in parallel to equilibrate with each other. Application to constant pH simulation in explicit water demonstrates the capability of this method. It is expected that the VMMS method will find more applications in biological problems related to the equilibrium of competing states.

Conformational Dynamics and Binding Free Energies of Inhibitors of BACE-1: From the Perspective of Protonation Equilibria

BACE-1 is the β-secretase responsible for the initial amyloidogenesis in Alzheimer’s disease, catalyzing hydrolytic cleavage of substrate in a pH-sensitive manner. The catalytic mechanism of BACE-1 requires water-mediated proton transfer from aspartyl dyad to the substrate, as well as structural flexibility in the flap region. Thus, the coupling of protonation and conformational equilibria is essential to a full in silico characterization of BACE-1. In this work, we perform constant pH replica exchange molecular dynamics simulations on both apo BACE-1 and five BACE-1-inhibitor complexes to examine the effect of pH on dynamics and inhibitor binding properties of BACE-1. In our simulations, we find that solution pH controls the conformational flexibility of apo BACE-1, whereas bound inhibitors largely limit the motions of the holo enzyme at all levels of pH. The microscopic pK_a values of titratable residues in BACE-1 including its aspartyl dyad are computed and compared between apo and inhibitor-bound states. Changes in protonation between the apo and holo forms suggest a thermodynamic linkage between binding of inhibitors and protons localized at the dyad. Utilizing our recently developed computational protocol applying the binding polynomial formalism to the constant pH molecular dynamics (CpHMD) framework, we are able to obtain the pH-dependent binding free energy profiles for various BACE-1-inhibitor complexes. Our results highlight the importance of correctly addressing the binding-induced protonation changes in protein-ligand systems where binding accompanies a net proton transfer. This work comprises the first application of our CpHMD-based free energy computational method to protein-ligand complexes and illustrates the value of CpHMD as an all-purpose tool for obtaining pH-dependent dynamics and binding free energies of biological systems.

Formation of insoluble amyloid plaques in the vascular and hippocampal areas of the brain characterizes Alzheimer’s disease, a devastating neurodegenerative disorder causing dementia. Site-specific hydrolytic catalysis of β-secretase, or BACE-1, is responsible for production of oligomerative amyloid β-peptide. As the catalytic activity of BACE-1 is pH-dependent and its structural dynamics are intrinsic to the catalysis, we examine the dependence of dynamics of BACE-1 on solution pH and its implications on the catalytic mechanism of BACE-1. Also, we highlight the importance of accurate description of protonation states of the titratable groups in computer-aided drug discovery targeting BACE-1. We hope the understanding of pH dependence of the dynamics and inhibitor binding properties of BACE-1 will aid the structure-based inhibitor design efforts against Alzheimer’s disease.