A method to construct the dynamic landscape of a bio-membrane with experiment and simulation

Biomolecular function is based on a complex hierarchy of molecular motions. While biophysical methods can reveal details of specific motions, a concept for the comprehensive description of molecular dynamics over a wide range of correlation times has been unattainable. Here, we report an approach to construct the dynamic landscape of biomolecules, which describes the aggregate influence of multiple motions acting on various timescales and on multiple positions in the molecule. To this end, we use 13C NMR relaxation and molecular dynamics simulation data for the characterization of fully hydrated palmitoyl-oleoyl-phosphatidylcholine bilayers. We combine dynamics detector methodology with a new frame analysis of motion that yields site-specific amplitudes of motion, separated both by type and timescale of motion. In this study, we show that this separation allows the detailed description of the dynamic landscape, which yields vast differences in motional amplitudes and correlation times depending on molecular position.

Supplementary Note 1: Experimental acquisition and data processing 37 1.1 Assignment and initial spectrum fits 38 POPC has been assigned previously, 1-3 with the 13 C assignment shown in Supplementary 39 Figure 1. In order to extract experimental relaxation rate constants as accurately as 40 possible, we simultaneously fit series of 1D spectra for each measurement to a decaying 41 (or recovering) exponential function. In order to achieve this, we use a reference fit of the 42 spectrum, also shown in Supplementary Figure 1, where the experimental, fitted, and 43 residual spectra are displayed. This is performed with the INFOS software, which has been 44 previously shown to accurately separate even heavily overlapped peaks, especially when 45 fitting a series of spectra simultaneously. 4

Experiments 62
Four types of experiments were used in this study to experimentally characterize dynamics. 63 These were 13 C T 1 relaxation, via saturation recovery, 1 H-13 C steady-state NOE, T 1ρ 64 relaxation ( 13 C transverse relaxation under spin-locking), and measurement of one-bond 65 1 H-13 C residual dipole couplings via DIPSHIFT. These are shown in Supplementary Figure  66 2, with experimental details given in Methods Table 1. 67 13 C T 1 experiments are performed by an initial saturation period of the 13 C 68 magnetization, followed by a variable length recovery period (τ). Following the recovery 69 period, a π/2-pulse is applied to the 13 C to read out the magnetization. Throughout the 13 C 70 saturation and recovery period, 1 H magnetization is saturated, allowing us to decouple the 71 evolution of 1 H magnetization from 13 C magnetization. Note that lipid samples are not 72 labeled (1% 13 C from natural abundance). This prevents a significant influence from 13 C 73 spin-diffusion. 74 Steady-state heteronuclear NOE was acquired by simply recycling the experiment 75 with a long delay (10 s, at least 3.4x 13 C T 1 ), once with 1 H saturation and once without 1 H 76 saturation. The ratio of these experiments yields the NOE enhancement ( 〈C z SS 〉 ), which is 77 related to the 13 C R 1 and 1 H-13 C NOE rate constant ( σ HC ) as follows: 78 (1) 13 C T 1ρ experiments were acquired approximately on-resonance (see Methods), in 79 order to avoid offset effects on the relaxation rate constants. 5 Experiments begin with a π/2-80 pulse on 13 C followed by variable length spin-lock, and acquisition directly after the spin-81 lock. 82 DIPSHIFT experiments were used to measure the one-bond 1 H-13 C dipole 83 couplings. 6 Experiments begin with a 13 C π/2-pulse followed by homonuclear decoupling 84 applied on 1 H for a variable fraction of one rotor period, followed by heteronuclear 85 decoupling for the remainder of the rotor period. A 13 C π-pulse is applied and the 13 C 86 chemical shift is allowed to refocus for one rotor period, followed by acquisition. 87 relaxation measurements, the full series of spectra for each data set are fitted to a decaying 103 exponential (decaying towards zero for R 1ρ and towards some constant for R 1 ), using the 104 INFOS "FitTrace" function, with variable amplitude and relaxation rate constant, but fixed 105 position and linewidth (taken from the initial fit). An example fit is found in Supplementary 106 Figure 3. NOE on and off signal intensities were found using the "FitSpec" function. 107 DIPSHIFT spectra were analyzed separately, with amplitudes extracted for each separate 108 spectrum. Then, for each peak, the extracted amplitudes are fitted against simulated 109 DIPSHIFT curves. Curves are simulated for 200 order parameters ranging between 0 and 1 110 (referenced to a 21.5 kHz H-C dipole coupling), where frequency-switched Lee-Goldburg 111 homonuclear decoupling is explicitly simulated, 7 and simulations include the correct number 112 of 1 H for the given resonance. This allows one to correctly capture the dipole coupling 113 scaling factor due to the decoupling, and to account for interference in dephasing due to 114

123
For all experiments, we obtain a standard deviation of the measured relaxation rate 124 constant or order parameter using a simple bootstrap approach. 8 For DIPSHIFT and R 1 /R 1ρ 125 measurements, we extract peak amplitudes from the series of spectra. Then, for each peak, 126 we refit the corresponding curve 200 times, where we randomly resample the original data 127 set. That is, suppose we have N data points for a given experiment and resonance 128 (including multiple time points and also repetition of those time points). Then, we select N 129 data points with replacement from the set, and refit the result. The standard deviation over 130 parameters resulting from the 200 bootstrapped data sets is reported. For the NOE 131 experiment, we extract error of the peak amplitudes using the "FitError" function of INFOS, 132 and combined with the bootstrapped error for 13 C R 1 , we apply the usual propagation-of-133 error rules to obtain an error for the NOE rate constant, although we note that this error 134 results almost entirely from 13 C R 1 error, as the NOE peak heights may be very accurately 135 extracted with INFOS. These errors are then used to determine the detector response error 136 via linear propagation-of-error. 137

Detector analysis 138
Detectors can be thought of as a linear recombination of experimental rate constants. This 139 is described in detail elsewhere, 9,10 but a brief review of the work and some updates to 140 computational methods is merited here. We first, assume that the correlation function of 141 motion (either experimental or MD-derived) is a sum of decaying exponential terms, such 142 that 143 (2) Then, the correlation function always begins at 1 and plateaus at S 2 . The rate of decay 144 may be described by one or more correlation times, where θ(z) indicates decay is 145 distributed as a function of correlation time (z is the log-correlation, z = log 10 (τ c / s) ). Then, 146 a relaxation rate constant (indexed by ζ) is equal to 147 where R ζ (z) is the sensitivity of experiment ζ. The sensitivity is the relaxation rate 148 constant's dependence on correlation time. That is, if a correlation function has mono-149 exponential decay with rate constant τ c and decays completely to 0 ( S 2 = 0 ), then the 150 resulting rate constant for experiment ζ is R ζ (z) for z = log 10 (τ c / s) . In this case, the integral 151 in (3) then adds up the scaled contributions for all correlation times in a multi-exponential 152 correlation function, where the scaled contribution is the product of the amplitude of decay 153 for that correlation time and the sensitivity of the experiment at that correlation time, R ζ (z) . 154 Examples of R ζ (z) are given in Supplementary Figure 4A (top). Then, the basic concept 155 behind detectors is that linear combinations of rate constants can be created to generate 156 detectors that are sensitive to specific timescale windows. Suppose we take a sum of 157 experimental rate constants, which we will define as 158 ρ n (θ ,S) = a n,ζ R ζ The linear combination yields a new parameter with a new relationship to the distribution of 159 correlation times of motion, (1− S 2 )θ(z) . The relationship is defined by the sensitivity of the 160 detector, ρ n (z) = Σ ζ a n,ζ R ζ (z) . Then, a detector response, ρ n (θ ,S) , quantifies the amplitude of 161 motion having correlation times in the window defined by the detector sensitivity.

206
A few additional considerations must be taken for analysis of MD-derived correlation 207 functions with detectors. First, we must deal with the offset term, S 2 . This term collects non-208 decaying components of the correlation function. Practically, however, there is no difference 209 between non-decaying components and very slowly decaying components (for example, 210 any motion at least ~5 times slower than the length of our MD simulations, so about 42 µs, 211 will not appear to decay during the simulation). Then, we simply neglect S 2 when analyzing 212 MD-derived correlation functions with detectors, and allow its contribution to be part of the 213 distribution, θ(z) , for z corresponding to arbitrarily long correlation times. Furthermore, we 214 require the relative standard deviations of the various time points in a correlation function 215 when optimizing detectors. We assume that the standard deviation is proportional to m −1/2 , 216 where m is the number of pairs of time points used to calculate that element of the 217 correlation function. m is given by N-n, where N is the total number of time points in the MD 218 simulation, and n is the index of the time point of the correlation function. Note that the first 219 element of the correlation function is exactly 1 in all cases. To deal with this numerically, we 220 set its standard deviation to 1x10 -7 , which is sufficiently small to force a nearly exact fit of 221 this time point, but sufficiently large to avoid dynamic range problems in our calculations. 222

Experimental and Simulated Detector Optimization 223
Detector optimization is performed as described previously: 10 M ! is an approximation of M, and in fact is the best approximation of M possible via linear 228 combination of t vectors. Columns of V t are then linearly recombined to form detectors. By 229 first selecting the largest singular values, we improve the data fit and minimize error of the 230 resulting detectors. Large choices of n will yield higher resolution detectors, but more error. 231 Linear recombination of the columns of V t may be optimized to either match some 232 target function, or to yield detectors with minimal width and overlap. The former case is 233 used when comparing MD results to NMR results: we already have NMR sensitivities, and 234 so we optimize MD sensitivities to match the NMR sensitivities. This may be achieved with 235 the 'lstsq' algorithm (linear least squares) in NumPy's linear algebra module. 13 Note that we 236 use the 15 largest singular values when using this approach. However, we must use the 237 latter case when optimizing the experimental detectors since we do not have an obvious 238 target function. In this case, we sweep through an array of 200 correlation times (log-239 spaced from 10 -14 to 10 -3 s), and at each correlation function, we attempt to optimize a 240 detector sensitivity that is equal to 1 at the current correlation time, and minimized, but non-241 negative elsewhere. This is achieved with the 'linprog' algorithm in SciPy's optimization 242 module. 14 Then, if n singular values are used, we find that for exactly n correlation times, 243 the maximum of the sensitivity is found at the same correlation time where the sensitivity is 244 forced to equal 1. We take the detectors from these n sensitivities. For experimental 245 analysis, we take n equal to 6. When analyzing MD data without matching sensitivities to 246 experiment, we take n equal to 6 or 9 (6: Supplementary Figure 10 which contains the detection vectors. 9 The corresponding matrices are given here 252 where the ρ n (θ ,S) are varied (bound between min(ρ n (z)) and max(ρ n (z)) = 1), using the 'lstsq' 257 algorithm of NumPy's linear algebra module. 13   to the lack of a bonded 1 H. For DIPSHIFT, we note that S 2 is already very small, and so we 261 simply set to the mean MD-simulated value (0.0365), with minimal impact on the results 262          detectors are compared to those obtained for other positions in Supplementary Figure 6. 296 One sees that ρ 0 (z) obtained for C' (black, dashed) extends to longer correlation times 297 than for other positions (color), and ρ 1 (z) is missing. ρ 2 (z) and ρ 3 (z) are shifted to slightly 298 shorter correlation times, whereas ρ 4 (z) and ρ 5 (z) are practically unchanged. lines. Missing NOE data requires us to omit one detector for analysis of C' relaxation, resulting in ρ 1 being 304 absent for C' analysis. Furthermore, ρ 0 (z) is shifted to longer correlation times, and ρ 2 (z) and ρ 3 (z) are 305 shifted to slightly shorter correlation times.

306
A second challenge in the experimental analysis is that ρ 0 (z) contains multiple 307 regions of non-zero sensitivity: one region capturing all motion faster than about 110 ps, 308 and a second window capturing motion slower than ~3.7 ns but faster than ~6.3 µs (a third 309 window appears around 1 ms, but the small experimental detector responses for ρ 4 and ρ 5 310 indicate that it is unlikely to contribute significantly to ρ 0 (θ ,S) ). Given the initial agreement 311 between experiment and simulation, we may use MD simulation to determine how 312 important those multiple windows are. In Supplementary Figure 7 frames. Note that if we assume that the individual correlation functions are multi-355 exponential, as in Supplementary Equation (2), then the product must also be multi-356 exponential. We demonstrate this for a product of two correlation functions, but the result 357 may be easily iterated over multiple correlation functions. 358 One sees that the total correlation function resulting from the product of two multi-359 exponential correlation functions is indeed itself multi-exponential. The resulting correlation 360 function has an offset term, S 2 = S 1 2 S 2 2 , includes the distribution corresponding to the first 361 correlation function, scaled by S 2 2 and the second correlation function, scaled by S 1 2 , and 362 finally a double integral containing effective correlation times resulting from the product of 363 two exponential terms coming from each of the correlation functions. Note that an explicit 364 expression for (1− S 2 )θ(z) may be obtained in analogy to the previous result from Smith et 365 al. 10 366 (1− S 2 )θ(z) = S 2 2 (1− S 1 2 )θ 1 (z) + S 1 2 (1− S 2 2 )θ 2 (z) + (1− S 1 2 )(1− S 2 2 ) θ 1 (z 1 )θ 2 (z + z 1 − log 10 (10 z 1 − 10 z )) 10  Figure 9A), and also take the chains as a whole (Supplementary 380 Figure 9B). We also attempt to separate internal motion into components. One approach is 381 to use a C-C bond as frame (one of the C should be bonded to the corresponding H), thus 382 separating rotation around the C-C bond from reorientation of the C-C bond 383 (Supplementary Figure 9C). In a second approach, we separate motion parallel to the MOI 384 from motion perpendicular to the MOI, by defining a frame that projects the bond direction 385 onto the plane perpendicular to the MOI (Supplementary Figure 9D). The procedure for 386 defining each frame is discussed below. 387

Definition of the frames of NMR interactions 398
Defining the frame of the NMR interaction (dipolar coupling) requires defining a full axis 399 system (x, y, and z-axes). While definition of the z-axis is usually obvious-for a one-bond 400 dipole coupling, the z-axis should simply be along the bond-the direction of the x-and y-401 axes is less so. In fact, the choice is somewhat arbitrary, with the relatively simple 402 requirement that the chosen axes move along with the NMR interaction. Then, for each H-403 C bond, we define the C-X (X being a heteronucleus, usually another C) leading towards 404 the g 2 carbon to be in the xz-plane (for g 2 , we take the g 1 carbon). Using the vector along z 405 (the H-C) bond, and a vector in the xz-plane, it is straightforward to obtain the x-and y-406 axes (y is the cross-product of the z and xz vectors, x is the cross product of the y and z 407 vectors). For the carbonyls, we assume the z-axis is along the C=O bond, and use the 408 single bonded oxygen to define the xz-plane. 409

Definition of RMS frames 410
RMS frames (root-mean square alignment frames) are defined by the rotation required to 411 rotate the current structure such that a selection of atoms are aligned to a reference 412 structure (acquiring a linear-least squares fit). The librational frames and head group RMS 413 frame use this approach. For the librational frames, for a given H-C bond (or C=O bond, for 414 the carbonyl), we select the central carbon and all atoms bonded to it as our reference 415 selection. The reference structure is simply taken as the positions of those atoms at the 416 beginning of the trajectory. For the alignment of the head group, we select the three 417 glycerol-backbone carbons (g 1 , g 2 , and g 3 ) and the three oxygen atoms bound to them. 418 To obtain vectors defining the frame at some time τ , we first subtract away the mean 419 position of all atoms in the reference structure and in the atoms at time τ . We then solve 420 for the rotation matrix that results in the best alignment of the reference atoms to the atoms Once H is obtained, we acquire its singular value decomposition (SciPy linear algebra 424 module 14 ), from which we may calculate the rotation matrix. 425 Then, the frame is defined by vectors ν X (τ ) , ν Y (τ ) , and ν Z (τ ) which are the columns of the 426 rotation matrix. 427

Definition of MOI frames 428
In each of the two chains, we separate overall motion of the moment of inertia (MOI) from 429 motion within the chain. Defining this frame is simply a matter of calculating the MOI of the 430 chain (we only take the C nuclei), and extracting the largest component, which defines the 431 z-axis of this frame. The MOI matrix (neglecting masses) is defined as 432 where r k is a vector to the point mass (center of mass). The direction of the largest 433 component (z) may be obtained by computing the eigenvalues and eigenvectors of I, and 434 then taking the vector corresponding to the largest eigenvalue. 435

Definition of bond frames 436
A bond frame is a frame defined simply by the direction of a bond. That is, we take the z-437 axis of the frame to be the bond direction. x-and y-axes are not required, and no further 438 computation is necessary (aside from vector normalization). 439

Separating parallel and perpendicular motion: MOIxy frame 440
We define a frame to separate internal motion in the chains into components parallel and 441 perpendicular to the MOI. This is achieved with a frame that projects the direction of the  Figure 9B). We also no longer separate the chains into upper and lower 466 parts when applying the MOI frame. The changes yield significant improvement in the 467 HG/BB, and minor improvements in the chains, so that we use RMS alignment of the 468 glycerol in the HG/BB and take the MOI of each chain, without separation into upper and 469 lower chains. 470 We continue to attempt to further separate internal motion into components. Our first 471 attempt is to separate internal motion by taking a C-C bond (one C should be in the H-C 472 bond), in order to separate rotation around the C-C bond from reorientation of the C-C 473 bond (we always take the bond leading towards the g 2 position). The result is severe 474 disagreement in Supplementary Figure 10D. This is the result of strong correlation between 475 these two motions, and so we cannot use this frame separation. Finally, we separate 476 motion within the two chains into components parallel and perpendicular to the MOI in 477 Supplementary Figure 10E 491 We use several short videos to help illustrate both detector responses and the behavior of 492 residual tensors due to motions. We also plot the separation of motions using frames, 493 juxtaposed both against residual tensors and against the resulting detector responses. Here 494 we briefly describe how the various pieces of information are plotted in the videos. 495

Time indicator 496
In each video, a time indicator is displayed. This shows how much time is elapsed in the 497 MD simulation within one second of the video. We log-space the MD frames shown in the 498 video so that this value will change, allowing one to view different timescales of motion as 499 the frame rate accelerates. For example, if the time indicator shows "1 s : 1 ns", this means 500 that from 0.5 s earlier (15 frames) in the video and 0.5 s later (15 frames), 1 ns of the MD 501 trajectory will be covered. Then, each video covers 125 ns of trajectory, but the indicator will 502 display a final value of 83 ns. Note that at the beginning of the video, we use "1 s : 2 ps", 503 but since we have saved frames only every 5 ps, resulting in some jumps at the beginning 504 of the movies. 505

Plotting detector responses 506
In plots of detector responses (e.g. det_overall.mov), we see detector responses as 3D 507 representations, with detector response encoded both as a color and atomic radius. The 508 information displayed changes as the motion accelerates in the video, according to the 509 following procedure: 510 1) We first take the correlation time corresponding to the current value on the time 511 indicator, so that the motions taking about 1 s to occur are then highlighted with the 512 detectors. We denote the log of this correlation time as z c 513 2) We calculate the radius from the detector responses and z c : 514 Note that this is just a weighted sum of the detector responses. If we are at the 515 maximum for detector n, then r ≈ 0.9 Å + (4 Å)ρ n (θ ,S) , since max(ρ n (z)) = 1. In practice, 516 it will be slightly larger since the other detector sensitivities are not quite zero. Then, the RGB value (on a scale 0-255) used is: 524 The normalized response, x n , is obtained by dividing by the maximum of the detector 525 responses for all bonds, and all detectors. Then, the RGB value is a weighted 526 average of the four colors above, where the weighting depends on the size of the 527 sensitivity of each detector at z c and the detector response, and the color tan, 528 corresponding to no response, where if the sum of detector responses times detector 529 sensitivities is 1, then tan makes no contribution, but if the sum is zero, then the color 530 is only tan. At the maximum of a given detector, for example ρ 2 , the color is 531 approximately a weighted average of green and tan, where if ρ 2 (θ ,S) = 1, then we will 532 have just green, but if ρ 2 (θ ,S) = 0 , then the color will be just tan.  Figure 12E), we use a Gaussian distribution, otherwise we fit to a 552 skewed Gaussian distribution (2x broader towards long correlation times). Note that 553 detector analysis of correlation functions does not distinguish between very slow motion 554 and no motion (i.e. S 2 . Practically, ρ 8 (θ ,S) is a good approximation for S 2 , but strictly 555 speaking, amplitude appearing in this detector may be from S 2 or from very slow motion). 556 Therefore, we simply insert S 2 into the distribution at the longest correlation time, i.e. z=-3 ( 557 τ c = 1 ms ). The resulting distributions, θ(z) , therefore always integrate to 1. 558 Gaussian Distribution:

Choice of functional form of the distribution 559
Nevzorov, Trouard, and Brown have shown that collective motions occurring over d 560 dimensions results in a power law where the spectral density is proportional to 16 561 We note that this power law results from a distribution of correlation times having the 562 following form 563 For d=3, increasing correlation times lead to a decreasing amplitude of the distribution, and 564 conversely, for d=1, increasing correlation times lead to an increasing amplitude, whereas 565 d=2 yields a uniform distribution. We may verify that θ(z) ∝ 10 z(1−d /2) leads to a spectral 566 density with the functional form in (21) by integrating the spectral density over the 567 distribution 568 The second term in the integral, 10 z ⋅ s / (1+ (ω ⋅10 z ⋅ s) 2 ) , is symmetric around 570 z 0 = − log 10 (ω ) , and so we assume θ(z) is approximately linear around this value 571 (reasonable for a broad, monotonic distribution as expected for a power law), such that 572 The term m(z − z 0 ) is antisymmetric about z 0 , whereas 10 z ⋅ s / (1+ (ω ⋅10 z ⋅ s) 2 ) is symmetric 573 about z 0 , such that the integral of the product of these terms is zero. Then, the integral 574 reduces to 575 A substitution of variables brings us to a familiar integral, and so we find that the result is 576 proportional to ω −(2−d /2) , which is the power-law from Brown and coworkers. 16 577 The original derivation by Brown and coworkers makes a few important 578 approximations. First, the distribution of correlation times results from there being many 579 modes of motion, each having a correlation time. However, one assumes these may span 580 from zero to infinitely long wavelengths, whereas in reality, the modes have a minimum 581 length (limited by the molecule size or bond lengths) and also a maximum length (limited at 582 least by vesicle size). One consequence of this assumption is that the distributions 583 themselves all integrate to infinity. A second assumption that may introduce some deviation 584 of the real distribution from the ideal power law is the calculation of the amplitude of a given 585 mode: the amplitude of a mode is given as 〈 δn(q,t) 2 〉 = kT / Kq 2 , as a result from the 586 classical equipartition theorem (q is the wave vector, K the elastic constant). 16 Then, one 587 assumes the resulting spectral density is thus proportional to the mode amplitude. 588 However, this neglects that faster motions result in an effective interaction tensor that has 589 been scaled, such that the influence of slower motions is less than would be estimated by 590 the power law (the distribution should decay faster than predicted by the power law). For 591 fast motion, this will have little effect, and if the total amplitude of motion is low (S 2 near to 592 1), the influence of such behavior will also be negligible, but in membranes, the total 593 amplitude may be quite large. We should be clear: these assumptions should not 594 significantly hinder the application of power-laws to experimental data, since one is typically 595 interested in the behavior of R 1 or R 1ρ at a relatively narrow range of correlation times. 596 However, if we are to characterize the distribution over a wider range of correlation times, 597 deviation of the integral of θ(z) , in particular, can be problematic. 598 While a thorough investigation into the best functional form to describe the 599 distribution of correlation times for collective motions is beyond the scope of this study, we 600 note that for d=3, a skewed Gaussian distribution is a good estimate, as given in 601 Supplementary Equation (15). At short correlation times, we know the distribution should 602 approach zero due to bounds on the minima mode length, although we do not have a 603 particular functional form. As one approaches long correlation times, our function should 604 decay slightly faster than predicted by the power law. In Supplementary Figure 11, we see 605 that the skewed Gaussian distribution has the desired property, with particularly close 606 correspondence to the power law for a full-width at half maximum of 2 (note that we later 607 find this value to be close to the fitted widths for internal motion, see Supplementary Figure  608 13). 609 For d=2, we simply take an un-skewed Gaussian distribution. Clearly, this 610 distribution cannot be actually uniform, since it would then integrate to infinity, but how it 611 should decay to zero for short or long correlation times is not predicted by the power-law, 612 so that we simply take a symmetric function. 613

Comparison of detector responses 619
Once the distribution is obtained, we may calculate the detector response, noting 620 that S 2 has been absorbed into θ(z) in (15), and furthermore, the integral may be 621 discretized. 622 While the results are very good, considering the multiple layers of complexity in this 695 analysis, we would like to finally refine the dynamic landscape using the NMR experimental 696 data. We cannot simply vary all parameters of the landscape because the number of free 697 parameters would outnumber the amount of experimental data obtained. On the other 698 hand, it is reasonable to vary a single free parameter for each measured resonance. In 699 Supplementary Figure 17, we see the largest disagreement for ρ 0 (θ ,S) and ρ 1 (θ ,S) , suggesting 700 that we should adjust the internal motions, since these have the largest detector responses 701 for ρ 0 and ρ 1 . Furthermore, we find that simulation usually overestimates ρ 0 (θ ,S) and 702 underestimates ρ 1 (θ ,S) , which points to a disagreement in correlation time. Therefore, at 703 each position, we vary τ max for the internal motion. Within chains, where motion is 704 separated into parallel and perpendicular components, we will vary the correlation time of 705 both motions using the same scaling factor, and for positions having experimentally 706 overlapping resonances, we will vary all correlation times also by a single scaling factor. 707 Therefore, we still only have a single free parameter per resonance. 708 The results are shown in Supplementary Figure 17