Collapse kinetics and chevron plots from simulations of denaturant-dependent folding of globular proteins

Additivity of Transfer Free Energies Enables the Description of Complex Protein Formulations in Implicit Solvent Molecular Dynamics Simulations

A complex 3D structure and the surrounding environment determine the function and stability of a protein. Various osmolytes can be added to a protein drug formulation to stabilize the native protein structure by preventing unfolding and aggregation. In this context, the concept of transfer free energy, which represents the change in chemical potential of a protein being transferred from water to an osmolyte solution, has emerged as a powerful tool to elucidate the energetics involved in the protein-osmolyte interaction. In the present work, we experimentally determine the transfer free energies for the excipients sodium chloride, arginine hydrochloride, and polysorbate 20, which are frequently used in pharmaceutical protein formulations. We show that these excipients display distinct patterns of exclusion or interaction toward different moieties on the protein surface. Furthermore, we report that the free energy cost for transferring a protein to a formulation composed of multiple components can be calculated by summing up the contributions of the individual components. This finding suggests that additivity applies to the transfer free energies. We demonstrate that this additive behavior can be leveraged to accurately and efficiently model complex protein formulations. Additionally, we discuss how transfer free energies can be incorporated within implicit solvent molecular dynamics calculations, providing a direct link between experiments and simulations. Our molecular dynamics results show good agreement with experimental data for lysozyme, interferon α-2a, and granulocyte colony-stimulating factor, for both single- and multicomponent matrices, demonstrating the validity of our approach.

Salts Influence IDP Properties by Modulating the Population of Conformational Clusters

Salts readily alter the physical properties of intrinsically disordered proteins (IDPs) rich in charged residues. Using a coarse-grained IDP model and computer simulations, we investigated how salts affect the heterogeneous conformational ensemble and segment-level structures of the IDP prothymosin-α, classified as a polyelectrolyte. We show that clusters of conformations with distinct structural features are present within the conformational ensemble of prothymosin-α by projecting it onto a two-dimensional latent space with the aid of autoencoders. Although prothymosin-α is inherently disordered, there are preferred transitions between these clusters of conformations. Changing the salt concentration led to the formation of new conformational clusters or/and the disappearance of existing conformational clusters, contributing to changes in IDP properties. Shuffling the Skopelitian domain (C-terminal sequence) of prothymosin-α, known for its anticancer activity, resulted in a different conformational cluster, indicating that clusters with specific structures are related to a particular IDP function. The multiple conformational clusters with distinct structural features could be correlated to different IDP functions, and salts aid or inhibit these functions by modulating the population of conformations in the clusters.

SOP-MULTI: A Self-Organized Polymer-Based Coarse-Grained Model for Multidomain and Intrinsically Disordered Proteins with Conformation Ensemble Consistent with Experimental Scattering Data

Multidomain proteins with long flexible linkers and full-length intrinsically disordered proteins (IDPs) are best defined as an ensemble of conformations rather than a single structure. Determining high-resolution ensemble structures of such proteins poses various challenges by using tools from experimental structural biophysics. Integrative approaches combining available low-resolution ensemble-averaged experimental data and in silico biomolecular reconstructions are now often used for the purpose. However, extensive Boltzmann weighted conformation sampling for large proteins, especially for ones where both the folded and disordered domains exist in the same polypeptide chain, remains a challenge. In this work, we present a 2-site per amino-acid resolution SOP-MULTI force field for simulating coarse-grained models of multidomain proteins. SOP-MULTI combines two well-established self-organized polymer models─: (i) SOP-SC models for folded systems and (ii) SOP-IDP for IDPs. For the SOP-MULTI, we introduce cross-interaction terms between the beads belonging to the folded and disordered regions to generate conformation ensembles for full-length multidomain proteins such as hnRNP A1, TDP-43, G3BP1, hGHR-ECD, TIA1, HIV-1 Gag, polyubiquitin, and FUS. When back-mapped to all-atom resolution, SOP-MULTI trajectories faithfully recapitulate the scattering data over the range of the reciprocal space. We also show that individual folded domains preserve native contacts with respect to solved folded structures, and root-mean-square fluctuations of residues in folded domains match those obtained from all-atom molecular dynamics simulation trajectories of the same folded systems. SOP-MULTI force field is made available as a LAMMPS-compatible user package along with setup codes for generating the required files for any full-length protein with folded and disordered regions.

OLIVES: A Go̅-like Model for Stabilizing Protein Structure via Hydrogen Bonding Native Contacts in the Martini 3 Coarse-Grained Force Field

Coarse-grained molecular dynamics simulations enable the modeling of increasingly complex systems at millisecond timescales. The transferable coarse-grained force field Martini 3 has shown great promise in modeling a wide range of biochemical processes, yet folded proteins in Martini 3 are not stable without the application of external bias potentials, such as elastic networks or Go̅-like models. We herein develop an algorithm, called OLIVES, which identifies native contacts with hydrogen bond capabilities in coarse-grained proteins and use it to implement a novel Go̅-like model for Martini 3. We show that the protein structure instability originates in part from the lack of hydrogen bond energy in the coarse-grained force field representation. By using realistic hydrogen bond energies obtained from literature ab initio calculations, it is demonstrated that protein stability can be recovered by the reintroduction of a coarse-grained hydrogen bond network and that OLIVES removes the need for secondary structure restraints. OLIVES is validated against known protein complexes and at the same time addresses the open question of whether there is a need for protein quaternary structure bias in Martini 3 simulations. It is shown that OLIVES can reduce the number of bias terms, hereby speeding up Martini 3 simulations of proteins by up to ≈30% on a GPU architecture compared to the established Go̅MARTINI Go̅-like model.

Brewing COFFEE: A Sequence-Specific Coarse-Grained Energy Function for Simulations of DNA-Protein Complexes

DNA-protein interactions are pervasive in a number of biophysical processes ranging from transcription and gene expression to chromosome folding. To describe the structural and dynamic properties underlying these processes accurately, it is important to create transferable computational models. Toward this end, we introduce Coarse-grained Force Field for Energy Estimation, COFFEE, a robust framework for simulating DNA-protein complexes. To brew COFFEE, we integrated the energy function in the self-organized polymer model with side-chains for proteins and the three interaction site model for DNA in a modular fashion, without recalibrating any of the parameters in the original force-fields. A unique feature of COFFEE is that it describes sequence-specific DNA-protein interactions using a statistical potential (SP) derived from a data set of high-resolution crystal structures. The only parameter in COFFEE is the strength (λDNAPRO) of the DNA-protein contact potential. For an optimal choice of λDNAPRO, the crystallographic B-factors for DNA-protein complexes with varying sizes and topologies are quantitatively reproduced. Without any further readjustments to the force-field parameters, COFFEE predicts scattering profiles that are in quantitative agreement with small-angle X-ray scattering experiments, as well as chemical shifts that are consistent with NMR. We also show that COFFEE accurately describes the salt-induced unraveling of nucleosomes. Strikingly, our nucleosome simulations explain the destabilization effect of ARG to LYS mutations, which do not alter the balance of electrostatic interactions but affect chemical interactions in subtle ways. The range of applications attests to the transferability of COFFEE, and we anticipate that it would be a promising framework for simulating DNA-protein complexes at the molecular length-scale.

Salt Effects on Caffeine across Concentration Regimes

Salts affect the solvation thermodynamics of molecules of all sizes; the Hofmeister series is a prime example in which different ions lead to salting-in or salting-out of aqueous proteins. Early work of Tanford led to the discovery that the solvation of molecular surface motifs is proportional to the solvent accessible surface area (SASA), and later studies have shown that the proportionality constant varies with the salt concentration and type. Using multiscale computer simulations combined with vapor-pressure osmometry on caffeine-salt solutions, we reveal that this SASA description captures a rich set of molecular driving forces in tertiary solutions at changing solute and osmolyte concentrations. Central to the theoretical work is a new potential energy function that depends on the instantaneous surface area, salt type, and concentration. Used in, e.g., Monte Carlo simulations, this allows for a highly efficient exploration of many-body interactions and the resulting thermodynamics at elevated solute and salt concentrations.

Folding kinetics of an entangled protein

The possibility of the protein backbone adopting lasso-like entangled motifs has attracted increasing attention. After discovering the surprising abundance of natively entangled protein domain structures, it was shown that misfolded entangled subpopulations might become thermosensitive or escape the homeostasis network just after translation. To investigate the role of entanglement in shaping folding kinetics, we introduce a novel indicator and analyze simulations of a coarse-grained, structure-based model for two small single-domain proteins. The model recapitulates the well-known two-state folding mechanism of a non-entangled SH3 domain. However, despite its small size, a natively entangled antifreeze RD1 protein displays a rich refolding behavior, populating two distinct kinetic intermediates: a short-lived, entangled, near-unfolded state and a longer-lived, non-entangled, near-native state. The former directs refolding along a fast pathway, whereas the latter is a kinetic trap, consistently with known experimental evidence of two different characteristic times. Upon trapping, the natively entangled loop folds without being threaded by the N-terminal residues. After trapping, the native entangled structure emerges by either backtracking to the unfolded state or threading through the already formed but not yet entangled loop. Along the fast pathway, trapping does not occur because the native contacts at the closure of the lasso-like loop fold after those involved in the N-terminal thread, confirming previous predictions. Despite this, entanglement may appear already in unfolded configurations. Remarkably, a longer-lived, near-native intermediate, with non-native entanglement properties, recalls what was observed in cotranslational folding.

Brewing COFFEE: A sequence-specific coarse-grained energy function for simulations of DNA-protein complexes

DNA-protein interactions are pervasive in a number of biophysical processes ranging from transcription, gene expression, to chromosome folding. To describe the structural and dynamic properties underlying these processes accurately, it is important to create transferable computational models. Toward this end, we introduce Coarse grained force field for energy estimation, COFFEE, a robust framework for simulating DNA-protein complexes. To brew COFFEE, we integrated the energy function in the Self-Organized Polymer model with Side Chains for proteins and the Three Interaction Site model for DNA in a modular fashion, without re-calibrating any of the parameters in the original force-fields. A unique feature of COFFEE is that it describes sequence-specific DNA-protein interactions using a statistical potential (SP) derived from a dataset of high-resolution crystal structures. The only parameter in COFFEE is the strength (λ D N A P R O ) of the DNA-protein contact potential. For an optimal choice of λ D N A P R O , the crystallographic B-factors for DNA-protein complexes, with varying sizes and topologies, are quantitatively reproduced. Without any further readjustments to the force-field parameters, COFFEE predicts the scattering profiles that are in quantitative agreement with SAXS experiments as well as chemical shifts that are consistent with NMR. We also show that COFFEE accurately describes the salt-induced unraveling of nucleosomes. Strikingly, our nucleosome simulations explain the destabilization effect of ARG to LYS mutations, which does not alter the balance of electrostatic interactions, but affects chemical interactions in subtle ways. The range of applications attests to the transferability of COFFEE, and we anticipate that it would be a promising framework for simulating DNA-protein complexes at the molecular length-scale.

Residue-Dependent Transition Temperatures and Denaturant Midpoints in the Folding of a Multidomain Protein

As a consequence of the finite size of globular proteins, it is expected that there should be dispersions in the global melting temperature (T_m) and the denaturation midpoint (C_m). Thermodynamic considerations dictate that the dispersions, ΔT_m in T_m, and ΔC_m in C_m, should decrease with N, the number of residues in the protein. We performed coarse-grained simulations of the self-organized polymer (SOP) model of the multidomain protein adenylate kinase (ADK) with N = 214 in order to calculate thermal and denaturation unfolding titration curves. The results show that ΔT_m/T_m and ΔC_m/C_m are nonzero and follow the previously established ( Phys. Rev. Lett. 2004, 93, 268107) thermodynamic 1/N scaling for proteins accurately. For ADK, the dispersions are small (≈0.004), which implies that the melting temperature is more or less unique, which is unlike in BBL (N = 40) where ΔT_m/T_m ≈ 0.03.

Metamorphic proteins under a computational microscope: Lessons from a fold-switching RfaH protein

Metamorphic proteins constitute unexpected paradigms of the protein folding problem, as their sequences encode two alternative folds, which reversibly interconvert within biologically relevant timescales to trigger different cellular responses. Once considered a rare aberration, metamorphism may be common among proteins that must respond to rapidly changing environments, exemplified by NusG-like proteins, the only transcription factors present in every domain of life. RfaH, a specialized paralog of bacterial NusG, undergoes an all-α to all-β domain switch to activate expression of virulence and conjugation genes in many animal and plant pathogens and is the quintessential example of a metamorphic protein. The dramatic nature of RfaH structural transformation and the richness of its evolutionary history makes for an excellent model for studying how metamorphic proteins switch folds. Here, we summarize the structural and functional evidence that sparked the discovery of RfaH as a metamorphic protein, the experimental and computational approaches that enabled the description of the molecular mechanism and refolding pathways of its structural interconversion, and the ongoing efforts to find signatures and general properties to ultimately describe the protein metamorphome.

A Transfer Free Energy Based Implicit Solvent Model for Protein Simulations in Solvent Mixtures: Urea-Induced Denaturation as a Case Study

We developed a method for implicit solvent molecular dynamics simulations of proteins in solvent mixtures (model with implicit solvation thermodynamics, MIST). The MIST method introduces experimental group transfer free energies to the generalized Born formulation for generating molecular trajectories without the need for developing rigorous explicit-solvent force fields for multicomponent solutions. As a test case, we studied the urea-induced denaturation of the Trp-cage miniprotein in water. We demonstrate that our method allows efficient exploration of the conformational space of the protein in only a few hundreds of nanoseconds of all-atom unbiased simulations. Furthermore, selective implementation of the transfer free energies of specific peptide groups, backbone, and side chains enables us to decouple their specific energetic contributions to the conformational changes of the protein. The approach herein developed can readily be extended to the investigation of complex matrices as well as to the characterization of protein aggregation. The MIST method is implemented in Plumed (ver. 2.8) as a separate module called SASA.

Cooperativity and Folding Kinetics in a Multidomain Protein with Interwoven Chain Topology

Although a large percentage of eukaryotic proteomes consist of proteins with multiple domains, not much is known about their assembly mechanism, especially those with intricate native state architectures. Some have a complex topology in which the structural elements along the sequence are interwoven in such a manner that the domains cannot be separated by cutting at any location along the sequence. Such proteins are multiply connected multidomain proteins (MMPs) with the three-domain (NMP, LID, and CORE) phosphotransferase enzyme adenylate kinase (ADK) being an example. We devised a coarse-grained model to simulate ADK folding initiated by changing either the temperature or guanidinium chloride (GdmCl) concentration. The simulations reproduce the experimentally measured melting temperatures (associated with two equilibrium transitions), FRET efficiency as a function of GdmCl concentration, and the folding times quantitatively. Although the NMP domain orders independently, cooperative interactions between the LID and the CORE domains are required for complete assembly of the enzyme. Kinetic simulations show that, on the collapse time scale, multiple interconnected metastable states are populated, attesting to the folding heterogeneity. The network of kinetically connected states reveals that the CORE domain folds only after the NMP and LID domains, reflecting the interwoven nature of the chain topology.

Temperature and Guanidine Hydrochloride Effects on the Folding Thermodynamics of WW Domain and Variants

We used simulations based on an all-atom Go model to calculate the folding temperatures (T_fs) and free energies (ΔGs) of two variants of the WW domain, which is a small all-β-sheet protein. The results, without adjusting any parameter, are in good agreement with experiments, thus validating the simulations. We then used the molecular transfer model to predict the changes in the ΔGs and T_fs as the guanidine hydrochloride concentration is varied. The predictions can be readily tested in experiments.

Shape changes and cooperativity in the folding of the central domain of the 16S ribosomal RNA

Both the small and large subunits of the ribosome, the molecular machine that synthesizes proteins, are complexes of ribosomal RNAs (rRNAs) and a number of proteins. In bacteria, the small subunit has a single 16S rRNA whose folding is the first step in its assembly. The central domain of the 16S rRNA folds independently, driven either by Mg²⁺ ions or by interaction with ribosomal proteins. To provide a quantitative description of ion-induced folding of the ∼350-nucleotide rRNA, we carried out extensive coarse-grained molecular simulations spanning Mg²⁺ concentration between 0 and 30 mM. The Mg²⁺ dependence of the radius of gyration shows that globally the rRNA folds cooperatively. Surprisingly, various structural elements order at different Mg²⁺ concentrations, indicative of the heterogeneous assembly even within a single domain of the rRNA. Binding of Mg²⁺ ions is highly specific, with successive ion condensation resulting in nucleation of tertiary structures. We also predict the Mg²⁺-dependent protection factors, measurable in hydroxyl radical footprinting experiments, which corroborate the specificity of Mg²⁺-induced folding. The simulations, which agree quantitatively with several experiments on the folding of a three-way junction, show that its folding is preceded by formation of other tertiary contacts in the central junction. Our work provides a starting point in simulating the early events in the assembly of the small subunit of the ribosome.

Energy Landscape of Ubiquitin Is Weakly Multidimensional

Single molecule pulling experiments report time-dependent changes in the extension (X) of a biomolecule as a function of the applied force (f). By fitting the data to one-dimensional analytical models of the energy landscape, we can extract the hopping rates between the folded and unfolded states in two-state folders as well as the height and the location of the transition state (TS). Although this approach is remarkably insightful, there are cases for which the energy landscape is multidimensional (catch bonds being the most prominent). To assess if the unfolding energy landscape in small single domain proteins could be one-dimensional, we simulated force-induced unfolding of ubiquitin (Ub) using the coarse-grained self-organized polymer-side chain (SOP-SC) model. Brownian dynamics simulations using the SOP-SC model reveal that the Ub energy landscape is weakly multidimensional (WMD), governed predominantly by a single barrier. The unfolding pathway is confined to a narrow reaction pathway that could be described as diffusion in a quasi-1D X-dependent free energy profile. However, a granular analysis using the P_fold analysis, which does not assume any form for the reaction coordinate, shows that X alone does not account for the height and, more importantly, the location of the TS. The f-dependent TS location moves toward the folded state as f increases, in accord with the Hammond postulate. Our study shows that, in addition to analyzing the f-dependent hopping rates, the transition state ensemble must also be determined without resorting to X as a reaction coordinate to describe the unfolding energy landscapes of single domain proteins, especially if they are only WMD.

Double Domain Swapping in Human γC and γD Crystallin Drives Early Stages of Aggregation

Human γD (HγD) and γC (HγC) are two-domain crystallin (Crys) proteins expressed in the nucleus of the eye lens. Structural perturbations in the protein often trigger aggregation, which eventually leads to cataract. To decipher the underlying molecular mechanism, it is important to characterize the partially unfolded conformations, which are aggregation-prone. Using a coarse grained protein model and molecular dynamics simulations, we studied the role of on-pathway folding intermediates in the early stages of aggregation. The multidimensional free energy surface revealed at least three different folding pathways with the population of partially structured intermediates. The two dominant pathways confirm sequential folding of the N-terminal [Ntd] and the C-terminal domains [Ctd], while the third, least favored, pathway involves intermediates where both the domains are partially folded. A native-like intermediate (I*), featuring the folded domains and disrupted interdomain contacts, gets populated in all three pathways. I* forms domain swapped dimers by swapping the entire Ntds and Ctds with other monomers. Population of such oligomers can explain the increased resistance to unfolding resulting in hysteresis observed in the folding experiments of HγD Crys. An ensemble of double domain swapped dimers are also formed during refolding, where intermediates consisting of partially folded Ntds and Ctds swap secondary structures with other monomers. The double domain swapping model presented in our study provides structural insights into the early events of aggregation in Crys proteins and identifies the key secondary structural swapping elements, where introducing mutations will aid in regulating the overall aggregation propensity.

Seeded Self-Assembly of Charge-Terminated Poly(3-hexylthiophene) Amphiphiles Based on the Energy Landscape

The creation of 1D π-conjugated nanofibers with precise control and optimized optoelectronic properties is of widespread interest for applications as nanowires. "Living" crystallization-driven self-assembly (CDSA) is a seeded growth method of growing importance for the preparation of uniform 1D fiber-like micelles from a range of crystallizable polymeric amphiphiles. However, in the case of polythiophenes, one of the most important classes of conjugated polymer, only limited success has been achieved to date using block copolymers as precursors. Herein, we describe studies of the living CDSA of phosphonium-terminated amphiphilic poly(3-hexylthiophene)s to prepare colloidally stable nanofibers. In depth studies of the relationship between the degree of polymerization and the self-assembly behavior permitted the unveiling of the energy landscape of the living CDSA process. On the basis of the kinetic and thermodynamic insight provided, we have been able to achieve an unprecedented level of control over the length of low dispersity fiber-like micelles from 40 nm to 2.8 μm.

Gō model revisited

This review discusses Gō models broadly used in biomolecular simulations. I start with a brief description of the original lattice model study by Nobuhiro Gō. Then, the theory of protein folding behind Gō model, free energy approaches, and off-lattice Gō models are reviewed. I also mention a stringent test for the assumption in Gō models given from all-atom molecular dynamics simulations. Subsequently, I move to application of Gō models to protein dynamical functions. Various extension of Gō models is also reviewed. Finally, some publicly available tools to use Gō models are listed.

Sequence Effects on Size, Shape, and Structural Heterogeneity in Intrinsically Disordered Proteins

Intrinsically disordered proteins (IDPs) lack well-defined three-dimensional structures, thus challenging the archetypal notion of structure-function relationships. Determining the ensemble of conformations that IDPs explore under physiological conditions is the first step toward understanding their diverse cellular functions. Here, we quantitatively characterize the structural features of IDPs as a function of sequence and length using coarse-grained simulations. For diverse IDP sequences, with the number of residues ( NT) ranging from 20 to 441, our simulations not only reproduce the radii of gyration ( Rg) obtained from experiments, but also predict the full scattering intensity profiles in excellent agreement with small-angle X-ray scattering experiments. The Rg values are well-described by the standard Flory scaling law, Rg = Rg0 NTν, with ν ≈ 0.588, making it tempting to assert that IDPs behave as polymers in a good solvent. However, clustering analysis reveals that the menagerie of structures explored by IDPs is diverse, with the extent of heterogeneity being highly sequence-dependent, even though ensemble-averaged properties, such as the dependence of Rg on chain length, may suggest synthetic polymer-like behavior in a good solvent. For example, we show that for the highly charged Prothymosin-α, a substantial fraction of conformations is highly compact. Even if the sequence compositions are similar, as is the case for α-Synuclein and a truncated construct from the Tau protein, there are substantial differences in the conformational heterogeneity. Taken together, these observations imply that metrics based on net charge or related quantities alone cannot be used to anticipate the phases of IDPs, either in isolation or in complex with partner IDPs or RNA. Our work sets the stage for probing the interactions of IDPs with each other, with folded protein domains, or with partner RNAs, which are critical for describing the structures of stress granules and biomolecular condensates with important cellular functions.

Role of Disulfide Bonds and Topological Frustration in the Kinetic Partitioning of Lysozyme Folding Pathways

Disulfide bonds in proteins can strongly influence the folding pathways by constraining the conformational space. Lysozyme has four disulfide bonds and is widely studied for its antibacterial properties. Experiments on lysozyme infer that the protein folds through a fast and a slow pathway. However, the reasons for the kinetic partitioning in the folding pathways are not completely clear. Using a coarse-grained protein model and simulations, we show that two out of the four disulfide bonds, which are present in the α-domain of lysozyme, are responsible for the slow folding pathway. In this pathway, a kinetically trapped intermediate state, which is close to the native state, is populated. In this state, the orientations of α-helices present in the α-domain are misaligned relative to each other. The protein in this state has to partially unfold by breaking down the interhelical contacts between the misaligned helices to fold to the native state. However, the topological constraints due to the two disulfide bonds present in the α-domain make the protein less flexible, and it is trapped in this conformation for hundreds of milliseconds. On disabling these disulfide bonds, we find that the kinetically trapped intermediate state and the slow folding pathway disappear. Simulations mimicking the folding of protein without disulfide bonds under oxidative conditions show that the native disulfide bonds are formed as the protein folds, indicating that folding guides the formation of disulfide bonds. The sequence of formation of the disulfide bonds is Cys64-Cys80 → Cys76-Cys94 → Cys30-Cys115 → Cys6-Cys127. Any disulfide bond that forms before its precursor in the sequence has to break and follow the sequence for the protein to fold. These results show that lysozyme also serves as a very good model system to probe the role of disulfide bonds and topological frustration in protein folding. The predictions from the simulations can be verified by single-molecule fluorescence resonance energy transfer or single-molecule pulling experiments, which can probe heterogeneity in the folding pathways.

Salt Effects on Protein Folding Thermodynamics

Salts differ in their ability to stabilize protein conformations, thereby affecting the thermodynamics and kinetics of protein folding. We developed a coarse-grained protein model that can predict salt-induced changes in protein properties by using the transfer free-energy data of various chemical groups from water to salt solutions. Using this model and molecular dynamics simulations, we probed the effect of seven different salts on the folding thermodynamics of the DNA binding domain of lac repressor protein ( lac-DBD) and N-terminal domain of ribosomal protein (NTL9). We show that a salt can act as a protein stabilizing or destabilizing agent depending on the protein sequence and folded state topology. The computed thermodynamic properties, especially the m values for various salts, which reveal the relative ability of a salt to stabilize the protein folded state, are in quantitative agreement with the experimentally measured values. The computations show that the degree of protein compaction in the denatured ensemble strongly depends on the salt identity, and for the same variation in salt concentration, the compaction in the protein dimensions varies from ∼4% to ∼30% depending on the salt. The transition-state ensemble (TSE) of lac-DBD is homogeneous and polarized, while the TSE of NTL9 is heterogeneous and diffusive. Salts induce subtle structural changes in the TSE that are in agreement with Hammond's postulate. The barrier to protein folding tends to disappear in the presence of moderate concentrations (∼3-4 m) of strongly stabilizing salts.

Denaturants Alter the Flux through Multiple Pathways in the Folding of PDZ Domain

Although we understand many aspects of how small proteins (number of residues less than about hundred) fold, it is a major challenge to quantitatively describe how large proteins self-assemble. To partially overcome this challenge, we performed simulations using the self-organized polymer model with side chains (SOP-SC) in guanidinium chloride (GdmCl), using the molecular transfer model (MTM), to describe the folding of the 110-residue PDZ3 domain. The simulations reproduce the folding thermodynamics accurately including the melting temperature (T_m), the stability of the folded state with respect to the unfolded state. We show that the calculated dependence of ln k_obs (k_obs is the relaxation rate) has the characteristic chevron shape. The slopes of the chevron plots are in good agreement with experiments. We show that PDZ3 folds by four major pathways populating two metastable intermediates, in accord with the kinetic partitioning mechanism. The structure of one of the intermediates, populated after polypeptide chain collapse, is structurally similar to an equilibrium intermediate. Surprisingly, the connectivities between the intermediates and hence, the fluxes through the pathways depend on the concentration of GdmCl. The results are used to predict possible outcomes for unfolding of PDZ domain subject to mechanical forces. Our study demonstrates that, irrespective of the size or topology, simulations based on MTM and SOP-SC offer a theoretical framework for describing the folding of proteins, mimicking precisely the conditions used in experiments.

Site-specific time-resolved FRET reveals local variations in the unfolding mechanism in an apparently two-state protein unfolding transition

Using multi-site time-resolved FRET, it is shown that equilibrium unfolding of monellin is not only heterogeneous, but that the degree of non-cooperativity differs between the sole α-helix and different parts of the β-sheet.

Modeling the effects of lattice defects on microtubule breaking and healing

Microtubule reorganization often results from the loss of polymer induced through breakage or active destruction by energy-using enzymes. Pre-existing defects in the microtubule lattice likely lower structural integrity and aid filament destruction. Using large-scale molecular simulations, we model diverse microtubule fragments under forces generated at specific positions to locally crush the filament. We show that lattices with 2% defects are crushed and severed by forces three times smaller than defect-free ones. We validate our results with direct comparisons of microtubule kinking angles during severing. We find a high statistical correlation between the angle distributions from experiments and simulations indicating that they sample the same population of structures. Our simulations also indicate that the mechanical environment of the filament affects breaking: local mechanical support inhibits healing after severing, especially in the case of filaments with defects. These results recall reports of microtubule healing after flow-induced bending and corroborate prior experimental studies that show severing is more likely at locations where microtubules crossover in networks. Our results shed new light on mechanisms underlying the ability of microtubules to be destroyed and healed in the cell, either by external forces or by severing enzymes wedging dimers apart. © 2016 Wiley Periodicals, Inc.

Linking computation and experiments to study the role of charge-charge interactions in protein folding and stability

Over the past two decades there has been an increase in appreciation for the role of surface charge-charge interactions in protein folding and stability. The perception shifted from the belief that charge-charge interactions are not important for protein folding and stability to the near quantitative understanding of how these interactions shape the folding energy landscape. This led to the ability of computational approaches to rationally redesign surface charge-charge interactions to modulate thermodynamic properties of proteins. Here we summarize our progress in understanding the role of charge-charge interactions for protein stability using examples drawn from my own laboratory and touch upon unanswered questions.

Molecular investigations into the mechanics of a muscle anchoring complex

The titin-telethonin complex, essential for anchoring filaments in the Z-disk of the sarcomere, is composed of immunoglobulin domains. Surprisingly, atomic force microscopy experiments showed that it resists forces much higher than the typical immunoglobulin domain and that the force distribution is unusually broad. To investigate the origin of this behavior, we developed a multiscale simulation approach, combining minimalist and atomistic models (SOP-AT). By following the mechanical response of the complex on experimental timescales, we found that the mechanical stability of titin-telethonin is modulated primarily by the strength of contacts between telethonin and the two titin chains, and secondarily by the timescales of conformational excursions inside telethonin and the pulled titin domains. Importantly, the conformational transitions executed by telethonin in simulations support its proposed role in mechanosensing. Our SOP-AT computational approach thus provides a powerful tool for the exploration of the link between conformational diversity and the broadness of the mechanical response, which can be applied to other multidomain complexes.

Folding of Protein L with Implications for Collapse in the Denatured State Ensemble

A fundamental question in protein folding is whether the coil to globule collapse transition occurs during the initial stages of folding (burst phase) or simultaneously with the protein folding transition. Single molecule fluorescence resonance energy transfer (FRET) and small-angle X-ray scattering (SAXS) experiments disagree on whether Protein L collapse transition occurs during the burst phase of folding. We study Protein L folding using a coarse-grained model and molecular dynamics simulations. The collapse transition in Protein L is found to be concomitant with the folding transition. In the burst phase of folding, we find that FRET experiments overestimate radius of gyration, Rg, of the protein due to the application of Gaussian polymer chain end-to-end distribution to extract Rg from the FRET efficiency. FRET experiments estimate ≈6 Å decrease in Rg when the actual decrease is ≈3 Å on guanidinium chloride denaturant dilution from 7.5 to 1 M, thereby suggesting pronounced compaction in the protein dimensions in the burst phase. The ≈3 Å decrease is close to the statistical uncertainties of the Rg data measured from SAXS experiments, which suggest no compaction, leading to a disagreement with the FRET experiments. The transition-state ensemble (TSE) structures in Protein L folding are globular and extensive in agreement with the Ψ-analysis experiments. The results support the hypothesis that the TSE of single domain proteins depends on protein topology and is not stabilized by local interactions alone.

Force-dependent switch in protein unfolding pathways and transition-state movements

Although it is known that single-domain proteins fold and unfold by parallel pathways, demonstration of this expectation has been difficult to establish in experiments. Unfolding rate, [Formula: see text], as a function of force f, obtained in single-molecule pulling experiments on src SH3 domain, exhibits upward curvature on a [Formula: see text] plot. Similar observations were reported for other proteins for the unfolding rate [Formula: see text]. These findings imply unfolding in these single-domain proteins involves a switch in the pathway as f or [Formula: see text] is increased from a low to a high value. We provide a unified theory demonstrating that if [Formula: see text] as a function of a perturbation (f or [Formula: see text]) exhibits upward curvature then the underlying energy landscape must be strongly multidimensional. Using molecular simulations we provide a structural basis for the switch in the pathways and dramatic shifts in the transition-state ensemble (TSE) in src SH3 domain as f is increased. We show that a single-point mutation shifts the upward curvature in [Formula: see text] to a lower force, thus establishing the malleability of the underlying folding landscape. Our theory, applicable to any perturbation that affects the free energy of the protein linearly, readily explains movement in the TSE in a β-sandwich (I27) protein and single-chain monellin as the denaturant concentration is varied. We predict that in the force range accessible in laser optical tweezer experiments there should be a switch in the unfolding pathways in I27 or its mutants.

Using the folding landscapes of proteins to understand protein function

Proteins fold on a biologically-relevant timescale because of a funnel-shaped energy landscape. This landscape is sculpted through evolution by selecting amino-acid sequences that stabilize native interactions while suppressing stable non-native interactions that occur during folding. However, there is strong evolutionary selection for functional residues and these cannot be chosen to optimize folding. Their presence impacts the folding energy landscape in a variety of ways. Here, we survey the effects of functional residues on folding by providing several examples. We then review how such effects can be detected computationally and be used as assays for protein function. Overall, an understanding of how functional residues modulate folding should provide insights into the design of natural proteins and their homeostasis.

Functionally Relevant Specific Packing Can Determine Protein Folding Routes

Functional residues can modulate the folding mechanisms of proteins. In some proteins, mutations to such residues can radically change the primary folding route. Is it possible then to learn more about the functional regions of a protein by investigating just its choice of folding route? The folding and the function of the protein Escherichia coli ribonuclease H (ecoRNase-H) have been extensively studied and its folding route is known to near-residue resolution. Here, we computationally study the folding of ecoRNase-H using molecular dynamics simulations of structure-based models of increasing complexity. The differences between a model that correctly predicts the experimentally determined folding route and a simpler model that does not can be attributed to a set of six aromatic residues clustered together in a region of the protein called CORE. This clustering, which we term "specific" packing, drives CORE to fold early and determines the folding route. Both the residues involved in specific packing and their packing are largely conserved across E. coli-like RNase-Hs from diverse species. Residue conservation is usually implicated in function. Here, the identified residues either are known to bind substrate in ecoRNase-H or pack against the substrate in the homologous human RNase-H where a substrate-bound crystal structure exists. Thus, the folding mechanism of ecoRNase-H is a byproduct of functional demands upon its sequence. Using our observations on specific packing, we suggest mutations to an engineered HIV RNase-H to make its function better. Our results show that understanding folding route choice in proteins can provide unexpected insights into their function.

Topological and sequence information predict that foldons organize a partially overlapped and hierarchical structure

It has been suggested that proteins have substructures, called foldons, which can cooperatively fold into the native structure. However, several prior investigations define foldons in various ways, citing different foldon characteristics, thereby making the concept of a foldon ambiguous. In this study, we perform a Gō model simulation and analyze the characteristics of substructures that cooperatively fold into the native-like structure. Although some results do not agree well with the experimental evidence due to the simplicity of our coarse-grained model, our results strongly suggest that cooperatively folding units sometimes organize a partially overlapped and hierarchical structure. This view makes us easy to interpret some different proposal about the foldon as a difference of the hierarchical structure. On the basis of this finding, we present a new method to assign foldons and their hierarchy, using structural and sequence information. The results show that the foldons assigned by our method correspond to the intermediate structures identified by some experimental techniques. The new method makes it easy to predict whether a protein folds sequentially into the native structure or whether some foldons fold into the native structure in parallel.

Nucleotides regulate the mechanical hierarchy between subdomains of the nucleotide binding domain of the Hsp70 chaperone DnaK

The regulation of protein function through ligand-induced conformational changes is crucial for many signal transduction processes. The binding of a ligand alters the delicate energy balance within the protein structure, eventually leading to such conformational changes. In this study, we elucidate the energetic and mechanical changes within the subdomains of the nucleotide binding domain (NBD) of the heat shock protein of 70 kDa (Hsp70) chaperone DnaK upon nucleotide binding. In an integrated approach using single molecule optical tweezer experiments, loop insertions, and steered coarse-grained molecular simulations, we find that the C-terminal helix of the NBD is the major determinant of mechanical stability, acting as a glue between the two lobes. After helix unraveling, the relative stability of the two separated lobes is regulated by ATP/ADP binding. We find that the nucleotide stays strongly bound to lobe II, thus reversing the mechanical hierarchy between the two lobes. Our results offer general insights into the nucleotide-induced signal transduction within members of the actin/sugar kinase superfamily.

Native contact density and nonnative hydrophobic effects in the folding of bacterial immunity proteins

The bacterial colicin-immunity proteins Im7 and Im9 fold by different mechanisms. Experimentally, at pH 7.0 and 10°C, Im7 folds in a three-state manner via an intermediate but Im9 folding is two-state-like. Accordingly, Im7 exhibits a chevron rollover, whereas the chevron arm for Im9 folding is linear. Here we address the biophysical basis of their different behaviors by using native-centric models with and without additional transferrable, sequence-dependent energies. The Im7 chevron rollover is not captured by either a pure native-centric model or a model augmented by nonnative hydrophobic interactions with a uniform strength irrespective of residue type. By contrast, a more realistic nonnative interaction scheme that accounts for the difference in hydrophobicity among residues leads simultaneously to a chevron rollover for Im7 and an essentially linear folding chevron arm for Im9. Hydrophobic residues identified by published experiments to be involved in nonnative interactions during Im7 folding are found to participate in the strongest nonnative contacts in this model. Thus our observations support the experimental perspective that the Im7 folding intermediate is largely underpinned by nonnative interactions involving large hydrophobics. Our simulation suggests further that nonnative effects in Im7 are facilitated by a lower local native contact density relative to that of Im9. In a one-dimensional diffusion picture of Im7 folding with a coordinate- and stability-dependent diffusion coefficient, a significant chevron rollover is consistent with a diffusion coefficient that depends strongly on native stability at the conformational position of the folding intermediate.

In order to fold correctly, a globular protein must avoid being trapped in wrong, i.e., nonnative conformations. Thus a biophysical account of how attractive nonnative interactions are bypassed by some amino acid sequences but not others is key to deciphering protein structure and function. We examine two closely related bacterial immunity proteins, Im7 and Im9, that are experimentally known to fold very differently: Whereas Im9 folds directly, Im7 folds through a mispacked conformational intermediate. A simple model we developed accounts for their intriguingly different folding kinetics in terms of a balance between the density of native-promoting contacts and the hydrophobicity of local amino acid sequences. This emergent principle is extensible to other biomolecular recognition processes.

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.

The role of high-dimensional diffusive search, stabilization, and frustration in protein folding

Proteins are polymeric molecules with many degrees of conformational freedom whose internal energetic interactions are typically screened to small distances. Therefore, in the high-dimensional conformation space of a protein, the energy landscape is locally relatively flat, in contrast to low-dimensional representations, where, because of the induced entropic contribution to the full free energy, it appears funnel-like. Proteins explore the conformation space by searching these flat subspaces to find a narrow energetic alley that we call a hypergutter and then explore the next, lower-dimensional, subspace. Such a framework provides an effective representation of the energy landscape and folding kinetics that does justice to the essential characteristic of high-dimensionality of the search-space. It also illuminates the important role of nonnative interactions in defining folding pathways. This principle is here illustrated using a coarse-grained model of a family of three-helix bundle proteins whose conformations, once secondary structure has formed, can be defined by six rotational degrees of freedom. Two folding mechanisms are possible, one of which involves an intermediate. The stabilization of intermediate subspaces (or states in low-dimensional projection) in protein folding can either speed up or slow down the folding rate depending on the amount of native and nonnative contacts made in those subspaces. The folding rate increases due to reduced-dimension pathways arising from the mere presence of intermediate states, but decreases if the contacts in the intermediate are very stable and introduce sizeable topological or energetic frustration that needs to be overcome. Remarkably, the hypergutter framework, although depending on just a few physically meaningful parameters, can reproduce all the types of experimentally observed curvature in chevron plots for realizations of this fold.

Benchmarking all-atom simulations using hydrogen exchange

Long-time molecular dynamics (MD) simulations are now able to fold small proteins reversibly to their native structures [Lindorff-Larsen K, Piana S, Dror RO, Shaw DE (2011) Science 334(6055):517-520]. These results indicate that modern force fields can reproduce the energy surface near the native structure. To test how well the force fields recapitulate the other regions of the energy surface, MD trajectories for a variant of protein G are compared with data from site-resolved hydrogen exchange (HX) and other biophysical measurements. Because HX monitors the breaking of individual H-bonds, this experimental technique identifies the stability and H-bond content of excited states, thus enabling quantitative comparison with the simulations. Contrary to experimental findings of a cooperative, all-or-none unfolding process, the simulated denatured state ensemble, on average, is highly collapsed with some transient or persistent native 2° structure. The MD trajectories of this protein G variant and other small proteins exhibit excessive intramolecular H-bonding even for the most expanded conformations, suggesting that the force fields require improvements in describing H-bonding and backbone hydration. Moreover, these comparisons provide a general protocol for validating the ability of simulations to accurately capture rare structural fluctuations.

Hamiltonian Mapping Revisited: Calibrating Minimalist Models to Capture Molecular Recognition by Intrinsically Disordered Proteins

Molecular recognition by intrinsically disordered proteins (IDPs) plays a central role in many critical cellular processes. Toward achieving detailed mechanistic understanding of IDP-target interactions, here we employ the "Hamiltonian mapping" methodology, which is rooted in the weighted histogram analysis method (WHAM), for the fast and efficient calibration of structure-based models in studies of IDPs. By performing reference simulations on a given Hamiltonian, we illustrate for two model IDPs how this method can extrapolate thermodynamic behavior under a range of modified Hamiltonians, in this case representing changes in the binding affinity (K_d) of the system. Given sufficient conformational sampling in a single trajectory, Hamiltonian mapping accurately reproduces K_d values from direct simulation. This method may be generally applied to systems beyond IDPs in force field optimization and in describing changes in thermodynamic behavior as a function of external conditions for connection with experiment.

Thermodynamic description of Hofmeister effects on the LCST of thermosensitive polymers

Cosolvent effects on protein or polymer collapse transitions are typically discussed in terms of a two-state free energy change that is strictly linear in cosolute concentration. Here we investigate in detail the nonlinear thermodynamic changes of the collapse transition occurring at the lower critical solution temperature (LCST) of the role-model polymer poly(N-isopropylacrylamide) [PNIPAM] induced by Hofmeister salts. First, we establish an equation, based on the second-order expansion of the two-state free energy in concentration and temperature space, which excellently fits the experimental LCST curves and enables us to directly extract the corresponding thermodynamic parameters. Linear free energy changes, grounded on generic excluded-volume mechanisms, are indeed found for strongly hydrated kosmotropes. In contrast, for weakly hydrated chaotropes, we find significant nonlinear changes related to higher order thermodynamic derivatives of the preferential interaction parameter between salts and polymer. The observed non-monotonic behavior of the LCST can then be understood from a not yet recognized sign change of the preferential interaction parameter with salt concentration. Finally, we find that solute partitioning models can possibly predict the linear free energy changes for the kosmotropes, but fail for chaotropes. Our findings cast strong doubt on their general applicability to protein unfolding transitions induced by chaotropes.

Propensity to form amyloid fibrils is encoded as excitations in the free energy landscape of monomeric proteins

Protein aggregation, linked to many of diseases, is initiated when monomers access rogue conformations that are poised to form amyloid fibrils. We show, using simulations of src SH3 domain, that mechanical force enhances the population of the aggregation-prone (N(⁎)) states, which are rarely populated under force free native conditions but are encoded in the spectrum of native fluctuations. The folding phase diagrams of SH3 as a function of denaturant concentration ([C]), mechanical force (f), and temperature exhibit an apparent two-state behavior, without revealing the presence of the elusive N(⁎) states. Interestingly, the phase boundaries separating the folded and unfolded states at all [C] and f fall on a master curve, which can be quantitatively described using an analogy to superconductors in a magnetic field. The free energy profiles as a function of the molecular extension (R), which are accessible in pulling experiments, (R), reveal the presence of a native-like N(⁎) with a disordered solvent-exposed amino-terminal β-strand. The structure of the N(⁎) state is identical with that found in Fyn SH3 by NMR dispersion experiments. We show that the timescale for fibril formation can be estimated from the population of the N(⁎) state, determined by the free energy gap separating the native structure and the N(⁎) state, a finding that can be used to assess fibril forming tendencies of proteins. The structures of the N(⁎) state are used to show that oligomer formation and likely route to fibrils occur by a domain-swap mechanism in SH3 domain.