# Explicit models of motions to analyze NMR relaxation data in proteins

Nicolas Bolik-Coulon  Laboratoire des Biomolécules, LBM, Département de
Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS,
75005 Paris, France  Fabien Ferrage  Laboratoire des Biomolécules, LBM,
Département de Chimie, École Normale Supérieure, PSL University, Sorbonne
Université, CNRS, 75005 Paris, France

###### Abstract

Nuclear Magnetic Resonance (NMR) is a tool of choice to characterize molecular
motions. In biological macromolecules, pico- to nano-second motions, in
particular, can be probed by nuclear spin relaxation rates which depend on the
time fluctuations of the orientations of spin interaction frames. For the past
40 years, relaxation rates have been successfully analyzed using the Model
Free (MF) approach which makes no assumption on the nature of motions and
reports on the effective amplitude and time-scale of the motions. However,
obtaining a mechanistic picture of motions from this type of analysis is
difficult at best, unless complemented with molecular dynamics (MD)
simulations. In spite of their limited accuracy, such simulations can be used
to obtain the information necessary to build explicit models of motions
designed to analyze NMR relaxation data. Here, we present how to build such
models, suited in particular to describe motions of methyl-bearing protein
side-chains and compare them with the MF approach. We show on synthetic data
that explicit models of motions are more robust in the presence of rotamer
jumps which dominate the relaxation in methyl groups of protein side-chains.
We expect this work to motivate the use of explicit models of motion to
analyze MD and NMR data.

##  1 Introduction

Proteins are not rigid molecules and are better represented by ensembles of
accessible conformations.[1] The distributions of conformations and timescales
of transitions between them report on the thermodynamics and kinetics of
biomolecular processes. Nuclear Magnetic Resonance (NMR) is particularly
suited for the characterization of dynamic properties at the atomic level with
a large toolset to probe timescales ranging from picoseconds to seconds and
more [2, 3]. Fast local motions of protein backbone and side-chains, on the
pico- and nanosecond range, are of special interest as they contribute to
allosteric pathways [4], catalytic reaction [5], the malleability of binding
sites to accomodate different ligands or, in the case of side-chains, act as
an entropy reservoir [6, 7].
A large variety of NMR experiments has been developped to measure the rates at
which ensembles of nuclear spin return to their equilibrium after a
perturbation [8], a process called nuclear spin relaxation. Relaxation is
driven by intra- and inter-atomic interactions which strengths are affected by
- or correlated to - the fluctuations of bond vectors orientations [9]. Thus,
relating relaxation rates to the time-dependent orientations of the
interaction frames allows to characterize the dynamics of bond vectors. This
can be quantified in the frame of the Bloch-Wangsness-Redfield (BWR)
relaxation theory which, in addition to quantum mechanic terms related to spin
interactions, uses spherical harmonics to describe the orientation of the
interaction frames in a fixed laboratory frame [10, 11, 12].
Relaxation properties depend on the correlation function for the orientation
of the principal axis frames of interactions contributing to relaxation. More
specifically, relaxation rates are expressed as linear combinations of its
Fourier transform, called the spectral density function, and evaluated at the
eigenfrequencies of the spin system. In order to obtain a quantitative
description of the motions from the analysis of NMR relaxation rates, an
analytical form of the correlation function has to be written with parameters
describing the amplitudes and timescales of the motions. This requires a
priori knowledge on the type of motions to build a model of correlation
function, and can lead to complex analytical expressions. In addition, when
analyzing a small set of experimental relaxation rates, several physically
different models may describe the experiments equally well.
These limitations have been overcome by the Model-Free (MF) approach which
aims at reproducing the correlation function for internal motions using a
decaying exponential characterized by an effective correlation-time and a
limit at infinite time called generalized order parameter [13]. This model has
later been complexified to the Extended Model-Free (EMF) approach [14] to
include two time-scales and accompagnying order parameters. Both have lead to
a large number of succesfull studies of protein dynamics by NMR [15, 16, 17,
18, 19, 20]. However, MF or EMF approaches are uninformative on the nature of
motions of protein backbone and side-chains, and relaxation analysis has to be
complemented with Molecular Dynamics (MD) simulations [21, 22, 23, 24, 25,
26].
In MF type of analysis, the global tumbling is often supposed to be isotropic,
which is the frame in which it was originially presented [13]. An approximated
correlation function for overall diffusion can be used when the diffusion
tensor is anisotropic [13, 27]. This limitation has also raised a number of
questions regarding the relevance of the analysis of relaxation data with MF-
type of correlation functions in the case of couplings between global (i.e
overall diffusion) and local motions [28, 29]. and has lead to the development
of the Slowly Relaxing Local Structure(SRLS) approach for the analysis of NMR
relaxation data [30, 31]. In addition, and as suggested by Freed and
coworkers, the assumptions of simple local geometry and motions with axial
symmetry have to be invoked in the frame of the MF and EMF approaches [32, 33,
34]. This aspect is particularly critical as MF order parameters have been
shown to be related to entropy, both for protein backbone [35] and side-chains
[36], and a mischaracterization of the motions will inevitably lead to a
mischaracterization of the thermodynamics [37].
Although perfectible, modern MD simulations inform us quantitatively on the
nature of molecular motions. It is now possible to obtain the a priori
knowledge required to develop explicit models of motions, building on
analytical models introduced between 1960 and 1980 [38, 39, 40]. Thus, the
complexity of the expressions of correlation functions is now outweight by the
wealth of information that can potentially be obtained by using such models.
Here, we review three models (the rotation on a cone, the rotamer jump and the
wobbling in a cone) and detail how to write complex correlation functions
adapted to protein backbone and side-chain dynamics. We discuss cases of
correlated motions and cases where the strengths of the nuclear spin
interactions include time-dependent amplitudes. We aim at presenting a
coherent framework to write correlation functions suitable for a wide range of
combinations of motions. We use synthetic data to compare explicit models of
motions and the MF correlation function, in particular in the presence of
asymmetric internal and overall motions. This work lays the theoretical
fundation for the analysis fo experimental NMR relaxation rates and molecular
dynamics simulations, which we present elsewhere [41].

##  2 Correlation function in the BWR theory

The description of the Bloch-Wangsness-Redfield (BWR) relaxation theory is
beyond the scope of this article and can be found elsewhere [42, 43, 9]. We
will only write the steps essential to reveal the correlation function in the
relaxation superoperator. The evolution of the density operator $\hat{\sigma}$
is given by the Liouville-von Neumann equation:

$\frac{\mathrm{d}\hat{\sigma}}{\mathrm{d}t}=-i\left[\hat{\mathcal{H}}(t),\hat{\sigma}(t)\right],$
(1)

where the Hamiltonian
$\hat{\mathcal{H}}(t)=\hat{\mathcal{H}}_{0}+\hat{\mathcal{H}}_{1}(t)$ is
expressed as the sum of stationnary $\hat{\mathcal{H}}_{0}$ and fluctuating
$\hat{\mathcal{H}}_{1}(t)$ Hamiltonians. This equation is transformed in the
interaction frame of the Hamiltonian $\hat{\mathcal{H}}_{0}$ (denoted by a
tilde) to yield, after considering the ensemble average (denoted by the
overbar) and hypotheses introduced by Bloch, Wangsness and Redfield:

$\frac{\overline{\mathrm{d}\tilde{\hat{\sigma}}(t)}}{\mathrm{d}t}=-\int\limits_{0}^{\infty}\overline{\left[\tilde{\hat{\mathcal{H}}}_{1}(t),[\tilde{\hat{\mathcal{H}}}_{1}(t+\tau),\tilde{\hat{\sigma}}(t)-\tilde{\hat{\sigma}}_{eq}]\right]}\mathrm{d}\tau,$
(2)

where $\tilde{\hat{\sigma}}_{eq}$ is the equilibrium density operator. The BWR
relaxation theory is a semi-classical formalism based on the separation of the
quantum spin part and the geometrical part when writting the Hamiltonian for
the interactions contributing to a relaxation mechanism:

$\hat{\mathcal{H}}_{1}(t)=\sum_{i}\zeta_{i}\sum^{2}\limits_{q=-2}(-1)^{q}V^{i}_{2,-q}(t)\hat{T}^{i}_{2,q},$
(3)

where $\zeta_{i}$ is the amplitude of the interaction $i$, $V^{i}_{2,-q}(t)$
is related to a rank-2 Wigner matrix and $\hat{T}^{i}_{2,q}$ is a rank-2
tensor spin operator with coherence order $q$ and which can further be
decomposed as a sum of eigenoperators ${\hat{A}_{2,q,p}^{i}}$ of the super-
operator $[\hat{\mathcal{H}}_{0},\cdot]$ with eigenvalues
$\omega_{2,q,p}^{(i)}$. After subsituting Eq. 3 in Eq. 2, we obtain:

$\displaystyle\frac{\overline{\mathrm{d}\tilde{\hat{\sigma}}(t)}}{\mathrm{d}t}=$
$\displaystyle-\sum\limits_{i,j}\zeta_{i}\zeta_{j}\sum^{2}\limits_{q,q^{\prime}=-2}\sum_{p,p^{\prime}}(-1)^{q+q^{\prime}}e^{i(\omega_{2,q,p}^{(i)}-\omega_{2,q^{\prime},p^{\prime}}^{(j)})t}\left[\hat{A}_{2,q,p}^{i},[\hat{A}_{2,q^{\prime},p^{\prime}}^{j,{\dagger}},\tilde{\hat{\sigma}}(t)]\right]\times$
(4) $\displaystyle\int\limits_{0}^{\infty}\langle
V^{i}_{2,-q}(t)V_{2,-q^{\prime}}^{j,\ast}(t+\tau)\rangle
e^{-i\omega_{2,q^{\prime},p^{\prime}}^{(j)}\tau}\mathrm{d}\tau,$

where $\langle\cdot\cdot\cdot\rangle$ denotes an ensemble average and
${\dagger}$ and $*$ stand for the conjugate transpose and complex conjugate,
respectively. The correlation function $C_{i,j}$ between interactions $i$ and
$j$ is now defined as:

$\langle
V^{i}_{2,-q}(0)V_{2,-q^{\prime}}^{j,\ast}(t)\rangle=\delta_{q,q^{\prime}}C_{i,j}(t).$
(5)

It can be written using the Wigner matrix notations:

$C_{i,j}(t)=\langle\mathcal{D}_{q,0}^{(2)*}(\Omega_{L,i},0)\mathcal{D}_{q,0}^{(2)}(\Omega_{L,j},t)\rangle,$
(6)

where $\Omega_{L,i}=\\{\alpha_{L,i},\beta_{L,i},\gamma_{L,i}\\}$ is the set of
three Euler angles orienting the interaction frame in the laboratory frame.

#### Useful properties of Wigner matrices.

Before going any further, we will review some mathematical properties of
Wigner rotation matrices that will be used when calculating correlation
functions.
Rank-2 spherical harmonics $\mathrm{Y}_{2m}$ are directly linked to Wigner
matrices $\mathcal{D}_{m0}^{(2)}$ as:

$\mathcal{D}_{m0}^{(2)}(\alpha,\beta,\gamma)=\sqrt{\frac{4\pi}{5}}\mathrm{Y}_{2m}^{*}(\beta,\alpha),$
(7)

with:

$\mathcal{D}_{mn}^{(2)}(\alpha,\beta,\gamma)=e^{-im\alpha}d_{mn}^{(2)}(\beta)e^{-in\gamma},$
(8)

where $d_{mn}$ is called the reduced Wigner matrix. From this definition (and
the one of the reduced Wigner matrices), it follows:

$\mathcal{D}_{mn}^{(2)}(\Omega)=(-1)^{-m-n}\mathcal{D}_{-m-n}^{(2)*}(\Omega).$
(9)

The addition theorem relates rank-2 spherical harmonics (or Wigner matrices)
to the second-order Legendre polynomial
$\mathcal{P}_{2}(x)=\frac{3x^{2}-1}{2}$:

$\mathcal{P}_{2}(\vec{x}\cdot\vec{y})=\sum_{m=-2}^{2}\mathcal{D}_{m0}^{(2)*}(\vec{y})\mathcal{D}_{m0}^{(2)}(\vec{x}).$
(10)

Wigner matrices can be decomposed into succesive rotations:

$\mathcal{D}_{mn}^{(2)}(\Omega_{A,B})=\sum_{x_{1}=-2}^{2}\sum_{x_{2}=-2}^{2}...\sum_{x_{k}=-2}^{2}\mathcal{D}_{mx_{1}}^{(2)}(\Omega_{A,X_{1}})\mathcal{D}_{x_{1}x_{2}}^{(2)}(\Omega_{X_{1},X_{2}})...\mathcal{D}_{x_{k}n}^{(2)}(\Omega_{X_{k},B}),$
(11)

where $\Omega_{X_{i},X_{j}}$ is the Euler angle for transformation between
frame $X_{i}$ and frame $X_{j}$.
Wigner matrices can be normalized:

$\int
d\Omega\mathcal{D}_{mn}^{L*}(\Omega)\mathcal{D}_{m^{\prime}n^{\prime}}^{L^{\prime}}(\Omega)=\frac{8\pi^{2}}{2L+1}\delta_{LL^{\prime}}\delta_{mm^{\prime}}\delta_{nn^{\prime}},$
(12)

where $\delta$ is the Kronecker function.
Finally, we have the following property:

$\sum_{m=-2}^{2}\mathcal{D}_{mn}^{(2)*}(\alpha,\beta,0)\mathcal{D}_{mn^{\prime}}^{(2)}(\alpha,\beta,0)=\delta_{nn^{\prime}}.$
(13)

It is worthwhile to specify the convention for Euler angle used throughout
this paper. If we write an Euler angle set $\Omega=\\{\alpha,\beta,\gamma\\}$,
for tranformation from a frame $\\{O,x,y,z\\}$ to
$\\{O^{\prime},x^{\prime},y^{\prime},z^{\prime}\\}$ (the two frame origins can
be different, which only yields a translation without affecting the
rotations), $\alpha$ is the rotation angle around the $Oz$ axis, $\beta$ is
the rotation angle about the new $Ox$ axis (axis rotated by an angle $\alpha$
around the $Oz$ axis), and $\gamma$ is the rotation about the new $Oz$ axis
(after applying the two previous rotations, only the last one affecting the
z-axis).

##  3 Explicit models of correlation function

###  3.1 Approach to calculate correlation functions

The approach to calculate the ensemble average in Eq. 6 is detailed in Ref.
[44]. If we write:

$C(t)=\langle A^{*}(\Omega(0))B(\Omega(t))\rangle,$ (14)

with $A$ and $B$ two Wigner matrices with potentially different indices, the
ensemble average expands into:

$C(t)=\int d\Omega_{0}\int d\Omega
P(\Omega_{0})P(\Omega,t|\Omega_{0},0)A^{*}(\Omega_{0})B(\Omega),$ (15)

where $P(\Omega)$ is the probability that the Euler angles for frame
transformation is $\Omega$, and $P(\Omega,t|\Omega_{0},0)$ is the conditional
probability that the Euler angle for frame transformation is $\Omega$ at time
$t$ when it was $\Omega_{0}$ at time $t=0$. In the following, we will
calculate the ensemble averages using a Master equation for the description of
relaxation-inducing processes:

$\frac{\partial}{\partial t}P(\Omega,t)=-KP(\Omega,t),$ (16)

which is a Fokker-Planck diffusion equation and where $P(\Omega,t)$ is the
probability of finding the Euler angles $\Omega$ for frame transformation, at
time $t$, and $K$ is a model-dependent operator describing the motions of
interest. The initial step is to diagonalize the operator $K$, with
eigenvalues $E_{n}$ and eigenvectors $\psi_{n}$:

$K\psi_{n}=E_{n}\psi_{n}.$ (17)

The conditional probability is then:

$P(\Omega,t|\Omega_{0},0)=\sum_{n}\psi_{n}^{*}(\Omega_{0})\psi_{n}(\Omega)e^{-E_{n}t},$
(18)

and:

$P(\Omega_{0})=\lim_{t\rightarrow\infty}P(\Omega,t|\Omega_{0},0).$ (19)

$P(\Omega_{0})$ can sometimes be evaluated by simple probabilistic and
geometric considerations. These probabilities are then inserted in Eq. 15 to
calculate the correlation function.

###  3.2 Global tumbling

In this section, we will only consider the rotational diffusion of the
protein: all interactions are fixed in the molecular frame. The correlation
function describes the motions of the principal axis frames of interactions in
the laboratory frame, which is equivalent to the motions of the protein in the
laboratory frame. We follow the thorough treatment of this situation presented
by Werbelow and Grant [45], but other approaches can be found elsewhere [46,
47, 48, 49],
The correlation function can be decomposed using Eq. 11:

$\displaystyle C_{i,j}(t)=$
$\displaystyle\langle\mathcal{D}_{q0}^{(2)*}(\Omega_{L,i},0)\mathcal{D}_{q0}^{(2)}(\Omega_{L,j},t)\rangle$
(20) $\displaystyle=$
$\displaystyle\sum_{a,a^{\prime}=-2}^{2}\langle\mathcal{D}_{qa}^{(2)*}(\Omega_{L,D},0)\mathcal{D}_{qa^{\prime}}^{(2)}(\Omega_{L,D},t)\rangle\mathcal{D}_{a0}^{(2)*}(\Omega_{D,i})\mathcal{D}_{a^{\prime}0}^{(2)}(\Omega_{D,j}),$

where $D$ is the frame associated to the diffusion tensor, fixed in the
molecular frame, $\Omega_{L,D}$ is the time-dependent Euler angle set for the
orientation of the diffusion tensor in the laboratory frame and
$\Omega_{D,k}$, $k=\\{i,j\\}$, are the time-independent Euler angle sets for
orientation of the interactions $i$ and $j$ in the diffusion frame, fixed in
the molecular frame. In the following, we define:

$C_{aa^{\prime}}(t)=\langle\mathcal{D}_{qa}^{(2)*}(\Omega_{L,D},0)\mathcal{D}_{qa^{\prime}}^{(2)}(\Omega_{L,D},t)\rangle.$
(21)

#### Master equation for overall diffusion.

The Master equation for global tumbling, written in a frame where the
rotational diffusion tensor is diagonal, is:

$\frac{\partial}{\partial
t}P(\Omega_{L,D},t)=-\sum_{j=x,y,z}D_{jj}L_{j}^{2}P(\Omega_{L,D},t),$ (22)

where $D_{jj}$ is the $j$th component of the diagonalized diffusion tensor and
$L_{j}$ is the associated angular momentum operators. We express the diffusion
operator in terms of raising and lowering operators:

$\mathcal{D}_{g}=\sum_{j=x,y,z}D_{jj}L_{j}^{2}=D_{+}L^{2}+\frac{1}{2}D_{-}(L_{+}^{2}+L_{-}^{2})+(D_{zz}-D_{+})L_{z}^{2},$
(23)

where:

$\displaystyle D_{\pm}=$ $\displaystyle\frac{1}{2}(D_{xx}\pm D_{yy}),$ (24)
$\displaystyle L_{\pm}=$ $\displaystyle L_{x}\pm iL_{y}.$

The eigen-equation is written:

$\mathcal{D}_{g}\Psi_{n}=E_{n}\Psi_{n},$ (25)

where the functions $\Psi_{n}$ are expanded in terms of the normalized Wigner
matrices:

$\Psi_{L,\mu,\kappa}(\Omega)=\sum_{k}a_{\kappa,k}\sqrt{\frac{2L+1}{8\pi^{2}}}\mathcal{D}_{\mu,k}^{L}(\Omega).$
(26)

#### Solving the eigen-equation.

The correlation function (Eq. 21) is expressed as:

$\displaystyle C_{aa^{\prime}}=$ $\displaystyle\int d\Omega_{0}\int d\Omega
P(\Omega_{0})\mathcal{D}_{qa}^{(2)*}(\Omega_{0})\mathcal{D}_{qa^{\prime}}^{(2)}(\Omega)\times$
(27)
$\displaystyle\sum_{L,\mu,k,k^{\prime},\kappa}\frac{2L+1}{8\pi^{2}}a_{\kappa,k}a_{\kappa,k^{\prime}}\mathcal{D}_{\mu,k}^{L*}(\Omega_{0})\mathcal{D}_{\mu,k^{\prime}}^{L}(\Omega)e^{-E_{L,\mu,\kappa}t}.$

The orthogonality condition of the Wigner matrices (Eq. 12) imposes $L=2$ when
calculating the integral. As a consequence, only 5 non-degenerate
eigenfunctions are needed to solve Eq. 22. With this simplification in hand,
the eigenvalues and eigenfunctions of the components of the angular momentum
operators are:

$\displaystyle L^{2}\mathcal{D}_{\mu,k}^{(2)}(\Omega)$
$\displaystyle=6\mathcal{D}_{\mu,k}^{(2)}(\Omega),$ (28) $\displaystyle
L_{z}^{2}\mathcal{D}_{\mu,k}^{(2)}(\Omega)$
$\displaystyle=k^{2}\mathcal{D}_{\mu,k}^{(2)}(\Omega).$

Table 1: Effect of the raising ($L_{+}$) and lowering ($L_{-}$) squared operators on the rank-2 Wigner matrices. k | $L_{+}^{2}\mathcal{D}_{\mu,k}^{(2)}(\Omega)$ | $L_{-}^{2}\mathcal{D}_{\mu,k}^{(2)}(\Omega)$
---|---|---
-2 | $\sqrt{24}\mathcal{D}_{\mu,0}^{(2)}(\Omega)$ | 0
-1 | $6\mathcal{D}_{\mu,1}^{(2)}(\Omega)$ | 0
0 | $\sqrt{24}\mathcal{D}_{\mu,2}^{(2)}(\Omega)$ | $\sqrt{24}\mathcal{D}_{\mu,-2}^{(2)}(\Omega)$
1 | 0 | $6\mathcal{D}_{\mu,-1}^{(2)}(\Omega)$
2 | 0 | $\sqrt{24}\mathcal{D}_{\mu,0}^{(2)}(\Omega)$

The effects of the raising and lowering squared operators are summed up in
Table 1. It is then straightforward to find 3 eigenvectors and associated
eigenvalue of $\mathcal{D}_{g}$ ($\kappa=1,2,3$ in Table 2), and the remaining
two must have the form:
$X(x_{1},x_{2})=x_{1}\sqrt{\frac{5}{8\pi^{2}}}\left(\mathcal{D}_{\mu,2}^{(2)}(\Omega)+\mathcal{D}_{\mu,-2}^{(2)}(\Omega)\right)+x_{2}\sqrt{\frac{5}{8\pi^{2}}}\mathcal{D}_{\mu,0}^{(2)}(\Omega)$,
where $x_{1},x_{2}\in\Re$. The effect of the operator $\mathcal{D}_{g}$ on
this eigenfunction is:

$\displaystyle\mathcal{D}_{g}X(x_{1},x_{2})=$ $\displaystyle
x_{1}\sqrt{\frac{5}{8\pi^{2}}}\left(\mathcal{D}_{\mu,2}^{(2)}(\Omega)+\mathcal{D}_{\mu,-2}^{(2)}(\Omega)\right)\left[6D_{+}+4(D_{zz}-D_{+})+\frac{\sqrt{24}}{2}D_{-}\frac{x_{2}}{x_{1}}\right]$
(29) $\displaystyle+$ $\displaystyle
x_{2}\sqrt{\frac{5}{8\pi^{2}}}\mathcal{D}_{\mu,0}^{(2)}(\Omega)\left[6D_{+}+\sqrt{24}D_{-}\frac{x_{1}}{x_{2}}\right],$

which yields to the following condition for $X(x_{1},x_{2})$ to be an
eigenvector:

$6D_{+}+4(D_{zz}-D_{+})+\frac{\sqrt{24}}{2}D_{-}\frac{x_{2}}{x_{1}}=6D_{+}+\sqrt{24}D_{-}\frac{x_{1}}{x_{2}}.$
(30)

In addition, the normalization condition for $X(x_{1},x_{2})$ gives:

$2x_{1}^{2}+x_{2}^{2}=1.$ (31)

We set $x_{1}^{\prime 2}=2x_{1}^{2}$ such that:

$x_{1}^{\prime 2}+x_{2}^{2}=1.$ (32)

We can find a number $x$ leading to:

$x_{1}^{\prime}=\sin(\frac{x}{2}),\,\,\,x_{2}=\cos(\frac{x}{2}).$ (33)

After inserting in Eq. 30 and defining
$\alpha=\frac{\sqrt{3}D_{-}}{D_{zz}-D_{+}}$, we obtain:

$x=-\tan^{-1}(\alpha),$ (34)

and one set of solution to Eq. 30:

$x_{1}=-\frac{1}{\sqrt{2}}\sin\left(\frac{\tan^{-1}(\alpha)}{2}\right),\,\,\,x_{2}=\cos\left(\frac{\tan^{-1}(\alpha)}{2}\right).$
(35)

A second set of solution is obtained by choosing:

$x_{1}^{\prime}=-\cos(\frac{x}{2}),\,\,\,x_{2}^{\prime}=\sin(\frac{x}{2}),$
(36)

which is a $\frac{\pi}{2}$-shift of the previous solution, and leads to:

$x_{1}^{\prime}=\frac{1}{\sqrt{2}}\cos\left(\frac{\tan^{-1}(\alpha)}{2}\right),\,\,\,x_{2}^{\prime}=\sin\left(\frac{\tan^{-1}(\alpha)}{2}\right).$
(37)

The eigenvalues associated to these two eigenvectors can be calculated easily
using Eq. 29. The five eigenvectors and associated eigenvalues of Eq. 25 are
gathered in Table 2 where we have defined
$\beta=\tan^{-1}\frac{\sqrt{3}D_{-}}{D_{zz}-D_{+}}$.

Table 2: Eigenvalues and eigenvectors for the operator associated to global tumbling. We define $\beta=\tan^{-1}\frac{\sqrt{3}D_{-}}{D_{zz}-D_{+}}$. $\kappa$ | $E_{\kappa}$ | $\Psi_{2,\mu,\kappa}(\Omega)/\sqrt{\frac{5}{8\pi^{2}}}=\sum_{k}a_{\kappa,k}\mathcal{D}_{\mu,k}^{(2)}(\Omega)$
---|---|---
1 | $5D_{+}+3D_{-}+D_{zz}$ | $\frac{1}{\sqrt{2}}\left(\mathcal{D}_{\mu,1}^{(2)}(\Omega)+\mathcal{D}_{\mu,-1}^{(2)}(\Omega)\right)$
2 | $5D_{+}-3D_{-}+D_{zz}$ | $\frac{1}{\sqrt{2}}\left(\mathcal{D}_{\mu,1}^{(2)}(\Omega)-\mathcal{D}_{\mu,-1}^{(2)}(\Omega)\right)$
3 | $2D_{+}+4D_{zz}$ | $\frac{1}{\sqrt{2}}\left(\mathcal{D}_{\mu,2}^{(2)}(\Omega)-\mathcal{D}_{\mu,-2}^{(2)}(\Omega)\right)$
4 | $4D_{+}+2D_{zz}+2\left[(D_{zz}-D_{+})^{2}+3D_{-}^{2}\right]^{\frac{1}{2}}$ | $\frac{1}{\sqrt{2}}\left(\mathcal{D}_{\mu,2}^{(2)}(\Omega)+\mathcal{D}_{\mu,-2}^{(2)}(\Omega)\right)\cos\frac{\beta}{2}+\mathcal{D}_{\mu,0}^{(2)}(\Omega)\sin\frac{\beta}{2}$
5 | $4D_{+}+2D_{zz}-2\left[(D_{zz}-D_{+})^{2}+3D_{-}^{2}\right]^{\frac{1}{2}}$ | $-\frac{1}{\sqrt{2}}\left(\mathcal{D}_{\mu,2}^{(2)}(\Omega)+\mathcal{D}_{\mu,-2}^{(2)}(\Omega)\right)\sin\frac{\beta}{2}+\mathcal{D}_{\mu,0}^{(2)}(\Omega)\cos\frac{\beta}{2}$

When the diffusion tensor shows some symmetry properties, the angular momentum
operator can have degenerate eigenfunctions. For example, if the molecule
behaves as a cylinder (axial symmetry), then $D_{xx}=D_{yy}$ (a frame
transformation can always lead to this situation), eliminating the $L_{+}$ and
$L_{-}$ parts of Eq. 23. The immediate consequence is that $D_{-}=0$ and two
of the three unique eigenvalues are doubly degenerate ($\kappa=1$ and 2, and
$\kappa=3$ and 4 for $D_{z}>D_{x}$ or $\kappa=3$ and 5 for $D_{z}<D_{x}$).
When the diffusion tensor is isotropic, $D_{xx}=D_{yy}=D_{zz}=D$, all
eigenvalues are degenerate and equal $6D$, usually refered to as the inverse
of the global tumbling correlation time: $\tau_{c}=1/(6D)$.

#### Static angle probability distribution.

It is clear from the form of the eigenvalues of the components of the angular
momentum that one only needs to calculate the eigenfunction associated with
the 0 eigenvalue in order to have a non-vanishing exponential term in Eq. 18
when calculating the limit with $t\rightarrow\infty$. It comes, from the
eigenvalues of $L^{2}$ (Eq. 28), that such an eigenfunction has $L=0$ and,
from the eigenvalues of $L_{z}$, $k=0$ such that it corresponds to
$D_{0,0}^{(0)}(\Omega)$. Then:

$P(\Omega_{0})=\frac{1}{8\pi^{2}}D_{0,0}^{(0)}(\Omega_{0})^{2}=\frac{1}{8\pi^{2}}.$
(38)

#### Integration of the correlation function.

Using the integration in Eq. 27:

$\displaystyle C_{aa^{\prime}}(t)=$ $\displaystyle\int d\Omega_{0}\int
d\Omega\sum_{\mu,k,k^{\prime}=-2}^{2}\sum_{\kappa=1}^{5}P(\Omega_{0})\frac{5}{8\pi^{2}}\times$
(39) $\displaystyle
e^{-E_{\kappa}t}a_{\kappa,k}a_{\kappa,k^{\prime}}\mathcal{D}_{qa}^{(2)*}(\Omega_{0})\mathcal{D}_{\mu
k}^{(2)*}(\Omega_{0})\mathcal{D}_{qa^{\prime}}^{(2)}(\Omega)\mathcal{D}_{\mu
k^{\prime}}^{(2)}(\Omega).$

We can use Eq. 9 to transform this equation into:

$\displaystyle C_{aa^{\prime}}(t)=$
$\displaystyle\frac{5}{64\pi^{4}}\sum_{\mu,k,k^{\prime}=-2}^{2}\sum_{\kappa=1}^{5}a_{\kappa,k}a_{\kappa,k^{\prime}}e^{-E_{\kappa}t}\times$
(40) $\displaystyle\times(-1)^{-k}\int
d\Omega_{0}\mathcal{D}_{qa}^{(2)*}(\Omega_{0})\mathcal{D}_{-\mu-k}^{(2)}(\Omega_{0})\times$
$\displaystyle\times(-1)^{-k^{\prime}}\int
d\Omega\mathcal{D}_{qa^{\prime}}^{(2)}(\Omega)\mathcal{D}_{-\mu-k^{\prime}}^{(2)*}(\Omega),$

According to Table 2, $k+k^{\prime}$ is always an even number such that
$(-1)^{-k-k^{\prime}}=1$ for all $k$ and $k^{\prime}$. It immediately follows
from the orthonormality condition (Eq. 12) that:

$C_{aa^{\prime}}(t)=\frac{5}{64\pi^{4}}\sum_{\mu,k,k^{\prime}=-2}^{2}\sum_{\kappa=1}^{5}\delta_{q-\mu}\delta_{a-k}\delta_{a^{\prime}-k^{\prime}}a_{\kappa,k}a_{\kappa,k^{\prime}}e^{-E_{\kappa}t}\frac{8\pi^{2}}{5}\times\frac{8\pi^{2}}{5},$
(41)

where $\delta$ is the Kronecker-delta function. It simplifies into:

$C_{aa^{\prime}}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}a_{\kappa,-a}a_{\kappa,-a^{\prime}}e^{-E_{\kappa}t}\,.$
(42)

From symmetry consideration, it follows that
$a_{\kappa,-a}a_{\kappa,-a^{\prime}}=a_{\kappa,a}a_{\kappa,a^{\prime}}$ so
that:

$C_{aa^{\prime}}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t},$

and the correlation function can now be written as:

$C_{i,j}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,i})\mathcal{D}_{a^{\prime},0}^{(2)}(\Omega_{D,j}).$
(43)

#### Special case 1: axial symmetry.

In this case, the Wigner matrices
$\mathcal{D}_{\mu,k}^{(2)}(\Omega),k\in[-2,2]$ are eigenfunctions of the
diffusion operator with associated eigenvalue
$6D_{\perp}+k^{2}\left(D_{\parallel}-D_{\perp}\right)$ where we define
$D_{\perp}=D_{xx}=D_{yy}$ and $D_{\parallel}=D_{zz}$. It follows:

$C_{aa^{\prime}}(t)=\delta_{aa^{\prime}}\frac{1}{5}e^{-(6D_{\perp}+a^{2}(D_{\parallel}-D_{\perp}))t},$
(44)

and the total correlation function is:

$C_{i,j}(t)=\frac{1}{5}\sum_{a=-2}^{2}e^{-(6D_{\perp}+a^{2}(D_{\parallel}-D_{\perp}))t}\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,i})\mathcal{D}_{a,0}^{(2)}(\Omega_{D,j}).$
(45)

#### Special case 2: isotropic tumbling.

In this case as well, the Wigner matrices
$\mathcal{D}_{\mu,k}^{(2)}(\Omega),k\in[-2,2]$ are eigenfunctions of the
diffusion operator with degenerate eigenvalue of $6D$. The well-known result
immediately follows:

$C_{aa^{\prime}}(t)=\delta_{aa^{\prime}}\frac{1}{5}e^{-t/\tau_{c}},$ (46)

where the global tumbling correlation time is defined as:

$\tau_{c}=\frac{1}{6D}.$ (47)

The total correlation function is:

$C_{i,j}(t)=\frac{1}{5}\sum_{a=-2}^{2}e^{-t/\tau_{c}}\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,i})\mathcal{D}_{a,0}^{(2)}(\Omega_{D,j}),$
(48)

and after using Eq. 10, we have:

$C_{i,j}(t)=\frac{1}{5}e^{-t/\tau_{c}}\mathcal{P}_{2}(\cos\theta_{i,j}),$ (49)

where $\theta_{i,j}$ is the angle between interactions $i$ and $j$. In the
case of small anisotropy of the overall diffusion tensor, this simpler model
of correlation function can be used with residue-specific correlation time for
global tumbling defined from the projection of the diffusion tensor onto the
interaction frame [50].

###  3.3 Rotamer jumps

In this section, we add one internal motion on top of the global tumbling. We
assume that the moiety in which the principal axis frame of the interactions
is anchored (e.g. a peptide plane, or methyl group) can adopt a finite number
of fixed conformations and that this moiety can exchange between all
conformations. Moreover, it is assumed that the transition events are
infinitely fast, that is, the system can always be found in one particular
conformation. This model has been used to describe methyl rotation as a three-
site jump[51, 52] and is well adapted to study side-chain motions which can
undergo fast rotamer transitions, as revealed by MD simulations [53, 24, 26,
37].
The treatment of rotameric jumps in this section is based on the work of
Wittebort and Szabo [39]. We will use the result from the previous section and
write the correlation function as:

$C_{i,j}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\langle\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,i},0)\mathcal{D}_{a^{\prime},0}^{(2)}(\Omega_{D,j},t)\rangle,$
(50)

where the ensemble average $\langle...\rangle$ accounts for the presence of
internal motions. This equation is valid only under the assumption that global
tumbling and internal motions are uncorrelated. This allows to introduce the
jump frame for a facilitated description of the jump. For example, in the case
of valines $\chi_{1}$ rotameric jumps (corresponding to different orientations
of the Cβ-Cγ1 axis), the jump frame main axis ($z_{J}$) points along the Cα-Cβ
bond, with origin at the Cβ. The orientation of $x_{J}$ axis is arbitrary. The
Wigner matrices can be split into successive frame transformations:

  * •

from the diffusion frame to the jump frame, with Euler angle $\Omega_{D,J}$.
This frame is fixed in the diffusion frame, and the transformation is time-
independent.

  * •

from the jump frame to the system frame (SF), with Euler angle
$\Omega_{J,SF}=\\{\alpha_{J,SF},\beta_{J,SF},0\\}$. The main axis of the
system frame points along the chemical bond defining the conformation. This
transformation is time-dependent.

  * •

from the system frame to the principal axis frame (PAF) of the interaction.
For instance, the PAF of a dipole-dipole interaction has its main axis
pointing along the internuclear vector. This is a time-independent
transformation, with Euler angle $\Omega_{SF,i}$.

These additional frame transformations lead to:

$\displaystyle C_{i,j}(t)=$
$\displaystyle\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime},c,c^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle\times$
(51)
$\displaystyle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,J})\mathcal{D}_{a^{\prime},b^{\prime}}^{(2)}(\Omega_{D,J})\mathcal{D}_{c,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{c^{\prime},0}^{(2)}(\Omega_{SF,j}).$

The conditional probability $P(\Omega,t|\Omega_{0},0)$ (Eq. 18) is found by
solving the Master equation:

$\frac{\partial}{\partial
t}p_{\alpha}(t)=\sum_{\beta=1}^{N}\mathcal{R}_{\alpha\beta}p_{\beta}(t),$ (52)

where $p_{\alpha}(t)$ is the population of state $\alpha$ at time $t$ and
$\mathcal{R}_{\alpha\beta}$ is an element of the exchange matrix $\mathcal{R}$
and corresponds to the exchange rate from state $\beta$ to $\alpha$. The
microscopic reversibility implies that:

$\mathcal{R}_{\alpha\beta}p_{\beta}^{eq}=\mathcal{R}_{\beta\alpha}p_{\alpha}^{eq},$
(53)

and diagonal elements of the exchange matrix are given by:

$\mathcal{R}_{\alpha\alpha}=-\sum_{\beta\neq\alpha}R_{\beta\alpha}.$ (54)

The boundary conditions associated to Eq. 52 are:

$\displaystyle(1)\,\,\,p(\beta,t=0|\alpha,0)=\delta_{\alpha\beta},$ (55)
$\displaystyle(2)\,\,\,\lim_{t\rightarrow\infty}p(\beta,t|\alpha,0)=p_{\beta}^{eq},$

Eq. 52 transforms into an eigen-equation as follows:

$\mathcal{R}X=\lambda X.$ (56)

The exchange matrix is in general not symmetric (it is only when equilibrium
populations are equal) but can be transformed into a symmetric matrix
$\tilde{\mathcal{R}}$ with coefficients:

$\displaystyle\widetilde{\mathcal{R}}_{\alpha\beta}$
$\displaystyle=\sqrt{\mathcal{R}_{\alpha\beta}\mathcal{R}_{\beta\alpha}},\,\,\,\alpha\neq\beta,$
(57) $\displaystyle\widetilde{\mathcal{R}}_{\alpha\alpha}$
$\displaystyle=\mathcal{R}_{\alpha\alpha}.$

The eigen-equation now reads:

$\tilde{\mathcal{R}}\tilde{X}=\lambda\tilde{X},$ (58)

where the $\alpha^{\mathrm{th}}$ coordinate of the eigenvector associated to
the eigenvalue $\lambda_{n}$ is:

$\tilde{X}_{\alpha}^{(n)}=\frac{1}{\sqrt{p_{\alpha}^{eq}}}X_{\alpha}^{(n)}.$
(59)

Since $\tilde{\mathcal{R}}$ is symmetric, its eigenvalues $\lambda_{n}$ are
all real and the associated eigenvectors $\tilde{X}^{(n)}$ are orthogonal. The
conditional probability is written as in Eq. 18 such that, for the second
condition in Eq. 55 to be met, we must have $\lambda_{0}=0$ (we order the
eigenvalues such that
$\left|\lambda_{n}\right|\leq\left|\lambda_{n+1}\right|$) and the associated
eigenvector is $\\{\tilde{X}_{n}^{(0)}\\}=\\{\sqrt{p_{n}^{eq}}\\}$.
The conditional probability can be written:

$p(\beta,t|\alpha,0)=\sqrt{\frac{p_{\beta}^{eq}}{p_{\alpha}^{eq}}}\sum_{n=0}^{N-1}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}e^{\lambda_{n}t},$
(60)

where the factor $\sqrt{\frac{p_{\beta}^{eq}}{p_{\alpha}^{eq}}}$ is introduced
to fulfill the second condition in Eq. 55 (the eigenvectors $\tilde{X}^{(n)}$
are orthogonal such that condition (1) is already met). It leads to the
following expression for the correlation function (note that integrals in Eq.
15 are replaced by discrete sums):

$\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle=\sum_{\alpha,\beta=1}^{N}\sum_{n=0}^{N-1}e^{\lambda_{n}t}\sqrt{p_{\alpha}^{eq}p_{\beta}^{eq}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF_{\alpha}})\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF_{\beta}}).$
(61)

The total correlation function accounting for overall rotational diffusion and
rotamer jumps is then written:

$C_{i,j}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},=-2}^{2}\sum_{\alpha,\beta=1}^{N}\sum_{n=0}^{N-1}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\sqrt{p_{\alpha}^{eq}p_{\beta}^{eq}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}e^{\lambda_{n}t}\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,i_{\alpha}})\mathcal{D}_{a^{\prime},0}^{(2)}(\Omega_{D,j_{\beta}}),$
(62)

where $\Omega_{D,i_{\alpha}}$ is the Euler angle set for transformation from
the diffusion frame to the interaction-$i$ frame in rotamer $\alpha$. Thus, in
addition to report on the kinetics of the exchange (populations and jump
rates), the correlation function for rotamer jump also depends on the
resulting distinct orientations of the interaction frames in the diffusion
frame. This can be particularly critical when the overal diffusion tensor is
anisotropic, as discussed below.

#### Order parameter for rotamer jumps.

The internal correlation function for rotamer jumps in Eq. 61 does not cancel
out only for $n=0$ when $t\rightarrow\infty$, so that the squared order
parameter for rotamer jumps is, after using Eq. 10 [54]:

$\mathcal{S}_{J}^{2}(i,j)=\sum_{\alpha=1}^{N}\sum_{\beta=1}^{N}p_{\alpha}^{eq}p_{\beta}^{eq}\mathcal{P}_{2}(\cos\theta_{i_{\alpha},j_{\beta}}),$
(63)

where $\theta_{i_{\alpha},j_{\beta}}$ is the angle between interaction $i$ in
rotamer $\alpha$ and interaction $j$ in rotamer $\beta$.

###  3.4 Rotation on a cone

The diffusion on a cone can be used, for example, to model the rotation of a
methyl group around its symmetry axis under the assumption that thermal energy
is much larger than the energy barrier for the rotation, so that the
probability distribution for the position of hydrogen atoms is uniform. Methyl
groups are now widely used in biomolecular NMR for their favorable relaxation
properties that precisely originate from the fast rotation around their
symmetry axis [55, 56, 57] .In particular, transverse-relaxation optimized
spectroscopy (TROSY) of methyl groups has made it possible to investigate
large proteins and protein complexes. This approach is based on the
cancelation between proton-proton dipolar auto- and cross-relaxation
mechanisms [55, 58] for which fast methyl rotation is essential [59]. In this
context, force fields used in MD simulations of proteins have been modified
recently using NMR experimental data to correctly account for the energy
barrier for the rotation [60, 61]. The solution of the Master equation (Eq.
16) has first been published by D. Wallach [38].
In order to describe the diffusion on a cone motion, an additional frame needs
to be introduced on top of the diffusion and interaction frames. Similarly to
the rotamer jumps, we call this frame the system frame (SF). Its main axis is
pointing along the axis of rotation, with the x- and y-axis rotating such that
the transformation from the SF to the interaction frame(s) is time-
independent. The Euler angle for transformation from the diffusion frame to
the system frame is
$\Omega_{D,SF}=\\{\alpha_{D,SF},\beta_{D,SF},\gamma_{D,SF}(t)\\}$ where only
the third angle is time-dependent. The correlation function is then:

$\displaystyle C_{i,j}(t)=$
$\displaystyle\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}e^{i\alpha_{D,SF}(a-a^{\prime})}d_{a,b}(\beta_{D,SF})d_{a^{\prime},b^{\prime}}(\beta_{D,SF})\times$
(64) $\displaystyle\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle\mathcal{D}_{b,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{b^{\prime},0}^{(2)}(\Omega_{SF,j}).$

The Master equation is written using the angular momentum operator
$L_{rot}^{2}$:

$\frac{\partial}{\partial
t}p(\gamma,t)=-D_{rot}L_{rot}^{2}p(\gamma,t)=D_{rot}\frac{\partial^{2}}{\partial\gamma^{2}}p(\gamma,t),$
(65)

where $D_{rot}$ is the rotational diffusion coefficient. The ensemble average
is calculated as:

$\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle=\int_{0}^{2\pi}\int_{0}^{2\pi}e^{ib\gamma_{0}}e^{-ib^{\prime}\gamma}p(\gamma_{0})p(\gamma,t|\gamma_{0},0)\mathrm{d}\gamma_{0}\mathrm{d}\gamma,$
(66)

where $p(\gamma_{0})=\frac{1}{2\pi}$ (all orientations are equi-probable), and
the conditional probability is expressed as:

$p(\gamma,t|\gamma_{0},0)=\sum_{n}\varphi_{n}^{*}(\gamma_{0})\varphi_{n}(\gamma)e^{-\lambda_{n}t},$
(67)

where $\varphi_{n}$ is the eigenfunction associated to the eigenvalue
$\lambda_{n}$ of the rotation operator $D_{rot}L_{rot}^{2}$. Eigenfunctions
and eigenvalues are given by:

$\displaystyle\varphi_{n}(\gamma)=\frac{1}{\sqrt{2\pi}}e^{in\gamma},$ (68)
$\displaystyle\lambda_{n}=D_{rot}n^{2}.$

Inserting in Eq. 66 leads to:

$\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle=\frac{1}{4\pi^{2}}\sum_{n}e^{-D_{rot}n^{2}t}\int_{0}^{2\pi}\int_{0}^{2\pi}e^{i\gamma_{0}(b-n)}e^{-i\gamma(b^{\prime}-n)}\mathrm{d}\gamma_{0}\mathrm{d}\gamma,$
(69)

which does not cancel out only for $n=b=b^{\prime}$, and leads to the
correlation function:

$\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle=\delta_{bb^{\prime}}e^{-D_{rot}b^{2}t}.$
(70)

The total correlation function is then:

$\displaystyle C_{i,j}(t)=$
$\displaystyle\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}e^{-D_{rot}b^{2}t}\times$
(71) $\displaystyle
e^{i\alpha_{D,SF}(a-a^{\prime})}d_{a,b}(\beta_{D,SF})d_{a^{\prime},b}(\beta_{D,SF})\mathcal{D}_{b,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{b,0}^{(2)}(\Omega_{SF,j}).$

#### Order parameter for rotation on a cone.

The only non-vanishing term in Eq. 70, when $t\rightarrow\infty$, is obtained
for $b=0$, so that the squared order-parameter for rotation on a cone is:

$\mathcal{S}_{R}^{2}(i,j)=\mathcal{P}_{2}(\cos\beta_{SF,i})\mathcal{P}_{2}(\cos\beta_{SF,j}).$
(72)

In the case of methyl-group rotation, with ideal tetrahedral geometry, it
leads to the well-known result
$\mathcal{S}_{met}^{2}(\mathrm{DD}_{\mathrm{CH}},\mathrm{DD}_{\mathrm{CH}})=1/9$,
where $\mathrm{DD}_{\mathrm{CH}}$ stands for the C-H dipole-dipole (DD)
interaction.

#### CSA/DD cross-correlation.

Figure 1: Correlation function for methyl group rotation. Auto-correlation
function for H-H (a) and C-H (b) motions. Calculated correlation functions are
shown in blue, and MF fits in dash orange. The correlation functions are shown
as a function of $D_{rot}\times t$ such that they are independent of the
diffusion constant $D_{rot}$. The expected and fitted MF parameters are shown
on each pannel.

In methyl groups, carbon-13 Chemical Shift Anisotropy (CSA) and DD cross-
correlated cross-relaxation rates can be measured [62]. In this case, the
carbon-13 CSA doesn’t undergo a rotation motion and the correlation function
in Eq. 64 is modified to:

$\displaystyle C_{\mathrm{CSA,DD}}(t)=\frac{1}{5}$
$\displaystyle\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,CSA})e^{ia^{\prime}\alpha_{D,SF}}d_{a^{\prime},b^{\prime}}(\beta_{D,SF})\times$
(73) $\displaystyle\langle
e^{ib^{\prime}\gamma_{D,SF}(t)}\rangle\mathcal{D}_{b^{\prime},0}^{(2)}(\Omega_{SF,DD}).$

Due to the methyl rotation, $\langle e^{ib^{\prime}\gamma_{D,SF}(t)}\rangle$
equals 0 on average unless $b^{\prime}=0$ such that:

$C_{\mathrm{CSA,DD}}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}e^{ia^{\prime}\alpha_{D,SF}}\mathcal{P}_{2}(\cos\beta_{SF,DD})\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,CSA})d_{a^{\prime},0}(\beta_{D,SF}),$
(74)

where $\mathcal{P}_{2}(\cos\beta_{SF,DD})=-1/3$ for methyl groups with perfect
tetrahedral geometry.

#### C-H and H-H auto-correlation.

The DD contribution to carbon-13 relaxation rates only reports on the C-H
correlation, while for protons (or deuterons in the case of specific labeling)
it reports on both C-H and H-H correlations. For the sake of simplicity, we
assume here an isotropic overall diffusion tensor and a perfect tetrahedral
geometry for the methyl group to expand the auto-correlation functions as:

$\begin{dcases}C_{CH}(t)=\frac{1}{5}e^{-t/\tau_{c}}\left(\frac{1}{9}+\frac{8}{27}e^{-D_{rot}t}+\frac{16}{27}e^{-4D_{rot}t}\right),\\\
C_{HH}(t)=\frac{1}{5}e^{-t/\tau_{c}}\left(\frac{1}{4}+\frac{3}{4}e^{-4D_{rot}t}\right).\end{dcases}$
(75)

It is interesting to note that the H-H auto-correlation function is written in
the form of the MF correlation function with generalized squared order
parameter $\mathcal{S}_{HH}^{2}=1/4$ and correlation time
$\tau_{HH}=1/(4D_{rot})$. The C-H auto-correlation function can be
approximated to a single-decay exponential with generalized squared order
parameter $\mathcal{S}_{CH}^{2}=1/9$ and correlation time defined to keep the
area of the correlation function constant:

$\tau_{CH}=\frac{1}{1-\mathcal{S}_{CH}^{2}}\sum_{b\neq
0}\frac{1}{b^{2}D_{rot}}\left[d_{b,0}(\beta_{met})\right]^{2}=2\tau_{HH},$
(76)

where $\beta_{met}=109.47^{\circ}$. Thus, C-H and H-H auto-correlation lead to
a methyl rotation correlation different, in theory, by a factor of 2 when
modeled using MF correlation functions. This is shown in Fig. 1 where the
calculated internal correlation function for methyl rotation is perfectly
fitted with the MF correlation function for H-H auto-correlation, but show
some deviations in the case of C-H auto-correlation with pseudo-correlation
time $D_{rot}\times\tau$ significantly different between the two types of
correlation despite the fact they report on the same motion. This has no
visible effects on the carbon longitudinal and transverse relaxation rates
(Fig. S1.a,b) and slight deviations in carbon-proton DD cross-relaxation rates
at high magnetic-fields (Fig. S1.c).
When the energy barrier for the rotation is significantly larger than thermal
energy [60, 61], a jump motion between three discrete positions might be more
physically relevant [39, 51, 52]. The rotation on a cone model if however
mathematically simpler, and for effective correlation time for methyl rotation
smaller than 50 ps, the two models are virtually indistinguishable in the
simple case evaluated in Fig. S2. When a three-site jump model is considered,
one should keep in mind that diffusion within each methyl rotamer energy well
is expected to cover a rather wide angle. This can be accounted for by the
restricted rotation in a cone model.

#### Restricted rotation on a cone

In the case where the rotation is restricted to a portion of the cone surface,
for example and without loss of generality for angles $\varphi_{D,SF}$ between
$\pm\gamma_{c}$, the Master equation Eq. 65 has reflecting boundary condition:

$\frac{\partial p(\gamma,t)}{\partial\gamma}|_{\gamma=\pm\gamma_{c}}=0.$ (77)

Solutions for the Master equation have been reported by Wittebort and Szabo
[39]:

$\displaystyle\varphi_{n}(\gamma)=\cos\frac{n\pi(\gamma-\gamma_{c})}{2\gamma_{c}},$
(78) $\displaystyle\lambda_{n}=D_{rot}\frac{n^{2}\pi^{2}}{4\gamma_{c}^{2}},$

leading to:

$\displaystyle\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle=\frac{1}{4\gamma_{c}^{2}}\sum_{n}e^{-D_{rot}\frac{n^{2}\pi^{2}}{4\gamma_{c}^{2}}t}\int_{-\gamma_{c}}^{+\gamma_{c}}e^{i\gamma_{0}b}\cos\frac{n\pi(\gamma_{0}-\gamma_{c})}{2\gamma_{c}}\mathrm{d}\gamma_{0}$
$\displaystyle\times$ (79)
$\displaystyle\int_{-\gamma_{c}}^{+\gamma_{c}}e^{-i\gamma
b^{\prime}}\cos\frac{n\pi(\gamma-\gamma_{c})}{2\gamma_{c}}\mathrm{d}\gamma$ ,

which integrates into:

$\displaystyle\langle
e^{i(b\gamma_{D,SF}(0)-b^{\prime}\gamma_{D,SF}(t))}\rangle=\frac{bb^{\prime}\gamma_{c}^{2}}{2(b^{2}\gamma_{c}^{2}-\frac{n^{2}\pi^{2}}{4})(b^{\prime
2}\gamma_{c}^{2}-\frac{n^{2}\pi^{2}}{4})}\times$ (80) $\displaystyle\left(\sin
b\gamma_{c}\sin b^{\prime}\gamma_{c}(1+(-1)^{n})+\cos b\gamma_{c}\cos
b^{\prime}\gamma_{c}(1-(-1)^{n})\right),$

and leads to the correlation function:

$\displaystyle
C_{i,j}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime}=-2}^{2}\sum_{n=0}^{+\infty}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,SF})\mathcal{D}_{a^{\prime},b^{\prime}}^{(2)}(\Omega_{D,SF})\mathcal{D}_{b,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{b^{\prime},0}^{(2)}(\Omega_{SF,j})\times$
(81)
$\displaystyle\frac{bb^{\prime}\gamma_{c}^{2}e^{-D_{rot}\frac{n^{2}\pi^{2}}{4\gamma_{c}^{2}}t}}{2(b^{2}\gamma_{c}^{2}-\frac{n^{2}\pi^{2}}{4})(b^{\prime
2}\gamma_{c}^{2}-\frac{n^{2}\pi^{2}}{4})}(\sin b\gamma_{c}\sin
b^{\prime}\gamma_{c}(1+(-1)^{n})+\cos b\gamma_{c}\cos
b^{\prime}\gamma_{c}(1-(-1)^{n})).$

The sum over $n$ has been reported to converge rapidly [39] such that only the
first 5 terms can be considered when calculating the correlation function.
Restricted rotation on a cone can also be treated using the Gaussian Axial
Fluctuation (GAF) model which assumes a Gaussian shape for the probability
distribution [63, 64, 65]. Such correlation function to analyse methyl
rotation could be combined with a jump-type of model to represent the dynamic
as a three-site hope with restricted diffusion within each site.

###  3.5 Wobbling in a cone

The last type of motion that we will consider here is the wobbling in a cone.
In this model, the bond vector undergoes restricted diffusion within a solid
angle that is constrained by a cone. It can be useful to model the motion of
any bonds such as a 15N-1H pair, or C-C bonds along protein side-chains. The
correlation function for this type of motion was initially introduced in
fluorescence [66] and dielectric spectroscopy [67]. As will be seen below, the
Master equation can be solved analytically, but the complexity of the solution
has led to the introduction of approximated forms for the correlation function
by the group of A. Ikegami [66], and later obtained similarly by G.Lipari and
A. Szabo [68, 69]. which guided the latter two towards the MF correlation
function [13]. However, the initial treatment of this model and the following
studies were performed under the assumption that the interaction of interest
was undergoing the wobbling motion, which makes it a priori not applicable to
the study of motions of a CC bond with 1H-13C relaxation. A. Kumar solved the
Master equation when the frame undergoing the diffusion in a cone is not the
interaction frame [70, 71]. We will follow his approach, as well as the one of
C. Wang and R. Pecora [40] who extended the work of M. Warchol and W. Vaughan
[67]. In the construction of this model, the diffusing bond vector also
undergoes rotation around itself, which is not likely to correctly reproduce
motions of side-chains C-C bond vectors. We introduce here a way to overcome
this difficulty.
For the description of this motion, we introduce the Wobbling Frame (WF) which
main axis points along the symmetry axis of the cone, and x-axis can point in
any direction. Similarly to what was presented in the previous sections, we
also include a System Frame (SF) which undergoes the diffusion motion.
Assuming the wobbling and global tumbling are uncorrelated, the correlation
function is written:

$\displaystyle C_{i,j}(t)=$
$\displaystyle\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime},c,c^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{WF,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{WF,SF},t)\rangle\times$
(82)
$\displaystyle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,WF})\mathcal{D}_{c,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{a^{\prime},b^{\prime}}^{(2)}(\Omega_{D,WF})\mathcal{D}_{c^{\prime},0}^{(2)}(\Omega_{SF,j}).$

The Master equation that solves the conditional probability is given by:

$\frac{\partial}{\partial t}p(\Omega,t)=-D_{W}L_{W}^{2}p(\Omega,t),$ (83)

where $D_{W}$ is the diffusion coefficient for the wobble motion and
$L_{W}^{2}$ is the angular momentum operator:

$L_{W}^{2}=-\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial}{\partial\theta}\right)-\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial\varphi^{2}},$
(84)

with reflecting boundary condition at $\theta=\beta_{c}$, with $\beta_{c}$ the
cone semi-angle opening:

$\frac{\partial}{\partial\theta}p(\Omega,t)|_{\theta=\beta_{c}}=0.$ (85)

#### Solving the Master equation.

The eigen-equation is written:

$D_{W}L_{W}^{2}c(\theta,\varphi)=Ec(\theta,\varphi),$ (86)

and we assume that the function $c$ can be separated into the product:

$c(\theta,\varphi)=\Theta(\theta)\Phi(\varphi).$ (87)

This leads to the following differential equation for $\Phi$:

$\frac{\partial^{2}}{\partial\varphi^{2}}\Phi(\varphi)\propto\Phi(\varphi).$
(88)

Following the work of Wang and Pecora [40], we write the complete set of
functions $\Phi$ as:

$\Phi_{m}(\varphi)=A_{m}\cos(m\varphi)+B_{m}\sin(m\varphi),$ (89)

where $m\in\ \mathbb{N}$, and $A_{m}$ and $B_{m}$ are coefficients associated
to the eigenvalue $m$. It follows:

$\frac{\partial^{2}}{\partial\varphi^{2}}\Phi_{m}(\varphi)=-m^{2}\Phi_{m}(\varphi),$
(90)

and we can write the differential equation for $\Theta$ after the variable
change $\mu\rightarrow\cos\theta$:

$2\mu\frac{\partial\Theta}{\partial\mu}-(1-\mu^{2})\frac{\partial^{2}\Theta}{\partial\mu^{2}}+\frac{m^{2}}{1-\mu^{2}}\Theta-\frac{E}{D_{W}}\Theta=0.$
(91)

We now define $E=D_{W}\nu_{m}(\nu_{m}+1)$ to obtain:

$(1-\mu^{2})\frac{\partial^{2}\Theta}{\partial\mu^{2}}-2\mu\frac{\partial\Theta}{\partial\mu}+\left(\nu_{m}(\nu_{m}+1)-\frac{m^{2}}{1-\mu^{2}}\right)\Theta=0.$
(92)

Solutions of the equations of the type Eq. 92 for $\mu\in[0,1]$ are given by
the Legendre associated functions $P_{\nu_{m}}^{m}$ of degree $\nu_{m}$ and
order $m$. There are an infinite number of solutions and the eigenfunctions
and eigenvalues are given by $P_{\nu_{m,n}}^{m}$ and
$D_{W}\nu_{m,n}(\nu_{m,n}+1)$ respectively, where the $\nu_{m,n}$ fulfill the
boundary condition Eq. 85. We choose the indices $n$ such that
$\nu_{m,n}<\nu_{m,n+1}$. It can be noted that $\nu_{0,0}=0$. The conditional
probability $p(\Omega,t|\Omega_{0},0)$ is now given by:

$p(\Omega,t|\Omega_{0},0)=\sum_{m=0}^{+\infty}\sum_{n=0}^{+\infty}e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}P_{\nu_{m,n}}^{m}(\cos\theta)(A_{m,n}\cos(m\varphi)+B_{m,n}\sin(m\varphi)),$
(93)

where we explicitely show the dependence in $m$ and $n$ of the coefficients
$A$ and $B$.

#### Legendre associated functions.

We give here the analytical expression of Legendre associated functions for
$m\geq 0$:

$P_{\nu_{m,n}}^{m}(z)=(-1)^{m}\frac{\Gamma(\nu_{m,n}+m+1)}{\Gamma(\nu_{m,n}-m+1)}\frac{(1-z)^{m/2}}{2^{m}m!}_{2}{F}{{}_{1}}(-\nu_{m,n}+m,\nu_{m,n}+m+1,m+1,\frac{1-z}{2}),$
(94)

where $\Gamma$ is the Gamma function defined for all $z\neq 0$:

$\Gamma(z)=\int_{0}^{\infty}e^{-t}t^{z-1}dt,$ (95)

and ${}_{2}{F}{{}_{1}}$ is the Hypergeometric function:

${}_{2}{F}{{}_{1}}(\alpha,\beta,\gamma,z)=1+\frac{\alpha\beta}{\gamma}\frac{z}{1!}+\frac{\alpha(\alpha+1)\beta(\beta+1)}{\gamma(\gamma+1)}\frac{z^{2}}{2!}+\cdot\cdot\cdot.$
(96)

The Hypergeometric function can be written using the Gamma function in a
compact form as:

$\begin{dcases}{}_{2}{F}{{}_{1}}(\alpha,\beta,\gamma,z)=1,\,\,\,\alpha=0\,\,\mathrm{or}\,\,\beta=0,\\\
{}_{2}{F}{{}_{1}}(\alpha,\beta,\gamma,z)=\sum_{k=0}^{\infty}\frac{\Gamma(\alpha+k)\Gamma(\beta+k)\Gamma(\gamma)}{\Gamma(\alpha)\Gamma(\beta)\Gamma(\gamma+k)}\frac{z^{k}}{k!},\alpha\neq
0\,\,\mathrm{and}\,\,\beta\neq 0.\end{dcases}$ (97)

The definition of legendre associated function can be extended to orders
having negative values:

$P_{\nu_{-m,n}}^{-m}(z)=(-1)^{m}\frac{\Gamma(\nu_{m,n}-m+1)}{\Gamma(\nu_{m,n}+m+1)}P_{\nu_{m,n}}^{m}(z),m\geq
0.$ (98)

#### Solving the boundary condition.

The conditional probability at time $t=0$ is:

$\displaystyle p(\Omega,t=0|\Omega_{0},0)=$
$\displaystyle\delta(\mu-\mu_{0})\delta(\varphi-\varphi_{0}),$ (99)
$\displaystyle=$
$\displaystyle\sum_{m=0}^{+\infty}\sum_{n=0}^{+\infty}P_{\nu_{m,n}}^{m}(\mu_{0})(A_{m,n}\cos(m\varphi_{0})+B_{m,n}\sin(m\varphi_{0})),$

where we set $\mu=\cos\theta$. Now, we can calculate the following two
integrals:

$\displaystyle I_{c}=\int_{\mu_{c}}^{1}d\mu\int_{0}^{2\pi}d\varphi
p(\Omega,t=0|\Omega_{0},0)P_{\nu_{m^{\prime},n^{\prime}}}^{m^{\prime}}(\mu)\cos
m^{\prime}\varphi,$ (100) $\displaystyle
I_{s}=\int_{\mu_{c}}^{1}d\mu\int_{0}^{2\pi}d\varphi
p(\Omega,t=0|\Omega_{0},0)P_{\nu_{m^{\prime},n^{\prime}}}^{m^{\prime}}(\mu)\sin
m^{\prime}\varphi.$

On one side, we have:

$\displaystyle I_{c}=P_{\nu_{m^{\prime},n^{\prime}}}^{m^{\prime}}(\mu_{0})\cos
m^{\prime}\varphi_{0},$ (101) $\displaystyle
I_{s}=P_{\nu_{m^{\prime},n^{\prime}}}^{m^{\prime}}(\mu_{0})\sin
m^{\prime}\varphi_{0}.$

On the other side, we can write:

$\displaystyle I_{c}=H_{m^{\prime},n^{\prime}}(\mu_{c})$
$\displaystyle\int_{0}^{2\pi}d\varphi\cos
m^{\prime}\varphi(A_{m^{\prime},n^{\prime}}\cos(m\varphi)+B_{m^{\prime},n^{\prime}}\sin(m^{\prime}\varphi)),$
(102) $\displaystyle I_{s}=H_{m^{\prime},n^{\prime}}(\mu_{c})$
$\displaystyle\int_{0}^{2\pi}d\varphi\sin
m^{\prime}\varphi(A_{m^{\prime},n^{\prime}}\cos(m\varphi)+B_{m^{\prime},n^{\prime}}\sin(m^{\prime}\varphi)),$

where $\mu_{c}=\cos\beta_{c}$ and $H_{m,n}$ is defined from the orthogonality
condition of Legendre associated functions [40]:

$\int_{\mu_{c}}^{1}P_{\nu_{m,n}}^{m}(\mu)P_{\nu_{m^{\prime},n^{\prime}}}^{m^{\prime}}(\mu)d\mu=\delta_{nn^{\prime}}\delta_{mm^{\prime}}H_{m,n}(\mu_{c}).$
(103)

After integration and identification, we can express the two coefficients
$A_{m,n}$ and $B_{m,n}$:

$\displaystyle A_{m,n}=\frac{\cos m\varphi_{0}}{\epsilon\pi
H_{m,n}(\mu_{c})}P_{\nu_{m,n}}^{m}(\mu_{0}),$ (104) $\displaystyle
B_{m,n}=\frac{\sin m\varphi_{0}}{\pi
H_{m,n}(\mu_{c})}P_{\nu_{m,n}}^{m}(\mu_{0}),$

where $\mu_{0}=\cos\theta_{0}$ and:

$\epsilon=\begin{dcases}2,\,\,\,m=0\\\ 1,\,\,\,m=\\{1,2\\}\end{dcases}.$ (105)

#### Conditional probability.

The conditional probability is now written:

$p(\Omega,t|\Omega_{0},0)=\sum_{m=0}^{+\infty}\sum_{n=0}^{+\infty}X_{m,n},$
(106)

where we define:

$X_{m,n}=e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}\frac{P_{\nu_{m,n}}^{m}(\mu)P_{\nu_{m,n}}^{m}(\mu_{0})}{\pi
H_{m,n}(\mu_{c})}\left(\frac{\cos m\varphi_{0}}{\epsilon}\cos m\varphi+\sin
m\varphi_{0}\sin m\varphi\right).$ (107)

We have:

$X_{0,n}=\frac{1}{2}\frac{P_{\nu_{0,n}}^{0}(\mu)P_{\nu_{0,n}}^{0}(\mu_{0})}{\pi
H_{0,n}(\mu_{c})}e^{-D_{W}\nu_{0,n}(\nu_{0,n}+1)t}$ (108)

and, for $m\neq 0$:

$X_{m,n}=e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}\frac{P_{\nu_{m,n}}^{m}(\mu)P_{\nu_{m,n}}^{m}(\mu_{0})}{\pi
H_{m,n}(\mu_{c})}\left(\cos m\varphi_{0}\cos m\varphi+\sin m\varphi_{0}\sin
m\varphi\right).$ (109)

We use the following properties for Legendre associated functions and their
roots:

$\displaystyle\nu_{-m,n}=\nu_{m,n},$ (110)
$\displaystyle\frac{P_{\nu_{-m,n}}^{-m}(\mu_{1})P_{\nu_{-m,n}}^{-m}(\mu_{2})}{H_{-m,n}(\mu_{c})}=\frac{P_{\nu_{m,n}}^{m}(\mu_{1})P_{\nu_{m,n}}^{m}(\mu_{2})}{H_{m,n}(\mu_{c})},$

to write:

$p(\Omega,t|\Omega_{0},0)=\frac{1}{2}\sum_{m=-\infty}^{+\infty}\sum_{n=0}^{+\infty}e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}\frac{P_{\nu_{m,n}}^{m}(\mu)P_{\nu_{m,n}}^{m}(\mu_{0})}{\pi
H_{m,n}(\mu_{c})}e^{im\varphi}e^{-im\varphi_{0}}.$ (111)

The factor $1/2$ comes from the equality $X_{-m,n}=X_{m,n}$ and the extension
of the sum on the index $m$ from $-\infty$ to $+\infty$, and the expression of
$X_{0,n}$ (Eq. 108). C. Wang and R. Pecora introduced the pseudo-spherical
harmonics [40]:

$Y_{\nu_{m,n}}^{m}(\Omega)=\frac{1}{\sqrt{2\pi
H_{m,n}(\mu_{c})}}P_{\nu_{m,n}}^{m}(\mu)e^{im\varphi},$ (112)

and wrote the conditional probability in a compact form:

$p(\Omega,t|\Omega_{0},0)=\sum_{m=-\infty}^{+\infty}\sum_{n=0}^{+\infty}e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}Y_{\nu_{m,n}}^{m}(\Omega)Y_{\nu_{m,n}}^{m*}(\Omega_{0}).$
(113)

#### Initial angle probability.

The conditional probability does not cancel out when $t\rightarrow\infty$ only
for $\nu_{m,n}=0$, that is $m=n=0$:

$p(\Omega_{0})=Y_{0}^{0}(\Omega)Y_{0}^{0}(\Omega_{0})=\frac{P_{0}^{0}(\mu)P_{0}^{0}(\mu_{0})}{2\pi
H_{0,0}(\mu_{c})}.$ (114)

It is then straightforward to show:

$\begin{dcases}p(\Omega_{0})=\frac{1}{2\pi(1-\mu_{c})},\,\,\,\theta_{0}\in[0,\beta_{c}],\\\
p(\Omega_{0})=0,\,\,\,\theta_{0}>\beta_{c}.\end{dcases}$ (115)

The correlation function for wobbling in a cone then reads:

$\displaystyle\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{WF,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{WF,SF},t)\rangle=\frac{1}{4\pi^{2}(1-\mu_{c})}\sum_{m=-\infty}^{+\infty}\sum_{n=0}^{+\infty}\frac{e^{-D_{W}\nu_{m,n}(\nu_{m,n}+1)t}}{H_{m,n}(\mu_{c})}\int_{0}^{2\pi}d\varphi_{0}\int_{0}^{2\pi}d\varphi\times$
(116)
$\displaystyle\int_{0}^{\beta_{c}}\sin\theta_{0}d\theta_{0}\int_{0}^{\beta_{c}}\sin\theta
d\theta
P_{\nu_{m,n}}^{m}(\cos\theta_{0})P_{\nu_{m,n}}^{m}(\cos\theta)\mathcal{D}_{b,c}^{(2)*}(\varphi_{0},\theta_{0},-\varphi_{0})\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\varphi,\theta,-\varphi)e^{im(\varphi-\varphi_{0})}.$

The third Euler angles in the Wigner matrices equal the opposite of the first
Euler angles in order for this correlation function to only account for the
motion of the C-C bond in a cone, and not any rotational motions around the
z-axis of the SF: if the third angle was set to 0, the orientation of the
x-axis of the system frame with respect to the wobbling frame main axis would
not change upon diffusion in the cone thus introducing an additional rotation
motion. When it was initially developped, this model was accounting for
motions of the interaction of interest [66, 67]. In this situation, the second
index of the Wigner matrices equals 0, and the third Euler angle does not play
any role. This is not the case when the interaction frame is not the one
undergoing the wobbling motion, and introducing this angle is essential to
best describe the motions.
Expanding the expression of the Wigner matrices and performing the integration
on $\varphi_{0}$ and $\varphi$ leads to the condition
$m=b-c=b^{\prime}-c^{\prime}$. Note that in the previous models, the condition
$m=b=b^{\prime}$ was obtained, which is correct when the interaction frame is
the diffusing frame since $c=c^{\prime}=0$ [40, 70, 71]. We finally can write
the total correlation function as:

$\displaystyle C_{i,j}(t)=\frac{1}{5}$
$\displaystyle\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime},c,c^{\prime}=-2}^{2}\delta_{b-c,b^{\prime}-c^{\prime}}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\sum_{n}\frac{e^{-D_{W}\nu_{b-c,n}(\nu_{b-c,n}+1)t}}{(1-\mu_{c})H_{b-c,n}(\mu_{c})}I_{b,c}^{n}(\beta_{c})I_{b^{\prime},c^{\prime}}^{n}(\beta_{c})\times$
(117)
$\displaystyle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,WF})\mathcal{D}_{c,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{a^{\prime},b^{\prime}}^{(2)}(\Omega_{D,WF})\mathcal{D}_{c^{\prime},0}^{(2)}(\Omega_{SF,j}),$

where:

$I_{b,c}^{n}(\beta_{c})=\int_{0}^{\beta_{c}}\sin\theta
d_{b,c}(\theta)P_{\nu_{b-c,n}}^{b-c}(\cos\theta)d\theta.$ (118)

#### Order parameter for wobbling in a cone.

The correlation function for wobbling is non-zero when $t\rightarrow+\infty$
for $b-c=n=0$. Thus, we have:

$\mathcal{S}_{W}^{2}(i,j)=\sum_{b=-2}^{2}\left(\frac{I_{b,b}^{0}(\beta_{c})}{1-\mu_{c}}\right)^{2}\mathcal{D}_{b,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{b,0}^{(2)}(\Omega_{SF,j}).$
(119)

In the particular case of auto-correlation and where the PAF of the
interactions are undergoing the wobbling motions, we have
$\Omega_{SF,i}=\Omega_{SF,j}=\\{0,0,0\\}$ and the order parameter simplifies
into the equation reported by K. Kinosita et al. [66]:

$\mathcal{S}_{W}^{2}(i_{W},i_{W})=\frac{1}{4}\cos^{2}\beta_{c}(1+\cos\beta_{c})^{2},$
(120)

where we write $i_{W}$ to indicate that the interaction $i$ is diffusing in a
cone.

###  3.6 Contribution of each motion to relaxation

Figure 2: Contribution of rotamer jump and wobbling in a cone motions to
relaxation. Carbon R1 (a), R2 (b) and carbon-proton $\sigma^{\mathrm{NOE}}$
(c) at 14.1 T in a hypothetical 13C1H2H2 valine methyl group, and associated
spectral density function for the Cβ-Cγ1 bond auto-correlations. Dashed lines
show the value of the relaxation rates and spectral density function without
wobbling. Relaxation rates (a-c) and spectral density function (d) calculated
in the abscence of rotamer jump are shown in green. In the presence of rotamer
jump, the case where all rotamer populations are equal (blue) and unequal
(orange) are distinguished. In this latter case, populations are $p_{1}=0.7$,
$p_{2}=0.2$ and $p_{1}=0.1$. Calculations are shown as a function of the
correlation time for wobbling, a function of the wobbling diffusion constant
and cone semi-angle opening, set here to $\beta_{c}=15$ deg, while the
diffusion constant $D_{W}$ is varied from $10^{6}$ to $10^{10}$ s-1. The blue
and orange vertical arrows indicate the values of the correlation time for
rotamer jump, respectively when populations are equal and unequal. The
spectral density functions are shown for $D_{W}=10^{8}$ s-1.

We evaluate here how each motion contributes to relaxation. Having an
exhaustive analysis is not feasible and we will focus on typical relaxation
rates: 13C-R1, 13C-R2 and 13C-1H-$\sigma^{\mathrm{NOE}}$ in a 13C1H2H2 methyl
group as encountered in a valine side-chain. The overall diffusion is supposed
to be isotropic with global tumbling correlation time $\tau_{c}=10$ ns. The
methyl rotation diffusion coefficient is set to $D_{rot}=5\times 10^{10}$ s-1,
which leads to a correlation time for rotation of $\tau_{rot}=10$ ps. When a
jump-type of motion is added, we suppose that the valine side-chain can jump
between three possible rotamer states, and we distinguish the situation where
each rotamers are equaly populated with equal exchange rate
$k_{ex}^{sym}=1/3\times 10^{9}$ s-1 between the different conformers (which
leads to a correlation time for exchange of $\tau_{ex}^{sym}=1$ ns), from the
situation where $p_{1}=0.7$, $p_{2}=0.2$ and $p_{3}=0.1$ (the state
numerotation is irrelevant since the global diffusion tensor is isotropic) and
all exchange rates $k_{ij}^{asym}=1/3\times 10^{9}$ s-1 with $i<j$ and
$k_{ij}^{asym}$ with $i>j$ is calculated to satisfy the microscropy
reversibility condition (Eq. 53) (which leads to a correlation time for
exchange of $\tau_{ex}^{asym}=1.8$ ns). The Euler angles associated to the
jump motion are $\Omega_{J}=\\{2n\pi/3,\beta_{CC},0\\}$ where $n=0,1,2$ and
$\beta_{CC}=76$ deg, which is typical for carbon side-chains. A wobbling in a
cone motion of the Cβ-Cγ1 bond is also considered, when mentioned. The cone
semi-angle opening is set to 15 deg (results for opening angles of 5, 30 and
60 deg are shown in Fig. S3), and the associated diffusion constant $D_{W}$ is
varied from $10^{6}$ to $10^{10}$ s-1, leading to variations of the wobbling
correlation time from 1 ps to 10 ns approximately for this cone angle opening.
Relaxation rates are calculated at a magnetic field of 14.1 T, and the
carbon-13 CSA is set to 20 ppm [72].
The rotamer exchange contributes significantly to the relaxation (compare
dashed lines in Fig. 2), with an R1 and $\sigma^{\mathrm{NOE}}$ increase by a
factor 1.5 approximately, while R2 decrease by a factor up to 3.8 when rotamer
jumps with equal populations are introduced. The effect of rotamer jumps is
reduced for R2 and $\sigma^{\mathrm{NOE}}$ when the populations are not equal,
and increased for R1. These effects can be rationalized by inspecting the
spectral density functions (Fig. 2d). In the presence of internal motions, the
spectral density function tends to be higher at frequencies in the range
$10^{9}$-$10^{10}$ rad.s-1, which are frequencies R1 and
$\sigma^{\mathrm{NOE}}$ are sensitive to, leading to an increase of these two
rates. The higher contribution of medium-frequency motions to the spectral
density function goes in pair with a decrease in the contribution of low-
frequency motions, which is the main determinant in the value of the
transverse relaxation rate $R_{2}$, explaining the large effects calculated in
the presence of rotamer exchange.
On the other hand, the wobbling of the Cβ-Cγ1 bond has smaller effects on
relaxation (Fig. 2) for the choosen opening of the cone. Wobbling reduces the
transverse relaxation rate $R_{2}$, but the effect on the longitudinal auto-
and cross-relaxation rates depends on the value of the diffusion constant
$D_{W}$: fast diffusion (low $\tau_{W}$) leads to a decrease in both R1 and
$\sigma^{\mathrm{NOE}}$; we obtain an increase for lower values of $D_{W}$
(higher values of $\tau_{W}$); and hardly any effects for very slow diffusion.
These results can be rationalized using the same arguments as for the rotamer
jump by analyzing the evolution of the spectral density function and the
relative contribution of each frequency range to the value of relaxation
rates. It is interesting to notice the similarities between the evolution of
the longitudinal relaxation rate R1 and $\sigma^{\mathrm{NOE}}$ as a function
of $\tau_{W}$ and the sensitivities calculated in the detector approach which
characterizes the amount of motions in a given range of correlation times [73,
74, 25, 75].
The contribution of the wobbling motions increases with increasing cone semi-
angle opening and can even become larger than the contribution of the rotamer
jump (Fig. S3). This can be explained by the value of order parameters, which
act as weights for the Lorentzian terms of the spectral density containing
contribution only from the global tumbling. The order parameters for rotamer
exchange when populations are equal and unequal are 0.17 and 0.43
respectively. In comparison, the wobbling order parameter equals 0.99, 0.90,
0.65 and 0.14 for values of cone semi-angle opening of 5, 15, 30 and 60 deg.
Thus, for $\beta_{c}=60$ deg, the decrease in $R_{2}$ from either a rotamer
exchange with equal population or the wobbling are similar for
$\tau_{W}\approx\tau_{ex}^{sym}$ (Fig. S3.h). The conclusions also apply when
an alternative model of motion is used where the wobbling of the Cα-Cβ is
allowed but the Cβ-Cγ is fixed in the rotamer frame (Fig. S4). These two
models may be difficult to discrimitate with experimental data.
As a general rule of thumb, one can reasonably expect rotamer jumps to be the
major source of relaxation for methyl groups in aliphatic side-chains, since
the amplitudes of wobbling are usually limited to small cone angle openings.
This has already been noted by Wand and coworker [76], but of course has to be
carefully investigated before neglecting one motion, for example with the use
of MD simulations.

##  4 Correlated motions

So far, correlation functions were written assuming statistical independence
between each motion. This assumption does not hold in some cases. Correlation
functions for correlated overal tumbling and jumps of entire domain have been
presented in the late 2000s [77, 78, 79, 80]. In the case of methyl groups in
aliphatic side-chains, the methyl rotation or C-C wobbling can potentially
depend on the rotameric state, introducing correlation of internal motions
[81, 41, 82]. We show here how to write the correlation function in the
presence of such correlated motions.

###  4.1 Correlated jumps and methyl rotation

Figure 3: Contribution of correlated rotamer jumps and methyl rotation to
relaxation. Magnetic field variations of carbon-13 R1 (a), R2 (b) and
carbon-13-proton DD cross-relaxation rate (c) in a 13C1H2H2 methyl group
exchanging between two rotamers states and calculated for three distributions
of populations. Relaxation rates calculated with a unique diffusion
coefficient for methyl rotation are shown in dash, and relaxation rates
calculated by considering the correlation between rotamer jumps and methyl
rotation motions are shown in plain lines. The overall diffusion is isotropic
with global tumbling correlation time $\tau_{c}=10$ ns. The population average
methyl rotation coefficient is $10^{11}$ s-1 and the difference between the
two state is $5\times 10^{10}$ s-1. The rotamer 1 has the highest diffusion
coefficient. Carbon-13 CSA is set to 20 ppm.

Using the same notations as above, the correlation funtion for correlated
rotamer jumps and methyl rotation is:

$\displaystyle C_{i,j}(t)=\frac{1}{5}$
$\displaystyle\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b,b^{\prime},c,c^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle\times$
(121)
$\displaystyle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,J})\mathcal{D}_{a^{\prime},b^{\prime}}^{(2)}(\Omega_{D,J})\mathcal{D}_{c,0}^{(2)*}(\Omega_{SF,i})\mathcal{D}_{c^{\prime},0}^{(2)}(\Omega_{SF,j}).$

In the Euler angle set
$\Omega_{J,SF}=\\{\varphi_{J,SF},\theta_{J,SF},\phi_{J,SF}\\}$,
$\varphi_{J,SF}$ rotates the jump frame to the direction of the populated
rotamer and is time dependent, $\theta_{J,SF}$ rotates the jump frame to align
it on the rotamer frame, and $\phi_{J,SF}$ rotates the resulting frame to
align its x-axis along the direction of one C-H bond and is time dependent. We
assume that both interactions $i$ and $j$ undergo rotation such that
$\Omega_{SF,i}=\Omega_{SF,j}$. Note that in the case where only one
interaction undergoes rotation, as in methyl-groups CSA/DD cross-correlated
cross-relaxation, the resulting correlation function is independent from the
rotation (Eq. 74) such that there is no possible correlation between jumps and
rotation. The ensemble average can be expressed as:

$\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle=\sum_{\alpha,\beta=1}^{N}p_{\alpha}^{eq}\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle_{\beta\alpha},$
(122)

where the sums run over all accessible states, and the notation
$\langle\cdot\cdot\cdot\rangle_{\beta\alpha}$ indicates that the ensemble
average is calculated with initial state $\alpha$ and final state $\beta$. It
can be calculated as follows:

$\displaystyle\langle\cdot\cdot\cdot\rangle_{\beta\alpha}=$
$\displaystyle\int_{0}^{2\pi}d\phi_{0}\int_{0}^{2\pi}d\phi
p(\phi_{0})p(\\{\varphi_{J,SF_{\beta}},\theta_{J,SF_{\beta}},\phi\\},t|\\{\varphi_{J,SF_{\alpha}},\theta_{J,SF_{\alpha}},\phi_{0}\\},0)\times$
(123)
$\displaystyle\mathcal{D}_{b,c}^{(2)*}(\\{\varphi_{J,SF_{\alpha}},\theta_{J,SF_{\alpha}},\phi_{0}\\})\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\\{\varphi_{J,SF_{\beta}},\theta_{J,SF_{\beta}},\phi\\}),$

where $\varphi_{J,SF_{k}}$ and $\theta_{J,SF_{k}}$ are the values of the Euler
angles to transform from the jump frame to the system frame in rotamer $k$.
The probability $p(\phi_{0})$ is:

$p(\phi_{0})=\frac{1}{2\pi},$ (124)

and in order to calculate the conditional probability, we will use the
notation:

$p(\\{\varphi_{J,SF_{\beta}},\theta_{J,SF_{\beta}},\phi\\},t|\\{\varphi_{J,SF_{\alpha}},\theta_{J,SF_{\alpha}},\phi_{0}\\},0)=p(\beta,\phi,t|\alpha,\phi_{0},0).$
(125)

The Master equation that solves the conditional probability is a combination
of the Master equations for rotamer exchange and diffusion on a cone:

$\frac{d}{dt}p(\beta,\phi,t|\alpha,\phi_{0},0)=\sum_{\gamma=1}^{N}R_{\beta\gamma}p(\gamma,\phi,t|\alpha,\phi_{0},0)-D_{rot,\beta}L_{rot}^{2}p(\beta,\phi,t|\alpha,\phi_{0},0),$
(126)

where $R_{\beta\gamma}$ are elements of the exchange matrix $\mathcal{R}$
(i.e. exchange rate from state $\gamma$ to state $\beta$) and $D_{rot,\beta}$
is the diffusion constant for diffusion on a cone in state $\beta$. We solve
Eq. 126 by writing the conditional probability in terms of eigenfunctions of
the angular momentum operator $L_{rot}^{2}$:

$p(\beta,\phi,t|\alpha,\phi_{0},0)=\sum_{n}\frac{1}{2\pi}e^{in(\phi-\phi_{0})}c_{n}^{\beta\alpha}(t),$
(127)

which, after insertion and identification in Eq. 126, leads to the following
differential equation for the functions $c_{n}^{\beta\alpha}$:

$\frac{d}{dt}c_{n}^{\beta\alpha}(t)=-D_{rot,\beta}n^{2}c_{n}^{\beta\alpha}(t)+\sum_{\gamma=1}^{N}R_{\beta\gamma}c_{n}^{\gamma\alpha}(t).$
(128)

It can be written in matrix form as:

$\frac{d}{dt}C_{n}^{\alpha}(t)=\left(\mathcal{R}-n^{2}\mathcal{D}_{rot}\right)C_{n}^{\alpha}(t),$
(129)

where $\mathcal{D}_{rot}$ is a diagonal matrix containing the methyl rotation
diffusion coefficients as diagonal elements and $C_{n}^{\alpha}(t)$ is a
column vector containing the elements $c_{n}^{\beta\alpha}(t)$ for all states
$\beta$. Similarly to the treatment of rotamer jumps, we define the
symmetrized pseudo-exchange matrix as:

$\tilde{\mathcal{R}}_{rot,n}=\tilde{\mathcal{R}}-n^{2}\mathcal{D}_{rot},$
(130)

where $\tilde{\mathcal{R}}$ is the symmetrized exchange matrix (Eq. 57). Then,
the functions $c_{n}^{\beta\alpha}$ can be explicitely written as:

$c_{n}^{\beta\alpha}(t)=\sqrt{\frac{p_{\beta}^{eq}}{p_{\alpha}^{eq}}}\sum_{m=1}^{N}\tilde{X}_{\alpha}^{(n,m)}\tilde{X}_{\beta}^{(n,m)}e^{\lambda_{n,m}t},$
(131)

where $N$ is the number of states and $\tilde{X}_{\zeta}^{(m,n)}$ is the value
of the eigenvector $\tilde{X}^{(n,m)}$ associated to the eigenvalue
$\lambda_{n,m}$ and at coordinate $\zeta$. After integration, we obtain the
correlation function for correlated rotamer jumps and methyl rotation:

$\displaystyle
C_{i,j}(t)=\frac{1}{5}\sum_{\kappa=1}^{5}\sum_{a,a^{\prime},b=-2}^{2}\sum_{\alpha,\beta}\sum_{n}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\sqrt{p_{\alpha}^{eq}p_{\beta}^{eq}}\tilde{X}_{\alpha}^{(b,n)}\tilde{X}_{\beta}^{(b,n)}e^{\lambda_{b,n}t}\times$
(132)
$\displaystyle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,SF_{\alpha}})\mathcal{D}_{a^{\prime},b}^{(2)}(\Omega_{D,SF_{\beta}})\left[d_{b,0}(\beta_{SF,i})\right]^{2}.$

In order to investigate how the correlation of motions affects relaxation
rates in a 13C1H2H2 methyl group, we consider a simple 2-state exchange, with
Euler angles for jumps of $\\{\pm\pi/2,\beta_{J,CC},0\\}$, where
$\beta_{JCC}=76$ deg is typical for aliphatic carbon chains and thus
corresponds to a $2\beta_{JCC}$-jump. The exchange rate $k_{12}$ is fixed to
$k_{12}=0.5\times 10^{9}$ s-1, and $k_{21}$ is calculated using Eq. 53 to
satisfy the miscroscopic reversibility condition. We impose an isotropic
overall diffusion tensor with global tumbling correlation time $\tau_{c}=10$
ns. Relaxation rates are calculated for three distributions of rotamers:
$p_{1}=0.3$, $p_{1}=0.5$ and $p_{1}=0.7$. Finally, the population average
methyl rotation diffusion coefficient is set to $10^{11}$ s-1, and the
difference of diffusion coefficients between the two states is $\Delta
D_{rot}=5\times 10^{10}$ s-1, with the highest diffusion coefficient being
associated to rotamer 1. We impose an axially symmetric CSA tensor for the 13C
nucleus with an anisotropy of 20 ppm.
Transverse relaxation rates calculated with the same properties for methyl
rotation in both rotamers are indistinguishable from the rates where the
correlation between motions are taken into account (Fig. 3b). However,
longitudinal auto-relaxation rates (Fig. 3a) and DD cross-relaxation rates
(Fig. 3c) show significant deviations depending on whether correlated motions
are considered or not.

###  4.2 Correlated jumps and wobbling in a cone

Here we treat the case where C-C wobbling depends on the rotamer state. The
correlation function can still be written as in Eq. 121 but the angles
$\varphi_{J,SF}$, $\theta_{J,SF}$ and $\phi_{J,SF}$ are all time-dependent,
which complicates the subsequent integration as in Eq. 123. In order to
overcome this difficulty, we decompose the Wigner matrices in Eq. 122 into:

$\displaystyle\langle\mathcal{D}_{b,d}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},d^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle=\sum_{\alpha,\beta=1}^{N}p_{\alpha}^{eq}\times$
(133)
$\displaystyle\sum_{c,c^{\prime}}\langle\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,W},0)\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,W},t)\mathcal{D}_{c,d}^{(2)*}(\Omega_{W,SF},0)\mathcal{D}_{c^{\prime},d^{\prime}}^{(2)}(\Omega_{W,SF},t)\rangle_{\beta\alpha},$

where the Wigner matrices with Euler angles $\Omega_{J,W}$ for transformation
from the jump frame to the wobbling frame can effectively be taken out of the
correlation function $\langle\cdot\cdot\cdot\rangle_{\beta\alpha}$ since they
account for the jump from state $\alpha$ to $\beta$:

$\displaystyle\langle\mathcal{D}_{b,d}^{(2)*}(\Omega_{J,SF},0)\mathcal{D}_{b^{\prime},d^{\prime}}^{(2)}(\Omega_{J,SF},t)\rangle=$
$\displaystyle\sum_{\alpha,\beta=1}^{N}p_{\alpha}^{eq}\sum_{c,c^{\prime}}\mathcal{D}_{b,c}^{(2)*}(\Omega_{J,W_{\alpha}})\mathcal{D}_{b^{\prime},c^{\prime}}^{(2)}(\Omega_{J,W_{\beta}})\times$
(134)
$\displaystyle\langle\mathcal{D}_{c,d}^{(2)*}(\Omega_{W,SF},0)\mathcal{D}_{c^{\prime},d^{\prime}}^{(2)}(\Omega_{W,SF},t)\rangle_{\beta\alpha},$

where $W_{\alpha}$ refers to the wobbling frame in state $\alpha$. We follow
the same approach introduced in the case of correlated jumps and rotation on a
cone to write the conditional probability
$p(\varphi,\theta,\beta,t|\varphi_{0},\theta_{0},\alpha,0)$:

$p(\varphi,\theta,\beta,t|\varphi_{0},\theta_{0},\alpha,0)=\sum_{m=-\infty}^{+\infty}\sum_{n=0}^{+\infty}\frac{e^{im(\varphi-\varphi_{0})}\mathcal{P}_{\nu_{m,n}^{\alpha}}^{m}(\mu_{0})\mathcal{P}_{\nu_{m,n}^{\beta}}^{m}(\mu)}{2\pi\sqrt{H_{m,n}(\mu_{c}^{\alpha})H_{m,n}(\mu_{c}^{\beta})}}c_{m,n}^{\beta\alpha}(t),$
(135)

where the Euler angle for transformation from the Wobbling frame (W) to the
System Frame (SF) is written $\Omega_{W,SF}=\\{\varphi,\theta,-\varphi\\}$,
$\mu=\cos\theta$ and $\mu_{c}^{\alpha}=\cos\beta_{c}^{\alpha}$ with
$\beta_{c}^{\alpha}$ the cone semi-angle opening in state $\alpha$. The index
$0$ refers to angles evaluated at time $t=0$. The time-dependent function
$c_{m,n}^{\beta\alpha}$ is found by solving:

$\frac{d}{dt}C_{m,n}^{\alpha}(t)=(\mathcal{R}-\mathcal{D}_{W}^{(m,n)})C_{m,n}^{\alpha}(t),$
(136)

where $C_{m,n}^{\alpha}(t)$ is a column vector containing the elements
$c_{m,n}^{\beta\alpha}(t)$ for all states $\beta$, $\mathcal{R}$ is the
exchange matrix and $\mathcal{D}_{W}^{(m,n)}$ is a diagonal matrix with
diagonal elements defined as
$D_{W,\beta}\nu_{m,n}^{\beta}(\nu_{m,n}^{\beta}+1)$ for all states $\beta$.
Using the same treatment as shown above, we obtain:

$c_{m,n}^{\beta\alpha}(t)=\sqrt{\frac{p_{\beta}^{eq}}{p_{\alpha}^{eq}}}\sum_{l=1}^{N}\tilde{X}_{\alpha}^{(n,m,l)}\tilde{X}_{\beta}^{(n,m,l)}e^{\lambda_{n,m,l}t},$
(137)

where the sum runs over all states $N$ and $\tilde{X}_{\alpha}^{(m,n,l)}$ is
the $\alpha^{\mathrm{th}}$ value of the eigenvector $\tilde{X}^{(n,m,l)}$
associated to the eigenvalue $\lambda_{n,m,l}$ for the symmetrized pseudo-
exchange matrix
$\tilde{\mathcal{R}}_{W,m,n}=\tilde{\mathcal{R}}-\mathcal{D}_{W}^{(m,n)}$.
Finally, the total correlation function can be calculated, and we obtain,
using the same notations as above:

$\displaystyle
C_{i,j}(t)=\frac{1}{5}\sum_{\kappa}\sum_{a,a^{\prime},b,^{\prime},c,c^{\prime}=-2}^{2}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}\sum_{\alpha,\beta}\sum_{m,n}\frac{\delta_{b-c,b^{\prime}-c^{\prime}}}{1-\mu_{c}^{\alpha}}\sqrt{p_{\alpha}^{eq}p_{\beta}^{eq}}\tilde{X}_{\alpha}^{(n,b-c,m)}\tilde{X}_{\beta}^{(n,b-c,m)}e^{\lambda_{n,b-c,m}t}\times$
(138)
$\displaystyle\frac{I_{b,c}^{n}(\theta_{c}^{\alpha})I_{b^{\prime},c^{\prime}}^{n}(\theta_{c}^{\beta})}{\sqrt{H_{b-c,n}(\mu_{c}^{\alpha})H_{b-c,n}(\mu_{c}^{\beta})}}\mathcal{D}_{a,b}^{*}(\Omega_{J,W_{\alpha}})\mathcal{D}_{a^{\prime},b^{\prime}}(\Omega_{J,W_{\beta}})\mathcal{D}_{c,0}^{*}(\Omega_{SF,i})\mathcal{D}_{c^{\prime},0}(\Omega_{SF,j}).$

Such a model can be computationally demanding and might be difficult to use on
experimental data where the contribution from wobbling motions is expected to
be small.

##  5 Time-dependent interaction strengths

When we introduced the BWR relaxation theory, we implicitly assumed that the
strengths of the interactions were time independent. Chemical bonds vibrate in
the femtosecond range by a few picometers, so that it is a valid approximation
in the case of the DD interactions between directly bonded nuclei, even if its
effective strength needs to be carefully set [83]. Similarly, bonds vibrations
lead to slight variations in CSA tensors for the nitrogen and carbonyl
carbon-13 of peptide planes [84]. S. Tang and D. Case indicate that the use of
a scaling factor for the CSA when analyzing relaxation data is a way to take
into account the motional averaging of the CSA tensor [84]. However, to the
best of our knowledge, no study focused on the variations of CSA tensors in
protein side-chains and their effect on relaxation. CSA tensors have been
measured by solid-state NMR [85] and shown to be conformation dependent,
particularly for alpha and beta carbon-13 nuclei [86, 87]. Rienstra and
coworkers showed that CSA tensors could be used as efficient constraints to
determine the structure of a protein [88, 89].
We restrict ourselves here to the simple case of a methyl group exchanging
between two rotamers in a protein undergoing isotropic overall diffusion with
correlation time $\tau_{c}=10$ ns. We assume perfect tetrahedral geometry for
the methyl group, and a diffusion coefficient for the rotation
$D_{rot}=5\times 10^{10}$ s-1. The set of Euler angles for transformation from
the jump frame to the rotamer frame is
$\Omega_{J,R}=\\{\pm\pi/2,\beta_{J},0\\}$ with $\beta_{J}=76$ deg, which is
typical for carbon chains. The carbon-13 CSA tensors are considered axially
symmetric and aligned along the C-C bond of the methyl group. The CSA of
proton is neglected [72]. We set the exchange rate from state 1 to 2 to
$k_{21}=0.5\times 10^{9}$ s-1 and the exchange rate $k_{21}$ is calculated to
satisfy the microreversibility condition (Eq. 53) depending on the three
situations we consider for the population of state 1: $p_{1}=0.3$, $p_{1}=0.5$
and $p_{1}=0.7$. Finally, the population average of the axially symmetric
carbon-13 CSA is set to $\sigma_{av}=18.2$ ppm which corresponds to the value
of the carbon-$\delta 1$ in isoleucine side-chains determined from cross-
correlated cross-relaxation rates [72]. We evaluate the effect of
$\Delta\sigma=\sigma_{2}-\sigma_{1}$, the difference in CSA between the two
states by focusing on carbon-13 longitudinal and transverse relaxation rates
and the 13C-CSA/DD cross-correlated cross-relaxation rates $\eta_{z}$ and
$\eta_{xy}$ at a magnetic field $B_{0}=14.1$ T [90]:

$\displaystyle R_{1}(^{13}\mathrm{C})=$
$\displaystyle\frac{1}{4}(\mathcal{J}_{CH}^{(I)}(\omega_{H}-\omega_{C})+3\mathcal{J}^{(I)}_{CH}(\omega_{C})+6\mathcal{J}^{(I)}_{CH}(\omega_{H}+\omega_{C}))$
(139)
$\displaystyle+\frac{4}{3}(\mathcal{J}^{(I)}_{CD}(\omega_{C}-\omega_{D})+3\mathcal{J}^{(I)}_{CD}(\omega_{C})+6\mathcal{J}^{(I)}_{CD}(\omega_{C}+\omega_{D}))+\mathcal{J}^{(I)}_{C}(\omega_{C}),$
$\displaystyle R_{2}(^{13}\mathrm{C})=$
$\displaystyle\frac{1}{8}(4\mathcal{J}^{(I)}_{CH}(0)+\mathcal{J}^{(I)}_{CH}(\omega_{H}-\omega_{C})+3\mathcal{J}^{(I)}_{CH}(\omega_{C})+6\mathcal{J}^{(I)}_{CH}(\omega_{H})+6\mathcal{J}^{(I)}_{CH}(\omega_{H}+\omega_{C}))$
$\displaystyle+\frac{2}{3}(4\mathcal{J}^{(I)}_{CD}(0)+\mathcal{J}^{(I)}_{CD}(\omega_{C}-\omega_{D})+3\mathcal{J}^{(I)}_{CD}(\omega_{C})+6\mathcal{J}^{(I)}_{CD}(\omega_{D})+6\mathcal{J}^{(I)}_{CD}(\omega_{C}+\omega_{D}))$
$\displaystyle+\frac{1}{6}(4\mathcal{J}^{(I)}_{C}(0)+3\mathcal{J}^{(I)}_{C}(\omega_{C})),$
$\displaystyle\eta_{z}(^{13}\mathrm{C})=$
$\displaystyle\sqrt{\frac{3}{2}}\mathcal{J}^{(I)}_{C,CH}(\omega_{C}),$
$\displaystyle\eta_{xy}(^{13}\mathrm{C})=$
$\displaystyle\frac{1}{2}\sqrt{\frac{1}{6}}(4\mathcal{J}^{(I)}_{C,CH}(0)+3\mathcal{J}^{(I)}_{C,CH}(\omega_{C})),$

where the superscript $(I)$ indicates that the strengths of the interactions
are included in the spectral density functions which are written:

$\displaystyle\mathcal{J}^{(I)}_{CH}(\omega)=d_{CH}^{2}\frac{2}{5}\sum_{a,b=-2}^{2}\sum_{\alpha,\beta=1}^{2}\sum_{n=0}^{1}\frac{\tau_{b,n}}{1+(\omega\tau_{b,n})^{2}}\sqrt{p_{\alpha}p_{\beta}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{D}_{a,b}^{(2)*}(\Omega_{J,R_{\alpha}})\mathcal{D}_{a,b}^{(2)}(\Omega_{J,R_{\beta}})d_{b,0}(\beta_{CCH})^{2},$
(140)
$\displaystyle\mathcal{J}^{(I)}_{CD}(\omega)=d_{CD}^{2}\frac{2}{5}\sum_{a,b=-2}^{2}\sum_{\alpha,\beta=1}^{2}\sum_{n=0}^{1}\frac{\tau_{b,n}}{1+(\omega\tau_{b,n})^{2}}\sqrt{p_{\alpha}p_{\beta}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{D}_{a,b}^{(2)*}(\Omega_{J,R_{\alpha}})\mathcal{D}_{a,b}^{(2)}(\Omega_{J,R_{\beta}})d_{b,0}(\beta_{CCH})^{2},$
$\displaystyle\mathcal{J}^{(I)}_{C}(\omega)=\frac{2}{3}\omega_{C}^{2}\frac{2}{5}\sum_{\alpha,\beta=1}^{2}\sum_{n=0}^{1}\sigma_{\alpha}\sigma_{\beta}\frac{\tau_{0,n}}{1+(\omega\tau_{0,n})^{2}}\sqrt{p_{\alpha}p_{\beta}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{P}_{2}(\cos\theta_{\alpha,\beta}),$
$\displaystyle\mathcal{J}^{(I)}_{C,CH}(\omega)=d_{CH}\sqrt{\frac{2}{3}}\omega_{C}\frac{2}{5}\mathcal{P}_{2}(\cos\beta_{CCH})\sum_{\alpha,\beta=1}^{2}\sum_{n=0}^{1}\sigma_{\alpha}\frac{\tau_{0,n}}{1+(\omega\tau_{0,n})^{2}}\sqrt{p_{\alpha}p_{\beta}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{P}_{2}(\cos\theta_{\alpha,\beta}),$

where $d_{CX}=-\mu_{0}\hbar\gamma_{C}\gamma_{X}/(4\pi r_{CX}),X={H,D}$ with
$\mu_{0}$ the permeability of free space, $\hbar$ the Planck’s constant
devided by $2\pi$, $\gamma_{A}$ the gyromagnetic ratio of nucleus $A$ and
$r_{CX}$ the distance between the 13C and nucleus X,
$\beta_{CCH}=180-109.47=70.53$ deg is the angle between the C-H bond and the
methyl group symmetry axis, $\theta_{\alpha,\beta}$ is the angle between the
vectors pointing along the directions of the C-C bonds in rotamers $\alpha$
and $\beta$ (that is $\theta_{\alpha,\beta}=0$ when $\alpha=\beta$ and
$\theta_{\alpha\beta}=2\beta_{CC}$ when $\alpha\neq\beta$), and:

$\tau_{b,n}^{-1}=\tau_{c}^{-1}-\lambda_{n}+b^{2}D_{rot},$ (141)

where $\lambda_{n}$ is the $n^{th}$ eigenvalue of the symmetrized exchange
matrix.

Figure 4: Evolution of the carbon-13 longitudinal R1 (a) and transverse R2 (b)
relaxation rates, as well as carbon longitudinal (c) and transverse (d) cross-
correlated cross-relaxation rate for a methyl group exchanging between two
rotamer positions as a function of the difference between the CSA of the two
rotamers. Calculations are performed for three equilibrium position for the
state 1 and by either considering a population-averaged CSA value (dash) or
distinct CSA tensors (solid).

The carbon auto-relaxation rates show very small deviations depending on
whether distinct CSA values are considered or not (Fig. 4a,b), which most
likely arises from the relatively small contribution of the CSA (about 1.5 %
for the carbon-R1 and 5 % for the carbon-R2) compared to the DD interactions.
The cross-correlated cross-relaxation rates calculated with distinct CSA
tensors show a linear variation with the difference in anisotropy between the
two states (Fig. 4c,d). This can be understood by expanding the spectral
density function for cross-correlation between the carbon-CSA and carbon-
proton DD interactions:

$\displaystyle\mathcal{J}^{(I)}_{C,CH}(\omega)=$ $\displaystyle
d_{CH}\sqrt{\frac{2}{3}}\omega_{C}\frac{2}{5}\mathcal{P}_{2}(\cos\beta_{CCH})\left[\sigma_{av}\left(\frac{\tau_{c}}{1+(\omega\tau_{c})^{2}}(p_{1}^{2}+p_{2}^{2}+2p_{1}p_{2}\mathcal{P}_{2}(\cos\beta_{J}))\right.\right.$
(142)
$\displaystyle\left.\left.+\frac{2\tau_{0,1}}{1+(\omega\tau_{0,1})^{2}}p_{1}p_{2}(1-\mathcal{P}_{2}(\cos\beta_{J}))\right)\right.$
$\displaystyle\left.+\Delta\sigma
p_{1}p_{2}(p_{1}-p_{2})(1-\mathcal{P}_{2}(\cos\beta_{J}))\left(\frac{\tau_{0,1}}{1+(\omega\tau_{0,1})^{2}}-\frac{\tau_{c}}{1+(\omega\tau_{c})^{2}}\right)\right],$

where the last line contains the $\Delta\sigma$-dependent part of the spectral
density function and highlights the linear variation of the cross-correlated
cross-relaxation rates when the difference of CSA is changed. For the chosen
geometry, we have:

$\displaystyle d_{CH}\mathcal{P}_{2}(\cos\beta_{CCH})>0,$ (143) $\displaystyle
1-\mathcal{P}_{2}(\cos\beta_{J})>0.$

In addition, we have chosen $\Delta\sigma=\sigma_{2}-\sigma_{1}>0$ and:

$\displaystyle\tau_{0,1}-\tau_{c}=-\tau_{c}^{2}\frac{k_{12}+k_{21}}{1+\tau_{c}(k_{12}+k_{21})}<0,$
(144)
$\displaystyle\frac{\tau_{0,1}}{1+(\omega_{C}\tau_{0,1})^{2}}-\frac{\tau_{c}}{1+(\omega_{C}\tau_{c})^{2}}>0,$

such that, when calculations are performed with distinct CSA tensors,
$\eta_{z}$ increases when $p_{1}>p_{2}$ and decreases otherwise, and since
$J^{(I)}_{C,CH}(0)\gg J^{(I)}_{C,CH}(\omega_{C})$, $\eta_{xy}$ decreases when
$p_{1}>p_{2}$ and increases otherwise. Finally, when $p_{1}=p_{2}$, the
$\Delta\sigma$-dependent term in Eq. 142 vanishes and the cross-correlated
cross-relaxation rates are independent from the difference in CSA value
between the two states.
Here, we used a simple model to highlight the contribution of internal
dynamics to relaxation when a spin system exchanges between discrete positions
with distinct CSA tensors. Calculations show that small effects on relaxation
rates can be expected from changes of the amplitude of the CSA of less than 10
ppm between conformers (Fig. 4). Investigations combining solid-state NMR,
relaxation and density functional theory calculations of CSA parameters may
offer experimental evidence of the contribution of such mechanisms to
relaxation.

##  6 Discussion

The models presented in the previous sections all share one property: they are
complex. Up to 3 correlation times might be necessary to model the global
tumbling; the number of decaying exponentials describing rotamer jumps
increases linearly with the number of accessible states; modeling the wobbling
in a cone motion involves evaluating difficult integrals. Out of the 4
considered motions, the rotation on a cone seems the most simple one. In
addition to their intrinsic complexity, another difficulty arises when
analyzing data, in particular NMR relaxation rates: how to choose one model
over the other without a priori knowledge of the type of motions involved?
These aspects naturally led to the popularization of the Model-Free
correlation function [13] to model internal motions:

$C_{MF}^{(i)}(t)=\mathcal{S}_{MF}^{2}+(1-\mathcal{S}_{MF}^{2})e^{-t/\tau_{MF}},$
(145)

where $\mathcal{S}_{MF}^{2}$ is the squared order parameter and $\tau_{MF}$ an
effective correlation time. The simplicity of this model makes it particularly
attractive, and has indeed led to a vast number of successful analyses of NMR
relaxation experiments [15, 91, 50, 27, 18, 92, 93, 94, 95, 96]. Over the past
few decades, the MF approach has been modified in order to include the effect
of cross correlation [97], or to include more than one correlation time and
study more complex systems [14, 98, 99, 100], in particular intrinsically
disordered proteins [101, 102, 103, 104, 105]. However, a number of questions
arises on this simple form of correlation function, and in particular
regarding the factorization of the correlation functions accounting for
internal motions (Eq. 145) and global tumbling. It is clear from section 3.2
that such factorization is mathematically exact only in the case of an
isotropic diffusion tensor, which is the frame in which the MF approach was
initially presented. An approximated form for the global tumbling correlation
function in the case of an axially symmetric diffusion tensor was proposed as
well, and consists in the sum of two decaying exponential, with decay
constants and relative weights that can be determined experimentally [13, 106,
107, 108]. Then, how is the factorization of the global tumbling and internal
motion correlation functions affecting the fitted values of
$\mathcal{S}_{MF}^{2}$ and $\tau_{MF}$? It must be noted that the SRLS model
of correlation function was introduced in part as an attempt to overcome this
difficulty [30, 28, 32]. Here, we evaluate the effect of the factorization of
the global tumbling correlation function on the value of the fitted MF
parameters. We focus on auto-correlation only:

$C(t)=\sum_{a,a^{\prime},b,b^{\prime}=-2}^{2}\langle\mathcal{D}_{q,a}^{(2)*}(\Omega_{L,D},0)\mathcal{D}_{q,a^{\prime}}^{(2)}(\Omega_{L,D},t)\rangle\mathcal{D}_{a,b}^{(2)*}(\Omega_{D,M})\mathcal{D}_{a^{\prime}b^{\prime}}^{(2)}(\Omega_{D,M}\langle\mathcal{D}_{b,0}^{(2)*}(\Omega_{M,i},0)\mathcal{D}_{b^{\prime},0}^{(2)}(\Omega_{M,i},t)\rangle,$
(146)

where $M$ denotes the frame in which the motion is best described,
$\langle\mathcal{D}_{q,a}^{(2)*}(\Omega_{L,D},0)\mathcal{D}_{q,a^{\prime}}^{(2)}(\Omega_{L,D},t)\rangle$
accounts for overall rotational diffusion, $\Omega_{D,SF}$ is the angle
orienting the diffusing frame in the global diffusion tensor frame, and
$\langle\mathcal{D}_{b,0}^{(2)*}(\Omega_{M,i},0)\mathcal{D}_{b^{\prime},0}^{(2)}(\Omega_{M,i},t)\rangle$
accounts for internal motions. The MF correlation function is written as:

$C_{MF}(t)=\sum_{a,a^{\prime}=-2}^{2}\langle\mathcal{D}_{q,a}^{(2)*}(\Omega_{L,D},0)\mathcal{D}_{q,a^{\prime}}^{(2)}(\Omega_{L,D},t)\rangle\mathcal{D}_{a,0}^{(2)*}(\Omega_{D,M})\mathcal{D}_{a^{\prime}0}^{(2)}(\Omega_{D,M})\left(\mathcal{S}^{2}_{MF}+(1-\mathcal{S}^{2}_{MF})e^{-t/\tau_{MF}}\right).$
(147)

In the case of an isotropic diffusion tensor, the factorization of the global
tumbling correlation function is exact and the correlation functions are
written:

$\displaystyle
C(t)=\frac{1}{5}e^{-t/\tau_{c}}\sum_{b=-2}^{2}\langle\mathcal{D}_{b,0}^{(2)*}(\Omega_{M,i},0)\mathcal{D}_{b,0}^{(2)}(\Omega_{M,i},t)\rangle,$
(148) $\displaystyle
C_{MF}(t)=\frac{1}{5}e^{-t/\tau_{c}}\left(\mathcal{S}^{2}_{MF}+(1-\mathcal{S}^{2}_{MF})e^{-t/\tau_{MF}}\right).$

Table 3: Values of the fitted MF parameters (correlation function Eq. 147) for 4 types of internal motions and the three possible symmetry properties of the overal diffusion tensor. In the case of the symmetric top and fully asymmetric overal diffusion tensor, we report the average and standard deviation over the possible orientations, i.e. the different values of $\theta_{D,M}$ and $\\{\varphi_{D,M},\theta_{D,M}\\}$ respectively. In these fits, the values and orientation of the diffusion tensor (angle $\Omega_{D,M}$) were fixed such that $\mathcal{S}_{MF}^{2}$ and $\tau_{MF}$ were the only adjustable parameters. The diffusion tensor for global tumbling has $D_{xx}=1.64\times 10^{7}$ s-1, $D_{yy}=1.82\times 10^{7}$ s-1 and $D_{zz}=2.35\times 10^{7}$ s-1 which corresponds to the diffusion tensor for the proximal Ubiquitin in a diubiquitin with linkage at lysine-11 [109]. When the symmetric model is used, eigenvalues are $D_{\parallel}=D_{zz}$ and $D_{\perp}=\frac{1}{2}(D_{xx}+D_{yy})$. For the isotropic diffusion, the eigenvalue equals $D=\frac{1}{3}(D_{xx}+D_{yy}+D_{zz})$. The orientation of the rotamer frames in the jump frame are $\Omega_{J,R}^{(2)}=\\{\pm\frac{\pi}{2},\beta_{CC}\\}$ when two states are considered, and $\Omega_{J,R}^{(3)}=\\{2n\frac{\pi}{3},\beta_{CC}\\}$ for the 3-state jump with $n=0,1,2$ and where $\beta_{CC}=76^{\circ}$. The 3-state jump model with unequal populations (bottom row) was simulated with $p_{1}=0.7$, $p_{2}=0.2$ and $p_{3}=0.1$. The rotation on a cone is modeled assuming a methyl group gegometry and with $D_{rot}=10^{11}$ s-1. For the wobbling motion, the cone semi-angle opening is $\beta_{c}=20^{\circ}$ and $D_{W}=10^{8}$ s-1. | Isotropic | Symmetric top | Fully asymmetric
---|---|---|---
| $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps) | $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps) | $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps)
rotation | $0.12$ | $4.8$ | $0.12\pm 0.00$ | $4.8\pm 0.0$ | $0.12\pm 0.00$ | $4.8\pm 0.0$
wobbling | $0.83$ | $339$ | $0.83\pm 0.00$ | $344\pm 5$ | $0.83\pm 0.00$ | $344\pm 5$
2-state jump | $0.83$ | $471$ | $0.77\pm 0.06$ | $1,290\pm 709$ | $0.77\pm 0.06$ | $1,319\pm 748$
$p_{1}=p_{2}$
3-state jump | $0.17$ | $336$ | $0.17\pm 0.00$ | $333\pm 3.0$ | $0.17\pm 0.00$ | $333\pm 3$
$p_{1}=p_{2}=p_{3}$
3-state jump | $0.43$ | $157$ | $0.43\pm 0.03$ | $161\pm 55$ | $0.43\pm 0.03$ | $164\pm 57$
$p_{1}\neq p_{2}\neq p_{3}$

We considered the diffusion tensor of diubiquitin with linkage at lysine-11
[109], which has anisotropy of 1.36 and rhombicity 0.44. The fitted parameters
for internal dynamics (Table 3) can be used as reference for the parameters
fitted in the case of an asymmetric overall diffusion tensor: if the
factorization does not affect the value of the fitted parameters, we would
expect to find the same values for different tensor symmetries.
For asymmetric diffusion tensor, the factorization of the global tumbling
correlation function is not mathematically exact. The two correlation
functions for the symmetrical top diffusion tensor read [13, 110]:

$\displaystyle
C(t)=\frac{1}{5}\sum_{a,b,b^{\prime}=-2}^{2}e^{-(6D_{\perp}+a^{2}(D_{\parallel}-D_{\perp}))t}d_{a,b}(\theta_{D,M})d_{a,b^{\prime}}(\theta_{D,M})\langle\mathcal{D}_{b,0}^{(2)*}(\Omega_{M,i},0)\mathcal{D}_{b,0}^{(2)}(\Omega_{M,i},t)\rangle,$
(149) $\displaystyle
C_{MF}(t)=\frac{1}{5}\left(\mathcal{S}^{2}_{MF}+(1-\mathcal{S}^{2}_{MF})e^{-t/\tau_{MF}}\right)\sum_{a=-2}^{2}e^{-(6D_{\perp}+a^{2}(D_{\parallel}-D_{\perp}))t}\left[d_{a,0}(\theta_{D,M})\right]^{2},$

and for the fully asymmetric tensor [27]:

$\displaystyle
C(t)=\frac{1}{5}\sum_{a,a^{\prime},b,b^{\prime}=-2}^{2}\sum_{\kappa=1}^{5}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}e^{i\varphi_{D,M}(a-a^{\prime})}d_{a,b}(\theta_{D,M})d_{a,b^{\prime}}(\theta_{D,M})\langle\mathcal{D}_{b,0}^{(2)*}(\Omega_{M,i},0)\mathcal{D}_{b,0}^{(2)}(\Omega_{M,i},t)\rangle,$
(150) $\displaystyle
C_{MF}(t)=\frac{1}{5}\left(\mathcal{S}^{2}_{MF}+(1-\mathcal{S}^{2}_{MF})e^{-t/\tau_{MF}}\right)\sum_{a,a^{\prime}=-2}^{2}\sum_{\kappa=1}^{5}a_{\kappa,a}a_{\kappa,a^{\prime}}e^{-E_{\kappa}t}e^{i\varphi_{D,M}(a-a^{\prime})}\left[d_{a,0}(\theta_{D,M})\right]^{2}.$

We can note here that the MF correlation function cannot distinguish the
orientation of the interaction frames in the diffusion frame, which is
particularly critical in the presence of a rotamer jump motion. The fitted MF
parameters agree well with the ones obtained when simulating an isotropic
global tumbling, except in the case of the two-state jump model where large
deviations of the fitted parameters are obtained when $\Omega_{D,M}$ is
changed (Table 3). Two- and three-state jump models have been studied
extensively by solid-state deuterium NMR, based on the investigation of
linewidths and relaxation rates [111, 112, 51, 81, 52]. J. Wand and co-workers
already distinguished the 2-state and 3-state jump models by considering the
resulting symmetry [76]. In our case of 3 exchanging states with equal
populations, the motion is azymutal symmetric, that is:

$\langle x^{2}\rangle=\langle y^{2}\rangle,\,\,\,\langle x\rangle=\langle
y\rangle=0,\,\,\,\mathrm{and}\,\,\,\langle xy\rangle=0,$ (151)

where $\langle\cdot\cdot\cdot\rangle$ denotes an ensemble average, and $x$ and
$y$ are the $x$\- and $y-$coordinates of the interaction frame main axis in
the jump frame (labelled $M$ in this section). For the 2-state jump model
presented here, the first condition in Eq. 151 is not fulfilled. It is not
fulfilled neither when populations in the 3-state jump model are not equal, in
which case the factorization of the asymmetric global tumbling correlation
functions does not accurately reproduces MF parameters obtained in the
isotropic tumbling case, and show large deviations depending on the
orientation in the diffusion tensor frame (Table 3). These conclusions also
apply for higher asymmetry of the diffusion tensor (Table S1).
These simulations suggest that the factorization of the global tumbling
correlation function does not alter the value of the fitted parameters in the
presence of a diffusing-type of internal motions (that is diffusion on a cone
or wobbling in a cone). When the sampling of the conformational space is
discrete (jump model), the factorization does not affect the value of the MF
fitted parameters only when the motion has azymutal symmetry, a strong
condition that may not be met in most cases when studying protein side-chains
dynamics. Overall, explicit models of motions are particularly powerful to
account for the orientation of interaction frames in the molecular frame. This
can prove itself especially adapted in the case of non-axially symmetric CSA
tensors. Their contribution to relaxation can be decomposed in contribution
from two orthogonal axially symmetric tensors which orientation in the
molecular frame are different and can easily be included in the explicit
models of motions.

##  7 Conclusion

We have reviewed explicit models of motions that can be relevant to study the
dynamic properties of biomolecules. For the past 40 years, the Model Free
approach has been prefered for the analysis of NMR relaxation data as it does
not require a priori knowledge on the nature of motions of the bond vectors.
The advances in MD simulations methods can now provide such information, for
instance with the help of the detector analysis which can easily distinguish
the contribution of motions originating from different times-scales [73, 74,
25, 75, 113]. Explicit models of motions as presented here can be used to
obtain a mechanistic picture of the motions from a combined NMR and MD
analysis [82, 41]. We finally discussed the use of MF when analyzing NMR
relaxation data when the overall diffusion tensor is anisotropic. In this
case, the MF fitted parameter are affected by the orientation of the
interaction tensor frames in the overall diffusion tensor frame, in particular
when the motion is not azymutal symmetric. We expect this theoretical work, in
pair with softwares for the calculation of the relaxation rates of nuclear
spins [114, 115, 90], to support the development of system-specific explicit
models of motions for the analysis of NMR relaxation.

## References

  * [1] H. Frauenfelder, S. G. Sligar, and P. G. Wolynes Science, vol. 254, p. 1598, 1991.
  * [2] C. Charlier, S. F. Cousin, and F. Ferrage Chem. Soc. Rev., vol. 45, p. 2410, 2016.
  * [3] T. R. Alderson and L. E. Kay Cell, vol. 184, p. 577, 2021.
  * [4] M. S. Choy, Y. Li, L. E. S. F. Machado, M. B. A. Kunze, C. R. Connors, X. Wei, K. Lindorff-Larsen, R. Page, and W. Peti Mol. Cell, vol. 65, p. 644, 2017\.
  * [5] R. Otten, R. A. P. Pádua, H. A. Bunzel, V. Nguyen, W. Pitsawong, M. Patterson, S. Sui, S. L. Perry, A. E. Cohen, D. Hilvert, and D. Kern Science, vol. 370, p. 1442, 2020.
  * [6] K. H. DuBay, G. R. Bowman, and P. L. Geissler Acc. Chem. Res., vol. 48, p. 1098, 2015.
  * [7] V. Kasinath, Y. Fu, K. A. Sharp, and A. J. Wand Angew. Chem. Int. Ed., vol. 54, p. 102, 2015.
  * [8] A. G. Palmer Chem. Rev., vol. 104, p. 3623, 2004.
  * [9] M. P. Nicholas, E. Eryilmaz, F. Ferrage, D. Cowburn, and R. Ghose Prog. Nucl. Magn. Reson. Spectrosc., vol. 57, p. 111, 2010.
  * [10] R. Wangsness and F. Bloch Phys. Rev., vol. 89, p. 728, 1953.
  * [11] A. G. Redfield IBM J. Res. Develop., vol. 1, p. 19, 1957.
  * [12] A. G. Redfield Advances in Magnetic and Optical Resonance, vol. 1, p. 1, 1965\.
  * [13] G. Lipari and A. Szabo J. Am. Chem. Soc., vol. 104, p. 4546, 1982.
  * [14] M. G. Clore, A. Szabo, A. Bax, L. E. Kay, P. C. Driscoll, and A. M. Gronenborn J. Am. Chem. Soc., vol. 112, p. 4989, 1990.
  * [15] G. Lipari and A. Szabo J. Am. Chem. Soc., vol. 104, p. 4559, 1982.
  * [16] L. K. Nicholson, L. E. Kay, D. M. Baldisseri, J. Arango, P. E. Young, A. Bax, and D. A. Torchia Biochemistry, vol. 31, p. 5253, 1992.
  * [17] N. A. Farrow, R. Muhandiram, A. U. Singer, S. M. Pascal, C. M. Kay, G. Gish, S. E. Schoelson, T. Pawson, J. D. Forman-Kay, and L. E. Kay Biochemistry, vol. 33, p. 5984, 1994.
  * [18] A. M. Mandel, M. Akke, and A. G. Palmer J. Mol. Biol, vol. 246, p. 144, 1995\.
  * [19] D. Fushman, N. Tjandra, and D. Cowburn J. Am. Chem. Soc., vol. 121, p. 8577, 1999.
  * [20] V. A. Jarymowycz and M. J. Stone Chem. Rev., vol. 106, p. 1624, 2006.
  * [21] K. Lindorff-Larsen, R. B. Best, M. A. DePristo, C. M. Dobson, and M. Vendruscolo Nature, vol. 433, p. 128, 2005.
  * [22] P. Robustelli, K. A. Stafford, and A. G. Palmer J. Am. Chem. Soc., vol. 134, p. 6365, 2012.
  * [23] N. Salvi, A. Abyzov, and M. Blackledge J. Phys. Chem. Lett., vol. 7, p. 2483, 2016.
  * [24] S. F. Cousin, P. Kadeřávek, N. Bolik-Coulon, Y. Gu, C. Charlier, L. Carlier, L. Bruschweiler-Li, T. Marquardsen, J.-M. Tyburn, R. Brüschweiler, and F. Ferrage J. Am. Chem. Soc., vol. 140, p. 13456, 2018.
  * [25] A. A. Smith, M. Ernst, S. Riniker, and B. H. Meier Angew. Chem. Int. Ed, vol. 58, p. 9383, 2019.
  * [26] F. Kümmerer, S. Orioli, D. Larsen-Harding, F. Hoffmann, Y. Gavrilov, K. Teilum, and K. Lindorff-Larsen J. Chem. Theory Comput., vol. 17, p. 5262, 2021.
  * [27] N. Tjandra, S. E. Feller, R. W. Pastor, and A. Bax J. Am. Chem. Soc., vol. 117, p. 12562, 1995.
  * [28] V. Tugarinov, Z. Liang, Y. E. Shapiro, J. H. Freed, and E. Meirovitch J. Am. Chem. Soc., vol. 123, p. 3055, 2001.
  * [29] E. Meirovitch, Y. E. Shapiro, V. Tugarinov, Z. Liang, and J. H. Freed J. Phys. Chem. B, vol. 107, p. 9883, 2003.
  * [30] A. Polimeno and J. H. Freed J. Phys. Chem., vol. 99, p. 10995, 1995.
  * [31] Z. Liang and J. H. Freed J. Phys. Chem. B, vol. 103, p. 6384, 1999.
  * [32] E. Meirovitch, Y. E. Shapiro, A. Polimeno, and J. H. Freed J. Phys. Chem. A, vol. 110, p. 8366, 2006.
  * [33] E. Meirovitch, A. Polimeno, and J. H. Freed J. Phys. Chem. B, vol. 110, p. 20615, 2006.
  * [34] E. Meirovitch, Y. E. Shapiro, A. Polimeno, and J. H. Freed J. Phys. Chem. B, vol. 111, p. 12865, 2007.
  * [35] D. Yang and L. E. Kay J. Mol. Biol., vol. 263, p. 369, 1996.
  * [36] D.-W. Li and R. Brüschweiler J. Am. Chem. Soc., vol. 131, p. 7226, 2009.
  * [37] F. Hoffmann, F. A. A. Mulder, and L. V. Schäfer J. Phys. Chem. B, vol. 126, p. 54, 2022.
  * [38] D. Wallach J. Chem. Phys., vol. 47, p. 5258, 1967.
  * [39] R. J. Wittebort and A. Szabo J. Chem. Phys., vol. 69, p. 1722, 1978.
  * [40] C. C. Wang and R. Pecora J. Chem. Phys., vol. 72, p. 5333, 1980.
  * [41] N. Bolik-Coulon, O. Languin-Cattoën, D. Carnevale, M. Zachrdla, D. Laage, F. Sterpone, G. Stirnemann, and F. Ferrage arXiv:2204.05813, 2022.
  * [42] M. Goldman J. Magn. Reson., vol. 60, p. 437, 1984.
  * [43] J. Cavanagh, W. J. Fairbrother, A. G. Palmer, M. Rance, and N. J. Skelton, Protein NMR Spectroscopy: Principles and Practice.  Elsevier Academic Press, 2007.
  * [44] P. Lunginbühl and K. Wüthrich Prog. Nucl. Magn. Reson. Spectrosc., vol. 40, p. 199, 2002.
  * [45] L. G. Werbelow and D. M. Grant Advances in Magnetic and Optical Resonance, vol. 9, p. 189, 1977.
  * [46] D. L. Favro Phys. Rev., vol. 119, p. 53, 1960.
  * [47] D. E. Woessner J. Chem. Phys., vol. 37, p. 647, 1962.
  * [48] W. A. Steele J. Chem. Phys., vol. 38, p. 2411, 1963.
  * [49] W. T. Huntress J. Chem. Phys., vol. 48, p. 3524, 1968.
  * [50] R. Brüschweiler, X. Liao, and P. E. Wright Science, vol. 268, p. 886, 1995.
  * [51] L. S. Batchelder, C. H. Niu, and D. A. Torchia J. Am. Chem. Soc., vol. 105, p. 2228, 1983.
  * [52] L. Vugmeyster and D. Ostrovsky Methods, vol. 148, p. 136, 2018.
  * [53] R. B. Best, J. Clarke, and M. Karplus J. Am. Chem. Soc., vol. 126, p. 7734, 2004.
  * [54] J. J. Chou, D. A. Case, and A. Bax J. Am. Chem. Soc., vol. 125, p. 8959, 2003.
  * [55] V. Tugarinov, P. M. Hwang, J. E. Ollerenshaw, and L. E. Kay J. Am. Chem. Soc., vol. 125, p. 10420, 2003.
  * [56] V. Tugarinov, V. Kanelis, and L. E. Kay Nat. Protoc., vol. 1, p. 749, 2006\.
  * [57] P. Gans, O. Hamelin, R. Sounier, I. Ayala, M. A. Durá, C. D. Amero, M. Noirclerc-Savoye, B. Franzetti, M. J. Plevin, and J. Boisbouvier Angew. Chem. Int. Ed, vol. 49, p. 1958, 2010.
  * [58] J. E. Ollerenshaw, V. Tugarinov, and L. E. Kay Magn. Reson. Chem., vol. 41, p. 843, 2003.
  * [59] N. Bolik-Coulon, S. F. Cousin, P. Kadeřávek, J.-N. Dumez, and F. Ferrage J. Chem. Phys., vol. 150, p. 224202, 2019.
  * [60] F. Hoffmann, F. A. A. Mulder, and L. V. Schäfer J. Phys. Chem. B, vol. 122, p. 5038, 2018.
  * [61] F. Hoffmann, F. A. A. Mulder, and L. V. Schäfer J. Chem. Phys., vol. 152, p. 084102, 2020.
  * [62] P. Pelupessy, F. Ferrage, and G. Bodenhausen J. Chem. Phys., vol. 126, p. 134, 2007.
  * [63] R. Brüschweiler and P. E. Wright J. Am. Chem. Soc., vol. 116, p. 8426, 1994.
  * [64] T. Bremi, R. Brüschweiler, and R. R. Ernst J. Am. Chem. Soc., vol. 119, p. 4272, 1997.
  * [65] T. Bremi and R. Brüschweiler J. Am. Chem. Soc., vol. 119, p. 6672, 1997.
  * [66] K. Kinosita, S. Kawato, and A. Ikegami Biophys. J., vol. 20, p. 289, 1977\.
  * [67] M. P. Warchol and W. E. Vaughan Advances in Molecular Relaxation and Interaction Processes, vol. 13, p. 317, 1978.
  * [68] G. Lipari and A. Szabo Biophys. J., vol. 30, p. 489, 1980.
  * [69] G. Lipari and A. Szabo J. Chem. Phys., vol. 75, p. 2971, 1981.
  * [70] A. Kumar and G. C. Levy J. Chem. Phys., vol. 85, p. 485, 1986.
  * [71] A. Kumar J. Chem. Phys., vol. 91, p. 1232, 1989.
  * [72] V. Tugarinov, C. Scheurer, R. Brüschweiler, and L. E. Kay J. Biomol. NMR, vol. 30, p. 397, 2004.
  * [73] A. A. Smith, E. Matthias, and B. H. Meier J. Chem. Phys., vol. 148, p. 045104, 2018.
  * [74] A. A. Smith, E. Matthias, B. H. Meier, and F. Ferrage J. Chem. Phys., vol. 151, p. 034102, 2019.
  * [75] A. A. Smith, N. Bolik-Coulon, E. Matthias, B. H. Meier, and F. Ferrage J. Biomol. NMR, vol. 75, p. 119, 2021.
  * [76] K. K. Frederick, K. A. Sharp, N. Warischalk, and J. A. Wand J. Chem. Phys. B, vol. 112, p. 12095, 2008.
  * [77] V. Wong, D. A. Case, and A. Szabo Proc. Natl. Acad. Sci. USA, vol. 106, p. 11016, 2009.
  * [78] Y. Ryabov, M. G. Clore, and C. D. Schwieters J. Chem. Phys., vol. 136, p. 034108, 2012.
  * [79] M. L. Gill and A. G. Palmer J. Phys. Chem. B, vol. 118, p. 11120, 2014\.
  * [80] Y. Ryabov Proteins, vol. 83, p. 1571, 2015.
  * [81] S. K. Straus, T. Bremi, and R. R. Ernst J. Biomol. NMR, vol. 10, p. 119, 1997.
  * [82] A. A. Smith J. Magn. Reson. Open, vol. 10-11, p. 100045, 2022.
  * [83] D. A. Case J. Biomol. NMR, vol. 15, p. 95, 1999.
  * [84] S. Tang and D. A. Case J. Biomol. NMR, vol. 38, p. 255, 2007.
  * [85] L. Liang, Y. Ji, K. Chen, P. Gao, Z. Zhao, and G. Hou Chem. Rev., vol. 122, p. 9880, 2022.
  * [86] M. Hong J. Am. Chem. Soc., vol. 122, p. 3762, 2000.
  * [87] B. J. Wylie, W. T. Franks, D. T. Graesser, and C. M. Rienstra J. Am. Chem. Soc., vol. 127, p. 11946, 2005.
  * [88] B. J. Wylie, C. D. Schwieters, E. Oldfield, and C. M. Rienstra J. Am. Chem. Soc., vol. 131, p. 985, 2009.
  * [89] B. J. Wylie, L. J. Sperling, A. J. Nieuwkoop, W. T. Franks, E. Oldfield, and C. M. Rienstra Proc. Natl. Acad. Sci. U.S.A., vol. 108, p. 16974, 2011\.
  * [90] N. Bolik-Coulon, P. Kadeřávek, P. Pelupessy, J.-N. Dumez, F. Ferrage, and S. F. Cousin J. Magn. Reson., vol. 313, p. 106718, 2020.
  * [91] D. M. Schneider, M. J. Dellwo, and A. J. Wand Biochemistry, vol. 31, p. 3645, 1992.
  * [92] N. A. Farrow, O. Zhang, J. D. Forman-Kay, and L. E. Kay Biochemistry, vol. 34, p. 868, 1995.
  * [93] L. K. Lee, M. Rance, W. J. Chazin, and A. G. Palmer J. Biomol. NMR., vol. 9, p. 287, 1997.
  * [94] J. B. Hall and D. Fushman J. Am. Chem. Soc., vol. 128, p. 7855, 2006.
  * [95] V. Z. Miloushev, F. Bahna, C. Ciatto, G. Ahlsen, B. Honig, L. Shapiro, and A. G. Palmer Structure, vol. 16, p. 1195, 2008.
  * [96] X. Liao, D. Long, D.-W. Li, R. Brüschweiler, and V. Tugarinov J. Phys. Chem. B, vol. 116, p. 606, 2012.
  * [97] L. E. Kay and D. A. Torchia J. Magn. Reson., vol. 95, p. 536, 1991.
  * [98] V. Calandrini, D. Abergel, and G. R. Kneller J. Chem. Phys., vol. 128, p. 145102, 2008.
  * [99] V. Calandrini, D. Abergel, and G. R. Kneller J. Chem. Phys., vol. 133, p. 145101, 2010.
  * [100] P. Calligari and D. Abergel J. Phys. Chem. B, vol. 116, p. 12955, 2012\.
  * [101] A. V. Buevich and J. Baum J. Am. Chem. Soc., vol. 121, p. 8671, 1999.
  * [102] A. V. Buevich, U. P. Shinde, M. Inouye, and J. Baum J. Biomol. NMR, vol. 121, p. 233, 2001.
  * [103] S. N. Khan, C. Charlier, R. Augustyniak, N. Salvi, V. Déjean, G. Bodenhausen, O. Lequin, P. Pelupessy, and F. Ferrage Biophys. J., vol. 109, p. 988, 2015.
  * [104] A. Hsu, F. Ferrage, and A. G. Palmer Biophys. J., vol. 115, p. 2301, 2018\.
  * [105] W. Adamski, N. Salvi, D. Maurin, J. Magnat, S. Milles, M. Ringkjobing Jensen, A. Abyzov, C. J. Moreau, and M. Blackledge J. Am. Chem. Soc., vol. 141, p. 17817, 2019.
  * [106] P. Dosset, J.-C. Hus, M. Blackledge, and D. Marion J. Biomol. NMR, vol. 16, p. 23, 2000.
  * [107] O. Walker, R. Varadan, and D. Fushman J. Magn. Reson., vol. 168, p. 336, 2004.
  * [108] E. J. d’Auvergne and P. R. Gooley J. Biomol. NMR, vol. 40, p. 121, 2008\.
  * [109] C. A. Castañeda, A. Chaturvedi, C. M. Camara, J. E. Curtis, S. Krueger, and D. Fushman Phys. Chem. Chem. Phys., vol. 18, p. 5771, 2016.
  * [110] G. Barbato, M. Ikura, L. E. Kay, R. W. Pastor, and A. Bax Biochemistry, vol. 31, p. 5269, 1992.
  * [111] C. M. Gall, J. A. DiVerdi, and S. J. Opella J. Am. Chem. Soc., vol. 103, p. 5039, 1981.
  * [112] D. M. Rice, R. J. Wittebort, R. G. Griffin, E. Meirovitch, E. R. Stimson, Y. C. Meinwald, J. H. Freed, and H. A. Scheraga J. Am. Chem. Soc., vol. 103, p. 7707, 1981.
  * [113] K. Zumpfe and A. A. Smith Front. Mol. Biosci., vol. 8, p. 727553, 2021.
  * [114] I. Kuprov, N. Wagner-Rundell, and P. J. Hore J. Magn. Reson., vol. 184, p. 196, 2007.
  * [115] C. Bengs and M. H. Levitt Magn. Reson. Chem., vol. 56, p. 374, 2018.

## Supplementary Material

Figure S1: Carbon longitudinal R1 (a) and transverse R2 relaxation rates (b,), and carbon-proton cross-relaxation rate $\sigma^{\mathrm{NOE}}$ (c) at 14.1 T in a 13C1H2H2 methy group for a correlation time for isotropic global tumbling $\tau_{c}=10$ ns shown as a function of the magnetic field. Two models of motions were used, each of them only considering the methyl rotation: rotation on a cone (blue) and model free (orange). The angle between the C-C and C-H bonds is fixed to $109.47^{\circ}$. The diffusion coefficient for methyl rotation is set to $D_{rot}=5\times 10^{10}$ s-1 leading to a correlation time of 10 ps. Figure S2: Carbon longitudinal R1 (a) and transverse R2 relaxation rates (b,), and carbon-proton cross-relaxation rate $\sigma^{\mathrm{NOE}}$ (c) at 14.1 T in a 13C1H2H2 methy group for a correlation time for isotropic global tumbling $\tau_{c}=10$ ns shown as a function of the correlation time for methyl rotation. Two models of motions were used, each of them only considering the methyl rotation: rotation on a cone (blue) and jumps (orange). The angle between the C-C and C-H bonds is fixed to $109.47^{\circ}$. In the case of jumps between discrete position, the angle between each of them is fixed to $2\pi/3$ and the populations are supposed to be equal for the three conformers. The correlation time for methyl rotation is calculated as $\tau_{met}=\frac{1}{1-\mathcal{S}^{2}_{m}}\sum_{a=-2}^{2}\sum_{\alpha,\beta=1}^{3}\sum_{n,\lambda_{n}!=0}-\frac{\sqrt{p_{\alpha}p_{\beta}}}{\lambda_{n}}\tilde{X}_{\alpha}^{(n)}\tilde{X}_{\beta}^{(n)}\mathcal{D}_{a,0}^{(2)*}(\Omega_{CCH})\mathcal{D}_{a,0}^{(2)}(\Omega_{CCH})$ where $\mathcal{S}^{2}_{m}$ is the order parameter for methyl rotation (identical for both models), $p_{\alpha}=1/3$ is the fractional population for conformer $\alpha$, $\tilde{X}_{\alpha}^{(n)}$ is the value $\alpha$ for the eigenvector of the exchange matrix associated to the eigenvalue $\lambda_{n}$ and $\Omega_{CCH}$ is the Euler angle for transformation from the C-C to C-H bond. The jump rate between each conformer is equal to $k_{ex}$ and the associated diffusion constant for methyl rotation used in the rotation on a cone model is calculated using Eq. 76 of the main text. Figure S3: Contribution of rotamer jump and wobbling in a cone motions to relaxation rates. Carbon longitudinal R1 (a, d, g) and transverse R2 relaxation rates (b, e, h), and carbon-proton cross-relaxation rate $\sigma^{\mathrm{NOE}}$ (c, f, i) at 14.1 T in a 13C1H2H2 methy group for a correlation time for isotropic global tumbling $\tau_{c}=10$ ns. Dashed horizontal lines show the value of the relaxation rates without wobbling. Relaxation rates calculated in the abscence of rotamer jumps are shown in green. In the presence of rotamer jumps, the case where all rotamer populations are equal (blue) and unequal (orange) are distinguished. In this latter case, populations are $p_{1}=0.7$, $p_{2}=0.2$ and $p_{1}=0.1$. Calculations are shown as a function of the correlation time for wobbling, a function of the wobbling diffusion constant and cone semi-angle opening $\beta_{c}$, while the diffusion constant is varied from $10^{6}$ to $10^{10}$ s-1. The blue and orange vertical arrows indicate the values of the correlation time for rotamer jumps, when populations are equal and unequal, respectively. Figure S4: Contribution of rotamer jump and wobbling in a cone motions to relaxation. Carbon R1 (a), R2 (b) and carbon-proton $\sigma^{\mathrm{NOE}}$ (c) at 14.1 T in a hypothetical 13C1H2H2 valine methyl group, and associated spectral density function for the Cβ-Cγ1 bond auto-correlations. The model of motion includes the methyl group rotation as well as (when mentioned) rotamer jump of the Cβ-Cγ bond and (when mentioned) the wobbling in a cone of the Cα-Cβ bond. Dashed lines show the value of the relaxation rates and spectral density function without wobbling. Relaxation rates (a-c) and spectral density function (d) calculated in the abscence of rotamer jump are shown in green. In the presence of rotamer jump, the case where all rotamer populations are equal (blue) and unequal (orange) are distinguished. In this latter case, populations are $p_{1}=0.7$, $p_{2}=0.2$ and $p_{1}=0.1$. Calculations are shown as a function of the correlation time for wobbling, a function of the wobbling diffusion constant and cone semi-angle opening, set here to $\beta_{c}=15$ deg, while the diffusion constant $D_{W}$ is varied from $10^{6}$ to $10^{10}$ s-1. The blue and orange vertical arrows indicate the values of the correlation time for rotamer jump, respectively when populations are equal and unequal. The spectral density functions are shown for $D_{W}=10^{8}$ s-1. Table S1: Values of the fitted MF parameters (correlation function Eq. 147) for 4 types of internal motions and the three possible symmetry properties of the overal diffusion tensor. In the case of the symmetric top and fully asymmetric overal diffusion tensor, we report the average and standard deviation over all possible orientations, i.e. the different values of $\theta_{D,M}$ and $\\{\varphi_{D,M},\theta_{D,M}\\}$ respectively. In these fits, the values and orientation of the diffusion tensor (angle $\Omega_{D,M}$) were fixed such that $\mathcal{S}_{MF}^{2}$ and $\tau_{MF}$ were the only adjustable parameters. The diffusion tensor for global tumbling has $D_{xx}=10^{7}$ s-1, $D_{yy}=5\times 10^{7}$ s-1 and $D_{zz}=10^{8}$ s-1. When the symmetric model is used, eigenvalues are $D_{\parallel}=D_{zz}$ and $D_{\perp}=\frac{1}{2}(D_{xx}+D_{yy})$. For the isotropic diffusion, the eigenvalue equals $D=\frac{1}{3}(D_{xx}+D_{yy}+D_{zz})$. The orientation of the rotamer frames in the jump frame are $\Omega_{J,R}^{(2)}=\\{\pm\frac{\pi}{2},\beta_{CC}\\}$ when two states are considered, and $\Omega_{J,R}^{(3)}=\\{2n\frac{\pi}{3},\beta_{CC}\\}$ for the 3-state jump with $n=0,1,2$ and where $\beta_{CC}=76^{\circ}$. The 3-state jump model with unequal populations (bottom row) was simulated with $p_{1}=0.7$, $p_{2}=0.2$ and $p_{3}=0.1$. The rotation on a cone is modeled assuming a methyl group gegometry and with $D_{rot}=10^{11}$ s-1. For the wobbling motion, the cone semi-angle opening is $\beta_{c}=20^{\circ}$ and $D_{W}=10^{8}$ s-1. | Isotropic | Symmetric top | Fully asymmetric
---|---|---|---
| $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps) | $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps) | $\mathcal{S}_{MF}^{2}$ |  $\tau_{MF}$ (ps)
rotation | $0.12$ | $4.5$ | $0.12\pm 0.00$ | $4.7\pm 0.0$ | $0.11\pm 0.00$ | $4.9\pm 0.0$
wobbling | $0.83$ | $340$ | $0.83\pm 0.00$ | $348\pm 27$ | $0.83\pm 0.00$ | $351\pm 27$
2-state jump | $0.82$ | $612$ | $0.56\pm 0.24$ | $1,368\pm 709$ | $0.52\pm 0.29$ | $1,840\pm 1,181$
$p_{1}=p_{2}$
3-state jump | $0.16$ | $319$ | $0.16\pm 0.01$ | $331\pm 21$ | $0.17\pm 0.00$ | $329\pm 20$
$p_{1}=p_{2}=p_{3}$
3-state jump | $0.43$ | $153$ | $0.43\pm 0.11$ | $188\pm 117$ | $0.42\pm 0.12$ | $186\pm 169$
$p_{1}\neq p_{2}\neq p_{3}$