Conformational relaxation and water penetration coupled to ionization of internal groups in proteins

Cytotoxic monastrol derivatives as adjective inhibitors of drug-resistant Eg5: a molecular dynamics perspective

The mitotic kinesin Eg5 is a motor protein involved in the formation of bipolar spindle and cell division. Eg5 is overexpressed in various cancer cells and Eg5 targeting agents are promising candidates for cancer therapy. Subsequent to the discovery of monastrol as a small-molecule Eg5 modulator, numerous inhibitors/modulators have been reported from which a few entered clinical trials. Mutagenic investigations specified declined sensitivity of Eg5 allosteric site to monastrol due to the occurrence of drug-resistant mutations in some cell cultures. Accordingly, identification of tight binders to the mutant Eg5 allosteric site is an invaluable strategy to devise more efficient Eg5 modulators. We have previously synthesized a few dihydropyrimidinethione (DHPMT)-based 5-carboxamide monastrol derivatives (1-5) with higher cytotoxicities against AGS (IC50 9.90-98.48 µM) and MCF-7 (IC50 15.20-149.13 µM) cancer cell lines than monastrol. Within a current study, a structural insight was offered into the binding mechanism of intended derivatives inside the mutant Eg5 loop5/α2/α3 allosteric pocket. Molecular docking of the DHPMT R and S-enantiomers unraveled top-scored Eg5 complexes. Molecular dynamics (MD) simulations were carried out on 5 superior complexes as (R)-2/D130V-Eg5, (R)-4/D130V-Eg5, (R)-5/D130V-Eg5, (R)-5/L214I-Eg5, (R)-5/R119L-Eg5, and the control groups monastrol/D130V-Eg5, monastrol/L214I-Eg5, monastrol/R119L-Eg5. Free energy calculations were conducted through conformational sampling of MD-driven binding trajectories. Our results provided structural details on probable interaction mechanism of the cytotoxic DHPMTs that are difficult to address experimentally. The outputs of the current study propose new monastrol derivatives as probable resistance-overwhelming Eg5 modulators.Communicated by Ramaswamy H. Sarma.

In silico target specific design of potential quinazoline-based anti-NSCLC agents

Non-small cell lung cancer (NSCLC) accounts for 85% of all lung cancers. In spite of great advances, treatment of the disease is a medical challenge. Epidermal-growth factor receptor (EGFR) has been taken as a promising cell surface target to develop anti-NSCLC therapies. The main bottleneck to attain clinical efficacy with current EGFR tyrosine kinase inhibitors (EGFR-TKIs) is the rapid spread of oncogenic mutations. Numerous efforts have been made for the synthesis of diverse EGFR-TKIs against resistance-conferring mutations. One of the best strategies to design potent agents would be to explore existing anti-NSCLC drugs at the nonclinical development stage and prioritize privileged structural patterns. Within current study, conformational stability of clinically frequent EGFR mutants (G719S, T790M, L858R and a double mutant form L858R/T790M) were validated via DynaMut and missense3D computational servers. Subsequently, structure activity relationship (SAR) and scaffold similarity inquiry were used to rationally propose a few erlotinib analogues. Intended molecules were subjected to molecular docking and top-scored binders were further analyzed through 50-ns all atom molecular dynamics (MD) simulations to infer the dynamic behavior. The aim was to offer potential binders to overwhelm clinically frequent EGFR-TK mutants. The linear interaction energy (LIE) method was applied to compute the binding free energies between EGFR and intended ligands. For this purpose, MD-based conformational sampling of ligand-enzyme complexes and ligand-water associations were used to acquire thermodynamic energy averages. Though mechanistic details are to be explored, results of the current study identify synthetically accessible quinazoline small molecules with potential affinity toward frequent EGFR-TK mutants.[Figure: see text]Communicated by Ramaswamy H. Sarma.

Gefitinib derivatives and drug-resistance: A perspective from molecular dynamics simulations

Epidermal-growth factor receptor (EGFR) is a transmembrane tyrosine kinase (TK) with a significant role in cell survival. EGFR is upregulated in various cancer cells and known as a druggable target. Gefitinib is a first-line TK inhibitor used against metastatic non-small cell lung cancer (NSCLC). Despite initial clinical response, a conserved therapeutic effect could not be achieved due to the occurrence of resistance mechanisms. Point mutations in EGFR genes are one of the major causes of rendered tumor sensitivity. To aid in the development of more efficient TKIs, chemical structures of prevailing drugs and their target binding patterns are very important. The aim of the present study was to propose synthetically-accessible gefitinib congeners with enhanced binding fitness to clinically frequent EGFR mutants. Docking simulations of intended molecules identified 1-(4-(3-chloro-4-fluorophenylamino)-7-methoxyquinazolin-6-yl)-3-(oxazolidin-2-ylmethyl) thiourea (23) as a top-binder structure inside G719S, T790 M, L858R and T790 M/L858R-EGFR active sites. Superior docked complexes were subjected to the entire 400 ns molecular dynamics (MD) simulations. Analysis of data revealed the stability of mutant enzymes upon binding to molecule 23. All mutant complexes with the exception of a T790 M/L858R-EGFR, were majorly stabilized through cooperative hydrophobic contacts. Pairwise analysis of H-bonds proved Met793 as the conserved residue with stable H-bond participations as hydrogen bond donor (Frequency 63-96%). Amino acid decomposition analysis confirmed the probable role of Met793 in complex stabilization. Estimated binding free energies indicated the proper accommodation of molecule 23 inside target active sites. Pairwise energy decompositions of stable binding modes revealed the energetic contribution of key residues. Although wet lab experiments are required to unravel the mechanistic details of mEGFR inhibition, MD results provide structural basis for those events that are difficult to address experimentally. The outputs of the current study may assist to design small molecules with high potency to mEGFRs.

Ion channel selectivity through ion-modulated changes of selectivity filter pKa values

In bacterial voltage-gated sodium channels, the passage of ions through the pore is controlled by a selectivity filter (SF) composed of four glutamate residues. The mechanism of selectivity has been the subject of intense research, with suggested mechanisms based on steric effects, and ion-triggered conformational change. Here, we propose an alternative mechanism based on ion-triggered shifts in pKa values of SF glutamates. We study the NavMs channel for which the open channel structure is available. Our free-energy calculations based on molecular dynamics simulations suggest that pKa values of the four glutamates are higher in solution of K+ ions than in solution of Na+ ions. Higher pKa in the presence of K+ stems primarily from the higher population of dunked conformations of the protonated Glu sidechain, which exhibit a higher pKa shift. Since pKa values are close to the physiological pH, this results in predominant population of the fully deprotonated state of glutamates in Na+ solution, while protonated states are predominantly populated in K+ solution. Through molecular dynamics simulations we calculate that the deprotonated state is the most conductive, the singly protonated state is less conductive, and the doubly protonated state has significantly reduced conductance. Thus, we propose that a significant component of selectivity is achieved through ion-triggered shifts in the protonation state, which favors more conductive states for Na+ ions and less conductive states for K+ ions. This mechanism also suggests a strong pH dependence of selectivity, which has been experimentally observed in structurally similar NaChBac channels.

Selenium Organic Content Prediction in Jengkol (Archidendron pauciflorum) and Its Molecular Interaction with Cardioprotection Receptors PPAR-γ, NF-κB, and PI3K

Selenium (Se) is a trace mineral found in plants with a distinct sulfuric odor that is cardioprotective and reported to have low toxicity. West Java, Indonesia, has a variety of plants with a distinct odor that are consumed raw, such as jengkol (Archidendron pauciflorum). This study is conducted to determine the Se content of jengkol using the fluorometric method, where the jengkol extract is separated, and the Se content is detected using high-pressure liquid chromatography (HPLC), combined with fluorometry. Two fractions with the highest Se concentration (A and B) are found and characterized using liquid chromatography mass spectrometry to predict the organic Se content by comparing the results with those in the external literature. The Se content of fraction (A) is found to be selenomethionine (m/z 198), gamma glutamyl-methyl-selenocysteine-(GluMetSeCys; m/z 313), and the Se-sulfur (S) conjugate of cysteine-selenoglutathione (m/z 475). Furthermore, these compounds are docked on receptors involved in cardioprotection. The receptors are peroxisome proliferator-activated receptor-γ (PPAR-γ), nuclear factor kappa-B (NF-κB), and phosphoinositide 3-kinase (PI3K/AKT). The interaction of receptor and ligan that has the lowest binding energy of the docking simulation is measured with molecular dynamic simulation. MD is performed to observe bond stability and conformation based on root mean square deviation, root mean square fluctuation, radius gyration, and MM-PBSA parameters. The results of the MD simulation show that the stability of the complex organic Se compounds tested with the receptors is lower than that of the native ligand, while the binding energy is lower than that of the native ligand based on the MM-PSBA parameter. This indicates that the predicted organic Se in jengkol, i.e., gamma-GluMetSeCys to PPAR-γ, gamma-GluMetSeCys AKT/PI3K, and Se-S conjugate of cysteine-selenoglutathione to NF-κB, has the best interaction results and provides a cardioprotection effect, compared to the molecular interaction of the test ligands with the receptors.

In silico repurposing of CNS drugs for multiple sclerosis

Multiple sclerosis (MS) is an autoimmune neurodegenerative disease affecting numerous people worldwide. While the relapsing subtypes of MS are to some extent treatable, the disease remains incurable leading to progressive disability. Limited efficacy of current small molecule drugs necessitates development of efficient and safe MS medications. Accordingly, drug repurposing is an invaluable strategy that recognizes new targets for known drugs especially in the field of poorly addressed therapeutic areas. Drug discovery largely depends on the identification of potential binding molecules to the intended biomolecular target(s). In this regard, current study was devoted to in silico repurposing of 263 small molecule CNS drugs to achieve superior binders to some MS-related targets. On the basis of molecular docking scores, thioxanthene and benzisothiazole-based antipsychotics could be identified as potential binders to sphingosine-1-phosphate lyase (S1PL) and cyclophilin D (CypD). Tightest interaction modes were observed for zuclopenthixol-S1PL (ΔG_b -7.96 kcal/mol) and lurasidone-CypD (ΔG_b -8.84 kcal/mol) complexes. Molecular dynamics (MD) simulations proved the appropriate and stable accommodation of top-ranked drugs inside enzyme binding sites during 100 ns. Hydroxyethyl piperazine of zuclopenthixol and benzisothiazole of lurasidone flipped inside the binding pocket to interact with adjacent polar and apolar residues. Solvent accessible surface area (SASA) fluctuations confirmed the results of binding trajectory analysis and showed that non-polar hydrophobic interactions played significant roles in acquired stabilities. Our results on lurasidone binding pattern were interestingly in accordance with previous reports on X-ray structures of other norbornane maleimide derivatives as CypD inhibitors. According to this, Asn144, Phe102 and Phe155 served as important residues in providing stable binding pose of lurasidone through both exo and endo conformations. Although experimental results are necessary to be achieved, the outcomes of this study proposed the potentiality of some thioxanthene and benzisothiazole-based antipsychotics for binding to S1PL and CypD, respectively, as MS-related targets.

Different pKa Shifts of Internal GLU8 in Human β-Endorphin Amyloid Revealing a Coupling of Internal Ionization and Stepwise Fibril Disassembly

As a functional amyloid, human β-endorphin amyloid fibril features a β-solenoid conformation and store peptide hormones within acidic secretory granules, which would be released into the blood through fibril disassembly when the cellular milieu pH increases from acidic to neutral level on exocytosis. To gain detailed atomic mechanism of β-endorphin amyloid fibrils' pH-responsive disassembly, we conduct constant pH molecular dynamics simulations to investigate the structural and dynamical properties of β-endorphin amyloid fibrils in experiencing the environmental pH changes. Our results demonstrate a clear pK_a shift of the internal ionizable residue of GLU8, and this shift becomes even more pronounced when it is buried more deeply in the amyloid fibrils. The unusual pK_a of GLU8 reveals that its protonation state changes from the protonated state in the acidic secretory granule to the deprotonated state in the neutral pH conditions in the blood, where the deprotonation of GLU8 leads to unfavorable interactions within the hydrophobic core of the amyloid and subsequent fibril disassembly. The different pK_a shifts of GLU8 relative to its positions in the amyloid fibril indicate that the β-endorphin amyloid fibril disassembly is a stepwise process, accounting for the experimental observation that the disassembly always initiates from the outermost layer. This study reveals the critical role of the protonation state of GLU8 in amyloid fibrils' pH-responsive disassembly, and provides clear insights for understanding the structural transitions of amyloids in hormone secretion. This study also provides theoretical basis for designing pH-sensitive biological tools for specific use with precise positioning of ionizable residues into the hydrophobic interior of proteins.

Protein pKa Prediction by Tree-Based Machine Learning

Protonation states of ionizable protein residues modulate many essential biological processes. For correct modeling and understanding of these processes, it is crucial to accurately determine their pKa values. Here, we present four tree-based machine learning models for protein pKa prediction. The four models, Random Forest, Extra Trees, eXtreme Gradient Boosting (XGBoost), and Light Gradient Boosting Machine (LightGBM), were trained on three experimental PDB and pKa datasets, two of which included a notable portion of internal residues. We observed similar performance among the four machine learning algorithms. The best model trained on the largest dataset performs 37% better than the widely used empirical pKa prediction tool PROPKA and 15% better than the published result from the pKa prediction method DelPhiPKa. The overall root-mean-square error (RMSE) for this model is 0.69, with surface and buried RMSE values being 0.56 and 0.78, respectively, considering six residue types (Asp, Glu, His, Lys, Cys, and Tyr), and 0.63 when considering Asp, Glu, His, and Lys only. We provide pKa predictions for proteins in human proteome from the AlphaFold Protein Structure Database and observed that 1% of Asp/Glu/Lys residues have highly shifted pKa values close to the physiological pH.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

Reverse Protonation of Buried Ion-Pairs in Staphylococcal Nuclease Mutants

Although buried titratable residues in protein cavities are often of major functional importance, it is generally challenging to understand their properties such as the ionization state and factors of stabilization based on experimental studies alone. A specific set of examples involve buried Glu-Lys pairs in a series of variants of Staphylococcal nuclease, for which recent structural and thermodynamic studies appeared to suggest that both the stability and the ionization state of the buried Glu-Lys pair are sensitive to its orientation (i.e., Glu23-Lys36 vs Lys23-Glu36). To further clarify the situation, especially ionization states of the buried Glu-Lys pairs, we have conducted extensive molecular dynamics simulations and free energy computations. Microsecond molecular dynamics simulations show that the hydration level of the cavity depends on the orientation of the buried ion-pair therein as well as its ionization state; free energy simulations recapitulate the relative stability of Glu23-Lys36 (EK) vs Lys23-Glu36 (KE) mutants measured experimentally, although the difference is similar in magnitude regardless of the ionization state of the Glu-Lys pair. A complementary set of free energy simulations strongly suggests that, in contrast to the original suggestion in the experimental analysis, the Glu and Lys residues prefer to adopt their charge-neutral rather than the ionized states. This result is consistent with the low dielectric constant computed for water in the cavity, which makes it difficult for the protein cavity to stabilize a pair of charged Glu-Lys residues, even with water penetration. The current study highlights the role of free energy simulations in understanding the ionization state of buried titratable residues and the relevant energetic contributions, forming the basis for the rational design of buried charge networks in proteins.

Determinants of conductance of a bacterial voltage-gated sodium channel

Through molecular dynamics (MD) and free energy simulations in electric fields, we examine the factors influencing conductance of bacterial voltage-gated sodium channel Na_vMs. The channel utilizes four glutamic acid residues in the selectivity filter (SF). Previously, we have shown, through constant pH and free energy calculations of pKa values, that fully deprotonated, singly protonated, and doubly protonated states are all feasible at physiological pH, depending on how many ions are bound in the SF. With 173 MD simulations of 450 or 500 ns and additional free energy simulations, we determine that the conductance is highest for the deprotonated state and decreases with each additional proton bound. We also determine that the pKa value of the four glutamic residues for the transition between deprotonated and singly protonated states is close to the physiological pH and that there is a small voltage dependence. The pKa value and conductance trends are in agreement with experimental work on bacterial Na_v channels, which show a decrease in maximal conductance with lowering of pH, with pKa in the physiological range. We examine binding sites for Na⁺ in the SF, compare with previous work, and note a dependence on starting structures. We find that narrowing of the gate backbone to values lower than the crystal structure's backbone radius reduces the conductance, whereas increasing the gate radius further does not affect the conductance. Simulations with some amount of negatively charged lipids as opposed to purely neutral lipids increases the conductance, as do simulations at higher voltages.

Electrostatic Environment of Proteorhodopsin Affects the pKa of Its Buried Primary Proton Acceptor

The protonation state of embedded charged residues in transmembrane proteins (TMPs) can control the onset of protein function. It is understood that interactions between an embedded charged residue and other charged or polar residues in the moiety would influence its pKa, but how the surrounding environment in which the TMP resides affects the pKa of these residues is unclear. Proteorhodopsin (PR), a light-responsive proton pump from marine bacteria, was used as a model to examine externally accessible factors that tune the pKa of its embedded charged residue, specifically its primary proton acceptor D97. The pKa of D97 was compared between PR reconstituted in liposomes with different net headgroup charges and equilibrated in buffer with different ion concentrations. For PR reconstituted in net positively charged compared to net negatively charged liposomes in low-salt buffer solutions, a drop of the apparent pKa from 7.6 to 5.6 was observed, whereas intrinsic pKa modeled with surface pH calculated from Gouy-Chapman predictions found an opposite trend for the pKa change, suggesting that surface pH does not account for the main changes observed in the apparent pKa. This difference in the pKa of D97 observed from PR reconstituted in oppositely charged liposome environments disappeared when the NaCl concentration was increased to 150 mM. We suggest that protein-intrinsic structural properties must play a role in adjusting the local microenvironment around D97 to affect its pKa, as corroborated with observations of changes in protein side-chain and hydration dynamics around the E-F loop of PR. Understanding the effect of externally controllable factors in tuning the pKa of TMP-embedded charged residues is important for bioengineering and biomedical applications relying on TMP systems, in which the onset of functions can be controlled by the protonation state of embedded residues.

Reservoir pH replica exchange

We present the reservoir pH replica exchange (R-pH-REM) method for constant pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a standard or enhanced molecular dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a constant pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from constant pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable molecules can sample two or more conformations separated by a large energy barrier, such that adequate sampling is not achieved on a time scale of standard constant pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small molecule with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calculations on this model system are in excellent agreement with a theoretical model.

Hydronium Ions Accompanying Buried Acidic Residues Lead to High Apparent Dielectric Constants in the Interior of Proteins

Internal ionizable groups are known to play important roles in protein functions. A mystery that has attracted decades of extensive experimental and theoretical studies is the apparent dielectric constants experienced by buried ionizable groups, which are much higher than values expected for protein interiors. Many interpretations have been proposed, such as water penetration, conformational relaxation, local unfolding, protein intrinsic backbone fluctuations, etc. However, these interpretations conflict with many experimental observations. The virtual mixture of multiple states (VMMS) simulation method developed in our lab provides a direct approach for studying the equilibrium of multiple chemical states and can monitor p K_a values along simulation trajectories. Through VMMS simulations of staphylococcal nuclease (SNase) variants with internal Asp or Glu residues, we discovered that cations were attracted to buried deprotonated acidic groups and the presence of the nearby cations were essential to reproduce experimentally measured p K_a values. This finding, combined with structural analysis and validation simulations, suggests that the proton released from a deprotonation process stays near the deprotonated group inside proteins, possibly in the form of a hydronium ion. The existence of a proton near a buried charge has many implications in our understanding of protein functions.

A Coupled Ionization-Conformational Equilibrium Is Required To Understand the Properties of Ionizable Residues in the Hydrophobic Interior of Staphylococcal Nuclease

Ionizable residues in the interior of proteins play essential roles, especially in biological energy transduction, but are relatively rare and seem incompatible with the complex and polar environment. We perform a comprehensive study of the internal ionizable residues on 21 variants of staphylococcal nuclease with internal Lys, Glu, or Asp residues. Using pH replica exchange molecular dynamics simulations, we find that, in most cases, the pK_a values of these internal ionizable residues are shifted significantly from their values in solution. Our calculated results are in excellent agreement with the experimental observations of the Garcia-Moreno group. We show that the interpretation of the experimental pK_a values requires the study of not only protonation changes but also conformational changes. The coupling between the protonation and conformational equilibria suggests a mechanism for efficient pH-sensing and regulation in proteins. This study provides new physical insights into how internal ionizable residues behave in the hydrophobic interior of proteins.

Anomalous Properties of Lys Residues Buried in the Hydrophobic Interior of a Protein Revealed with 15N-Detect NMR Spectroscopy

Ionizable residues buried in hydrophobic environments in proteins are essential for many fundamental biochemical processes. These residues titrate with anomalous pK_a values that are challenging to reproduce with structure-based calculations owing to the conformational reorganization coupled to their ionization. Detailed characterization of this conformational reorganization is of interest; unfortunately, the properties of buried Lys residues are difficult to study experimentally. Here we demonstrate the utility of ¹⁵N NMR spectroscopy to gain insight into the protonation state, state of hydration and conformational dynamics of the Nζ amino group of buried Lys residues. The experiments were applied to five variants of staphylococcal nuclease, with internal Lys residues that titrate with pK_a values ranging from 6.2 to 8.1. Direct detection of buried Lys residues with these NMR spectroscopy methods will enable correlation between thermodynamic and structural data as well as unprecedented examination of how conformational transitions coupled to their ionization affect their pK_a values.

Influence of glutamic acid residues and pH on the properties of transmembrane helices

Negatively charged side chains are important for the function of particular ion channels and certain other membrane proteins. To investigate the influence of single glutamic acid side chains on helices that span lipid-bilayer membranes, we have employed GWALP23 (acetyl-GGALW⁵LALALALALALALW¹⁹LAGA-amide) as a favorable host peptide framework. We substituted individual Leu residues with Glu residues (L12E or L14E or L16E) and incorporated specific ²H-labeled alanine residues within the core helical region or near the ends of the sequence. Solid-state ²H NMR spectra reveal little change for the core labels in GWALP23-E12, -E14 and -E16 over a pH range of 4 to 12.5, with the spectra being broader for samples in DOPC compared to DLPC bilayers. The spectra for samples with deuterium labels near the helix ends on alanines 3 and 21 show modest pH-dependent changes in the extent of unwinding of the helix terminals in DLPC and DOPC bilayers. The combined results indicate minor overall responses of these transmembrane helices to changes in pH, with the most buried residue E12 showing no pH dependence. While the Glu residues E14 and E16 may have high pK_a values in the lipid bilayer environment, it is also possible that a paucity of helix response is masking the pK_a values. Interestingly, when E16 is present, spectral changes at high pH report significant local unwinding of the core helix. Our results are consistent with the expectation that buried carboxyl groups aggressively hold their protons and/or waters of hydration.

Microscopic mechanisms that govern the titration response and pKa values of buried residues in staphylococcal nuclease mutants

To probe the microscopic mechanisms that govern the titration behavior of buried ionizable groups, microsecond explicit solvent molecular dynamics simulations are carried out for several mutants of Staphylococcal nuclease using both fixed charge and polarizable force fields. While the ionization of Asp 66, Glu 66, and Lys 125 lead to enhanced structural fluctuations and partial unfolding of adjacent α-helical regions, the ionization of Lys 25 causes local unfolding of adjacent β sheets. Using the sampled conformational ensembles, good agreement with experimental pK_a values is obtained with Poisson-Boltzmann calculations using a protein dielectric constant of 2-4 for V66D/E; slightly larger dielectric constants are needed for Lys mutants especially L25K, suggesting that structural responses beyond microseconds are involved in ionization of Lys 25. Overall, the set of unbiased simulations provides insights into the spatial and temporal scales of protein and solvent motions that dictate the diverse titration behaviors of buried protein residues. Proteins 2017; 85:268-281. © 2016 Wiley Periodicals, Inc.

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.

Origin of pKa Shifts of Internal Lysine Residues in SNase Studied Via Equal-Molar VMMS Simulations in Explicit Water

Protein internal ionizable groups can exhibit large shifts in pK_a values. Although the environment and interaction changes have been extensively studied both experimentally and computationally, direct calculation of pK_a values of these internal ionizable groups in explicit water is challenging due to energy barriers in solvent interaction and in conformational transition. The virtual mixture of multiple states (VMMS) method is a new approach designed to study chemical state equilibrium. This method constructs a virtual mixture of multiple chemical states in order to sample the conformational space of all states simultaneously and to avoid crossing energy barriers related to state transition. By applying VMMS to 25 variants of staphylococcal nuclease with lysine residues at internal positions, we obtained the pK_a values of these lysine residues and investigated the physics underlining the pK_a shifts. Our calculation results agree reasonably well with experimental measurements, validating the VMMS method for pK_a calculation and providing molecular details of the protonation equilibrium for protein internal ionizable groups. Based on our analyses of protein conformation relaxation, lysine side chain flexibility, water penetration, and the microenvironment, we conclude that the hydrophobicity of the microenvironment around the lysine side chain (which affects water penetration differently for different protonation states) plays an important role in the pK_a shifts.

Mechanism of pH-dependent activation of the sodium-proton antiporter NhaA

Escherichia coli NhaA is a prototype sodium-proton antiporter, which has been extensively characterized by X-ray crystallography, biochemical and biophysical experiments. However, the identities of proton carriers and details of pH-regulated mechanism remain controversial. Here we report constant pH molecular dynamics data, which reveal that NhaA activation involves a net charge switch of a pH sensor at the entrance of the cytoplasmic funnel and opening of a hydrophobic gate at the end of the funnel. The latter is triggered by charging of Asp164, the first proton carrier. The second proton carrier Lys300 forms a salt bridge with Asp163 in the inactive state, and releases a proton when a sodium ion binds Asp163. These data reconcile current models and illustrate the power of state-of-the-art molecular dynamics simulations in providing atomic details of proton-coupled transport across membrane which is challenging to elucidate by experimental techniques.

The pH dependence of the activity of Escherichia coli main sodium-proton antiporter NhaA is still not fully understood. Here, the authors use continuous constant pH molecular dynamics simulations to identify NhaA proton carrier residues and elucidate its gating and ion transport processes.

Charge-Neutral Constant pH Molecular Dynamics Simulations Using a Parsimonious Proton Buffer

In constant pH molecular dynamics simulations, the protonation states of titratable sites can respond to changes of the pH and of their electrostatic environment. Consequently, the number of protons bound to the biomolecule, and therefore the overall charge of the system, fluctuates during the simulation. To avoid artifacts associated with a non-neutral simulation system, we introduce an approach to maintain neutrality of the simulation box in constant pH molecular dynamics simulations, while maintaining an accurate description of all protonation fluctuations. Specifically, we introduce a proton buffer that, like a buffer in experiment, can exchange protons with the biomolecule enabling its charge to fluctuate. To keep the total charge of the system constant, the uptake and release of protons by the buffer are coupled to the titration of the biomolecule with a constraint. We find that, because the fluctuation of the total charge (number of protons) of a typical biomolecule is much smaller than the number of titratable sites of the biomolecule, the number of buffer sites required to maintain overall charge neutrality without compromising the charge fluctuations of the biomolecule, is typically much smaller than the number of titratable sites, implying markedly enhanced simulation and sampling efficiency.

pH-Dependent Population Shift Regulates BACE1 Activity and Inhibition

BACE1, a major therapeutic target for treatment of Alzheimer's disease, functions within a narrow pH range. Despite tremendous effort and progress in the development of BACE1 inhibitors, details of the underlying pH-dependent regulatory mechanism remain unclear. Here we elucidate the pH-dependent conformational mechanism that regulates BACE1 activity using continuous constant-pH molecular dynamics (MD). The simulations reveal that BACE1 mainly occupies three conformational states and that the relative populations of the states shift according to pH. At intermediate pH, when the catalytic dyad is monoprotonated, a binding-competent state is highly populated, while at low and high pH a Tyr-inhibited state is dominant. Our data provide strong evidence supporting conformational selection as a major mechanism for substrate and peptide-inhibitor binding. These new insights, while consistent with experiment, greatly extend the knowledge of BACE1 and have implications for further optimization of inhibitors and understanding potential side effects of targeting BACE1. Finally, the work highlights the importance of properly modeling protonation states in MD simulations.

Enhancing constant-pH simulation in explicit solvent with a two-dimensional replica exchange method

We present a new method for enhanced sampling for constant-pH simulations in explicit water based on a two-dimensional (2D) replica exchange scheme. The new method is a significant extension of our previously developed constant-pH simulation method, which is based on enveloping distribution sampling (EDS) coupled with a one-dimensional (1D) Hamiltonian exchange method (HREM). EDS constructs a hybrid Hamiltonian from multiple discrete end state Hamiltonians that, in this case, represent different protonation states of the system. The ruggedness and heights of the hybrid Hamiltonian's energy barriers can be tuned by the smoothness parameter. Within the context of the 1D EDS-HREM method, exchanges are performed between replicas with different smoothness parameters, allowing frequent protonation-state transitions and sampling of conformations that are favored by the end-state Hamiltonians. In this work, the 1D method is extended to 2D with an additional dimension, external pH. Within the context of the 2D method (2D EDS-HREM), exchanges are performed on a lattice of Hamiltonians with different pH conditions and smoothness parameters. We demonstrate that both the 1D and 2D methods exactly reproduce the thermodynamic properties of the semigrand canonical (SGC) ensemble of a system at a given pH. We have tested our new 2D method on aspartic acid, glutamic acid, lysine, a four residue peptide (sequence KAAE), and snake cardiotoxin. In all cases, the 2D method converges faster and without loss of precision; the only limitation is a loss of flexibility in how CPU time is employed. The results for snake cardiotoxin demonstrate that the 2D method enhances protonation-state transitions, samples a wider conformational space with the same amount of computational resources, and converges significantly faster overall than the original 1D method.

Coupling of remote alternating-access transport mechanisms for protons and substrates in the multidrug efflux pump AcrB

Membrane transporters of the RND superfamily confer multidrug resistance to pathogenic bacteria, and are essential for cholesterol metabolism and embryonic development in humans. We use high-resolution X-ray crystallography and computational methods to delineate the mechanism of the homotrimeric RND-type proton/drug antiporter AcrB, the active component of the major efflux system AcrAB-TolC in Escherichia coli, and one most complex and intriguing membrane transporters known to date. Analysis of wildtype AcrB and four functionally-inactive variants reveals an unprecedented mechanism that involves two remote alternating-access conformational cycles within each protomer, namely one for protons in the transmembrane region and another for drugs in the periplasmic domain, 50 Å apart. Each of these cycles entails two distinct types of collective motions of two structural repeats, coupled by flanking α-helices that project from the membrane. Moreover, we rationalize how the cross-talk among protomers across the trimerization interface might lead to a more kinetically efficient efflux system.

DOI:http://dx.doi.org/10.7554/eLife.03145.001

The interior of living cells is separated from their external environment by an enveloping membrane that serves as a protective barrier. To regulate the chemical composition of their interior, cells are equipped with specialized proteins in their membranes that move substances in and out of cells. Membrane proteins that expel molecules from the inside to the outside of the cell are called efflux pumps.

In Escherichia coli bacteria, an efflux pump known as AcrB is part of a system that removes toxic substances from the bacterial cell—such as the antibiotics used to treat bacterial infections. AcrB and other closely related efflux pumps in pathogenic bacteria are often polyspecific transporters—they can transport a large number of different toxic molecules. These efflux pump systems are also more active in bacteria that have been targeted by antibiotics, and therefore they help bacteria to evolve resistance to multiple drugs. The emergence of bacterial multi-drug resistance is a global threat to human health; to combat this phenomenon, it is essential to understand its molecular basis.

Each AcrB protein has three main parts or domains. The periplasmic domain, which is located between the two membranes that surround E. coli, works via an ‘alternating-access cycle’; that is, the shape of the periplasmic domain changes between three different forms in such a way that antibiotic molecules are first captured and subsequently squeezed through the protein towards the outside of the cell. However, the mechanism of the transmembrane domain—which is embedded in the innermost membrane of the bacterium and is the source of energy for the transport process—was not understood.

Here, Eicher et al. use X-ray crystallography to examine the three-dimensional structures of the AcrB efflux pump—and several inactive variants—in high detail. Combining these results with computer simulations reveals the mechanism used by the transmembrane domain to take up protons from the exterior and transport them into the cell. Proton transport also proceeds according to an alternating-access mechanism—and, although the transmembrane and periplasmic domains are far apart, their movements are tightly linked. Thus, because proton uptake releases energy, the transmembrane domain effectively powers the periplasmic domain to expel drugs and other molecules. Eicher et al. note that a similar mechanism has not been seen before in other efflux pumps or transporter proteins.

Understanding how AcrB works opens up new avenues that could be exploited to develop new drugs against multidrug resistant bacteria. Furthermore, Eicher et al. suggest that efflux pumps in humans closely related to AcrB might function in a similar way—including those required for regulation of cellular cholesterol, and for the correct development of embryos.

DOI:http://dx.doi.org/10.7554/eLife.03145.002

Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange

We present a new computational approach for constant pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike constant pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The analysis is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the minimum energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estimated with the EDS-HREX method agree well with the experimental pKa values. The mean absolute errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theoretical framework, which gives the correct description of multiple protonation states and good calculated pKa values.

Uncovering pH-dependent transient states of proteins with buried ionizable residues

The role of pH in regulating biological activity is ubiquitous, and understanding pH-mediated activity has traditionally relied on analyzing static biomolecular structures of highly populated ground states solved near physiological pH. However, recent advances have shown the increasing importance of transiently populated states, the characterization of which is extremely challenging but made plausible with the development of techniques such as relaxation dispersion NMR spectroscopy. To unlock the pH dependence of these transient states with atomistic-level details, we applied the recently developed explicit solvent constant pH molecular dynamics (CPHMD(MSλD)) framework to a series of staphylococcal nuclease (SNase) mutants with buried ionizable residues and probed their dynamics in different pH environments. Among our key findings is the existence of open states in all SNase mutants containing "buried" residues with highly shifted pKa's, where local solvation around the protonation site was observed. The calculated pKa demonstrated good agreement with experimental pKa's, with a low average unsigned error of 1.3 pKa units and correlation coefficient R(2) = 0.78. Sampling both open and closed states in their respective pH range, where they are expected to be dominant, was necessary to reproduce experimental pKa's, and in the most extreme examples of pKa shifts measured, it can be interpreted that the open-state structures are transient at physiological pH, contributing a small population of 1-2%. This suggests that buried ionizable residues can trigger conformational fluctuations that may be observed as transient-state structures at physiological pH. Furthermore, the coupled relationship of both open and closed states and their role in recapitulating macroscopic experimental observables suggest that structural analysis of buried residues may benefit from looking at structural pairs, as opposed to the conventional approach of looking at a single static ground-state conformation.

Identifying residues that cause pH-dependent reduction potentials

The pH dependence of the reduction potential E° for a metalloprotein indicates that the protonation state of at least one residue near the redox site changes and may be important for its activity. The responsible residue is usually identified by site-specific mutagenesis, which may be time-consuming. Here, the titration of E° for Chromatium vinosum high-potential iron-sulfur protein is predicted to be in good agreement with experiment using density functional theory and Poisson-Boltzmann calculations if only the sole histidine undergoes changes in protonation. The implementation of this approach into CHARMMing, a user-friendly web-based portal, allows users to identify residues in other proteins causing similar pH dependence.

Assessing the chemical accuracy of protein structures via peptide acidity

Although the protein native state is a Boltzmann conformational ensemble, practical applications often require a representative model from the most populated region of that distribution. The acidity of the backbone amides, as reflected in hydrogen exchange rates, is exquisitely sensitive to the surrounding charge and dielectric volume distribution. For each of four proteins, three independently determined X-ray structures of differing crystallographic resolution were used to predict exchange for the static solvent-exposed amide hydrogens. The average correlation coefficients range from 0.74 for ubiquitin to 0.93 for Pyrococcus furiosus rubredoxin, reflecting the larger range of experimental exchange rates exhibited by the latter protein. The exchange prediction errors modestly correlate with the crystallographic resolution. MODELLER 9v6-derived homology models at ~60% sequence identity (36% identity for chymotrypsin inhibitor CI2) yielded correlation coefficients that are ~0.1 smaller than for the cognate X-ray structures. The most recently deposited NOE-based ubiquitin structure and the original NMR structure of CI2 fail to provide statistically significant predictions of hydrogen exchange. However, the more recent RECOORD refinement study of CI2 yielded predictions comparable to the X-ray and homology model-based analyses.

Thermodynamic coupling of protonation and conformational equilibria in proteins: theory and simulation

Ionization-coupled conformational phenomena are ubiquitous in biology. However, quantitative characterization of the underlying thermodynamic cycle comprised of protonation and conformational equilibria has remained an elusive goal. Here we use theory and continuous constant pH molecular dynamics (CpHMD) simulations to provide a thermodynamic description for the coupling of proton titration and conformational exchange between two distinct states of a protein. CpHMD simulations with a hybrid-solvent scheme and the pH-based replica-exchange (REX) protocol are applied to obtain the equilibrium constants and atomic details of the ionization-coupled conformational exchange between open and closed states of an engineered mutant of staphylococcal nuclease. Although the coupling of protonation and conformational equilibria is not exact in the simulation, the results are encouraging. They demonstrate that REX-CpHMD simulations can be used to study thermodynamics of ionization-coupled conformational processes--which has not possible using present experimental techniques or traditional simulations based on fixed protonation states.

The pKa Cooperative: a collaborative effort to advance structure-based calculations of pKa values and electrostatic effects in proteins

The pK(a) Cooperative (http://www.pkacoop.org) was organized to advance development of accurate and useful computational methods for structure-based calculation of pK(a) values and electrostatic energies in proteins. The Cooperative brings together laboratories with expertise and interest in theoretical, computational, and experimental studies of protein electrostatics. To improve structure-based energy calculations, it is necessary to better understand the physical character and molecular determinants of electrostatic effects. Thus, the Cooperative intends to foment experimental research into fundamental aspects of proteins that depend on electrostatic interactions. It will maintain a depository for experimental data useful for critical assessment of methods for structure-based electrostatics calculations. To help guide the development of computational methods, the Cooperative will organize blind prediction exercises. As a first step, computational laboratories were invited to reproduce an unpublished set of experimental pK(a) values of acidic and basic residues introduced in the interior of staphylococcal nuclease by site-directed mutagenesis. The pK(a) values of these groups are unique and challenging to simulate owing to the large magnitude of their shifts relative to normal pK(a) values in water. Many computational methods were tested in this first Blind Prediction Challenge and critical assessment exercise. A workshop was organized in the Telluride Science Research Center to objectively assess the performance of many computational methods tested on this one extensive data set. This volume of Proteins: Structure, Function, and Bioinformatics introduces the pK(a) Cooperative, presents reports submitted by participants in the Blind Prediction Challenge, and highlights some of the problems in structure-based calculations identified during this exercise.

Arginine residues at internal positions in a protein are always charged

Many functionally essential ionizable groups are buried in the hydrophobic interior of proteins. A systematic study of Lys, Asp, and Glu residues at 25 internal positions in staphylococcal nuclease showed that their pK(a) values can be highly anomalous, some shifted by as many as 5.7 pH units relative to normal pK(a) values in water. Here we show that, in contrast, Arg residues at the same internal positions exhibit no detectable shifts in pK(a); they are all charged at pH ≤ 10. Twenty-three of these 25 variants with Arg are folded at both pH 7 and 10. The mean decrease in thermodynamic stability from substitution with Arg was 6.2 kcal/mol at this pH, comparable to that for substitution with Lys, Asp, or Glu at pH 7. The physical basis behind the remarkable ability of Arg residues to remain protonated in environments otherwise incompatible with charges is suggested by crystal structures of three variants showing how the guanidinium moiety of the Arg side chain is effectively neutralized through multiple hydrogen bonds to protein polar atoms and to site-bound water molecules. The length of the Arg side chain, and slight deformations of the protein, facilitate placement of the guanidinium moieties near polar groups or bulk water. This unique capacity of Arg side chains to retain their charge in dehydrated environments likely contributes toward the important functional roles of internal Arg residues in situations where a charge is needed in the interior of a protein, in a lipid bilayer, or in similarly hydrophobic environments.

pH replica-exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.

Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics

pH is an important parameter in condensed-phase systems, because it determines the protonation state of titratable groups and thus influences the structure, dynamics, and function of molecules in solution. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fixed and cannot adapt to changes of the local environment. Here, we present a method, implemented within the MD package GROMACS, for constant pH molecular dynamics simulations in explicit solvent that is based on the λ-dynamics approach. In the latter, the dynamics of the titration coordinate λ, which interpolates between the protonated and deprotonated states, is driven by generalized forces between the protonated and deprotonated states. The hydration free energy, as a function of pH, is included to facilitate constant pH simulations. The protonation states of titratable groups are allowed to change dynamically during a simulation, thus reproducing average protonation probabilities at a certain pH. The accuracy of the method is tested against titration curves of single amino acids and a dipeptide in explicit solvent.