Modelling the structure and interactions of intrinsically disordered proteins with multiple-replica, metadynamics-based sampling methods and force-field combinations
Intrinsically disordered proteins (IDPs) play a key role in many biological processes, including the formation of biomolecular condensates within cells. A detailed characterization of their configurational ensemble and structure-function paradigm is crucial for understanding their biological activity and for exploiting them as building blocks in material sciences. In this work, we incorporate bias-exchange metadynamics and parallel-tempering well-tempered metadynamics with CHARMM36m and CHARMM22* to explore the structural and thermodynamic characteristics of a short archetypal disordered sequence derived from a DEAD-box protein. The conformational landscapes emerging from our simulations are largely congruent across methods and forcefields. Nevertheless, differences in fine details emerge from varying forcefield/sampling method combinations. For this protein, our analysis identifies features that help to explain the low propensity of this sequence to undergo selfassociation in vitro, which can be common to all force-field/sampling method combinations. Overall, our work demonstrates the importance of using multiple forcefield/enhanced sampling method combinations for accurate structural and thermodynamic information in the study of general disordered proteins.
INTRODUCTION<!>RESULTS AND DISCUSSION
<p>Intrinsically disordered proteins (IDPs) are abundantly present in nature and play important roles in carrying out diverse biological functions, including cellular signaling, cytoplasmic compartment and regulation of gene expression [1][2][3][4][5][6][7][8] . A class of IDPs, defined as low complexity domains (LCDs), have been recently discovered in association with the dynamic formation of open compartments in cells [4][5][6][7][8] . These compartments, which are defined as condensates, are formed by liquid-liquid phase separation (LLPS) of proteins and nucleic acids and could underlie important functions and dysfunctions in biology. In addition to their role in biology, these LCDs represent promising building blocks for the design of synthesizing organelles to encode novel biochemical functionalities in bioengineering in a controllable and programmable manner, either alone or conjugated to soluble globular domains [9][10][11][12][13][14][15] .</p><p>IDPs can be characterized by the lack of stable well-defined native structures of folded proteins 16 , and amino acid compositions biased towards a high fraction of charged and polar residues, secondary-structure disrupters such as proline and glycine, and a low amount of bulky hydrophobic amino acids 1,[15][16][17][18][19][20] . It is believed the multivalent attractive interaction between side chains may give rise to favorable energetic gain that is responsible for counteracting the entropic loss during LLPS 15,[21][22] , and the phase behaviours of IDPs can be specifically encoded in their protein sequences [22][23][24][25][26][27][28][29][30][31][32] .</p><p>Nevertheless, our understanding of IDPs and their key role in mediating phase separations of multicomponent structures into coacervate assemblies is still largely limited. In the literature, mean-field theories such as Flory-Huggins theory and a recently emerging sticker-and-spacer model provide highly generalizable frameworks for describing the thermodynamics of the phase behavior of associate polymers in solution [33][34] , with emerging applications to complex biological systems both experimentally and computationally [34][35][36][37] . The thermodynamic driving force towards phase separation is likely to be inherently determined by the conformational free energy (FE) landscape of individual IDPs (intramolecular factor) and their molecular interactions with other sequences in solution including their own replicates (intermolecular factor). Thus, uncovering protein dynamics and their structure-function relationship can be pivotal for establishing versatile and sensitive protocols to design, program and predict the bottom-up assembly of multifunctional bio-inspired proteinbased materials based on IDPs.</p><p>Molecular dynamics (MD) simulations can in principle offer the opportunity of resolving structural, kinetic and thermodynamic characteristics of biological systems such as IDPs and IDP-rich bodies at the length and time scales that are not easily accessible experimentally 23,[38][39][40][41][42][43][44][45][46][47][48][49][50][51][52] . Recently, many protein folding problems were successfully solved by the advance of the deep learning method AlphaFold 53 . However, IDPs generally do not adopt a stably folded structure, featuring free energy landscapes containing a number of minima for different competing low-energy structures [16][17] . Thus, exploring the conformations of IDPs still largely rely on molecular dynamic simulations 16 . MD simulations complement state-of-the-art experimental methods such as NMR spectroscopy [54][55] , single-molecule Förster resonance energy transfer (sm-FRET) [56][57] , and small-angle X-ray and neutron scattering (SAXS and SANS) [58][59] , which may exhibit challenges in measuring molecular motions at atomic resolution and conformational heterogeneity associated with structural disorder 38,[42][43]50 . Nevertheless, conformational sampling in classical molecular dynamics (MD) simulations can be restrained to local minima of free energy surfaces (FES), and accessing the full and complex landscape of disordered proteins can be nontrivial with standard unbiased simulations. To mitigate the problem, various enhanced sampling schemes have been developed to help probing regions of the configurational space that could be rarely explored otherwise. Unbiased enhanced samplings are largely based upon the concept of parallel-tempering (PT) [51][52] . In temperature PT/replica-exchange algorithm, multiple replicas are simulated at different temperatures, in which conformations are exchanged at regular intervals, based on the Metropolis acceptance criterion. The stochastic nature of the method ensures generation of Boltzmann-weighted ensemble from which thermodynamic averages can be straightforwardly obtained. However, temperature differences between neighboring replicas must be moderate in order to yield practically large acceptance rates. For very large biomolecular systems simulated in explicit solvent, the number of replicas needed increases as O (f 1/2 ) with f as the system degrees of freedom, so reaching high temperature ensembles requires challenging computational costs 52,60 . In Replica Exchange with Solute Tempering (REST2) 38,41,46,[61][62] , such issue can be overcome by scaling the protein intramolecular potential energy, in which case the acceptance probability depends only on the differences of intramolecular solute energy and intermolecular solute-water energy. Recently, Replica exchange with hybrid tempering (REHT) 42 was developed to enhance exploring complex free energy landscape with large free energy barriers in the replica-exchange framework. At the same time, metadynamics-based techniques are powerful tools to accelerate the exploration of rare events such as protein folding and conformational transitions [39][40]23,[63][64][65][66][67][68][69][70][71][72] . These methods promote the exploration of free energy minima, by iteratively building a history-dependent bias potential as a sum of Gaussians defined as a function of a chosen set of collective variables (CVs). Boltzmann-weighted configuration ensemble can then be obtained via appropriate reweighting techniques 73-5 metadynamics (MetaD-WTE) allows the Gaussian height to be decreased with time so that the bias potential smoothly converges to the exact FE surface in the long-time limit [68][69] , while Metadynamics with Adaptive Gaussians can adapt the Gaussian width on the fly to the local features of FES in order to improve sampling efficiency 63 .</p><p>Here, we critically compare the application combinations of two metadynamicsbased, multiple-replica sampling methods and two protein force-fields to develop a systematic understanding of the conformational ensemble of an intrinsically disordered protein. The methods chosen are bias-exchange metadynamics (BEMD) [66][67] and parallel-tempering well-tempered metadynamics (PTMetaD-WTE) 39,47,[68][69][70] , which merge the advantages of PT method and metadynamics-based techniques. PTMetaD-WTE can significantly reduce the number of replicas required thanks to the potentialenergy bias introduced in the MetaD-WTE part; at the same time, the high temperature replicas from the PT part may compensate for the limited number of CVs directly biased to explore a high-dimensional phase space. On the other hand, in the BEMD method, a larger set of structural CVs is separately biased, and the system is often simulated at one temperature. Overall, the efficiency of the CV-based metadynamics methods depends on a suitable set of collective variables, which usually requires a priori knowledge of topological, chemical, and physical properties of the protein of interest.</p><p>In this context, we have selected structural CVs commonly used for exploring folded and disordered proteins 39,47,66 .</p><p>It is well-known that the results of MD simulations can strongly depend strongly on the accuracy of the applied protein-water force-fields, which may lead to large discrepancies with various experimental measurements 43,46,50,[76][77][78][79][80][81][82][83][84] . Examples have been observed with state-of-the-art Amber force-fields a99SB*-ILDN [85][86] with TIP3P   76 , the requirement of fourpoint water model significantly increases computational costs; the force-field is also reported to be too soluble for studying aggregations of some disordered proteins [91][92] .</p><p>Other important IDP-specific AMBER force fields include ff99IDPs [93][94] , ff14IDPs [94][95] and ff14IDPSFF 94,96 with the TIP3P water model 87 , which were optimized by incorporating different backbone torsion parameters from different amino acid groups including certain disorder-promoting amino acids. Here, we propose and demonstrate the use of two protein force-fields compatible with three-point CHARMM-modified TIP3P water, namely, CHARMM36m and CHARMM22*, in combination with our chosen sampling methods. In particular, CHARMM36m showed improved accuracy in generating polypeptide backbone conformational ensembles for intrinsically disordered proteins 79,97 , despite the issue of over-compact structures and over-stabilized helices 76,98 . The "helix-coil-balanced" CHARMM22*, instead, is not only able to reproduce experimental native-state structure and proteins folding rate, but also shows good agreement with experimental secondary structure propensities and NMR chemical shifts for many disordered proteins 76,89,98 .</p><p>In summary, we use BEMD and PTMetaD-WTE with CHARMM22* and CHARMM36m protein force-fields to effectively sample IDP conformational landscapes within a microsecond timescale. We apply this approach to elucidate the structural and thermodynamic properties of an archetypal IDP sequence derived from N-terminus of DEAD-box protein DHH1 [99][100] . DHH1N is a 46 amino-acid sequence that contains a low fraction of hydrophobic residues [Supplementary Note 1].</p><p>The peptide is enriched in polar residue such as Asparagine (Asn) and Threonine (Thr), negatively charged Aspartic acid (Asp), and positively charged large-sized Arginine (Arg) and Lysine (Lys), and with a few Prolines (Pro) distributed along the C-terminal half of the peptide. Experimental observations indicate that in vitro DHH1N does not undergo LLPS on its own at physiological pH, despite its composition is characteristic of low complexity domains involved in LLPS, and the fact that it participates in the formation of molecular adhesives to promote LLPS of chimera proteins 15 . In this work, we investigate a relatively short IDP that does not undergo phase separation, demonstrating how information from multiple sampling algorithms and forcefields reveal a detailed and congruent multi-dimensional FES within microsecond timescales.</p><p>This method represents a basis to next investigate the behavior of IDPs undergoing phase separation and connect structural properties with phase transition.</p><!><p>PTMetaD-WTE. To obtain a general understanding of the configurational ensemble of DHH1N using PTMetaD-WTE, we first calculated one-dimensional free energy (1D-FE) profiles as a function of individual collective variables at 300 K for both   To study more accurately the structural motifs emerging from the extensive sampling of DHH1N conformational ensemble, we also evaluate key secondarystructure content for CHARMM36m and CHARMM22* PTMetaD-WTE by means of the DSSP algorithm 106 . DHH1N is predominantly disordered, as demonstrated by the high number of loops/irregular elements, bends and turns for both CHARMM36m and CHARMM22*. When modelled with CHARMM36m, DHH1N displays three domains with a relatively high propensity to nucleate α-helices, located approximately along the sequences of Asn6-to-Asp16, Asp18-to-Asn26 and Thr36-to-Thr42 (Fig. 3c).</p><p>Interestingly, the propensity of α-helices along the N-terminal Asn6-to-Asp16 is comparatively low compared with those of the other two regions, which becomes negligible in the case of CHARMM22* PTMetaD-WTE (Fig. 3f) and unbiased simulation (Fig. 3i), suggesting that DHH1N simulated under CHARMM22* demonstrates negligible probability to form α-helices in the Asn-rich N-terminal domain.</p><p>Moreover, the sequence between Pro28 and Pro34, which separates the α-  Bias-exchange metadynamics. To complement the insights obtained with PTMetaD-WTE, we adopted BEMD method to further explore the conformational ensemble of DHH1N. We constructed 1D-FE profiles for CHARMM36m and CHARMM22* at 300  To further examine the conformational states of DHH1N, we construct 2D-FES for CHARMM36m and CHARMM22* BEMD and compare them with the 2D-FES from PTMetaD-WTE. Similarly, every simulation reveals a global minimum at low αcontent in the disordered cluster 1, accompanied by multiple shallow local minima, indicating a large number of substates as a result of conformational heterogeneity (Fig. 5). Cluster 1 of CHARMM36m and CHARMM22* BEMD comprises approximately 70 % of the equilibrium conformational population, similar to the result of CHARMM36m PTMetaD-WTE, while cluster 1 population of CHARMM22* PTMetaD-WTE is much lower due to relatively shorter simulation time (Supplementary Note 5, Table S3). The largest disparities between the 2D-FES are primarily due to different choices of force-fields. In detail, the CHARMM36m results</p><p>show that the nucleation of local α-helical domains is an activated process for both  Based upon the 2D free energy analysis, we proceeded to contact-maps and secondary-structures analysis for CHARMM36m and CHARMM22* of BEMD, in comparison with results from CHARMM36m and CHARMM22* of PTMetaD-WTE (Fig. 3). By focusing on residue-residue contacts with ~ 5 % frequencies, our analysis again reveals a sparse contact-map in which key contacts are mostly along the diagonal regions for all simulations. In respect of the off-diagonal regions, the four simulations show a lack of common residue-contacts. One exception is the Asn14-Ala24/Ser16-Lys22 domain of CHARMM36m BEMD. The contact-map of CHARMM36m BEMD displays prominent Asn14-Ala24 and Ser16-Lys22 contacts in this region, which is noticeably different from the other simulations within microsecond simulation length (Fig. 3i). At the same time, CHARMM36m BEMD explores a relatively lower number of α-helical structures, and a comparatively higher amount of β-strands in the domain of Asp18-to-Asn26, compared to the other method/force-field combinations (Fig. 3l).</p><p>To investigate possible relationship between the two features, we extracted structures with concurring Asn14-Ala24 and Ser16-Lys22 contacts from the CHARMM36m BEMD configurational ensemble and found that these residue-contacts are correlated with the presence of β-strands within Asn14-to-Asn26 (Supplementary Note 6, Figs.</p><p>S8, S9). The subgroup occupies ~ 18 % of the CHARMM36m BEMD population and spans both cluster 1 and 2 (Supplementary Note 6, Fig. S9b) [Fig. 5]. The relevant conformations, highlighted by the cyan circles in Fig. 5i, are less explored in the other force-field/method combinations (Figs. 5j, k and l).</p><p>We also constructed Cα -Cα contact-maps and secondary-structure analysis for clusters 1, 2 and 3 of all simulations (Fig. 6). For CHARMM36m BEMD, the concurring Asn14-Ala24 and Ser16-Lys22 contacts, which have been proven to be highly associated with the β-strands observed along Asn14-to-Asn26, emerge primarily in the less-ordered cluster 1 and 2 (Figs. 6j, k). In all cases, adjacent diagonal residuecontacts are focused within the two highest-order clusters (2 and 3), reflecting the presence of locally-ordered domains for the generally-disordered DHH1N. In contrast, the four force-field/method combinations share very few off-diagonal contacts, especially in the largest and most-disordered cluster 1, which contains about 70 % of overall population in most cases. A potential region of common contacts can be approximately between Pro28-Thr32 and Pro34-Leu40 (Fig. 6). The first sequence is capped by Pro28 and Pro34 and contains two adjacent units of positively-charged and bulky Lys29-Lys30), while the second sequence contains two adjacent units of negatively-charged Asp37-Asp38. Thus, electrostatic interactions may contribute to forming these residue-contacts, and it is also possible that participating in α-helical motifs reduces the likelihood of Pro34-Leu40 to interact with Pro28-Thr32 in cluster 2 and 3.</p><p>Overall, CHARMM36m BEMD shows a tendency to sample more β-sheetdisplaying conformations than the other force-field/method combinations. As shown earlier, a large fraction of the relevant structures corresponds to forming antiparallel βsheets along Asp18-to-Asn26 and α-helices along Thr36-to-Thr42, with concurring Asn14-Ala24 and Ser16-Lys22 residue contacts (Supplementary Note 5, Fig. S7; Supplementary Note 6, Figs. S8, S9). Such trend is not observed for CHARMM22* simulated under the BEMD method and CHARMM36m simulated under PTMetaD-WTE. One hypothesis could be that the feature could be specifically related with using the backbone-optimized force-field in combination with the sampling method. In the literature, it is also reported that BEMD explores a larger portions of phase space than simple temperature PT simulations 43 . Our evidence also shows that BEMD is able to sample a slightly larger conformational ensemble than PTMetaD-WTE (Figs. S7a and  b). Thus, it is possible that the observed discrepancy can be due to the fact that both αcontent (α-content) and β-content (β-content) were used as collective variables in BEMD, while the bias potential of PTMetaD-WTE (a function of Cα -Cα and Cγ -Cγ contacts) can only implicitly explore along the two secondary-structure reaction coordinates. Nevertheless, the difference may be relatively trivial in presence of the dominating disordered nature of DHH1N -over 70 % of the total population is contained in the low β-content (β-content ≤ 4) region for both CHARMM36m BEMD and PTWTE. Furthermore, all simulations provide a similar insight into a possible nucleation pathway for the partially disordered regions of DHH1N (schematically represented by the green dashed lines in Fig. 5). While the 2D-FES cannot reveal an accurate multi-dimensional pathway for the full conformational transition kinetics 109 , which could be a common challenge for representing protein free energies in lowdimensional space of popular CVs with physical meanings, our evidence indicates the dynamic disorder à a transition in the local domains of Asp18-to-Asn26 and Thr36to-Thr42 is likely to occur within a relatively compact DHH1N (Rg ~ 1.5 nm) with moderate sphericity (Fig. 5, Supplementary Note 5, Figs. S7g, h and i). Nevertheless, the high conformational heterogeneity of IDPs means that the exact locations and energies of the metastable states and barriers would require a more accurate free energy representation of reduced dimensionality. As a result of the designs and advantages of various protein force-fields and sampling methods, different portions of high-dimensional phase space may be preferably explored within simulation time, despite the high possibility that all results should converge in the long-time limit. Relying on the results from a single forcefield/sampling method simulation may restrain our understanding of the protein. For example, the temperature-dependence of certain protein properties can be easily recovered by analysing the sampling obtained from PTMetaD-WTE, which can provide key structure-property information to help understanding protein structures at different temperatures 110 ; on the other hand, BEMD simulations explore more collective variables and possibly higher portions of phase space within limited time, which can be preferred if certain metastable states can be better explored along certain CV directions.</p><p>Hence, in order to develop a comprehensive picture, it may be beneficial to include multiple FF/sampling methods combinations in the study of IDPs. For example, PTMetaD-WTE of multiple CV-combinations can also be conducted to compare the conformational space sampled within the same timeframe.</p><p>Our evidence shows that DHH1N has an average Rg between 1.54 and 1.64 nm (Table I), which is comparable with that of the other phase-separating IDPs that are much longer in length 23 . Such feature indicates a relatively low level of compactness.</p><p>In the literature, it is reported that the phase-separating behavior of disordered proteins can be generally associated with their single-chain compactness, and the sequencedeterminant of the compaction of disordered proteins is not only related to the overall protein charge, but also affected by the organization of charged and aromatic residues along the length 24,32,[111][112][113][114][115] . At physiological pH, DHH1N is electrostatically neutral, but its charged residues, i.e., negatively-charged Asp and positively-charged Lys, Arg, are mostly focused within the C-terminal portion of the sequence from Asp16 onwards.</p><p>However, DHH1N does not contain high-probability off-diagonal residue-contacts with that are shared by all force-field/method combinations, implying an absence of salt bridges or favorable electrostatic attractions between oppositely-charged patches of DHH1N, which could be key factors in inducing complex coacervation of many phaseseparating IDPs 21,23 . Scarce presence of aromatic amino acids in DHH1N also means that the protein is not like to form condensates stabilized by cation-π, and π-π interactions.</p><p>The comprehensive approach implemented in this work, which we have initially applied to analyse a relatively short IDP that does not undergo phase separation in vitro, represents a prospecting platform to investigate the conformational ensemble and thermodynamic driving force of other IDPs that exhibit liquid-liquid phase separation under a broad range of solution conditions. Such analysis will provide crucial information at the atomistic level which will be key to unravel biological phase separation in general and assist the design of new building blocks for advanced proteinbased materials and microreactors.  the public repository of the PLUMED consortium 126 . We used an exchange frequency of every 500 steps and the acceptance probabilities are between 8-20 %, with the higher acceptance rates mainly for exchanges between high-temperature replicas.</p><p>Bias-exchange metadynamics. BEMD 66-67 simulations were carried out adopting seven collective variables (the number of Cα -Cα contacts, Cγ -Cγ hydrophobic contacts and backbone H-bonds, dihedral correlation, α-content, antiparallel β-content and parallel β-content, as defined in Supplementary Note 2) and eight replicas, with one collective variable per replica plus the unbiased replica. One of the initial guess configurations was first equilibrated in NpT ensemble at 1 atm and 300 K for about 5 ns per replica in NpT ensemble according to the discussed protocol (vide supra); BEMD simulations were performed in NVT ensemble at 300 K, collecting 1000 ns of data for each replica; exchange of conformations between two randomly selected replicas were periodically attempted every 10000 simulation steps. CV were biased according to ordinary metadynamics scheme, adding the bias potential every 2500 simulation steps using a height value equal to 0.3 kJ mol -1 . After 160 ns, the system explored a wide region for each collective variable and we introduced loose lower and upper boundaries to improve convergence (Supplementary Note 8). The PLUMED input files required to reproduce the results of BEMD are available on PLUMED-NEST (www.plumednest.org), the public repository of the PLUMED consortium 126 .</p><p>Unbiased simulations. Unbiased simulations were performed adopting only CHARMM22* 88 force-field. Nine of the input guess structures were solvated and equilibrated at 1 atm and 300 K according to the discussed protocol (vide supra); 1000 ns MD simulations were subsequently performed in NVT ensemble at 300 K for each system. The nine set of data were concatenated for data analysis.</p>