<?xml-model href='http://www.tei-c.org/release/xml/tei/custom/schema/relaxng/tei_all.rng' schematypens='http://relaxng.org/ns/structure/1.0'?><TEI xmlns="http://www.tei-c.org/ns/1.0">
	<teiHeader>
		<fileDesc>
			<titleStmt><title level='a'>Quantitative prediction of ensemble dynamics, shapes and contact propensities of intrinsically disordered proteins</title></titleStmt>
			<publicationStmt>
				<publisher></publisher>
				<date>09/09/2022</date>
			</publicationStmt>
			<sourceDesc>
				<bibl> 
					<idno type="par_id">10414931</idno>
					<idno type="doi">10.1371/journal.pcbi.1010036</idno>
					<title level='j'>PLOS Computational Biology</title>
<idno>1553-7358</idno>
<biblScope unit="volume">18</biblScope>
<biblScope unit="issue">9</biblScope>					

					<author>Lei Yu</author><author>Rafael Brüschweiler</author><author>Bert L. de Groot</author>
				</bibl>
			</sourceDesc>
		</fileDesc>
		<profileDesc>
			<abstract><ab><![CDATA[Intrinsically disordered proteins (IDPs) are highly dynamic systems that play an important role in cell signaling processes and their misfunction often causes human disease. Proper understanding of IDP function not only requires the realistic characterization of their three-dimensional conformational ensembles at atomic-level resolution but also of the time scales of interconversion between their conformational substates. Large sets of experimental data are often used in combination with molecular modeling to restrain or bias models to improve agreement with experiment. It is shown here for the N-terminal transactivation domain of p53 (p53TAD) and Pup, which are two IDPs that fold upon binding to their targets, how the latest advancements in molecular dynamics (MD) simulations methodology produces native conformational ensembles by combining replica exchange with series of microsecond MD simulations. They closely reproduce experimental data at the global conformational ensemble level, in terms of the distribution properties of the radius of gyration tensor, and at the local level, in terms of NMR properties including              15              N spin relaxation, without the need for reweighting. Further inspection revealed that 10–20% of the individual MD trajectories display the formation of secondary structures not observed in the experimental NMR data. The IDP ensembles were analyzed by graph theory to identify dominant inter-residue contact clusters and characteristic amino-acid contact propensities. These findings indicate that modern MD force fields with residue-specific backbone potentials can produce highly realistic IDP ensembles sampling a hierarchy of nano- and picosecond time scales providing new insights into their biological function.]]></ab></abstract>
		</profileDesc>
	</teiHeader>
	<text><body xmlns="http://www.tei-c.org/ns/1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xlink="http://www.w3.org/1999/xlink">
<div xmlns="http://www.tei-c.org/ns/1.0"><head>INTRODUCTION</head><p>Intrinsically disordered proteins (IDPs) and protein regions (IDRs) are an integral part of the proteomes of many different organisms with more than 30% of all eukaryotic proteins possessing 40 or more consecutive disordered residues.(1, 2) While IDPs and IDRs in isolation do not adopt well-defined three-dimensional (3D) structures, they often play important biological roles in molecular recognition processes by interacting in specific ways with binding partners that are typically well-ordered. <ref type="bibr">(3)</ref><ref type="bibr">(4)</ref><ref type="bibr">(5)</ref> For instance, the human oncoprotein protein p53 possesses the Nterminal transactivation domain (p53TAD) that binds to the N-terminal domain of human MDM2 protein adopting a stable a-helix.( <ref type="formula">6</ref>) Prokaryotic ubiquitin-like protein (Pup) is another IDP that is directly linked to protein degradation folding into an a-helix when binding to Mpa protein. <ref type="bibr">(7)</ref> In addition to binding to their target protein(s), IDPs can also be involved in liquidliquid phase separation (LLPS). <ref type="bibr">(8)</ref><ref type="bibr">(9)</ref><ref type="bibr">(10)</ref><ref type="bibr">(11)</ref> LLPS is the segregation of molecules in solution into a condensed phase and a dilute phase with high and low biomolecular concentrations. These membraneless droplet-like compartments formed by IDPs and other biomolecules are important for cellular function. Knowledge of the structural and dynamic propensities of IDPs both in isolation and in complex biological environments is essential for understanding these processes and their role in human diseases.</p><p>In order to relate IDP sequences to biological function, detailed knowledge of IDP conformational ensembles is needed. The description of conformational ensembles can range from local secondary structure populations to explicit ensembles in 3D space with atomic resolution.(12) Some of the earliest approaches generate random coil conformational ensembles that are subsequently refined against a host of experimental data reflecting both local and global structural features. <ref type="bibr">(13)</ref><ref type="bibr">(14)</ref><ref type="bibr">(15)</ref> These approaches continue to be successfully applied through integrative modeling provided that a large amount of high quality experimental data is available for each system under investigation. <ref type="bibr">(16,</ref><ref type="bibr">17)</ref> Even when data from various complementary experimental techniques are being used, the amount of experimental information obtainable is still sparse when compared to the information needed to uniquely characterize large, highly heterogeneous structural ensembles that are the hallmark of IDPs. As a consequence, the amount of information that can be gained and that is not directly reflected in the experimental data used to refine the ensemble is restricted to robust descriptors ranging from coarse-grained to global that can be compared with predictions by polymer theory under various assumptions. <ref type="bibr">(16)</ref> In addition, site-specific interaction information, such as transient inter-residue contacts, can be obtained at medium to low resolution from paramagnetic relaxation enhancement <ref type="bibr">(PRE)</ref> experiments by attaching electron spin labels to selected sites. <ref type="bibr">(15,</ref><ref type="bibr">17)</ref> Because empirical ensembles generated based on such data lack a time axis, they do not include dynamics time scales of IDPs associated with interconversion rates between substates and, hence, they do not inform about an essential part of the energy landscape.</p><p>From a theoretical and computational perspective, all-atom molecular dynamics <ref type="bibr">(MD)</ref> simulations are an attractive alternative to empirical approaches for the generation of IDP conformational ensembles, including dynamic time scale information, for the comprehensive interpretation of experimental results. <ref type="bibr">(18)</ref> However, for many years limitations in computer power precluded the generation of statistically well-converged results and MD force fields primarily developed for ordered proteins turned out to be unsatisfactory for applications to IDPs.</p><p>With the continuing increase in computer power, the quality of sampling has reached a level that allows rigorous validation by quantitative comparison with a rich body of experimental data. In cases where discrepancies are observed between simulation and experiment, as is commonly the case, approaches have been developed that use restraining or reweighting that bias the original simulation to obtain results that agree better with experimental data. <ref type="bibr">(19)</ref><ref type="bibr">(20)</ref><ref type="bibr">(21)</ref><ref type="bibr">(22)</ref><ref type="bibr">(23)</ref><ref type="bibr">(24)</ref><ref type="bibr">(25)</ref><ref type="bibr">(26)</ref> When not only the conformational ensemble but also the underlying dynamics time scales are of interest, suitable rescaling of the MD time step or correlation times of the dominant motional modes can be applied to improve agreement with experiment. <ref type="bibr">(27)</ref><ref type="bibr">(28)</ref><ref type="bibr">(29)</ref><ref type="bibr">(30)</ref> Because these methods can often improve the unaltered simulations only within certain boundaries, they are best suited when the original predictions are fairly close to experimental data. <ref type="bibr">(31)</ref> Although these methods rarely fail to produce better agreement, at least on average for those experimental parameters directly used as restraints or for reweighting, they naturally depend on large amounts of experimental data of good quality as input for each protein system studied. This amounts to a laborious experimental effort that needs to be repeated for each new protein system as the experimental data are proteinspecific rendering them non-transferrable between systems.</p><p>An alternative and more principled approach is to improve the MD force fields themselves enabling them to increasingly accurately predict experimental data in a way that is fully transferrable between protein systems, both ordered and disordered. This premise has led to a recent proliferation of protein force field developments <ref type="bibr">(32)</ref><ref type="bibr">(33)</ref><ref type="bibr">(34)</ref><ref type="bibr">(35)</ref><ref type="bibr">(36)</ref><ref type="bibr">(37)</ref> and new explicit water models <ref type="bibr">(38)</ref><ref type="bibr">(39)</ref><ref type="bibr">(40)</ref> specifically geared toward the improved representation of disordered proteins. In a significant development, residue-specific force fields have been introduced.(41) These force fields use in addition coil library information from the Protein Data Bank (PDB) by incorporating the individual backbone j,y propensities of each residue type. <ref type="bibr">(41)</ref><ref type="bibr">(42)</ref><ref type="bibr">(43)</ref><ref type="bibr">(44)</ref><ref type="bibr">(45)</ref><ref type="bibr">(46)</ref><ref type="bibr">(47)</ref> Such residue-specific force fields, in combination with suitable water models, can provide an improved representation of disordered states while retaining the properties of ordered proteins.</p><p>With respect to water models, TIP4P-D and closely related derivatives have been notably successful in preventing overly compact conformations by favoring more extended IDP structures showing improved agreement with experiment. <ref type="bibr">(38)</ref> Besides global properties, such as the radii of gyration and asphericities, IDP ensembles and trajectories should also accurately reproduce local dihedral angle distributions and secondary structure propensities. Moreover, they should also replicate dynamic and kinetic IDP properties, such as librational motions and time scales of interconversion between conformational substates. Such information is important for understanding recognition events between IDPs and their binding targets, including IDP interactions with other disordered biomolecules, for example, during the formation of LLPS condensates. Experimental IDP dynamics information can be gained from fluorescence depolarization spectroscopy, <ref type="bibr">(48)</ref> F&#246;rster resonance energy transfer (FRET), <ref type="bibr">(16)</ref> and nuclear magnetic resonance (NMR) relaxation.(15) NMR 15 N longitudinal R1 and transverse R2 spin relaxation rates are exquisitely sensitive to the dynamics of disordered proteins and the underlying time scales.(49-51) R2 relaxation rates, for example, have been linked to residual intramolecular interactions in chemically unfolded proteins.(51-53) 15 N R1 and R2 rates can be experimentally determined for each protein residue and therefore they are valuable for validating MD simulations with respect to amplitudes and time scales of IDP dynamics. <ref type="bibr">(29,</ref><ref type="bibr">(54)</ref><ref type="bibr">(55)</ref><ref type="bibr">(56)</ref> We recently developed the AMBER ff99SBnmr2 force field by modifying the backbone dihedral angle potentials of each amino-acid residue type to reproduce the j,y dihedral angle distributions found in a random coil library. <ref type="bibr">(57)</ref> The ff99SBnmr2 force field has been validated against experimental nuclear magnetic resonance (NMR) scalar 3 J-couplings of a-synuclein and b-amyloid IDPs demonstrating that this force field accurately reproduces their sequencedependent local backbone structural propensities. <ref type="bibr">(58)</ref> The primary goal of this work is to learn whether state-of-the-art replica exchange and extended MD simulations of IDPs can also realistically reproduce NMR R1, R2 relaxation rates with their strong and unique dependence on motional time scales without the need of any additional corrections such as constraints or reweighting. Moreover, in-depth analysis of the MD trajectories generated yields a wealth of information about the radius of gyration tensor distribution and dominant dynamics modes allowing graph-theory based identification of specific inter-residue interaction propensities and residue clusters for the better understanding of IDP behavior.</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>RESULTS</head><p>Ensemble properties of radius of gyration tensor. The radius of gyration Rg(t) is shown as a function of time for representative 1-&#181;s MD trajectories of p53TAD and Pup in Fig. <ref type="figure">1A,</ref><ref type="figure">B</ref> (see also Fig. <ref type="figure">S1</ref>). The trajectories exhibit predominantly stationary stochastic behavior reflecting random expansion and contraction of the overall IDP size with the mean value (blue horizontal lines) in good agreement with the experimentally determined &lt;Rg&gt; (black line) or the predicted &lt;Rg&gt; from polymer theory (Eq. 6). The MD-distributions of Rg of all 10 MD trajectories are shown as histograms in Fig. <ref type="figure">1C,</ref><ref type="figure">D</ref>. The Flory exponent n of the polymer scaling law was determined from the REMD ensembles at 298 K. Using r0 = 1.927 &#197;, we obtain a value of n = 0.601 for Pup, which closely matches the theoretical value ntheory = 0.588 of a fully disordered, self-avoiding random coil. <ref type="bibr">(59,</ref><ref type="bibr">60)</ref> For p53TAD, the REMD &lt;Rg&gt; value of 28.1 &#197; is in almost perfect agreement with experiment(61) (28.0 &#197;) corresponding to n = 0.624, which clearly exceeds ntheory. The characteristic time scales of Rg(t) fluctuations can be obtained from the timecorrelation functions CRg(t) (Eq. 5), which are well-converged over the course of the 1-&#181;s trajectories (Fig. <ref type="figure">1E,</ref><ref type="figure">F</ref>). CRg(t) of both proteins decay in good approximation biexponentially with reconfigurational correlation times ta &#8773; 10 ns and tb &#8773; 55 ns. The normalized variance of the Rg(t) fluctuations, given by <ref type="bibr">(1)</ref> is almost the same for p53TAD (0.03) and Pup (0.04). The ensemble distribution of the gyration tensor S (Eq. 2) contains information about the deviation of individual MD snapshots from spherical shape, which can be directly compared with a random Gaussian chain serving as a perfect random coil (Fig. <ref type="figure">2</ref>).(62) Both proteins show unimodal asphericity distributions (Eq. 3) with maxima around A &#8773; 0.18, which qualitatively differ from the Gaussian chain model (Fig. <ref type="figure">2C</ref>) peaking at A = 0. Compared to p53TAD, Pup has a higher tendency to adopt a more spherical conformation. Another useful measure of the overall shape of individual snapshots is the prolateness P (Eq. 4). The distribution of P is bimodal for both proteins with the global maximum corresponding to prolate-shaped (cigar-like) structures (P = 1) and a second (local) maximum corresponding to disk-like structures (P = -1). The distribution of the prolateness of Pup is more balanced between positive and negative values with &lt;P&gt; = 0.2 than for p53TAD, which has a higher tendency to adopt prolate-shaped conformers (&lt;P&gt; = 0.35), whereas the Gaussian chain distribution (&lt;P&gt; = 0.3) lies between the two IDP distributions. The distinct asphericity distribution and increased prolateness of p53TAD is at the origin of its increased &lt;Rg&gt; over the Gaussian random coil model. agreement with experiment(63) evidenced by small RMSEs (0.10 s -1 for p53TAD and 0.12 s -1 for Pup) and Pearson correlation coefficients R of 0.78 for p53TAD and 0.86 for Pup (Fig. <ref type="figure">3A,</ref><ref type="figure">B</ref>).</p><p>R2 relaxation rates determined from the simulations are also in good agreement with experiment with correlation coefficients R of 0.88 for p53TAD and 0.70 for Pup and RMSEs of 0.84 s -1 for p53TAD and 0.81 s -1 for Pup and (Fig. <ref type="figure">3C,</ref><ref type="figure">D</ref>). It can be seen that the simulations tend to underestimate R1 and overestimate R2 rates, although only slightly, in a manner that is notably uniform for the R1 values of both proteins and for the R2 values of p53TAD. The 10 N-terminal residues of p53TAD are very flexible with small R2's, which closely follow the experiment. For Pup, differences in R2 between MD and experiment display the same trend and are most pronounced for residues 30-48. The error bars of the computed relaxation rates, which represent the root-mean-square deviations over all 10 MD trajectories, are fairly uniform along the polypeptide chains and systematically larger for R2 than for R1, again with the exception of the 10 N-terminal residues of p53TAD. For both proteins, not all 10 1-&#181;s MD trajectories individually reproduce the experimental data equally well. Either 1 (p53TAD) or 2 (Pup)</p><p>trajectories have more compact average IDP structures, which quantitatively affect the agreement with experiment (Fig. <ref type="figure">S2</ref>).</p><p>Correlation times of backbone N-H bond vectors in both proteins fitted from the average correlation functions range from picoseconds to about 20 ns (Fig. <ref type="figure">3E,</ref><ref type="figure">F</ref>). Consistent with the finding for other IDPs, <ref type="bibr">(55,</ref><ref type="bibr">64)</ref> the dominant contribution to the time correlation functions stems from dynamics on the intermediate time scale around 1 ns reporting about backbone j,y jumps.</p><p>Fast dynamics on the time scale of 100 ps or faster report on local 15 N-1 H bond librations, similar to those observed in secondary structures of folded proteins,( <ref type="formula">65</ref>) and slower dynamics on the time scale between 3 and 20 ns reports on collective IDP chain motions. The presence of slower modes correlate with increased R2 values most pronounced for residues 30-48 in Pup. This is consistent with relaxation theory (Eq. 12), which predicts that in solution transverse spin relaxation rates R2 are in good approximation proportional to the effective overall correlation time experienced by the 15 N-1 H spin pairs.   Inter-residue contact propensities. Different residues along the polypeptide chain display different tendencies to form contacts with other residues. Fig. <ref type="figure">5A</ref>,B shows the average number of contacts per snapshot for each residue, which was calculated as the total number of contacts formed by a residue divided by the total number of MD snapshots. To better visualize the different behaviors, the residues were divided into four distinct groups: the majority of residues that form 0.5-1.5 contacts per snapshot (colored in black), residues that form an unusually small number of contacts (&lt; 0.5) (colored in blue), residues that form a moderately large number of contacts (1.5-2) (colored in yellow), and residues that form a relatively large number of contacts (&gt; 2) are colored in red. For Pup, there are three distinct regions that form the largest numbers of contacts (red) comprising residues (1) Lys7, Arg8, (2) Arg28, Arg29, and (3) Arg56. They perfectly align with the three centers of Fig. <ref type="figure">3</ref> with elevated R2 values, namely (1) Arg8, (</p><p>Arg29, and (3) Arg56. For p53TAD, the residue that forms the largest number of contacts is Arg65, which is surrounded by residues with a number of contacts below average between 0.5 and 1.0. This rationalizes why R2 of Arg65 shows a local maximum that is still lower than R2 in other regions of p53TAD, such as residues 19-26 forming a residue cluster with an intermediate number of contacts. Notably, the 11 N-terminal residues of p53TAD display a lower-thanaverage amount of contacts, which is consistent with low R2 values observed across all 10 individual MD trajectories. When the same type of contact analysis is performed with side-chain atoms only, a similar behavior is observed with only a small, systematic reduction in contacts (Fig. <ref type="figure">S3</ref>) reflecting that the majority of medium-to long-range inter-residue contacts are made by side-chain atoms.</p><p>We also grouped the number of contacts per snapshot formed by each residue according to residue type and normalized them by the number of residues of the same type. The resulting value for each amino acid residue type present in p53TAD and Pup reflects their inherent contact propensity (Fig. <ref type="figure">5C,</ref><ref type="figure">D</ref>). These profiles display the following trends: positively charged residues arginine and lysine are on average most prone to form contacts, followed by hydrophobic residues isoleucine and leucine as well as aromatic residues tryptophan and phenylalanine.</p><p>Negatively charged residues aspartate and glutamate, however, are least disposed to form contacts. This may be also a consequence that both IDPs are overall negatively charged (-14e for p53TAD and -12e for Pup). When acidic residues outnumber basic residues, the former tend to repulse each other, thereby increasing Rg, while the latter have more options to interact with an acidic residue than vice versa leading to an increase of the contact propensity of basic over acidic residues. </p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>Contact analysis by graph theory.</head><p>To investigate the nature of some of the most frequent pairwise contacts in these IDPs, the MD snapshots were analyzed by graph theory where each snapshot is represented as an undirected graph with each residue corresponding to an edge and an inter-residue contact corresponds to an edge connecting the two residues (nodes). The resulting graphs were then analyzed in terms of clusters, which are disconnected graph components that do not have any edges to nodes outside of the cluster. On average 6.0 clusters per snapshot are found for p53TAD and 5.4 clusters for Pup. The probabilities of a cluster to have a given size are represented for both IDPs by the histograms of cluster sizes (Fig. <ref type="figure">6A</ref>), which reveal that clusters consisting of 2 nodes are most abundantly present (around 40%) in both p53TAD and Pup. Moreover, the cluster size probability decreases rapidly with increasing size. For instance, the fraction of clusters with 10 or more nodes (residues) is only 2-3%. Despite their sequence independence and different lengths, the two IDPs have strikingly similar cluster size distributions. The number of edges grows on average linearly with the number nodes (straight solid line), which is much slower than the quadratic behavior of complete graphs (dashed line, Fig. <ref type="figure">6B</ref>). In fact, most of the clusters formed during MD simulations are sparse graphs with a relatively small average edge-to-node ratio of 1.54, which is indicative of tree-like graphs consisting mostly of linear branches with few cross-links.  <ref type="formula">48</ref>) These collective motions involve medium to longer-range interactions between residues that can be elucidated by graph theoretical analysis of the MD trajectories described here. For Pup, many of these slower motional modes have correlation times around 3-4 ns whereas for p53TAD they are on average twice as large.</p><p>For both proteins the three distinct bands of time scales are pervasive across their polypeptide sequence (Fig. <ref type="figure">3E,</ref><ref type="figure">F</ref>).</p><p>MD methodology has made great strides in recent years to toward an increasingly realistic representation of disordered proteins. <ref type="bibr">(26)</ref> Besides experimental scattering data, quantitative NMR has played a key role for the independent validation of MD ensembles.</p><p>Because NMR spin relaxation parameters fully quantitatively reflect IDP dynamics at atomiclevel resolution both in terms of motional amplitudes and time scales, their accurate reproduction by MD has been an important but also very challenging task. A recent comparison of commonly used MD force fields that do not use residue-specific backbone potentials showed for several IDPs significant force-field dependences with the best results obtained when the analysis was restricted to average correlation functions of chunks of 10-ns subtrajectories.( <ref type="formula">56</ref>) The need to exclude slower time-scale motions, which are prominent in both experimental data and simulations (see for example Fig. <ref type="figure">3</ref>), may reflect the lack of convergence due to limited sampling. Beneficial for all simulations was the improvement of the TIP4P-D water model over TIP3P preventing overly collapsed IDP ensembles, which is consistent with other computational studies. <ref type="bibr">(38,</ref><ref type="bibr">57)</ref> Because of the observed discrepancies between experiments and MD simulations, some studies applied post factum adjustments to the MD simulations in order to improve agreement, which include uniform or selective scaling of the MD time scale or correlation times <ref type="bibr">(27)</ref><ref type="bibr">(28)</ref><ref type="bibr">(29)</ref><ref type="bibr">(30)</ref> or the reweighting of sub-trajectories.( <ref type="formula">64</ref>) Here, we chose a different approach:</p><p>rather than relying on post factum modifications, we use the residue-specific ff99SBnmr2 force field, which was specifically designed for the improved representation of IDPs without the need of any corrections.(57, 58) A correction-free MD approach has recently been reported for the intrinsically disordered SH4UD protein with the Amber ff03ws force field, which does not use residue-type independent backbone dihedral angle potentials, and no time-scale dependent data, such as NMR spin relaxation, were used for validation.( <ref type="formula">79</ref>) NMR chemical shifts were backcalculated using SHIFTX2,(80) which, besides 3D structural information, makes extensive use of protein sequence data. Here, we back-calculated NMR chemical shifts using PPM(81) (Fig. <ref type="figure">S4</ref>), which only uses the physical parametrization of chemical shifts with respect to 3D protein structure of each snapshot, achieving very good agreement with experiment. <ref type="bibr">(73)</ref> The close correspondence observed between experimental and computed 15 N relaxation R1 and R2 relaxation rates for both IDPs studied here (Fig. <ref type="figure">3</ref>), without the need for post factum corrections, attests to the accuracy and robustness of the computational protocol used. It applies REMD for the generation of conformational ensembles belonging to different temperatures from which 10 representative structures at 300 K were randomly selected as starting structures for 1&#181;s MD trajectories whereby all simulations made use of the ff99SBnmr2 force field and the TIP4P-D water model. MD-derived longitudinal 15 N R1 follow the shapes of the experimental R1</p><p>profiles with a small tendency to underestimate the experimental 15 N R1 rates by 4-6% whereas</p><p>15 N R2 relaxation rates overestimate the experimental values on average by 26% for Pup and 34% for p53TAD. This level of agreement is significantly better than for previously reported comparisons of this type.</p><p>Few individual trajectories (10-20%) show systematically larger differences with respect to experiment than the rest. For the proteins studied here, they are trajectories #4 of p53TAD and #5 and #7 of Pup (Fig. <ref type="figure">S2</ref>, Table <ref type="table">S4</ref>). These trajectories are characterized by the persistent formation of secondary structure (#4 of p53TAD and #5 of Pup) (Fig. <ref type="figure">4C,</ref><ref type="figure">D</ref>) or by a collapsed overall geometry with a reduced &lt;Rg&gt; compared to the other trajectories (trajectory #7 of Pup) (Table <ref type="table">S5</ref>). At the individual trajectory level, these outlier trajectories are in poorer agreement with experimental data and their removal from the set of 10 trajectories during the backcalculation of relaxation rates further improves the agreement with experiment (Fig. <ref type="figure">S3</ref>). From such diagnostic analysis it follows that these outlier trajectories are either overrepresented in the original simulations or the result of simulation artifacts, for example, caused by inaccuracies of the underlying force field. Removal of individual trajectories based on comparison with experiment should be applied with great care and be reserved primarily for diagnostic purposes, such as the analysis of shortcomings of the simulations. While post factum trajectory selection or reweighting can provide better agreement with experiment, it is generally unclear whether the altered ensembles are in fact consistent with an alternative, physics-based force field, thereby complicating the physical interpretation of such ensembles.</p><p>Although it is difficult to identify individual force field terms responsible for the IDP behaviour observed in the outlier trajectories, these results can nonetheless provide useful input to guide future force field improvements. With more computer power, it will be possible to gain better statistics by generating a larger number of trajectories for the improved sampling of conformational space allowing the more rigorous assessment of the underlying force field, the water model, and other aspects of the computational methods used. Conversely, such insights may allow the further improvement of force fields and methods for applications also to other proteins. In fact, the ff99SBnmr1 force field, which is the parent force field of ff99SBnmr2, was developed and optimized using this strategy by the systematic reweighting of MD snapshots based on many trial force fields using experimental NMR data of intact proteins. <ref type="bibr">(82)</ref> The good agreement of the MD simulation with experimental observables both motivates and justifies the analysis of other protein properties observed in the MD trajectories that are difficult to measure. This includes the analysis of transient inter-residue interactions. The molecular driving forces of these interactions are fundamentally similar to those of ordered proteins although average hydration properties may differ. <ref type="bibr">(79)</ref> In contrast to ordered proteins, inter-residue interactions between non-sequential amino acids are short-lived. Therefore, the time-averaged interaction maps (Fig. <ref type="figure">4A,</ref><ref type="figure">B</ref>) offer only partial insights as they conceal the compositions and distributions of instantaneous interaction clusters. In fact, the relatively large network reflected by the average contact map contrasts the much smaller size of graphs that exist at any given time, which attests to the very heterogeneous and transient nature of instantaneous contact clusters. The highest occupancy of pairwise contacts found is around 0.5, which mostly belong to (i,i+3) contacts. For a list of the most frequent pairwise contacts, see Tables <ref type="table">S2,</ref><ref type="table">S3</ref>.</p><p>Snapshot by snapshot analysis revealed the dominance of small cluster sizes over larger ones (Fig. <ref type="figure">6</ref>). For both p53TAD and Pup, clusters with 2 or 3 residues make up more than 50% of all clusters and clusters with more than 10 residues have notably low occurrence, although their formation could be functionally relevant during molecular recognition events. Because clusters consisting of residue pairs dominate intra-residual interactions in both IDPs, further analysis of the interaction network was performed based on pairwise contacts. Contact maps were generated for p53TAD and Pup averaged over all MD trajectories and pairwise contacts that have occupancies larger than 0.2 visualized as separate graphs (Fig. <ref type="figure">4E,</ref><ref type="figure">F</ref> The majority of clusters are linear graphs with few circular sub-graphs leading to the linear relationship between the number of nodes and number of edges (Fig. <ref type="figure">6B</ref>). Acidic residues tend to have low cluster participation whereas arginine residues have the highest participation in both proteins (Fig. <ref type="figure">5A,</ref><ref type="figure">B</ref>). This difference in cluster participation between cationic and anionic residues is also evident in Fig. <ref type="figure">5C</ref>,D. Among the neutral amino acids, those with larger sidechains are more prone to interactions with non-neighboring residues due to their intrinsically larger distance range. In fact, Pro, Val, Ser, Ala, Gly have the lowest interaction propensities among neutral residues and among pairs of chemically similar residues, such as Gln vs. Asn and Leu vs. Val, the larger residue (Gln, Leu) dominates the smaller one (Asn, Val).</p><p>A primary biological function of p53TAD is to negatively regulate p53 by interacting with the ubiquitin ligases MDM2 and MDMX for the degradation of p53. This interaction is one of the earliest and best studied interactions between an IDP and a folded protein both by experiment(67-69) and computation. <ref type="bibr">(90)</ref> In order to better understand the molecular recognition mechanism underlying the formation of this complex, a realistic and accurate description of the free state of p53TAD is of central importance. For MD studies, the choice of the protocol, especially of the force field and water model, is consequential. A recent unbiased REMD study of free p53TAD reported the detailed comparison using five different MD force fields all without residue-specific backbone potentials. Based on 1-&#181;s long replicas major differences were revealed in terms of the structural propensities among them and also with respect to experimental data.( <ref type="formula">91</ref>) An even longer simulation of residues 10-39 of p53TAD for a total length of 1.4 ms analyzed by Markov state models identified substantial populations of b-sheets across the sequence,(92) a behavior that is at variance with the above mentioned REMD ensembles(91) as well as with experimental solution NMR data.( <ref type="formula">67</ref>) These together with many other studies show that force fields need to be chosen following extensive testing to ensure that long trajectories, generated with considerable computational effort, offer the most realistic biophysical insights about these highly complex, heterogeneous systems.</p><p>In </p><p>where is cartesian coordinate a (b) (= x, y, z) of atom i in the coordinate system that has its origin in the center of mass of the molecule. Diagonalization of S yields three non-negative eigenvalues from which the radius of gyration Rg is obtained, , the asphericity A,(98, 99)</p><p>and the prolateness P,(</p><p>The asphericity measures the degree to which the three axis lengths of the ellipsoid of inertia (eigenvalues) are equal, whereas the prolateness P indicates whether the largest or smallest axis length is closer to the middle axis length. P takes values between -1 and 1, quantifying the transition from oblate to prolate shapes. Normalized time-correlation functions of Rg(t), made offset-free, were computed according to</p><p>as an average over all 1-&#181;s MD trajectories.</p><p>According to polymer theory, for an unfolded polymer the ensemble-averaged Rg scales with the number of residues N as(62, 77)</p><p>where r0 is a constant reflecting the average size of a residue and the Flory exponent n determines the overall compactness of the polymer serving as a reference.</p><p>A = (&#955; 3 -&#955; 2 ) 2 + (&#955; 3 -&#955; 1 ) 2 + (&#955; 2 -&#955; 1 ) 2 2(&#955; 3 + &#955; 2 + &#955; 1 ) 2</p><p>Back-calculation of R1, R2 relaxation rates. For IDPs, the normalized time-autocorrelation function C(t) of the lattice part of the spin-relaxation active magnetic dipole-dipole interaction cannot be factorized into an overall tumbling part and an internal dynamics part. Rather, we compute the full C(t) directly from an MD trajectory using the second-order Legendre polynomial:</p><p>&#9001;3[&#119942;(&#120591;)&#119942;(&#120591; + &#119905;)] " -1&#9002;</p><p>where e(t) is the unit vector defining the 15 N-1 H bond orientation whereby snapshots were not aligned with respect to a reference snapshot. The angular brackets indicate averaging from time t = 0 to TMD -t, where TMD is the total trajectory length. The calculation of C(t) was efficiently performed by the fast Fourier transform (FFT) using the Wiener-Khinchin theorem. For acceptable statistical convergence, the analysis of C(t) was limited to its initial portion from t = 0 -TMD /3. Next, a multiexpoential decay function was fitted to C(t):( <ref type="formula">101</ref>)</p><p>where Ai and ti are the best fitting parameters subject to the conditions:</p><p>The spectral density function J(w) can be then analytically obtained via </p><p>where</p></div>
<div xmlns="http://www.tei-c.org/ns/1.0"><head>45</head><p>) " ( 6 "5</p><p>) " &#120574; 2 " &#120574; 1 " &#9001;&#119903; 12 $7 &#9002; " and &#119888; ++ = ! !8</p><p>&#8710;&#120590; " . &#181;0 is the permeability of vacuum, h is Plank's constant, gH and gN are the gyromagnetic ratios of 1 H and 15 N, and rNH = 1.02 &#197; is the backbone N-H bond length. The 15 N chemical shift anisotropy was set to Ds = -160 ppm.</p><p>Analysis of inter-residue contacts and residue clusters by graph theory. Contact analysis was performed on all snapshots of the MD simulations of both p53TAD and Pup. A contact is considered formed when the nearest distance between atoms from two different residues is smaller than 4 &#197;. First-neighbor contacts (between residues i,i+1), and second-neighbor contacts (between residues i,i+2) were excluded since they are present for most residues. For each residue in p53TAD and Pup, the total number of contacts formed by a particular residue is determined and normalized by the number of MD snapshots. Each snapshot was converted to a graph where residues are represented as nodes and contacts between two residues are represented as edges between them. The initial graph was then decomposed into a maximal number of disconnected graph components called clusters, i.e. there is no edge between any node in the cluster and any node outside the cluster. The size of a cluster corresponds to the number of its nodes.   Table <ref type="table">S1</ref>. MD and REMD simulation details for p53TAD and Pup.</p><p>Table <ref type="table">S2</ref>. Most frequent pairwise residue contacts in p53TAD from MD simulations.</p><p>Table <ref type="table">S3</ref>. Most frequent pairwise residue contacts in Pup from MD simulations.</p><p>Table <ref type="table">S4</ref>. Chemical shift comparisons for p53TAD.</p><p>Table <ref type="table">S5</ref>. Radius of gyration pf p53TAD and Pup.</p></div></body>
		</text>
</TEI>
