Folding intermediate in the villin headpiece domain arises from disruption of a N-terminal hydrogen-bonded network

A Guide to the Continuous Constant pH Molecular Dynamics Methods in Amber and CHARMM [Article v1.0]

Like temperature and pressure, solution pH is an important environmental variable in biomolecular simulations. Virtually all proteins depend on pH to maintain their structure and function. In conventional molecular dynamics (MD) simulations of proteins, pH is implicitly accounted for by assigning and fixing protonation states of titratable sidechains. This is a significant limitation, as the assigned protonation states may be wrong and they may change during dynamics. In this tutorial, we guide the reader in learning and using the various continuous constant pH MD methods in Amber and CHARMM packages, which have been applied to predict pK a values and elucidate proton-coupled conformational dynamics of a variety of proteins including enzymes and membrane transporters.

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.

Common mechanism of thermostability in small α- and β-proteins studied by molecular dynamics

Protein thermostability is important to evolution, diseases, and industrial applications. Proteins use diverse molecular strategies to achieve stability at high temperature, yet reducing the entropy of unfolding seems required. We investigated five small α-proteins and five β-proteins with known, distinct structures and thermostability (T_m ) using multi-seed molecular dynamics simulations at 300, 350, and 400 K. The proteins displayed diverse changes in hydrogen bonding, solvent exposure, and secondary structure with no simple relationship to T_m . Our dynamics were in good agreement with experimental B-factors at 300 K and insensitive to force-field choice. Despite the very distinct structures, the native-state (300 + 350 K) free-energy landscapes (FELs) were significantly broader for the two most thermostable proteins and smallest for the three least stable proteins in both the α- and β-group and with both force fields studied independently (tailed t-test, 95% confidence level). Our results suggest that entropic ensembles stabilize proteins at high temperature due to reduced entropy of unfolding, viz., ΔG = ΔH - TΔS. Supporting this mechanism, the most thermostable proteins were also the least kinetically stable, consistent with broader FELs, typified by villin headpiece and confirmed by specific comparison to a mesophilic ortholog of Thermus thermophilus apo-pyrophosphate phosphohydrolase. We propose that molecular strategies of protein thermostabilization, although diverse, tend to converge toward highest possible entropy in the native state consistent with the functional requirements. We speculate that this tendency may explain why many proteins are not optimally structured and why molten-globule states resemble native proteins so much.

Generalized Born Implicit Solvent Models for Biomolecules

It would often be useful in computer simulations to use an implicit description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation and can be very efficient compared to the explicit treatment of the solvent. Here, we review a particular class of so-called fast implicit solvent models, generalized Born (GB) models, which are widely used for molecular dynamics (MD) simulations of proteins and nucleic acids. These approaches model hydration effects and provide solvent-dependent forces with efficiencies comparable to molecular-mechanics calculations on the solute alone; as such, they can be incorporated into MD or other conformational searching strategies in a straightforward manner. The foundations of the GB model are reviewed, followed by examples of newer, emerging models and examples of important applications. We discuss their strengths and weaknesses, both for fidelity to the underlying continuum model and for the ability to replace explicit consideration of solvent molecules in macromolecular simulations.

Atomistic Peptide Folding Simulations Reveal Interplay of Entropy and Long-Range Interactions in Folding Cooperativity

Understanding how proteins fold has remained a problem of great interest in biophysical research. Atomistic computer simulations using physics-based force fields can provide important insights on the interplay of different interactions and energetics and their roles in governing the folding thermodynamics and mechanism. In particular, generalized Born (GB)-based implicit solvent force fields can be optimized to provide an appropriate balance between solvation and intramolecular interactions and successfully recapitulate experimental conformational equilibria for a set of helical and β-hairpin peptides. Here, we further demonstrate that key thermodynamic properties and their temperature dependence obtained from replica exchange molecular dynamics simulations of these peptides are in quantitative agreement with experimental results. Useful lessons can be learned on how the interplay of entropy and sequentially long-range interactions governs the mechanism and cooperativity of folding. These results highlight the great potential of high-quality implicit solvent force fields for studying protein folding and large-scale conformational transitions.

Optimization of the GBMV2 implicit solvent force field for accurate simulation of protein conformational equilibria

Accurate treatment of solvent environment is critical for reliable simulations of protein conformational equilibria. Implicit treatment of solvation, such as using the generalized Born (GB) class of models arguably provides an optimal balance between computational efficiency and physical accuracy. Yet, GB models are frequently plagued by a tendency to generate overly compact structures. The physical origins of this drawback are relatively well understood, and the key to a balanced implicit solvent protein force field is careful optimization of physical parameters to achieve a sufficient level of cancellation of errors. The latter has been hampered by the difficulty of generating converged conformational ensembles of non-trivial model proteins using the popular replica exchange sampling technique. Here, we leverage improved sampling efficiency of a newly developed multi-scale enhanced sampling technique to re-optimize the generalized-Born with molecular volume (GBMV2) implicit solvent model with the CHARMM36 protein force field. Recursive optimization of key GBMV2 parameters (such as input radii) and protein torsion profiles (via the CMAP torsion cross terms) has led to a more balanced GBMV2 protein force field that recapitulates the structures and stabilities of both helical and β-hairpin model peptides. Importantly, this force field appears to be free of the over-compaction bias, and can generate structural ensembles of several intrinsically disordered proteins of various lengths that seem highly consistent with available experimental data. © 2017 Wiley Periodicals, Inc.

Effective Approximation of Molecular Volume Using Atom-Centered Dielectric Functions in Generalized Born Models

The generalized Born (GB) theory is a prime choice for implicit treatment of solvent that provides a favorable balance between efficiency and accuracy for reliable simulation of protein conformational equilibria. In GB, the dielectric boundary is a key physical property that needs to be properly described. While it is widely accepted that the molecular surface (MS) should provide the most physical description, most existing GB models are based on van der Waals (vdW)-like surfaces for computational simplicity and efficiency. A simple and effective approximation to molecular volume is explored here using atom-centered dielectric functions within the context of a generalized Born model with simple switching (GBSW). The new model, termed GBSW/MS2, is as efficient as the original vdW-like-surface-based GBSW model, but is able to reproduce the Born radii calculated from the "exact" Poisson-Boltzmann theory with a correlation of 0.95. More importantly, examination of the potentials of mean force of hydrogen-bonding and charge-charge interactions demonstrates that GBSW/MS2 correctly captures the first desolvation peaks, a key signature of true MS. Physical parameters including atomic input radii and peptide backbone torsion were subsequently optimized on the basis of solvation free energies of model compounds, potentials of mean force of their interactions, and conformational equilibria of a set of helical and β-hairpin model peptides. The resulting GBSW/MS2 protein force field reasonably recapitulates the structures and stabilities of these model peptides. Several remaining limitations and possible future developments are also discussed.

Modulation of the disordered conformational ensembles of the p53 transactivation domain by cancer-associated mutations

Intrinsically disordered proteins (IDPs) are frequently associated with human diseases such as cancers, and about one-fourth of disease-associated missense mutations have been mapped into predicted disordered regions. Understanding how these mutations affect the structure-function relationship of IDPs is a formidable task that requires detailed characterization of the disordered conformational ensembles. Implicit solvent coupled with enhanced sampling has been proposed to provide a balance between accuracy and efficiency necessary for systematic and comparative assessments of the effects of mutations as well as post-translational modifications on IDP structure and interaction. Here, we utilize a recently developed replica exchange with guided annealing enhanced sampling technique to calculate well-converged atomistic conformational ensembles of the intrinsically disordered transactivation domain (TAD) of tumor suppressor p53 and several cancer-associated mutants in implicit solvent. The simulations are critically assessed by quantitative comparisons with several types of experimental data that provide structural information on both secondary and tertiary levels. The results show that the calculated ensembles reproduce local structural features of wild-type p53-TAD and the effects of K24N mutation quantitatively. On the tertiary level, the simulated ensembles are overly compact, even though they appear to recapitulate the overall features of transient long-range contacts qualitatively. A key finding is that, while p53-TAD and its cancer mutants sample a similar set of conformational states, cancer mutants could introduce both local and long-range structural modulations to potentially perturb the balance of p53 binding to various regulatory proteins and further alter how this balance is regulated by multisite phosphorylation of p53-TAD. The current study clearly demonstrates the promise of atomistic simulations for detailed characterization of IDP conformations, and at the same time reveals important limitations in the current implicit solvent protein force field that must be sufficiently addressed for reliable description of long-range structural features of the disordered ensembles.

Tumor suppressor p53 is the most frequently mutated protein in human cancers. Clinical studies have suggested that the type of p53 mutation can be linked to cancer prognosis, response to drug treatment, and patient survival. It is thus crucial to understand the molecular basis of p53 inactivation by various types of mutations, so as to understand the biological outcomes and assess potential cancer intervention strategies. Here, we utilize a recently developed replica exchange with guided annealing enhanced sampling technique to calculate well-converged atomistic conformational ensembles of the intrinsically disordered transactivation domain (TAD) of tumor suppressor p53 and several cancer-associated mutants in an implicit solvent protein force field. The calculated ensembles are in quantitative agreement with several types of existing NMR data on the wild-type protein and the K24N mutant. The results suggest that, while all sequences sample a similar set of conformational substates, cancer mutants could introduce both local and long-range structural modulations and in turn perturb the balance of p53 binding to various regulatory proteins and further alter how this balance is regulated by multisite phosphorylation of p53-TAD. The study also reveals important limitations in implicit solvent for simulations of disordered proteins like p53-TAD.

Quantitative assessments of the distinct contributions of polypeptide backbone amides versus side chain groups to chain expansion via chemical denaturation

In aqueous solutions with high concentrations of chemical denaturants such as urea and guanidinium chloride (GdmCl) proteins expand to populate heterogeneous conformational ensembles. These denaturing environments are thought to be good solvents for generic protein sequences because properties of conformational distributions align with those of canonical random coils. Previous studies showed that water is a poor solvent for polypeptide backbones, and therefore, backbones form collapsed globular structures in aqueous solvents. Here, we ask if polypeptide backbones can intrinsically undergo the requisite chain expansion in aqueous solutions with high concentrations of urea and GdmCl. We answer this question using a combination of molecular dynamics simulations and fluorescence correlation spectroscopy. We find that the degree of backbone expansion is minimal in aqueous solutions with high concentrations of denaturants. Instead, polypeptide backbones sample conformations that are denaturant-specific mixtures of coils and globules, with a persistent preference for globules. Therefore, typical denaturing environments cannot be classified as good solvents for polypeptide backbones. How then do generic protein sequences expand in denaturing environments? To answer this question, we investigated the effects of side chains using simulations of two archetypal sequences with amino acid compositions that are mixtures of charged, hydrophobic, and polar groups. We find that side chains lower the effective concentration of backbone amides in water leading to an intrinsic expansion of polypeptide backbones in the absence of denaturants. Additional dilution of the effective concentration of backbone amides is achieved through preferential interactions with denaturants. These effects lead to conformational statistics in denaturing environments that are congruent with those of canonical random coils. Our results highlight the role of side chain-mediated interactions as determinants of the conformational properties of unfolded states in water and in influencing chain expansion upon denaturation.

Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism

pH is a ubiquitous regulator of biological activity, including protein-folding, protein-protein interactions, and enzymatic activity. Existing constant pH molecular dynamics (CPHMD) models that were developed to address questions related to the pH-dependent properties of proteins are largely based on implicit solvent models. However, implicit solvent models are known to underestimate the desolvation energy of buried charged residues, increasing the error associated with predictions that involve internal ionizable residue that are important in processes like hydrogen transport and electron transfer. Furthermore, discrete water and ions cannot be modeled in implicit solvent, which are important in systems like membrane proteins and ion channels. We report on an explicit solvent constant pH molecular dynamics framework based on multi-site λ-dynamics (CPHMD(MSλD)). In the CPHMD(MSλD) framework, we performed seamless alchemical transitions between protonation and tautomeric states using multi-site λ-dynamics, and designed novel biasing potentials to ensure that the physical end-states are predominantly sampled. We show that explicit solvent CPHMD(MSλD) simulations model realistic pH-dependent properties of proteins such as the Hen-Egg White Lysozyme (HEWL), binding domain of 2-oxoglutarate dehydrogenase (BBL) and N-terminal domain of ribosomal protein L9 (NTL9), and the pKa predictions are in excellent agreement with experimental values, with a RMSE ranging from 0.72 to 0.84 pKa units. With the recent development of the explicit solvent CPHMD(MSλD) framework for nucleic acids, accurate modeling of pH-dependent properties of both major class of biomolecules-proteins and nucleic acids is now possible.

Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH

Recent development of constant pH molecular dynamics (CpHMD) methods has offered promise for adding pH-stat in molecular dynamics simulations. However, until now the working pH molecular dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are critical for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titration coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titration with simultaneous ionization or neutralization of a co-ion in solution. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliphatic dicarboxylic acids with varying carbon chain length. The average absolute deviation from the experimental pK(a) values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titration coordinates, while maintaining charge neutrality of the system improves the accuracy in the calculated electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom molecular dynamics in the near future.

Towards Accurate Prediction of Protonation Equilibrium of Nucleic Acids

The role of protonated nucleotides in modulating the pH-dependent properties of nucleic acids is one of the emerging frontiers in the field of nucleic acid biology. The recent development of a constant pH molecular dynamics simulation (CPHMD^MSλD) framework for simulating nucleic acids has provided a tool for realistic simulations of pH-dependent dynamics. We enhanced the CPHMD^MSλD framework with pH-based replica exchange (pH-REX), which significantly improves the sampling of both titration and spatial coordinates. The results from our pK_a calculations for the GAAA tetraloop, which was predicted with lower accuracy previously due to sampling challenges, demonstrates that pH-REX reduces the average unsigned error (AUE) to 0.7 pK_a units, and the error of the most poorly predicted residue A17 was drastically reduced from 2.9 to 1.2 pK_a unit. Lastly, we show that pH-REX CPHMD^MSλD simulations can be used to identify the dominant conformation of nucleic acid structures in alternate pH environments. This work suggests that pH-REX CPHMD^MSλD simulations provide a practical tool for predicting nucleic acid protonation equilibrium from first-principles, and offering structural and mechanistic insight into the study of pH-dependent properties of nucleic acids.

pH-dependent dynamics of complex RNA macromolecules

The role of pH-dependent protonation equilibrium in modulating RNA dynamics and function is one of the key unanswered questions in RNA biology. Molecular dynamics (MD) simulations can provide insight into the mechanistic roles of protonated nucleotides, but it is only capable of modeling fixed protonation states and requires prior knowledge of the key residue's protonation state. Recently, we developed a framework for constant pH molecular dynamics simulations (CPHMDMSλD) of nucleic acids, where the nucleotides' protonation states are modeled as dynamic variables that are coupled to the structural dynamics of the RNA. In the present study, we demonstrate the application of CPHMDMSλD to the lead-dependent ribozyme; establishing the validity of this approach for modeling complex RNA structures. We show that CPHMDMSλD accurately predicts the direction of the pKa shifts and reproduces experimentally-measured microscopic pKa values with an average unsigned error of 1.3 pKa units. The effects of coupled titration states in RNA structures are modeled, and the importance of conformation sampling is highlighted. The general accuracy of CPHMDMSλD simulations in reproducing pH-dependent observables reported in this work demonstrates that constant pH simulations provides a powerful tool to investigate pH-dependent processes in nucleic acids.

pH-replica exchange molecular dynamics in proteins using a discrete protonation method

Protonation equilibria in biological molecules modulates structure, dynamics, and function. A pH-replica exchange molecular dynamics (pH-REMD) method is described here to improve the coupling between conformational and protonation sampling. Under a Hamiltonian replica exchange setup, conformations are swapped between two neighboring replicas, which themselves are at different pHs. The method has been validated on a series of biological systems. We applied pH-REMD to a series of model compounds, to an terminally charged ADFDA pentapeptide, and to a heptapeptide derived from the ovomucoid third domain (OMTKY3). In all of those systems, the predicted pK(a) by pH-REMD is very close to the experimental value and almost identical to the ones obtained by constant pH molecular dynamics (CpH MD). The method presented here, pH-REMD, has the advantage of faster convergence properties due to enhanced sampling of both conformation and protonation spaces.

Residual structures, conformational fluctuations, and electrostatic interactions in the synergistic folding of two intrinsically disordered proteins

To understand the interplay of residual structures and conformational fluctuations in the interaction of intrinsically disordered proteins (IDPs), we first combined implicit solvent and replica exchange sampling to calculate atomistic disordered ensembles of the nuclear co-activator binding domain (NCBD) of transcription coactivator CBP and the activation domain of the p160 steroid receptor coactivator ACTR. The calculated ensembles are in quantitative agreement with NMR-derived residue helicity and recapitulate the experimental observation that, while free ACTR largely lacks residual secondary structures, free NCBD is a molten globule with a helical content similar to that in the folded complex. Detailed conformational analysis reveals that free NCBD has an inherent ability to substantially sample all the helix configurations that have been previously observed either unbound or in complexes. Intriguingly, further high-temperature unbinding and unfolding simulations in implicit and explicit solvents emphasize the importance of conformational fluctuations in synergistic folding of NCBD with ACTR. A balance between preformed elements and conformational fluctuations appears necessary to allow NCBD to interact with different targets and fold into alternative conformations. Together with previous topology-based modeling and existing experimental data, the current simulations strongly support an “extended conformational selection” synergistic folding mechanism that involves a key intermediate state stabilized by interaction between the C-terminal helices of NCBD and ACTR. In addition, the atomistic simulations reveal the role of long-range as well as short-range electrostatic interactions in cooperating with readily fluctuating residual structures, which might enhance the encounter rate and promote efficient folding upon encounter for facile binding and folding interactions of IDPs. Thus, the current study not only provides a consistent mechanistic understanding of the NCBD/ACTR interaction, but also helps establish a multi-scale molecular modeling framework for understanding the structure, interaction, and regulation of IDPs in general.

Intrinsically disordered proteins (IDPs) are now widely recognized to play fundamental roles in biology and to be frequently associated with human diseases. Although the potential advantages of intrinsic disorder in cellular signaling and regulation have been widely discussed, the physical basis for these proposed phenomena remains sketchy at best. An integration of multi-scale molecular modeling and experimental characterization is necessary to uncover the molecular principles that govern the structure, interaction, and regulation of IDPs. In this work, we characterize the conformational properties of two IDPs involved in transcription regulation at the atomistic level and further examine the roles of these properties in their coupled binding and folding interactions. Our simulations suggest interplay among residual structures, conformational fluctuations, and electrostatic interactions that allows efficient synergistic folding of these two IDPs. In particular, we propose that electrostatic interactions might play an important role in facilitating rapid folding and binding recognition of IDPs, by enhancing the encounter rate and promoting efficient folding upon encounter.

Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent

The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, respectively, and are assumed to be unprotonated under physiological conditions. However, evidence from recent NMR and X-Ray crystallography studies has revealed the prevalence of protonated adenine and cytosine in RNA macromolecules. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first constant pH MD simulations (CPHMD) for nucleic acids in explicit solvent in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λ(Nexp) for λ that was recently developed for Multi-Site λ-Dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calculated pKa values of simple nucleotides are in a good agreement with experimentally measured values, with a mean absolute error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biological properties of RNA macromolecules.

Sequence-dependent pKa shift induced by molecular self-assembly: insights from computer simulation

The control of catalytic activity using molecular self-assembly is of fundamental interest. Recent experiments (Muller et al., Angew. Chem., Int. Ed., 2009, 48, 922-925) have demonstrated that two sequence isomers of β-peptides show remarkably different activity as an amine catalyst for a retro-aldol cleavage reaction, a difference attributed to the ability of one of the sequences to form large aggregates. The self-assembly and catalytic activity of these two isomers are investigated using constant pH molecular dynamics (CPHMD), for an atomistic model of β-peptides in implicit solvent. Simulations show that the globally amphiphilic (GA) isomer, which experimentally has high activity, forms large aggregates, while the non-GA isomer forms aggregates that are at most three or four molecules in size. The pK(a) shift of the βK-residues is significantly higher in the GA isomers that make a large aggregate. Since the decrease in pK(a) of the side-chain ammonium group is the main driving force for amine catalysis, the calculations are consistent with experiment. We find that the buried βK residues become entirely deprotonated, and the pK(a) shift for other titratable βK residues is accompanied mainly by a clustering of solvent exposed βK residues. We conclude that simulations can be used to understand catalytic activity due to self-assembly.

Predicting pKa values with continuous constant pH molecular dynamics

Knowledge of pK(a) values is important for understanding structure and function relationships in proteins. Over the past two decades, theoretical methods for pK(a) calculations have been mainly based on macroscopic models, in which the protein is considered as a low-dielectric cavity embedded in a high-dielectric continuum. In recent years, constant pH molecular dynamics methods have been developed based on a microscopic description of the protein. We describe here the methodology of continuous constant pH molecular dynamics (CPHMD), which has emerged as one of the most robust and accurate tools for predicting protein pK(a)s and for the study of pH-modulated conformational dynamics. We illustrate the utility of CPHMD by the calculation of pK(a)s for surface residues in ribonuclease A, buried residues in staphylococcal nuclease, and titratable groups in the intrinsically flexible protein α-lactalbumin. We will compare the CPHMD results with experimental data as well as calculations from PB-based and empirical methods. These examples demonstrate the accuracy and robustness of the CPHMD method and its ability to capture the correlation between ionization equilibria and conformational dynamics as well as the local dielectric response to structural rearrangement. Finally, we discuss future improvement of the CPHMD method.

Force-induced unfolding simulations of the human Notch1 negative regulatory region: possible roles of the heterodimerization domain in mechanosensing

Notch receptors are core components of the Notch signaling pathway and play a central role in cell fate decisions during development as well as tissue homeostasis. Upon ligand binding, Notch is sequentially cleaved at the S2 site by an ADAM protease and at the S3 site by the γ-secretase complex. Recent X-ray structures of the negative regulatory region (NRR) of the Notch receptor reveal an auto-inhibited fold where three protective Lin12/Notch repeats (LNR) of the NRR shield the S2 cleavage site housed in the heterodimerization (HD) domain. One of the models explaining how ligand binding drives the NRR conformation from a protease-resistant state to a protease-sensitive one invokes a mechanical force exerted on the NRR upon ligand endocytosis. Here, we combined physics-based atomistic simulations and topology-based coarse-grained modeling to investigate the intrinsic and force-induced folding and unfolding mechanisms of the human Notch1 NRR. The simulations support that external force applied to the termini of the NRR disengages the LNR modules from the heterodimerization (HD) domain in a well-defined, largely sequential manner. Importantly, the mechanical force can further drive local unfolding of the HD domain in a functionally relevant fashion that would provide full proteolytic access to the S2 site prior to heterodimer disassociation. We further analyzed local structural features, intrinsic folding free energy surfaces, and correlated motions of the HD domain. The results are consistent with a model in which the HD domain possesses inherent mechanosensing characteristics that could be utilized during Notch activation. This potential role of the HD domain in ligand-dependent Notch activation may have implications for understanding normal and aberrant Notch signaling.

Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange

A computational tool that offers accurate pKa values and atomically detailed knowledge of protonation-coupled conformational dynamics is valuable for elucidating mechanisms of energy transduction processes in biology, such as enzyme catalysis and electron transfer as well as proton and drug transport. Toward this goal we present a new technique of embedding continuous constant pH molecular dynamics within an explicit-solvent representation. In this technique we make use of the efficiency of the generalized-Born (GB) implicit-solvent model for estimating the free energy of protein solvation while propagating conformational dynamics using the more accurate explicit-solvent model. Also, we employ a pH-based replica exchange scheme to significantly enhance both protonation and conformational state sampling. Benchmark data of five proteins including HP36, NTL9, BBL, HEWL, and SNase yield an average absolute deviation of 0.53 and a root mean squared deviation of 0.74 from experimental data. This level of accuracy is obtained with 1 ns simulations per replica. Detailed analysis reveals that explicit-solvent sampling provides increased accuracy relative to the previous GB-based method by preserving the native structure, providing a more realistic description of conformational flexibility of the hydrophobic cluster, and correctly modeling solvent mediated ion-pair interactions. Thus, we anticipate that the new technique will emerge as a practical tool to capture ionization equilibria while enabling an intimate view of ionization coupled conformational dynamics that is difficult to delineate with experimental techniques alone.

Intrinsically disordered proteins in a physics-based world

Intrinsically disordered proteins (IDPs) are a newly recognized class of functional proteins that rely on a lack of stable structure for function. They are highly prevalent in biology, play fundamental roles, and are extensively involved in human diseases. For signaling and regulation, IDPs often fold into stable structures upon binding to specific targets. The mechanisms of these coupled binding and folding processes are of significant importance because they underlie the organization of regulatory networks that dictate various aspects of cellular decision-making. This review first discusses the challenge in detailed experimental characterization of these heterogeneous and dynamics proteins and the unique and exciting opportunity for physics-based modeling to make crucial contributions, and then summarizes key lessons from recent de novo simulations of the structure and interactions of several regulatory IDPs.

Constant pH replica exchange molecular dynamics in biomolecules using a discrete protonation model

A constant pH replica exchange molecular dynamics (REMD) method is proposed and implemented to improve coupled protonation and conformational state sampling. By mixing conformational sampling at constant pH (with discrete protonation states) with a temperature ladder, this method avoids conformational trapping. Our method was tested and applied to seven different biological systems. The constant pH REMD not only predicted pKa correctly for small, model compounds but also converged faster than constant pH molecular dynamics (MD). We further tested our constant pH REMD on a heptapeptide from ovomucoid third domain (OMTKY3). Although constant pH REMD and MD produced very close pKa values, the constant pH REMD showed its advantage in the efficiency of conformational and protonation state samplings.

Lambda-dynamics free energy simulation methods

Free energy calculations are fundamental to obtaining accurate theoretical estimates of many important biological phenomena including hydration energies, protein-ligand binding affinities and energetics of conformational changes. Unlike traditional free energy perturbation and thermodynamic integration methods, lambda-dynamics treats the conventional "lambda" as a dynamic variable in free energy simulations and simultaneously evaluates thermodynamic properties for multiple states in a single simulation. In the present article, we provide an overview of the theory of lambda-dynamics, including the use of biasing and restraining potentials to facilitate conformational sampling. We review how lambda-dynamics has been used to rapidly and reliably compute relative hydration free energies and binding affinities for series of ligands, to accurately identify crystallographically observed binding modes starting from incorrect orientations, and to model the effects of mutations upon protein stability. Finally, we suggest how lambda-dynamics may be extended to facilitate modeling efforts in structure-based drug design.

Partially folded equilibrium intermediate of the villin headpiece HP67 defined by 13C relaxation dispersion

Identification and characterization of ensembles of intermediate states remains an important objective in describing protein folding in atomic detail. The 67-residue villin headpiece, HP67, consists of an N-terminal subdomain (residues 10-42) that transiently unfolds at equilibrium under native-like conditions and a highly stable C-terminal subdomain (residues 43-76). The transition between folded and unfolded states of the N-terminal domain has been characterized previously by (15)N NMR relaxation dispersion measurements (Grey et al. in J Mol Biol 355:1078, 2006). In the present work, (13)C spin relaxation was used to further characterize backbone and hydrophobic core contributions to the unfolding process. Relaxation of (13)C(alpha) spins was measured using the Hahn echo technique at five static magnetic fields (11.7, 14.1, 16.4, 18.8, and 21.1 T) and the Carr-Purcell-Meiboom-Gill (CPMG) relaxation dispersion method at a static magnetic field of 14.1 T. Relaxation of methyl (13)C spins was measured using CPMG relaxation dispersion experiments at static magnetic fields of 14.1 and 18.8 T. Results for (13)C and (15)N spins yielded a consistent model in which the partially unfolded intermediate state of the N-terminal subdomain maintains residual structure for residues near the unprotonated His41 imidazole ring and in the interface between the N- and C-terminal subdomains. In addition, a second faster process was detected that appears to represent local dynamics within the folded state of the molecule and is largely confined to the hydrophobic interface between the N- and C-terminal subdomains.

Quantitative determination of site-specific conformational distributions in an unfolded protein by solid-state nuclear magnetic resonance

Solid-state nuclear magnetic resonance (NMR) techniques are used to investigate the structure of the 35-residue villin headpiece subdomain (HP35) in folded, partially denatured, and fully denatured states. Experiments are carried out in frozen glycerol/water solutions, with chemical denaturation by guanidine hydrochloride (GdnHCl). Without GdnHCl, two-dimensional solid-state (13)C NMR spectra of samples prepared with uniform (13)C labeling of selected residues show relatively sharp cross-peaks at chemical shifts that are consistent with the known three-helix bundle structure of folded HP35. At high GdnHCl concentrations, most cross-peaks broaden and shift, qualitatively indicating disruption of the folded structure and development of static conformational disorder in the frozen denatured state. Conformational distributions at one residue in each helical segment are probed quantitatively with three solid-state NMR techniques that provide independent constraints on backbone varphi and psi torsion angles in samples with sequential pairs of carbonyl (13)C labels. Without GdnHCl, the combined data are well fit by alpha-helical conformations. At [GdnHCl]=4.5 M, corresponding to the approximate denaturation midpoint, the combined data are well fit by a combination of alpha-helical and partially extended conformations at each site, but with a site-dependent population ratio. At [GdnHCl]=7.0 M, corresponding to the fully denatured state, the combined data are well fit by a combination of partially extended and polyproline II conformations, again with a site-dependent population ratio. Two entirely different models for conformational distributions lead to nearly the same best-fit distributions, demonstrating the robustness of these conclusions. This work represents the first quantitative investigation of site-specific conformational distributions in partially folded and unfolded states of a protein by solid-state NMR.

Atomistic details of the disordered states of KID and pKID. Implications in coupled binding and folding

Intrinsically disordered proteins (IDPs) are a newly recognized class of functional proteins for which a lack of stable tertiary fold is required for function. Because of the heterogeneous and dynamical nature, molecular modeling is necessary to provide the missing details of disordered states of IDP that are crucial for understanding their functions. In particular, generalized Born (GB) implicit solvent, combined with replica exchange (REX), might offer an optimal balance between accuracy and efficiency for modeling IDPs. We carried out extensive REX simulations in an optimized GB force field to characterize the disordered states of a regulatory IDP, KID domain of transcription factor CREB, and its phosphorylated form, pKID. The results revealed that both KID and pKID, though highly disordered on the tertiary level, are compact and mainly occupy a small number of helical substates. Interestingly, although phosphorylation of KID Ser133 leads only to marginal changes in average helicities on the ensemble level, underlying conformational substates differ significantly. In particular, pSer133 appears to restrict the accessible conformational space of the loop region and thus reduces the entropic cost of KID folding upon binding to the KIX domain of CREB-binding protein. Such an expanded role of phosphorylation in the KID:KIX recognition was not previously recognized because of a lack of substantial conformational changes on the ensemble level and inaccessibility of the structural details from experiments. The results also suggest that an implicit solvent-based modeling framework, despite various existing limitations, might be feasible for accurate atomistic simulation of small IDPs in general.

Heterogeneity and dynamics in villin headpiece crystal structures

The villin headpiece domain (HP67) is the C-terminal F-actin-binding motif that confers F-actin-bundling activity to villin, a component of the actin bundles that support the brush-border microvilli. It has been investigated extensively by both experimental and theoretical measurements. Our laboratory, for example, has determined both its NMR and its crystal structures. This study presents the structures of HP67 and its pH-stabilized mutant (H41Y) in a different crystal form and space group. For both constructs, two molecules are found in each asymmetric unit in the new space group P6(1). While one of the two structures (Mol A) is structurally similar to our previously determined structure (Mol X), the other (Mol B) has significant deviations, especially in the N-terminal subdomain, where lattice contacts do not appear to contribute to the difference. In addition, the structurally most different crystal structure, Mol B, is actually closer to the averaged NMR structure. Harmonic motions, as suggested by the B-factor profiles, differ between these crystal structures; crystal structures from the same space group share a similar pattern. Thus, heterogeneity and dynamics are observed in different crystal structures of the same protein even for a protein as small as villin headpiece.

Dielectric relaxation in proteins: the computational perspective

In photoexcitation and electron transfer, a new dipole or charge is introduced, and the structure is adjusted. This adjustment represents dielectric relaxation, which is the focus of this review. We concentrate on a few selected topics. We discuss linear response theory, as a unifying framework and a tool to describe non-equilibrium states. We review recent, molecular dynamics simulation studies that illustrate the calculation of dynamic and thermodynamic properties, such as Stokes shifts or reorganization free energies. We then turn to the macroscopic, continuum electrostatic view. We recall the physical definition of a dielectric constant and revisit the decomposition of the free energy into a reorganization and a static term. We review some illustrative continuum studies and discuss some difficulties that can arise with the continuum approach. In conclusion, we consider recent developments that will increase the accuracy and broaden the scope of all these methods.

Recent advances in implicit solvent-based methods for biomolecular simulations

Implicit solvent-based methods play an increasingly important role in molecular modeling of biomolecular structure and dynamics. Recent methodological developments have mainly focused on the extension of the generalized Born (GB) formalism for variable dielectric environments and accurate treatment of nonpolar solvation. Extensive efforts in parameterization of GB models and implicit solvent force fields have enabled ab initio simulation of protein folding to native or near-native structures. Another exciting area that has benefited from the advances in implicit solvent models is the development of constant pH molecular dynamics methods, which have recently been applied to the calculations of protein pK(a) values and the studies of pH-dependent peptide and protein folding.

Implicit modeling of nonpolar solvation for simulating protein folding and conformational transitions

Accurate description of the solvent environment is critical in computer simulations of protein structure and dynamics. An implicit treatment of solvent aims to capture the mean influence of water molecules on the solute via direct estimation of the solvation free energy. It has emerged as a powerful alternative to explicit solvent, and provides a favorable compromise between computational cost and level of detail. We review the current theory and techniques for implicit modeling of nonpolar solvation in the context of simulating protein folding and conformational transitions, and discuss the main directions for further development. It is demonstrated that the current surface area based nonpolar models have severe limitations, including insufficient description of the conformational dependence of solvation, over-estimation of the strength of pair-wise nonpolar interactions, and incorrect prediction of anti-cooperativity for three-body hydrophobic associations. We argue that, to improve beyond current level of accuracy of implicit solvent models, two important aspects of nonpolar solvation need to be incorporated, namely, the length-scale dependence of hydrophobic association and solvent screening of solute-solute dispersion interactions. We recognize that substantial challenges exist in constructing a sufficiently balanced, yet reasonably efficient, implicit solvent protein force field. Nonetheless, most of the fundamental problems are understood, and exciting progress has been made over the last few years. We believe that continual work along the frontiers outlined will greatly improve one's ability to study protein folding and large conformational transitions at atomistic detail.

Linking folding with aggregation in Alzheimer's beta-amyloid peptides

Growing evidence suggests that the beta-amyloid (Abeta) peptides of Alzheimer's disease are generated in early endosomes and that small oligomers are the principal toxic species. We sought to understand whether and how the solution pH, which is more acidic in endosomes than the extracellular environment, affects the conformational processes of Abeta. Using constant pH molecular dynamics simulations of two model peptides, Abeta(1-28) and Abeta(10-42), we found that the folding landscape of Abeta is strongly modulated by pH and is most favorable for hydrophobically driven aggregation at pH 6. Thus, our theoretical findings substantiate the possibility that Abeta oligomers develop intracellularly before secretion into the extracellular milieu, where they may disrupt synaptic activity or act as seeds for plaque formation.