Looking at the Disordered Proteins through the Computational Microscope.

Intrinsically disordered proteins (IDPs) have attracted wide interest over the past decade due to their surprising prevalence in the proteome and versatile roles in cell physiology and pathology. A large selection of IDPs has been identified as potential targets for therapeutic intervention. Characterizing the structure-function relationship of disordered proteins is therefore an essential but daunting task, as these proteins can adapt transient structure, necessitating a new paradigm for connecting structural disorder to function. Molecular simulation has emerged as a natural complement to experiments for atomic-level characterizations and mechanistic investigations of this intriguing class of proteins. The diverse range of length and time scales involved in IDP function requires performing simulations at multiple levels of resolution. In this Outlook, we focus on summarizing available simulation methods, along with a few interesting example applications. We also provide an outlook on how these simulation methods can be further improved in order to provide a more accurate description of IDP structure, binding, and assembly.

Introduction

Intrinsically disordered proteins (IDPs) are a class of proteins that are an exception to the widely accepted Anfinsen’s sequence-structure–function dogma, in which the amino acid sequence of a protein defines its unique three-dimensional fold that is inherently linked to the function. Instead, IDPs exist as a diverse ensemble of rapidly interconverting structures.[1] The structural disorder, unstructuredness, or fuzziness is abundant in eukaryotic proteome. About 25% of the eukaryotic proteins are predicted to be IDPs, and as high as 55% of eukaryotic proteins are predicted to consist of an intrinsically disordered region (IDR) longer than 30 residues.[2] IDPs, typically composed of 30 to 1000 residues, exhibit significant chain compaction, contrary to the typical random-coil-like behavior of chemically denatured proteins.[3] The function of an IDP often involves interaction with partners resulting in “fuzzy complexes”,[4] which encompasses molecular recognition, membraneless organelle formation, amyloidogenesis, cellular signaling, and regulation.[5] In addition, a number of IDPs are pathological targets in the age-associated amyloid disorders.[6] The known examples include Amyloid beta peptide (Aβ) and tau peptide in Alzheimer’s disease (AD), α-synuclein in Parkinson’s disease (PD) and synucleopathies, huntingtin protein in Huntington’s disease (HD), and islet amyloid polypeptide (hIAPP) in type-2 diabetes. Accumulation of the intracellular and extracellular protein aggregates, often with amyloid-like characteristics, is a pathological hallmark of those disorders.[7,8] The vast literature supports higher toxicity of the prefibrillar, early aggregate species compared to the mature fibril.[9−14] In healthy individuals, formation of those potentially toxic aggregates is prevented via cellular protection mechanism(s), such as molecular chaperones and protein degradation machinery.[15]

The detailed structural characterization of disordered proteins using conventional experimental techniques remains extremely challenging, due to their structural heterogeneity and dynamic nature. IDPs undergo conformational transitions at a time scale spanning several orders of magnitude (ps-ms), which is often difficult to address using a single experimental approach. Small-angle X-ray scattering (SAXS)[16] and single molecule Forster resonance energy transfer (FRET)[17] experiments have provided valuable information on the chain dimension, while NMR[18] has been useful in extracting secondary structure. However, the interpretation of those ensemble-averaged experimental observations has been hindered by the lack of a molecular model. Additionally, identification of transient structural epitope(s) within an IDP ensemble appears to be crucial, as those might have functional implication or disease relevance. Molecular simulations play an essential complementary role by allowing generation of sufficiently accurate structural ensembles that can directly interpret and/or predict experimental results. Simulations give mechanistic and dynamic insights into the IDP assembly as well as binding with partners. Additionally, simulations provide an in silico platform for designing and testing therapeutic strategies.

In this Outlook, we will summarize the existing computational approaches for IDP structure characterization along with some example applications, accentuate their strengths and limitations, and comment on future directions for improving those approaches. Computational approaches for generating IDP conformational ensembles can be essentially classified into two primary categories: (1) knowledge-based methods that use experimental data, mostly NMR, and (2) first principle or de novo approaches that combine empirical force fields with Molecular Dynamics or Monte Carlo simulations. Before getting into the details, it must be mentioned that these methods have been used for studying folded globular proteins for years. However, an extensive amount of research shows the need for modifications and improvements of existing tools, as well as calls for development of new methods.

Knowledge-Based Approaches

In the knowledge-based approaches, random coil ensembles are first generated using methods such as TraDES[19] and Flexible-Meccano.[20] The generated ensemble is then reweighted to match experimentally derived conformational constraints. Examples include the Monte Carlo optimization approach (e.g., ENSEMBLE[21] and ASTEROIDS[22]), Bayesian weighting approach,[23] and maximum entropy approach.[24] Albeit useful, lowering the generation bias toward statistical coil ensembles, improvements in back-calculation algorithms, and simultaneous optimization with respect to local, global, and dynamic experimental observables are some of the improvements needed for accurate identification of existing structures in an IDP ensemble.[25]

De Novo Simulations

• Explicit Solvent Simulations

The most accurate way to model a system is by including both protein and water at their atomistic detail. Different combinations of popular protein force fields (GROMOS,[26] OPLS/AA,[27] CHARMM,[28] AMBER[29]) and water models (TIP3P,[30] SPC,[31] SPC/E[32]) are commonly used for this purpose. However, the large spatiotemporal scales associated with IDP dynamics and aggregation often make the conventional atomistic explicit solvent simulations computationally unaffordable. Until recently, full-blown atomistic simulations in explicit water have been successfully used to study the structural landscape of a number of IDP sequences in their monomeric and small oligomeric form (for a review of recent work on this topic see ref (33) and references cited hereafter). The central idea behind conducting such simulations is to identify the intrinsic propensities of the IDP sequence and changes thereof due to disease mutations, environmental factors (solvent quality, proximity to lipid membranes, post-translational modification, etc.), which can affect function, aggregation, and toxicity. Thermodynamics and kinetics of the higher-order assembly of shorter amyloidogenic peptide fragments have been also investigated.[34−36] Simulations have been extensively used to investigate IDP interactions with various partners, for example, with small molecule,[37,38] peptide inhibitor,[39] membrane,[40,41] and chaperone.[42] Use of a higher than experimental concentration of the peptide and/or enhanced sampling protocols is becoming routine to ensure adequate sampling of an IDP’s conformational landscape. While standard molecular dynamics simulations can provide both structural and dynamical information on IDPs, sufficient convergence of the equilibrated structural ensembles can be more efficiently accomplished using enhanced sampling techniques, such as replica exchange molecular dynamics,[43,44] metadynamics,[45] and related methods, e.g. bias-exchange metadynamics,[46] solute tempering,[47−49] etc. Recently, parallel tempering in the well-tempered ensemble (PT-WTE) method has become available, which allows for temperature based sampling enhancement with reduced computational cost, as compared to standard REMD.[50,51] In the following, we will discuss examples of how simulations can help in characterizing IDP ensembles, both native and modified.

Transient Structure Identification

Transient local structure formation within monomeric IDP has been postulated to impact aggregation, loss-of-function, and toxicity, which can be efficiently investigated using de novo simulations. For example, atomistic simulations of the monomeric TDP-43 protein revealed local helix stabilization at the C-terminus in the presence of ALS mutations, consistent with NMR experiments.[52] Another simulation study showed that the known aggregation propensity of four amylin sequences, implicated in type-2 diabetes, can potentially be explained by their local helical propensity differences.[53]

A computationally well-studied IDP is Aβ peptide that is attractive due to its length (40–43 residues), disease-relevance (Alzheimer’s), and wealth of available experimental data on the wild-type (WT) and mutated sequences (Figure ). A number of familial mutations reported in Aβ sequence result in early onset of AD. Simulations from several groups have shed light onto the structural alterations resulting from those point mutations in the full-length and fragment Aβ, particularly within the monomeric and dimeric landscape. In the following, we will review in detail the differential effect of two recently reported N-terminal mutations, one AD-causative and one protective, on the monomeric and dimeric Aβ landscapes (Figure ). The A2V mutation is associated with the early onset of AD in homozygous carriers, whereas it offers some protection in the heterozygous state.[55] The A2T mutation reportedly has an AD-protective effect.[56] Correspondingly, a number of atomistic REMD studies using different combination of protein force fields and water models have reported differences in the monomeric configuration between the wild-type, A2V, and A2T variants of full-length Aβ peptide as well as of Aβ fragments.[54,57,58] Although all three variant monomers were found to form collapsed, disordered ensembles in water (Figure , bottom), they were strikingly different in terms of transient secondary and tertiary structure formation (Figure ). Simulations indicate the transient presence of a β-hairpin motif within the WT Aβ monomer, which is comprised of the central and C-terminal hydrophobic residues (Figure , magenta rectangle),[54,59−61] similar to the β-hairpin structure found in the complex of Aβ40 with an affibody protein ZAβ3.[62] The N-terminus forms limited interaction with the rest of the peptide in those short-lived β-hairpin structures. Interestingly, an enhanced population of a double β-hairpin motif involving the NTR has been observed in simulations of a A2V monomer,[54,57] which contains a strong resemblance to a solid-state NMR model of neurotoxic Aβ42 oligomers[63] (Figure , red rectangle). In contrast, the A2T mutation enhances disorder within the Aβ42 ensemble due to atypical long-range interactions between the NTR, CHC, and 22–28 turn region (Figure , green rectangle).[54,58] The effect of homotypic and heterotypic interactions between WT, A2V, and A2T monomers on the Aβ40/42 dimeric landscape has been also studied using atomistic simulations.[64−68] While aggregation-prone WT/WT and A2V/A2V homodimers are found to populate highly stabilized interpeptide structures with a high content of interpeptide β-sheets or β-hairpin conformations, the AD-protective WT/A2V and WT/A2T Aβ42 heterotypical interactions favored formation of weak dimers that are highly disordered and collapsed in nature.[67−69] These results suggest that A2V and A2T mutations differentially alter the misfolding into the transient monomeric β-hairpin structures that are thought to play a crucial role in aggregation, as suggested by earlier simulations[70,71] and experiments.[72] The overly stabilized double β-hairpin structures might explain the longer lag phase as well as the atypical spherical oligomer population during A2V Aβ42 aggregation,[73] which might be linked to its disease-causing effect in homozygous individuals. β-Hairpin lowering in the A2T Aβ42 ensemble is interesting, given its relatively lower long-term potentiation (LTP) inhibition ability observed in rat CA1 hippocampal CA1 neurons[73] and the reported toxicity relevance of a monomeric β-hairpin structure.[72] The weak dimer formation observed in protective WT/A2V and WT/A2T mixture simulations agrees with the oligomeric properties reported in FT-IR and ANS fluorescence experiments.[74] Additionally, the enhanced disorder in the heterodimeric population might explain the absence of toxic dodecamer in the protective WT/A2T(A2V) mixture, as revealed in IM-MS experiments.[75] Taken together, existing simulations imply that the mutated Aβ N-terminus can modulate an aggregation-prone, toxic β-hairpin population in a sequence-dependent manner.

Figure 1

Sequence (top) and simulated ensemble (bottom) of the wild-type (WT), causative A2V, and protective A2T Aβ42 variant monomer.[54]

Figure 2

Changes in the Aβ42 conformational landscape due to A2V and A2T mutation.[54]

Toward Structure-Based Therapeutic Design

Designing Aβ-targeting small molecules,[77−79] peptide fragments,[80,39,81] or antibodies[82] is an active area of AD research due to the widely accepted amyloid cascade hypothesis. It has been suggested that a better understanding of structure-toxicity relationships, much earlier detection of AD, and drugs that consider multiple targets (Aβ, tau) and can pass blood-brain barrier are needed for effective drug design.[83] Atomistic simulations are constantly being employed in designing mechanism-based antiamyloid molecules against monomeric, oligomeric, and fibrillated forms of Aβ. The atypical interactions triggered by an altered N-terminus, as seen within monomeric and dimeric Aβ, call for further simulations toward Aβ ΝΤR-based therapeutic design.

In fact, recent experiments have reported aggregation and/or toxicity inhibiting property of Aβ N-terminal fragments.[84−86] Experiments also showed that linking NTR Aβ fragments with an arginine-rich TAT sequence can facilitate their delivery to brain and cell membrane translocation.[84] The 1-6A2VAβ-TAT molecule exhibited robust antiamyloidogenic properties in vitro and protected transgenic C. elegans and human neuroblastoma cells from Aβ toxicity.[84] Preclinical studies in AD mouse models indicated that short-term treatment with 1-6A2VAβ-TAT impedes Aβ aggregation and cerebral amyloid deposition.[84] Inspired by those experiments, the interaction of WT, 1-6A2V, and 1-6A2T Aβ hexapeptide with the full-length Aβ42 monomer and its effect on the monomeric Aβ42 landscape have been recently investigated using atomistic simulations[76] (Figure ). While all three hexapeptide bindings led to a robust lowering of transient WT hairpin population in simulations, interesting alternative structural features emerged. A2V-bound Aβ42 preferably visited a disordered population (Figure a), while a central helix was more frequently populated within the A2T-bound monomer (Figure b). Such redirection of an IDP to atypical disordered conformations induced by an inhibitor/drug interaction may have implications in the lowering of IDP aggregation and toxicity in general, as found in experiments for α-synuclein, tau protein, and Aβ.[79,77,87,88] Taken together, it is becoming more evident that synergistic simulation-experimental investigations are going to play a crucial role in rational design of mechanism-based therapeutics for IDP-associated diseases in the forthcoming days.

Figure 3

Aβ N-terminal hexapeptide binding stabilizes alternative Aβ42 structures.[76] (a) Unstructured Aβ42 bound to 1-6A2V and (b) helical Aβ42 bound to 1-6A2T.

• Implicit Solvent Models

A cheaper alternative way to simulate IDP conformational dynamics is by representing solvent in an implicit manner. In implicit solvent models, the solute–solvent interaction is presented using a mean field approach, in which the solvation free energy is represented as a function of the atomic coordinates of the solute. Approximating the solvation free energy relies on the solvent-accessible surface area (SASA) models or the solvent contact models. The accurate and efficient formulation of the solvent effect remains a nontrivial theoretical challenge, and detailed discussion of this topic is beyond the scope of this Outlook. Apparently, many implicit solvent models suffer from the same problems as the explicit solvent models do, such as overly collapsed disordered states.[89] Implicit solvent models that appear promising for IDP simulations are based on the EEF1 model originally developed by Lazaridis and Karplus.[90] The EEF1 model decouples and independently estimates the contributions of the charged and nonpolar groups toward the solvation free energy of the solute. A conceptually similar, but much superior, model is the ABSINTH force field of Vitalis and Pappu[91] that additionally leverages the strengths of the Poisson framework and better represents polyelectrolyte and polyampholyte IDPs. The ABSINTH model has been extensively and successfully used in recent years to delineate the sequence-structure–function paradigm of IDPs.[92,93] For example, ABSINTH simulations combined with experiments have recently uncovered how charge patterning in the “RAM” disordered region of the Notch receptor can modulate its transcriptional activity.[94]

• Coarse-Grained Models

Coarse-grained (CG) models further decrease computational requirements (usually by ∼100× compared to fully atomistic explicit solvent simulations) by eliminating few degrees of freedom of the peptide and thus allow simulations of larger length- or time scales. These models are sometimes simulated with the Discrete Molecular Dynamics (DMD) method that involves discontinuous potentials, in order to further gain computational efficiency.[95] CG protein models can assume various levels of reduced representation, from a single monomer (peptide unit), to an intermediate level (multiple beads per residue), to a single residue (single bead per residue) resolution. The most coarse-grained are the lattice models in which a single residue occupies a single site of the lattice. While the lattice models allow efficient sampling of the complete phase diagram, due to their inability to capture realistic kinetics and structures they can only be used to answer very generic questions. Therefore, coarse-grained protein models are mostly off-lattice in nature. The interaction parameters in those models are derived either by matching with the atomistic simulations and/or with experimental data, as in the MARTINI model,[96] or by using a knowledge-based approach, as in PRIME20.[97] The lack of granularity in CG models leads to obvious limitations, such as failure to correctly estimate secondary structure. For example, one limitation of the widely used MARTINI model is the use of secondary and tertiary structure constraints, which prevents exploration of protein conformational fluctuations. It has been shown that inclusion of peptide backbone polarization effects within the MARTINI model framework can alleviate that limitation.[98,99] A combination of physics-based and knowledge-based potentials has been also used, such as in OPEP.[100] AWSEM (associated memory, water mediated, structure and energy model) developed by Wolynes, Papoian, and co-workers is another CG model that is comprised of both physics-based and knowledge-based potentials.[101] Over the years, these CG models have provided valuable insights onto the mechanism of peptide assembly. For example, an AWSEM simulation study of Aβ40 aggregation has characterized the “amyloid funnel” that reveals a prefibrillar to fibrillar interconversion bottleneck at the pentamer level.[102] Many coarse-grained models are built on adaptations of the so-called Go-like model that uses a native-centric potential.[103] Further reparametrization and inclusion of additional terms in the energy function were performed to better suit the task of IDP structure modeling, e.g. fine-tuning the local structural propensity and accounting for the non-native and electrostatic interactions. Inclusion of those essential factors in the CG model has allowed dissecting the mechanism of coupled folding and binding of disordered proteins, which is highly valuable, as experimental characterization as well as atomistic simulation of IDP binding is nontrivial. One important observation is formation of a transient encounter complex at the initial stage of IDP binding with a partner, which is driven by non-native, intermolecular interactions, often electrostatic in nature.[104−106]

While CG models allow exploration of dimensions, time scales, and rare events that are inaccessible to atomistic simulations, they often lack transferability and fail to accurately reproduce experimental observables (structure, thermodynamics, kinetics, solvation effect), when compared to all-atom models. Further development, refinement, and optimization of CG models to make them more reliable and predictive in the context of IDP simulations is an active area of research. For example, a CG modeling framework at a residue-level resolution has been recently introduced to study the long time scale liquid–liquid phase separation process of IDPs.[107]

Outlook

The emerging synergy between simulations of various resolution and experiments has established simulations as a promising tool for delineating the sequence-structure-(mis)function relationship of disordered proteins. The structural transience and/or heterogeneity of soluble oligomers revealed in simulations further encourages use of computational modeling toward developing strategies for rational manipulation of IDP assembly. Nevertheless, the need for reassessment and reparametrization of existing all-atom protein force fields and water models is becoming increasingly evident, as those parameters were originally purposed for reproducing the properties of the folded, globular proteins. For example, two separate studies (refs (58) and (108)) compared multiple modern force fields and reported comparable performance in characterizing Aβ conformations. All force fields, despite their fundamental differences, yielded similar tertiary structure sampling as well as adequate agreement with NMR data such as J-coupling. While convergence among different force fields at accurately representing the IDP conformational ensemble is encouraging, more work needs to be done to further fine-tune force fields or update solvent models[109] for maintaining accurate secondary structure balance and chain dimension. The crucial role of hydration water and hydrodynamic interactions in IDP misfolding and aggregation is evident.[110,111] One effort in this direction involves designing force field(s) optimized for right secondary structure balance (AMBER ff03*)[112] to capturing temperature cooperativity in protein folding (AMBER03w)[113] to reproducing unfolded protein dimensions consistent with experiment (AMBERws).[114] The key difference between ff03* and ff03w is the use of the TIP4P/2005 water model in the latter, which is expected to provide a better representation of solvation effects such as temperature and pressure dependence. In ff0ws, protein–water interactions are strengthened by 10% from standard mixing rules to achieve the right balance of cross-interactions. The improvements in the resulting force field highlight the important role of capturing water-mediated interactions appropriately in a biomolecular force field. Recent development of the CHARMM36m[115] force field is another notable effort in this direction, which involves refinement of the CHARMM36 protein force field, mainly by using a refined backbone CMAP potential and improved description of salt-bridge interaction, in order to correct the overestimation of left-handed α-helix population. This study also highlights the role of protein–water dispersion interactions in correcting the general bias of physics-based atomistic models toward overly compact IDP ensembles.

Future force field development research should also focus on accurate modeling of the highly dynamic metal-IDP interactions, due to their pivotal roles in neurodegenerative diseases.[116] For example, metal coordination is known to alter the aggregation kinetics and morphology of Aβ aggregates.[117] Metal ions in unpolarized atomistic simulations are treated as simple van der Waals spheres, nonbonded models with atoms, or with bonded models by employing artificial bonds.[118] Nevertheless, the polarization effects are crucial, as the polarization energy can range from very small values to values around 10–20% for the total interaction energy.[119] Therefore, a number of polarizable modeling approaches are continuously being developed, such as fluctuating charge, Drude oscillator, and the induced dipole models.[118] However, the usage of polarizable force fields increases the computational requirements compared to unpolarized approaches.

The development of more accurate force fields will benefit from an extensive, open-source database comprising simulated ensembles and experimental data of various functional and toxic IDPs. To achieve such a platform, a crowd-sourcing approach is likely needed. This will also help tackling the “Big Data” problem, as with increasing computational power,[120,121] emerging efficient sampling schemes,[122] and wide interest in protein simulations, the cumulative amount of simulation data is becoming massive. How to merge, analyze, and build predictive models from this massive amount of data is becoming a research challenge.[123] Application of machine learning and statistical methods will likely play a critical role in tackling this big data challenge.

The open simulation data platform can also be helpful in developing novel coarse-grained models that are transferable (e.g., globular proteins to disordered proteins, in solution to on membrane) and more accurate. Such CG models can bias and guide all-atom simulations for further validation, establish complementary tools to experiments, and improve our current understanding of IDP behavior in cellular environments. For example, explicit representation of protein hydration water should be considered in CG models,[98,124] as hydration water surrounding IDPs appears crucial for determining their structure and function.[125,126] In addition, both atomistic and CG models should be calibrated as well as combined with various, often fragmentary, experimental data, for enabling deeper interpretation.

In order to obtain the full mechanistic picture of IDP function, the peptides need to be simulated in a crowded, cell-like environment. In addition, the short time scale structural dynamics should be combined with the assembly/binding occurring at the larger spatiotemporal scale (Figure ). As shown in Figure , such investigations may shed light onto the relation between functional droplets and amyloid-related toxicity. While the computational efficiency of coarse-grained models makes them an attractive choice for this purpose, they need to be carefully combined with higher resolution models, highlighting the need for constructing novel multiscale modeling strategies. One way to achieve that goal is by allowing smooth, on-the-fly exchange of particles between different resolutions.[127,128]

Figure 4

A schematic diagram depicting potential states populated on the IDP assembly landscape.