{'fraction_non_alphanumeric': 0.03221593382673052, 'fraction_numerical': 0.013849586417065738, 'mean_word_length': 4.705580998214701, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theoretical and Computational Framework for the Analysis of the Relaxation Properties of Arbitrary Spin Systems. Application to High-Resolution Relaxometry

A wide variety of nuclear magnetic resonance experiments rely on the prediction and analysis of relaxation processes. Recently, innovative approaches have been introduced where the sample travels through a broad range of magnetic fields in the course of the experiment, such as dissolution dynamic nuclear polarization or high-resolution relaxometry. Understanding the relaxation properties of nuclear spin systems over orders of magnitude of magnetic fields is essential to rationalize the results of these experiments. For example, during a high-resolution relaxometry experiment, the absence of control of nuclear spin relaxation pathways during the sample transfers and relaxation delays leads to systematic deviations of polarization decays from an ideal mono-exponential decay with the pure longitudinal relaxation rate. These deviations have to be taken into account to describe quantitatively the dynamics of the system. Here, we present computational tools to (1) calculate analytical expressions of relaxation rates for a broad variety of spin systems and (2) use these analytical expressions to correct the deviations arising in high-resolution relaxometry experiments. These tools lead to a better understanding of nuclear spin relaxation, which is required to improve the sensitivity of many pulse sequences, and to better characterize motions in macromolecules.


Introduction
The development of most Nuclear Magnetic Resonance (NMR) experiments requires the understanding of relaxation properties. Improvement in the sensitivity and resolution have been obtained, ranging from the use of an optimum excitation angle with respect to the longitudinal relaxation, known as the Ernst angle [1], to the development of Transverse Relaxation Optimized SpectroscopY (TROSY) experiments [2,3] that exploit relaxation interferences [4][5][6][7]. An in depth investigation of relaxation processes is particularly critical to design and interpret several classes of experiments which are based on moving the sample through a broad range of magnetic fields. A variety of such experiments have been designed recently: (1) The existence of Long-Lived States (LLS) [8,9] was revealed by the combination of high-field coherent evolution and low-field relaxation; (2) In dissolution Dynamic Nuclear Polarization (dDNP) [10][11][12], the hyperpolarized sample is transferred back and forth between the polarizing magnetic center and the high-field spectrometer through magnetic fields that can be as low as the earth magnetic field; (3) Multi-scale dynamics can be characterized with Fast-Field Cycling (FFC) relaxometry [13] where the magnetic field is switched from ca. 1 T down to ca. 100 µT; (4) A sample-shuttle apparatus can be used to combine relaxometry experiments with high-field NMR [14][15][16][17] to gain atomic resolution description of molecular dynamics; (5) This kind of device can also be used to investigate relaxation properties of spin terms that are only relevant at low fields [18]; (6) A sample shuttle may couple two magnetic centers in a two-field NMR spectrometer [19] to record multi-dimensional experiments where spins are manipulated at two vastly different fields [19][20][21][22]. Sample-shuttling experiments have been used to measure longitudinal relaxation rates over orders of magnitude of magnetic fields and characterize the dynamics of membrane vesicules [23], protein backbone [17,24] and side-chains [25]. This type of experiments, called High-Resolution Relaxometry (HRR), consists in the measurement of relaxation rates over a broad range of magnetic field while preserving the high resolution of conventional high-field magnets (i.e. higher than 9 T) [14,15]. This approach relies on moving the NMR sample in the stray field of a commercial magnet to measure longitudinal relaxation rates over orders of magnitude of magnetic field. The sample is transfered back in the high-field magnetic center for detection, thus ensuring high sensitivity and resolution. During a high-resolution relaxometry experiment, the sample is moved outside of the magnetic center where no radiofrequency pulse can be applied. Thus, relaxation decays acquired using HRR suffer from two types of systematic errors. First, the effective density operator at the beginning of the relaxation delay is usually different from the desired longitudinal operator due to cross-relaxation during the sample transfers. Second, cross-relaxation pathways during the relaxation delay may lead to multi-exponential polarization decays. Therefore, the analysis of experimental HRR rates requires to account for these systematic deviations in order to accurately determine the motional parameters of the system under study. We introduced an iterative correction procedure called Iterative Correction for the Analysis of Relaxation Under Shuttling (ICARUS) [17,26] for the correction of HRR relaxation rates. Using symbolic expressions of magnetic-field dependent relaxation matrices, the HRR experiments are simulated and measured relaxometry relaxation rates are corrected so that a reliable analysis of the dynamic properties of the system under study can be performed. Thus, the development of tools to simulate spin relaxation for a broad variety of field trajectories is of great interest, in several areas of magnetic resonance [27][28][29]. Here, we present a toolbox that combines two programs. The first one, RedKite, provides analytical expressions of relaxation rates and relaxation matrices for arbitrary spin system. The second one, ICARUS, is used to retrieve accurate estimates of longitudinal relaxometry relaxation rates that are further used to determine the parameters describing the dynamics of the system. ICARUS simulates the HRR experiments using analytical expressions obtained from RedKite.
RedKite has been written in Mathematica (version 12.0) [30] to perform efficiently analytical calculations using the SpinDynamica (version 2.15.1b10) [29] package and the so called "BRW engine" to simplify the computation of relaxation rates [28]. This version of ICARUS has been written in Python (version 3.5). This language has the advantage of being free and easy to install, allowing for relatively fast numerical evaluations, and being easy to customize by the user. ICARUS is written as a framework so that users can define the spin systems, relaxation matrices and spectral density functions relevant for their applications. In this paper, we first describe succinctly our approach to calculate relaxation rates efficiently and apply this method on an isolated 15 N 1 H spin system using RedKite. We illustrate the power of these tools with a detailed presentation of the recently published analysis of carbon-13 HRR in { 13 C 1 H 2 H 2 }-methyl groups in the protein Ubiquitin [25] and test the validity of key hypotheses made during the analysis. In particular, we use two-field NMR to determine the relevant interactions to describe the relaxation properties of { 13 C 1 H 2 H 2 }-methyl groups, and verify the validy of the correction at 0.33 T.

Calculation of relaxation superoperators with RedKite
The full description of the Bloch-Wangsness-Redfield (BWR) relaxation theory in liquid-state NMR is beyond the scope of this article and can be found elsewhere [5,[31][32][33][34]. A condensed version is presented here. The evolution of the density operatorσ(t) is described by the Liouville-von Neumann equation, in units of : The HamiltonianĤ of the system can be expressed as the sum of a stationary partĤ 0 and a fluctuating partĤ 1 (t):Ĥ (t) =Ĥ 0 +Ĥ 1 (t).
This equation can be transformed in the interaction frame of the stationary Hamiltonian H 0 . An operatorÔ transformed into the interaction frame is labeled with a tilde: The frame transformation of the full HamiltonienĤ requires the subtraction of the Zeeman HamiltonienĤ 0 , so that the Liouville-von Neumann equation now reads: After developing a second-order time-dependent perturbation, the Liouville-von Neumann equation in the interaction frame can be written as: In the frame of the BWR theory, the following hypotheses are made to calculate the ensemble average of the evolution of the density operator: i) for an ensemble average, denoted by the horizontal bar, σ (0),H 1 (t) averages to zero, and ii) a time t can be found that is short enough such that the evolution of the spin system is negligible on the interval [0, t] but that is much larger than the typical correlation times for the fluctuations ofH 1 (t). The evolution of the density matrixσ(t) over time for an ensemble average, under a perturbation HamiltonianH 1 (t), can now be expressed as: This equation can be further simplified using the irreducible tensor representation in order to separate the angular and spin parts of the Hamiltonian. The perturbation Hamiltoniañ H 1 (t) may include several interactions, identified by the label i. Each of them can be written as the sum of the product of time-dependent spatial variables V l,−q (t) and tensor spin operatorsT l,q of rank l and coherence order q (which is usually simply called order): where ζ i is the amplitude of the interaction i. The irreducible tensorT i l,q can be expressed as a linear combination of eigenoperators {Â i l,q,p } of the superoperator [Ĥ 0 , ·], with eigenvalues ω (i) l,q,p :T i l,q = pÂ i l,q,p .
These eigenoperators can be written in the interaction frame as: A i l,q,p (t) = exp (iĤ 0 t)Â i l,q,p exp (−iĤ 0 t) = e iω (i) l,q,p tÂi l,q,p .
In the interaction frame, we now have: SinceH 1 is Hermitian, we can also write: where ( †) denotes the hermitian conjugate of the operator, and ( * ) the complex conjugate. Substituting Eq. 10 and 11 into Eq. 6 gives: The correlation function C i,j between the interations i and j is defined as: where δ is the Kronecker delta. Oscillating terms are neglected as they average to zero much faster than the evolution of the density operator (secular approximation) under relaxation. Thus, only secular terms for which ω l ,q ,p contribute to Eq. 12. Only rank-2 (l = 2) tensors are relevant to describe dipole-dipole and quadrupolar interactions. For the CSA interaction, the rank-1 tensor part (antisymmetric) is usually neglected. Note that, in the presence of highly anisotropic motions, the contribution of the antisymmetric CSA (rank-1 tensors) may account for up to 10 % of the contribution of the CSA rank-2 tensors to auto-relaxation [35,36]. In the following, only rank-2 tensors are considered. The spectral density function is defined as the Fourier tranform of the correlation function: Inserting the spectral density function in Eq. 12 and applying the above approximations leads to the following expression of the Master equation: The final step consists in transforming Eq. 6 from the interaction representation back to the Schrödinger representation given in Eq. 1. For this, we invert Eq. 3: with time-derivative: Inserting Eq. 15 into Eq. 17 leads to: We now define the relaxation super-operatorR as: The relaxation rate between operatorsÂ andB is: IfÂ =B, we speak of an auto-relaxation rate, whileÂ =B refers to a cross-relaxation rate, if i = j, it is an auto-correlated relaxation rate, and if i = j a cross-correlated relaxation rate. These rates can easily be calculated analytically using the BRW engine [28]. It consists in calculating the double commutator for each pair of spin tensors with identical eigenfrequencies and multiplying them by the spectral density function evaluated at this frequency. The implementation of this algorithm in Mathematica [30] is detailed for an isolated 15 N 1 H spin pair (Section 3.1) and a 13 C 1 H 2 H 2 methyl group with a vicinal deuterium (Supplementary Materials).

Expectation value of spin operators
The expectation value of a specific operator after an evolution period t is obtained from the calculation of the propagator:P withL the Liouvillian. Eq. 21 assumes a constant Liouvillian over the interval t, including a constant Hamiltonian. This assumption does not hold when pulses are applied, or in field-varying experiments, such as in dDNP or relaxometry. In dDNP, the sample is polarized using a microwave source at a specific field outside the NMR spectrometer, dissolved and pushed into the spectrometer, so that the sample experiences successively: the static field of the polarizer, the fields of the trajectory between the polarizer and the spectrometer, and the static field of the NMR spectrometer [10]. In a relaxometry experiment, the fields during the polarization, relaxation and detection periods are potentially all different [13,15]. In these cases, the evolution time t is decomposed in periods that are small enough so that the field can be considered constant, and the propagator equals: whered P i is the propagator during the interval δt i for which the magnetic field equals B i . When pulses are applied, which is typically the case in standard pulse sequences for the measurement of relaxation rates [37], cross-relaxation pathways may no longer be active and Eq. 22 can be simplified using averaged Liouvillian theory [38,39]. For example, for the measurement of longitudinal relaxation rates of nitrogen-15 in a 15 N-1 H spin pair, proton π-pulses are applied during the relaxation delay. In the abscence of such pulses, the Liouvillian reads:L where the relaxation matrix has been written in the basis formed by the spin operators {N z ,Ĥ z , 2N zĤz } and R N 1 (respectively R H 1 ) refers to nitrogen-15 (respectively proton) longitudinal relaxation rate, R NH to the two-spin order relaxation rate, σ NH to the dipole-dipole (DD) cross-relaxation rate between the nitrogen-15 and proton, and δ N (respectively δ H ) to the CSA-DD cross-correlated cross-relaxation rate involving the nitrogen-15 (respectively proton) CSA. After applying a proton π-pulse, the Liouvillian is transformed according to:L whereP π is the propagator for an ideal proton π-pulse: When the evolution delay before and after the pulse are equal, the proton inversion pulse After definition of the operator basis, Hamiltonian operators are defined. After indicating the operator of the basis studied during the experiment, a reduction of the size of the basis is performed. Rates of interest are defined as well. c) Calculations produce analytical expressions for the relaxation rates and the relaxation matrix. Blue rectangles: user inputs. Yellow rectangles: calculated outputs. Pink rectangle and purple triangles: tasks performed by RedKite.
leads to the following average Liouvillian over the whole relaxation period: Over this time period, the spin operatorN z is an eigenvector of the relaxation matrix, and the time-evolution of its expectation value is given by: which is the usual mono-exponential decay used for the analysis of relaxation rates measurements (note that the evolution towards an effective saturated state is obtained from the averaging of consecutive scans [38,40]). By constrast, an accurate analysis of relaxation properties in the abscence of radio-frequency pulses, or in field-varying experiments, requires the full relaxation matrix.

RedKite in Mathematica
The computation of the relaxation rates is highly efficient with the formalism of the BRW engine [28] which does not require an explicit expression of the Wigner matrices defining the correlation function (Eq. 13). Relaxation rates are first expressed as a function of the spectral density function J (ω, θ i , θ j ) where θ k is the orientation of the interaction k in the system frame (SF) of the chemical moiety. This frame corresponds to an arbitrary frame in which the orientation of the interactions are calculated. The different steps of RedKite are presented in the flowchart shown in Fig. 1. We will illustrate the use of RedKite on an isolated pair of spin-1/2 nuclei: a 15 N-1 H pair. We have used RedKite to analyse HRR data recorded on 13 C 1 H 2 H 2 specifically labelled isoleucine-δ1 methyl groups of the protein Ubiquitin [25], and to study the relaxation properties of 13 C 1 H 3 methyl groups during a HZQC experiments [41].

Definition of the spin system
The first step is to define the spin system by specifying for each nuclear spin the nucleus type with its isotopic number, and a unique label for each spin which is used for identification. We present as an illustration the example of a simple spin system composed of an isolated 15 N-1 H pair. The spin system is therefore defined as: Nuclei = {{"15N","NA"}, {"1H", "HA"}}; where "NA" and "HA" refer to the Nitrogen-15 and Proton respectively, before running the SpinDynamica [29] SetSpinSystem command: The NuclearSpinQuantumNumber command is implemented in SpinDynamica [29] and defines the quantum spin number of the considered nucleus. The geometry of the spin system is defined next. We define an array of size n × 3 (where n is the number of nuclei in the spin system, in our case 2) containing the position of each atom in a Cartesian axis system. In our example, we set the nitrogen nucleus at the origin of the axis system and the proton 1.02 Å away from the nitrogen in the z-direction: To complete the definition of the spin system, the Chemical Shift Anisotropy (CSA) and quadrupolar properties have to be defined. The nuclei for which the CSA will be considered must be defined as such. In our example, we will only consider the nitrogen CSA: It is possible to give a numerical value to the CSA or keep its value as an analytical parameter. We will consider this latter case here: Note that defining δ csa [2] is not necessary since the proton CSA is neglected. Similarly, the strength of the quadrupolar interaction does not need to be defined (see in the Supplementary Materials for an example that includes quadrupolar interactions).

9
The orientations of the CSA tensor have to be given (either numerically or analytically). For the sake of simplicity, we choose an alignement along the N-H axis: vectorNum "CSA" The index 1 refers to the first spin in the spin system (i.e. the nitrogen-15). There is also a possibility to consider asymmetric CSA tensors. In this case, the asymmetric CSA tensor is decomposed in two axially symmetric components.

Definition of spin tensors and Hamiltonian
Three different types of interactions are considered in RedKite: the dipolar couplings, the CSA (in the case where at least one spin has a CSA) and the quadrupolar couplings (in the case where spins with m s > 1/2 are present in the spin system). Analytical forms of these Hamiltonian operators are calculated automatically. Other Hamiltonian operators can be defined and added if other interactions or effects are considered. Calculation of Hamiltonian operators requires the definition of spin-tensor operators.
SpinDynamica already contains their definition, but each tensor of coherence order-q is given as a linear combination of eigentensors [29]. Consequently, SpinDynamica tensors can be linear combinations of eigenvectors with different eigenfrequencies, which is an inappropriate basis to perform the secular approximation (based on the equality of eigenfrequencies of two eigenvectors). The secular approximation is better performed with complete separation of the tensor operators. The definition of each tensor has already been reported for each considered interactions (dipole-dipole, CSA and quadrupolar) [42] and their definition in Mathematica can be found in Tables S3-S5. In the case of non-equivalent homonuclear spin systems, performing the secular approximation is more complex, especially at low fields, where the oscillation frequency in Eq. 12 can be comparable to the relaxation rates. Numerical tools, such as Spinach [28], are available to study such systems. The Hamiltonian, as written in RedKite, can be found in the Supplementary Materials. In the definition of the Hamiltonian, we introduce the function M, similarly to the BRW engine [28], which depends on the operator coherence order m being considered, its associated eigenfrequency, a time t at which the Hamiltonian is calculated, and the orientation of the interaction. The function M is useful when calculating the double commutators to obtain relaxation rates (as detailed in Section 2.1). Products of the function M appear, which are simplified according to: where KroneckerDelta[x, y] = 1 if x = y and 0 otherwise, l and k are associated to tensor coherence order, f1 and f2 to the tensor eigenfrequencies, t the time at which the Hamiltonian is calculated, and i and j are the orientation of the interactions in the molecular frame. G[t, f1, i, j] is the correlation function evaluated at time t and is further replaced by the spectral density function evaluated at frequency f 1. For auto-correlation, i = j, while cross-correlation is obtained when i = j.

Operator of interest
We define the operator of interest as the initial state where the polarization has been stored. In HRR, it is the longitudinal Zeeman term. In our case, we are interested in the nitrogen-15 longitudinal relaxation rates, which is defined by: where opI is a SpinDynamica [29] command to define operators, here the N z operator.

Analytical and numerical spin state restriction
The number of terms in the basis is equal to 4 n for n spin-1/2 nuclear spins. Hence, in this two-spin system there are 16 terms, which is still a workable number. For more complex spin systems, reducing the size of the basis is essential. We only keep the terms contributing to the relaxation of the operator of interest following the scheme of Fig. 2. First, only terms with the same coherence order as the operator of interest are selected (indicated in blue in Fig. 2a). Then, the secular approximation removes all non-secular terms in the interaction frame (Fig. 2b). Cross-relaxation rates with the operator of interest in this reduced basis are calculated (Fig. 2c) and the operators with no crossrelaxation with the operator of interest are discarded from the basis (here this last step only removes the identity operator E, Fig. 2d). This step is basis-dependent and some indirect cross-relaxation pathways affecting the operator of interest may be suppressed. An additional step can be applied for large spin systems to sort and select only major cross-relaxation pathways. In our example of an isolated 15 N-1 H spin pair with a CSA on the nitrogen-15, only 3 terms remain in the basis:

Calculations
Once the basis has been defined, the relaxation matrix can be calculated:  Figure 2: Reduction of the matrix size for our case example of a 15 N-1 H spin system. a) A 15 N-1 H isolated spin pair has 16 operators in its basis. b) The first step of the matrix reduction size consists in keeping only terms that have the same coherence order as the spin-term of interest, leading to 6 terms in the basis. c) The secular approximation allows another level of size reduction: only terms that are secular with the Zeeman Hamiltonian are kept in the basis. Two terms are removed at this stage. d) In the absence of cross-relaxation with the spin term of interest N z , the identity operator is removed from the basis and the final basis contains 3 operators. In this graphical representation of the relaxation matrices, a red square indicates a non-zero value for the corresponding relaxation rate. The blue rectangles contain the selected part of the relaxation matrix after each steps of the size reduction. Normalization factors for the spin operators have been omitted for clarity.
where R N 1 and R H 1 refer to the nitrogen-15 and proton longitudinal relaxation rates respectively, R NH to the auto-relaxation rate of the two-spin order, σ NH to the dipole-dipole cross-relaxation rate between nitrogen-15 and proton and δ N to the CSA-(dipole-dipole) cross-relaxation rate due to the cross-correlation of the nitrogen-15 CSA and the dipoledipole coupling: the dipolar coefficient between the proton and the nitrogen-15, r NH the distance separating the two nuclei, γ X the gyromagnetic ratio of nucleus X, the Plank constant devided by 2π, µ 0 the permeability of free space, and ∆σ N = σ zz − σxx+σyy 2 the CSA of the nitrogen-15 with σ kk the k th diagonal element of the chemical shift tensor. J is the spectral density function and is expressed as a function of the proton (ω H ) and nitrogen-15 (ω N ) Larmor frequencies.
All types of relaxation rates in this spin system can be calculated. In such a spin system, it is relatively easy to record longitudinal and transverse relaxation rates for the nitrogen-15 nucleus, as well as the cross-relaxation rate with the proton. These rates are calculated by: where Rate is the implemented command to calculate relaxation rates as described in the previous section. This leads to the expression of transverse relaxation rate for nitrogen-15:

Model selection and formating
The user has to provide at least one definition of spectral density function in order to have a model for the dynamics of the system. In our case, we can use a model-free approach [43] with a correlation time for global tumbling τ c , one order parameter S 2 and 13 an effective correlation time for internal motions τ int : This function is implemented in RedKite as: At this point, the relaxation rates seen above can be expressed as a function of the parameters of the dynamics of the system (order parameter and correlation times). Numerical calculations can be performed if values for the parameters of the spectral density function are provided.

Preparing for ICARUS
In order to use the results obtained in RedKite for the analysis of HRR, symbolic expressions have to be exported. Exporting to ICARUS requires that all variables have Latin-only characters as the interpretation of non-Latin characters is not implemented in ICARUS. During the export process, the spectral density function is provided by the user as: JofInterest = JNH; The user can export the first derivatives of the relaxation rates with respect to all the variables (magnetic field excluded as it is not useful in the following analysis). All the expressions of the relaxation matrix and the relaxation rates (and first derivatives if required) are saved in separate files named respectively RelaxationMatrix.txt for the entire relaxation matrix, Rate.txt for the relaxation rates defined in the RatesOfInterest array, and Ratederiv_Variable.txt where Rate refers to the considered relaxation rate and Variable to the variable name by which the rate is derivated. The first derivatives of the relaxation rates can be used in minimization procedures. An additional file named PositionOfInterest.txt is also created and contains the position of the operator of interest in the relaxation matrix (N z in our case example).

ICARUS implementation
In this paper, we show as an example how RedKite can be used for the analysis of HRR experiments. Other applications of RedKite have been published elsewhere [41,44], and can be envisioned, as relaxation rates can be obtained for any spin system. We detail here the analysis of HRR relaxation rates.

Accurate estimation of relaxation rates from high-resolution relaxometry measurements
High-resolution relaxometry can be used to obtain a precise description of the dynamics of spin systems over orders of magnitude of timescales [17,25,26]. The analysis is based on the measurement of longitudinal relaxation rates over a broad range of magnetic fields (typically from a few tenths of Tesla up to about 20 T). A reliable description of the motions requires accurate estimates of the relaxation rates. During each high-resolution relaxometry experiment, the sample is transferred outside of the magnetic center to a defined position z relax in the stray field above the magnet (Fig. 3). During the two transfers (from high to low field, and back) and the relaxation delay, all relaxation pathways are active. In contrast to the example presented in Section 2.2, measured polarization decays can be affected by cross-relaxation and therefore cannot be used as is to determine longitudinal relaxation rates accurately (this is true for any relaxation experiment where pulses can not be applied during the relaxation period). Doing so would lead to systematic deviations in the parameters used to describe the dynamics of the system. Simulating the experiment including the time when the sample is outside the superconducting coil allows one to take into account cross-relaxation pathways and to estimate accurate relaxation rates. The complete relaxation period in a high-resolution relaxometry experiment includes three delays at constant fields and two transfers through a strong gradient of magnetic field. The simulation of the experiment is performed by calculating the propagator during the highlighted part of the pulse sequence in Fig. 3b. For convenience, it is written as a product of individual propagators: whereP HF,1 andP HF,2 are the propagators calculated at high field, respectively before and after shuttling,P LF is the propagator calculated at the low field position andP up (re-spectivelyP down ) is the propagator calculated during the motion up (respectively down) from the high-field to the low-field position (respectively from the low-field to the highfield position). This decomposition allows the calculation of the segmental propagators using either Eq. 21 (for constant Liouvillian superoporators) or Eq. 22 (for time-dependent Liouvillian superoporators). The propagators for constant-field positions (i.e.P HF,1 ,P LF andP HF,2 ) are calculated using Eq. 21 and the relaxation matrix calculated at high field (R HF ) and low field (R LF ):P HF,i (t HF,i ) = e −t HF,iRHF , t a t a t a t a GARP g 1 g 1 g 2 g 3 g 3 g 4 g 5 g 5 Figure 3: Description of an HRR scheme. a) The position of the sample is changed during the course of the experiment. It is first polarized at high field, and transfered to a chosen position in the stray field of the superconducting magnet, characterized by a lower magnetic field, for relaxation. The sample is then moved back to the high-field position for detection. Pannel adapted from [25]. b) A typical pulse sequence used to record HRR experiment. During the analysis of HRR rates, the highlighted part of the pulse sequence (blue) is simulated. Black narrow (respectively wide empty) rectangles represent π/2-pulses (respectively π-pulses). Pulses are applied along the x-axis if not otherwise stated (by the ϕ i ). The amplitude of pulse field gradients are labeled g i . Additional experimental details can be found in [26].
The simulation of the transfers through the magnetic field gradient is performed by subdividing the experiment into intervals of few milli-seconds δt that still fulfill the conditions of Redfield theory. In order to stay in the Redfield hypothesis, δt must be large compared to the correlation time of the system to extend the integration to infinity in Eq. 5. In addition, δt must be sufficiently small in order to perform a discretization of the integral over the full sample trajectory. In the case of high-resolution relaxometry with a sample traveling at ≈10 m.s −1 over at most 1 m, we considered a δt of 1 ms, which corresponds, at most, to a change of about 10 % of the magnetic field between two consecutive steps. The propagators dP(δt, z(t)) for these small steps are obtained following Eq. 21: whereR(z(t)) is the relaxation matrix evaluated at the position z(t) along the bore of the magnet and characterized by its magnetic field (note: the field profile can be mapped using a gaussmeter). The experimental field profile is fitted to a polynomial expansion in ICARUS. Each propagator dP(δt, z(t)) is field dependent due to the field dependence of the relaxation matrix. The propagator for the motions up to and down from the position z relax are defined as the products of the infinitesimal propagators dP: dP up (δt, (z(n × δt))), where n up max (respectively n down max ) is defined by t up transfer = n max × δt (respectively t down transfer = n down max × δt) with t up transfer (respectively t down transfer ) the delay of transfer to the top (respectively down) position. In these calculations, the relaxation matrix is derived using the analytical expression obtained from RedKite, a model of motions and a set of parameters of dynamics. The expectation value for the operator of interest at the end of the full relaxation period (delays at high field and low field as well as the two transfers in between) can then be extracted from the calculated propagator for each relaxation delay. The simulated decay as a function of the relaxation time is fitted with a mono-exponential decay function with an effective longitudinal relaxation rate R sim (Table 1 sums up our nomenclature for the different calculated and measured relaxometry relaxation rates). All relaxation pathways are active during the transfers between high and low-field positions. The initial density operator is partially projected onto the eigenvectors of the relaxation matrix (relaxation modes) of lowest eigenvalues. Thus, the simulated decay rate R sim is a priori lower than the pure longitudinal relaxation rate R calc calculated using the parameters of dynamics. We define the correction factor for each relaxometry experiment as the ratio between these two rates for an experiment j (corresponding to a specific low field B (j) LF ) and a LF is the low field chosen for relaxation and D i are the parameters of the spectral density function used to describe the dynamics of residue i.

Label
Parameters Description where E j are the experimental parameters (shuttling times and relaxation delays), and D i are the parameters of dynamics. The correction factor is applied to each corresponding measured relaxometry data R meas (E j , B (j) LF ): The correction is performed iteratively (Fig. 4). The set of parameters D i for the first iteration is obtained from the analysis of the accurate relaxation rates, i.e measured with the use of pulses, typically on high-field magnets. Then corrected relaxometry relaxation rates are analyzed alongside high-field relaxation rates. A new set of parameters of dynamics is extracted from this ensemble of relaxation rates. In the next iteration, these parameters of dynamics are used to simulate the experiment and compute improved corrections of experimental rates to estimate the accurate low-field relaxation rates. This is repeated until the correction factors converge. The final set of high-field and corrected relaxometry relaxation rates can then be used to extract the distribution of the parameters of local motions in a Markov-Chain Monte-Carlo (MCMC) procedure and thus evaluate the median value and uncertainty of these parameters (see below).

Compiling expressions in the FunctionsFile.py script
Information about the relaxation properties of the spin system are contained in an independent script called FunctionsFile with expressions of the relaxation rates (and their  a) After a FunctionsFile has been obtained from RedKite, ICARUS can be run, using, among other inputs, relaxation rates recorded on standard high-field spectrometers and the high-resolution relaxometry data. Accurate relaxometry relaxation rates are obtained, and a Markov-Chain Monte-Carlo (MCMC) analysis of these corrected rates and highfield relaxation rates leads to values of parameters describing the dynamics of the system and their distribution. b) Flow chart of the ICARUS procedure. Accurate high-field (HF) relaxation rates are used to obtain an initial set of parameters for the dynamics of the system. These parameters are used to simulate the high-resolution relaxometry experiments (using the same experimental set up, i.e. shuttling time, delays, magnetic fields) from which biased simulated R 1 are extracted, and also to calculate the accurate expected R 1 . The ratios of these two calculated rates are called correction factors. The product of experimental decay rates and correction factors are corrected experimental low field (LF) relaxometry relaxation rates. Together with the high-field relaxation data, the corrected rates are used to determine a new set of parameters of dynamics, further used in the next correction iteration. Convergence is not evaluated within ICARUS and the number of iterations remains a choice of the user. However, we recommend to verify the convergence of the correction factors, as these ones are essential in the determination of the final parameters of the dynamics. Typically three or four iterations are sufficient. 19 derivatives if required) and the relaxation matrix in the considered basis. The Functions-File can be edited and adapted to the spin system under investigation. RedKite outputs first need to be converted from Mathematica to Python format and compiled in this FunctionsFile.py script. This task is performed by the RedKite2ICARUS.py program. Briefly, it takes as input all the output files from RedKite (Section 3.1) and asks for variables names (the ones that have to be fitted, usually parameters defining the spectral density function) and the ones that characterize the system and are not fitted (e.g., CSA tensors). It is also possible to set the CSA as a fitted variable. In the case where the overall diffusion frame is asymmetric, ICARUS requires a file containing the orientations of internuclear vectors in the anisotropic diffusion frame. Creating such a file has been implemented in RedKite2ICARUS.

Fitting parameters of the model of motion to relaxation rates
The program ICARUS (Iterative Correction for the Analysis of Relaxation Under Shuttling) [17,26] has been entirely written in Python (version 3.5). The detailed description on how to use ICARUS has been already published elsewhere [26]. The key parts of the code are the fitting of parameters of a user-defined model of motion using accurate (generally high field) relaxation rates and corrected relaxometry rates as experimental constraints, as well as the simulation of the experiments (as detailed in Section 3.2.1).
Fitting the parameters of the model relies on the basin-hopping function implemented in the scipy.optimize Python library with the L-BFGS-B method [45] for χ 2 minimization: where R model,i are the calculated relaxation rates and R exp,i are the measured relaxation rates with experimental error σ exp,i . Bounds of the dynamics parameters are provided by the user in the GUI. The basinhopping function allows the use of first derivatives of the relaxation rates in the fitting (provided in the FunctionsFile as explained above), usually leading to faster minimization. An additional minimization based on a grid search has been implemented in order to avoid local minimum traps. This step is time-consuming and optional. The core of the code does not contain information about a particular spin system nor experimental set up, such that the usage of ICARUS can be extended to any situation (spin system or model of motion). Data and experimental set up are loaded as separate text files using the GUI, and all analytical expressions of relaxation rates are contained in the independent FunctionsFile script.

ICARUS output and MCMC
The outputs of ICARUS have already been described [26]. Fig. 5 shows selected figures created by ICARUS. Briefly, output figures consist of the fit of the stray field gradient (Fig. 5a), profiles of the relaxation rates (accurate, calculated, and corrected in the case of relaxometry data) at each field throughout the protein sequence (Fig. 5b, c, d), fits of all the relaxation rates for each residue (Fig. 5e, f), bar plots of fitted parameters (Fig. 5g). Several text files are created which contain corrected relaxometry relaxation rates, the set of fitted parameters after the fit of the accurate relaxation rates only, and of the whole data set (accurate and corrected relaxometry data) as well as the correction factors after each iteration of ICARUS. Finally, scripts are also created. One allows the user to calculate all the defined relaxation rates and the relaxation matrix using the final fitted parameters with the use of a GUI where the magnetic field and residue number of interest have to be set. The other scripts are created only if a PDB ID for the protein of interest has been provided in the GUI, and are meant to be run in PyMOL in order to color the structure according to the final set of fitted parameters (order parameters, correlation times, etc...) and the final χ 2 (Eq. 34) to facilitate the visualisation of the results over the protein structure (one file is created for each of these parameters). An example is shown in Fig. 5h. In order to provide a better analysis of the dynamics, a Markov-Chain Monte-Carlo (MCMC) analysis of accurate and corrected relaxometry relaxation rates should be performed. We have written a script that directly reads ICARUS output folders to perform an MCMC using the emcee Python library [46]. The MCMC analysis provides a better error evaluation of the parameters of dynamics as well as potential correlations between them. A README file explaining how to use the MCMC program is provided with the script. Overall, the RedKite-ICARUS suite is intended to allow for an efficient (a complete analysis of 15 N relaxometry data on Ubiquitin can be obtained within two hours on a standard laptop computer) and highly flexible (it can be extended to broad range of spin systems, with all types of model of motions and for most commonly measured relaxation rates) analysis of high-resolution relaxometry data. RedKite and RedKite2ICARUS create the scaffold (FunctionsFile) that is used by ICARUS. The use of ICARUS is convenient with a simple graphical user interface. After the correction of the relaxometry relaxation rate, a final Markov-Chain Monte-Carlo analysis is performed by a script that reads directly ICARUS output folders (Fig. 4a). 22 4 Application to { 13 C 1 H 2 H 2 }-methyl groups using HRR and 2F-NMR Motions of protein side-chains are important for their function. These motions have been investigated thanks to NMR methodological development and selective labeling strategies based on the clever use of metabolic pathways [47][48][49]. The averaging of the dipolar interactions arising from their fast rotation confers favourable relaxation properties to methyl groups. They make good candidates for the study of side-chain motions, in particular in the hydrophobic core of proteins where they constitute an entropy reservoir [50,51], or at protein-protein and protein-ligand binding interfaces where their motions can allow a re-modeling for a better complementary interaction with the binding partner. In this context, we have recently performed a detailed analysis of the motions of isoleucine-δ 1 methylgroup on the selectively labeled protein U-[ 2 H, 15

Model of correlation function
Different models of correlation function for a wide variety of molecular systems have been suggested in the past [43,[52][53][54][55][56][57][58]. In our analysis of high field and relaxometry relaxation rates on { 13 C 1 H 2 H 2 }-methyl group of Ubiquitin, the data recorded at low fields (lower than 5 T) allowed a better characterization of the complexity of motions that can occur in a methyl-bearing side-chain, in particular χ 1 /χ 2 rotameric transitions in isoleucine residues on nanosecond timescales [25]. The analysis was based on the Extended Model Free (EMF) description of the CC bond motions. Assuming (i) isotropic tumbling of the protein characterized by a correlation time τ c , (ii) EMF for CC bonds motions, (iii) perfect tetrahedral symmetry for the methyl group with a characteristic correlation time for the methyl group rotation τ met associated to an order parameter S 2 met (θ i,j ) [59] and (iv) statistical independence between methyl group rotation, motions of the methyl group axis and overall rotational diffusion, the correlation function can be modeled by: where: with S 2 met (θ i,j ) = P 2 (cos θ i )P 2 (cos θ j ) and P 2 is the second order Legendre polynomial function, P 2 (x) = (3x 2 − 1)/2, θ k is the angle between the principal axis of an axially symmetric interaction k vector and the CC-axis (methyl group symmetry axis) and θ i,j the angle between the principal axes of two (possibly identical) axially symmetric interactions i and j. The order parameters S 2 f and S 2 s characterize motions of the system frame and are associated with the correlation times τ f and τ s , respectively. The overall order parameter is defined as S 2 = S 2 f S 2 s . The value of the angles θ k and θ i,j is constrained by the geometry of the spin system. The corresponding spectral density function is: where τ a −1 = τ −1 a + τ −1 c and τ a −1 = τ −1 a + τ −1 c + τ −1 met . In the following, J AB will be used to denote the dipole-dipole auto-correlation between nuclei A and B, J A for the CSA auto-correlation of nucleus A, J AB,CD for the dipoledipole/dipole-dipole cross-correlation between the spin pairs AB and CD, J A,BC for the cross-correlation between the CSA of nucleus A and the dipole-dipole interaction between nuclei B and C. Finally, the index Q will be used to denote the quadrupolar interactions. These notations follow conventions proposed by Werbelow and Grant [60]. As detailed bellow, in our treatment of the relaxometry data, the effects of the surrounding deuterium nuclei arising from the labelling of the protein have to be considered. These have been taken into account by adding a single additional deuterium nucleus in the spin system. For simplicity, while we consider the additional dipolar contributions to relaxation rates of the { 13 C 1 H 2 H 2 } spin system, we do not include this additional nucleus in our basis. We approximated the spectral density function for the correlations involving this vicinal deuterium D vic to be described by Eq. 37, although it is not part of the methyl group.

Relaxation rates
In our analysis of high-field and relaxometry relaxation rates on { 13 C 1 H 2 H 2 }-methyl groups of Ubiquitin, longitudinal and transverse carbon-13 autorelaxation rates, longitudinal proton autorelaxation rates and dipolar cross-relaxation rates were used. Dipolar relaxation with an effective vicinal deuterium was considered. The set-up of RedKite for such a spin system is detailed in Supplementary Materials. The contribution of the proton CSA to relaxation is expected to be negligible [61], and is not considered in the following. The CSA tensor of the carbon-13 nucleus is assumed to be symmetric and aligned with the CC bond. Expressions of the relaxation rates are given in the following equations: where d AB is the dipolar coefficient between atoms A and B and equals −(µ 0 γ A γ B )/(4πr 3 AB ) with µ 0 the permeability of free space, the Planck's constant divided by 2π, γ X the gyromagnetic ratio of nucleus X and r AB the internuclear distance between nuclei A and B, ∆σ C is the chemical shift anisotropy of the carbon-13 nucleus and ω X = −γ X B 0 is the Larmor frequency for the nuclei X at a magnetic field B 0 . The geometry of the methyl group was assumed to be tetrahedral with r CH = r CD = 111.5 pm leading to r HD = 182 pm. The distance r CD vic is determined during the ICARUS analysis as described below. 25

Relaxation matrix
The secularized basis for the subspace that includesĈ z in a { 13 C 1 H 2 H 2 }-methyl group contains 14 terms: where C, H, D 1 and D 2 refer to the carbon, proton, deuterium 1 and deuterium 2, respectively, as defined in the spin system in RedKite. The deuterium 1 and 2 are considered magnetically equivalent and can be exchanged by symmetry (see Fig. 7c for a visualisization of the geometry of the system).
As shown below, the analysis of the relaxation properties of the { 13 C 1 H 2 H 2 }-methyl groups of Ubiquitin during a relaxometry experiment can be performed with satisfactory accuracy in the subspace spanned by the three operators: leading to the following relaxation matrix: where R 1 ( 13 C), R 1 ( 1 H) and σ CH are defined above and: The expression of the secularized relaxation matrix can be found in the Supplementary Materials.

Size of the relaxation matrix
The ICARUS protocol aims at obtaining accurate estimates of low-field relaxation rates by accounting for the effects of cross-relaxation on the longitudinal relaxation decays during a high-resolution relaxometry experiment. This estimate is based on the simulation of the relaxometry experiments, where the sample travels through a broad range of magnetic fields. In order to obtain a reliable description of relaxation over orders of magnitude of magnetic fields, simulations must use appropriate relaxation matrices as well as expressions of relaxation rates, with accurate parameters for the amplitudes of interactions and the description of the spectral density function. The full Liouville space for a { 13 C 1 H 2 H 2 } spin system is spanned by a large basis of (2 × 1 2 + 1) 2×n 1/2 × (2 × 1 + 1) 2×n 1 = 1296 spin terms, with n 1/2 and n 1 the number of spin-half and spin-one respectively (Fig. 6a). An efficient calculation requires to minimize the size of the Liouville space where the evolution of the density operator is calculated. We have reduced the size of the subspace using the steps described in Section 3.1 for 15 N-1 H spin systems. First, we have considered the subspace only spanned by zero-quantum coherences and population operators (Fig. 6b). We then applied the secular approximation, and calculated all cross-relaxation terms with theĈ z operator, in order to keep only non zero terms, i.e. terms that cross-relax withĈ z , reducing the size of the basis to 14 terms (Fig. 6d). Cross-relaxation and autorelaxation rates in this 14-element basis have been calculated at the lowest and highest magnetic fields used during our HRR experiments, i.e 0.33 T and 14.1 T, using parameters obtained after a preliminary ICARUS analysis (for Ile-3) performed using B reduced,3 (Eq. 40, Fig. 6e). The inspection of these two relaxation matrices justifies the use of a basis containing only 3 operators as cross-relaxation rates involving other operators are either negligible (cross relaxation fromĈ z to another operator can be neglected if the ratio of this cross-relaxation rate to the auto-relaxation rate ofĈ z is small) or involve an operator with an auto-relaxation rate much larger than the auto-relaxation rate ofĈ z and the cross-relaxation rate withĈ z (see the Supplementary Materials for the proof that cross-relaxation with fast relaxing operator do not contribute to the polarization decay of slowly relaxing operators). At both magnetic fields, the largest cross-relaxation rate with the carbon-13 longitudinal polarization is the dipolar cross-relaxation with the proton longitudinal polarization. At low magnetic field (0.33 T), even a 2-operator basis { 1 3Ĉ z , 1 3Ĥ z } would be sufficient to describe the relaxation properties of a { 13 C 1 H 2 H 2 }methyl group as cross-relaxation towards other terms is either very small or towards fast-relaxing terms. However, the subspace should include the two-spin order 2Ĉ zĤz at high field (14.1 T). Thus, high-resolution relaxometry experiments in { 13 C 1 H 2 H 2 }-methyl groups have been simulated in the small subspace spanned by the three operators (Ĉ z , H z and 2Ĉ zĤz ). This subspace was used throughout our analysis of carbon-13 HRR in

Proton relaxation and surrounding deuterium
Proton longitudinal relaxation rates R 1 ( 1 H) were measured at three magnetic fields (0.33, 14.1 and 18.8 T) using standard high-field magnets (18.8 T and 14.1 T) and a 2F-NMR spectrometer operating at 14.1 T and 0.33 T [21]. These rates were also calculated after an ICARUS analysis of high-field and HRR rates considering intra-methyl group interactions only. The predicted relaxation rates are systematically lower than those measured at 0.33 T, 14.1 T and 18.8 T (Fig. 7a, b). Thus, even if relaxation rates in a { 13 C 1 H 2 H 2 }methyl group are dominated by the contributions of internal interactions, another contribution to relaxation has to be taken into account to describe proton relaxation. The differences between the measured and calculated R 1 ( 1 H) rates were assigned to the effect of the neighbouring deuterium nuclei. Adding the dipolar interactions with surrounding deuterium nuclei leads to non-negligible contributions to relaxation to both the proton and the carbon-13. The closest neighbouring deuterium nuclei are the 2 Hγ 1 and 2 Hγ 2 sites of the isoleucine side-chain, but other deuterium nuclei may also be in close proximity to the methyl group especially within the hydrophobic core of the protein. The correlation function for the fluctuations of the corresponding internuclear vectors are expected to vary. In particular, these interactions are expected to be affected in different ways by the fast rotation of the methyl group. We modeled the surrounding deuterium nuclei by a single deuterium at an effective distance (Fig. 7c). The interaction of the proton and carbon-13 nuclei of the methyl group with this deuterium accounts for the interaction with all the other deuterium nuclei of the protein. We used two adjustable parameters to describe its position, defining its coordinates in the Cartesian axis system: the y-and z-coordinate were fitted while the x-coordinate was fixed to 0. The position of the effective surrounding deuterium nucleus is determined independently for each residue using proton relaxation rates as well as all relaxation rates used in the ICARUS iterations (accurate and corrected) and keeping the other parameters constant (i.e. the parameters describing the dynamics). When fitting the parameters of the model during further ICARUS analysis, the effective position of the surrounding deuterium is kept constant. Introducing the contribution of the surrounding deuterium and performing the whole ICARUS analysis again preserves the agreement between the measured and calculated proton longitudinal relaxation rates (Fig. 7a, b). The surrounding deuterium has an effect on the correction factors (Fig. 7d) which leads to differences of corrected HRR rates between 0 and 4 % (Fig. 7e). Correction factors depend on the magnetic field and generally increase with decreasing magnetic. It must be pointed out that non-monotonous changes in the correction factors profiles in Fig. 7d are due to differences in shuttling and waiting delays at low magnetic fields (Fig. S2). The effective distances with the surrounding deuterium nucleus are close to extracted distances from the NMR structure of Ubiquitin (Fig. 7f, PDB 1D3Z). The dipolar inter-d) Reduced    The black line is shown as a guide for perfect equality between the two rates. c) Geometry of the methyl group and position of the effective neighbouring deuterium. The distance r C−D vic = r 2 y,D vic + r 2 z,D vic is determined using additional relaxation rates as explained in the main text. d) Correction factors as a function of the magnetic field for Ile-30 and Ile-44 with and without an effective vicinal deuterium nucleus. e) Corrected relaxometry relaxation rates for Ile-30 and Ile-44 with and without including an effective vicinal deuterium nucleus. f ) Comparison of the distance of the vicinal deuterium with the carbon-13 nucleus obtained from the analysis of proton relaxation (red, ICARUS) to the calculated distance to an effective deuterium nucleus that accounts for either only the 2 Hγ 1 and 2 Hγ 2 nuclei of the isoleucine residue (green) or all the hydrogens (blue) in the structure of Ubiquitin (PDB ID: 1D3Z). In these NMR derived structures, the distances were averaged over the 10 models present in the PDB file. In each model, the distance equals r C−D vic =  action between the methyl group and the effective deuterium is included in the following iterations of the ICARUS analysis.

Convergence of the iterative correction
The number of iteration steps is expected to be dependent on the spin system under study. In the case of the { 13 C 1 H 2 H 2 }-spin system, the convergence was reached after 2 iterations (Fig. 8a) for all residues except residue 44. Some slight instability in the convergence of the correction at low field is observed for this residue (Fig. 8b) but the amplitude of change (1-2 % at most) has a negligible effect on the values of the corrected relaxation rates (Fig. 8c).

Influence of the model of spectral density function on the correction
Different models can be used to describe the motions in a methyl group. Eq. 37 gives a rather complex description of the motion, but a simpler model can be tested by reducing the number of internal dynamics parameters to 3 by only considering the global tumbling, the methyl-group rotation with one fitted correlation time and C-C axis motions with only one fitted correlation time and one order parameter. The spectral density function for this model is: with the same definitions as above and where τ int is an internal correlation time for the motion of the C-C axis. Correction factors obtained for the two spectral density functions are shown in Fig. 8d and e. They are identical for Ile-30 where both models fit the experimental data well. In contrast, the correction is slightly different for the two models of motion for Ile-44 (Fig. 8e), where the 5-parameters model is in better agreement with the experiments. Yet, the variation on the corrected rates is small (between 1 and 2 %, Fig. 8f) with equally small effects on the analysis. The ICARUS analysis requires a model that accounts for the overall changes of the spectral density function on the range of frequencies probed during the experiments but it does not require that the used model reproduces all subtle details of the spectral density function: small variations of the value of the spectral density function at a specific frequency have negligible effects on the correction.

Scaling of the CSA/dipole-dipole cross-correlated cross-relaxation rates
Our combined analysis of low-field longitudinal and high-field transverse relaxation has allowed us to obtain the value of the CSA for each residue in addition to parameters of internal motions, except for Ile-44 for which chemical exchange prevented the analysis of the carbon-13 transverse relaxation rates [25]. In order to validate our analysis, a series of relaxation rates were measured as detailed hereafter: accurate low field carbon longitudinal relaxation rates [21] as well as high-field longitudinal CSA/dipole-dipole (CSA/DD) cross-correlated cross-relaxation rates (cross-relaxation betweenĈ z and 2Ĉ zĤz refered to as η C z ). These relaxation rates were not used during the analysis of the relaxometry relaxation rates, but calculated using the set of motional parameters obtained after correction of the relaxometry data. The calculated longitudinal CSA/DD cross-relaxation rates were strongly correlated to measurements at 14.1 T and 18.8 T but significantely overestimated (Fig. 9a). In order to have a better description of the CSA/DD cross-correlation, a scaling factor was applied directly to this term in the relaxation matrix. The scaling factor was calculated as the averaged inverse correlation coefficient between the unscaled and measured η C z at 14.1 T and 18.8 T and equals 0.505. A number of hypothesis can be made to explain the origin of the scaling factor: i) the carbon-13 CSA may be overestimated since it is determined essentially from transverse relaxation rates, which may suffer from small chemical exchange contributions; ii) the carbon-13 CSA may not be perfectly alligned with the C-C bond; iii) the form of the spectral density function may not describe correctly the motions of the methyl group; iv) the amplitude of the carbon-13 CSA may be rotamer-dependent. To understand the origin of this scaling factor, we also measured the carbon transverse CSA/DD cross-correlated cross-relaxation rates (η C xy ). The calculated relaxation rates correlate with the measurement, with an averaged inverse correlation coefficient between the calculated and measured η C xy at 14.1 T and 18.8 T of 0.629 (Supplementary Materials  Fig. S3). The discrepency between the scaling factors of the longitudinal and transverse CSA/DD cross-correlated cross-relaxation rates can not be accounted for only from a miss-evaluation of the carbon-13 CSA (under our assumptions of axially symmetry and perfect alignment allong the CC bond). Thus, it is likely that the model of correlation function does not describe entirely the complexity of the motions in the methyl group, and additional work toward this direction has to be done. For example, transitions between rotamers may be better modelled with instantateous jumps. The analysis of the relaxometry relaxation data was performed again after applying the scaling factor to longitudinal CSA/DD relaxation rates. As expected, the agreement between calculated and measured CSA/DD cross-relaxation rates is significantly improved by the use of a scaling factor (Fig. 9a). Low-field correction factors are not sensitive to the scaling of a CSA-dependent relaxation rate (Fig. 9b). At moderate and high field, the effect is larger with a reduction of the correction by about 2 % which has limited impact on the analysis.

Validation of the correction with the suppression of crossrelaxation pathways
Using the recently developed 2F-NMR spectrometer [19,20], we measured, among other relaxation rates, the longitudinal carbon-13 relaxation rates at 0.33 T with suppression of cross-relaxation pathways [21]. The rates of the seven isoleucines acquired at 0.33 T have been compared to measured and corrected relaxometry relaxation rates at the same magnetic field (Fig. 10a). The uncorrected relaxometry rates R 1 ( 13 C) are systematically lower than the accurate relaxation rates. This stresses the fact that the relaxometry relaxation rates have to be corrected in order to reach a reliable analysis of the properties the dynamics of the system. Corrected rates are in excellent agreement with the accurate R 1 ( 13 C) rates measured with the two-field system. This comparison validates the ICARUS approach on this spin system. In addition, experiments have been recorded at 14.1 T with and without pulses during the relaxation delay. Corresponding relaxation rates are displayed in Fig. 10b. The high-field experiment recorded without control of crossrelaxation pathways is similar to a shuttling experiments. Correction factors seem to be   Figure 10: Validation of the correction protocol. a) Correlation plot between the relaxometry uncorrected (blue) and corrected (orange) carbon R 1 at 0.33 T with the measured two-field R 1 ( 13 C). b) Correlation plot between the pseudo-relaxometry uncorrected (blue) and corrected (orange) R 1 ( 13 C) with the accurate relaxation rates measured at 14.1 T.
The black line is shown as a guide for perfect equality between the two rates.

Conclusion
In this paper, we have presented a general framework for the analysis of high-resolution relaxometry data. First, RedKite is a powerful Mathematica notebook to calculate relaxation rates and entire relaxation matrices in any nuclear spin system. We have shown how it can be used for the analysis of HRR, but it can also be applied more generally for the study of relaxation properties. Second, ICARUS is a Python-based program designed to analyze relaxometry datasets accounting for the effects of multiple cross-relaxation pathways. The two toolkits have been developed in order to be easily adapted to other spin systems, diffusion tensors and models of motions. Conclusions drawn here in the case of a { 13 C 1 H 2 H 2 }-methyl group with respect to the effect of the size of the relaxation matrix, the number of iteration of ICARUS or the model for the spectral density function may be different in other systems. Overall, a complete analysis by RedKite and ICARUS can be performed quickly, allowing one to evaluate these effects efficiently. Our approach to correct high-resolution relaxometry data has been cross-validated by the measurements of accurate low-field relaxation rates.

Materials and methods
Methods to obtain carbon-13 and proton longitudinal relaxation rates at 0.33 T were previously described [21] and are based on the use of a two-field spectrometer operating at 14.1 T and 0.33 T [19,20]. Proton longitudinal relaxation rates at 14. The longitudinal and transverse cross-correlated cross-relaxation rates (η C z and η C xy ) were measured using the symmetrical reconversion principle [40,62]. For enhanced sensitivity, cross-relaxation experiments were accumulated with 8-times more scans than autorelaxation experiments. The longitudinal cross-correlated cross-relaxation rate at 18.8 T was determined with a relaxation delay of 1.5 s, while at 14.1 T the experiment was performed with the relaxation delays of 1.0, 1.5, and 2.0 s. The measurement of the transverse cross-correlated cross relaxation rate was done using a spin lock irradiation with amplitudes of 2031 and 2062 Hz at 14.1 and 18.8 T, respectively. The alignment of the spins into the direction of the spin-lock field and back to z-direction was achieved using adiabatic half passage pulses. The calibration of the spin lock rf amplitude was done by measuring the scaling of scalar couplings under off-resonance continuous wave irradiation. The transverse cross-correlated cross relaxation rate at 18.8 T was determined from a single experiment performed with the relaxation delay 250 ms, while the experiment was repeated twice with the relaxation delays 175 and 250 ms at 14.1 T. The measurement of the "relaxometry-like" relaxation rate at 14.1 T was performed with the standard pulse program to measure longitudinal relaxation rates [37], but all pulses usually applied during the relaxation period were omitted. The experiment was measured twice, first with the relaxation delays 0.06*, 0. 18  1 Size-reduction of relaxation matrices by removing

fast-relaxing operators
Here, we will show that fast-relaxing terms of a relaxation matrix can be discarded (as done in Section 4.2.5 of the main text) in order to reduce the size of the relaxation matrix and save computational time. For the sake of simplicity, we consider a 2x2 Liouvillian: The characteristic polynomial of L is: with I the identity matrix. The roots are given by: with: Let's assume R 1 R 1 , σ. A first order approximation in R 1 and σ of √ ∆ leads to: such that the eigenvalues of L are R 1 and R 1 . The associated eigenvectors approximate to {1, 0} and {0, 1} and the autorelaxation of the operator of interest can be considered mono-exponential with decay rate of R 1 . The fast relaxing operator does not contribute to the relaxation of the slowly relaxing operator. This can be verified by simulating the polarization decay. We will set R 1 = 1 s −1 , σ = 0.5 s −1 and vary R 1 . We can compute the polarization decay (associated with the operator of interest with autorelaxation rate R 1 ) following Section 2.2 of the main text (Fig. S1). The polarization decay can be fitted to a mono-exponential decay, and fitted relaxation rates are reported in Table S1. It is clear that the fast relaxing operator has negligeable effects on the polarization decay when R 1 R 1 .

Correlation functions and spectral density functions
The choice of the model of motions is a key step in the analysis of relaxation rates to characterize quantitatively protein dynamics. The description of models of correlation functions can be found elsewhere [43,[52][53][54][55][56][57]. Any analytical form of the spectral density function can be used in RedKite and ICARUS. Assuming that different types of motions are statistically independent, the overall correlation function C i,j associated to auto-or  Figure S1: Simulated polarization decay (plain) and exponential fit (dash) for different values of R 1 relaxation rates. cross-correlation of interaction(s) (i, j ) can be written as the product of the correlation function of overall rotation C g , assumed here to be isotropic, and of the individual motions C n i,j , all supposed to be independent and isotropic: In model-free approaches, the overall rotation correlation function C i,j is described by a single exponential decay for isotropic diffusion, or a sum of exponentials for axially symmetric or fully anisotropic rotational diffusion [63]. The correlation function used for the model-free C MF i,j and extended model-free C EMF i,j approaches are: where θ i,j is the angle between the principal axes of the two interactions, P 2 (x) is the second order Legendre polynomial P 2 (x) = (3x 2 − 1)/2, τ g the correlation time for the global tumbling. The correlation function for the model-free approach is defined by the effective correlation time τ int and the order parameter S 2 . In the extended modelfree correlation function, τ s (respectively τ f ) is the correlation time associated with the order parameter S 2 s (respectively S 2 f ) for the slower (respectively faster) motion. The corresponding spectral density functions J MF i,j (ω) and J EMF i,j (ω) can be used for both auto-and cross-correlation of interactions: where τ a is the effective correlation time defined as τ a −1 = τ −1 a + τ −1 g . Other correlation functions can be used depending on the system under study. For example, the correlation function can be written as a sum of exponential functions: The corresponding spectral density is: In the case of relaxation in a methyl group, assuming the statistical independence of the methyl group rotation, the motions of the methyl group axis and the overall rotational diffusion, the correlation function C met i,j can be expressed as the product of the three corresponding correlation functions: C g for the global tumbling, C i,j rot for the methyl group rotation, C axis for the complex motions of the methyl group. The correlation function was given in the main text (Eq. 36). The rotation of the methyl group is an anisotropic motion characterized by the correlation time τ met and the order parameter S 2 met (θ i,j ) imposed by the geometry of the methyl group (supposed to be a tetrahedron, three corners of which are occupied by the proton and the two deuterium nuclei and the center by the carbon-13) and the relative orientations of the principal axes of interactions i and j with respect to the methyl axis. Motions of the methyl group axis are described by an extended model-free correlation function, with the parameters S 2 f , τ f , S 2 s , and τ s , as is detailed in the main text. 44 3 Set up of RedKite for the { 13 C 1 H 2 H 2 }-methyl groups of Ubiquitin with a vicinal deuterium Here, we show the most important command lines used to calculate relaxation rates and relaxation matrix of a { 13 C 1 H 2 H 2 }-methyl group with a vicinal deuterium nucleus. The carbon-13 is set at the origin of the Cartesian axis system, the 1 H is in the Oyz plan, as is the vicinal deuterium, which position is determined by two unknown (later optimized) variables describing its position along axes Oy and Oz (ryCD and rzCD, respectively). The two deuterium nuclei of the methyl group are mirror image of one another with respect to the Oyz plane. We define a System Frame with z-axis along the symmetry axis of the methyl group, i.e. the Oz axis:

Definition of the spin system
The orientation of the interactions relative to the System Frame is important when studyng the dynamics of the methyl groups, in particular their rotation around the symmetry axis, and are used in the definition of the spectral density function (see main text). We only consider the CSA for the carbon-13 nucleus, assumed to be axially symmetric: with value CSAValue which will be a variable optimized during the analysis of relaxation data: and we define the orientations of the considered quadrupolar interactions: vectorNum "Quad" where the command Vec extracts the vector between the two entries (the two nuclei). In the following analytical expressions of relaxation rates, the intensity of the quadrupolar interaction will be labelled ζ Q .

Spectral density function
We used the same spectral density function written in Eq. 37 of the main text. We assumed the vicinal deuterium nucleus follows the same model of spectral density function, even if it is not sensitive to the rotation of the methyl group as the 13 C, 1 H and deuterium nuclei are. Note that the two parameters used to position the effective vicinal deuterium nucleus change the effect of the methyl group rotation on relative correlation functions.

Experiment number
Magnetic field (T) Figure S2: Experimental delays for the 25 experiments used in the analysis of the dynamics of isoleucine-δ1-methyl groups of Ubiquitin, and ordered from the highest magnetic field at which relaxation takes place to the lowest. The time labels refer to the decomposition of the free-relaxation part of the pulse-sequence, as shown in Fig. 3 of the main text. The blue curve (right y-axis) shows the variation of the magnetic field for each experiment (associated with an increase of shuttling height). Experiments 1, 2 and 4 were performed on high-field spectrometers, with no shuttle.

Hamiltonian in RedKite
We report here the definition of the Hamiltonian as written in RedKite. Constants are defined in Table S2