Origin of orientation‐dependent R1 (=1/T1) relaxation in white matter

In a recent MRI study, it was shown that the longitudinal relaxation rate, R1, in white matter (WM) is influenced by the relative orientation of nerve fibers with respect to the main magnetic field (B0). Even though the exact nature of this R1 orientation dependency is still unclear, it can be assumed that the origin of the phenomenon can be attributed to the anisotropic and unique molecular environment within the myelin sheath surrounding the axons. The current work investigates the contribution of dipolar induced R1 relaxation of the myelin associated hydrogen nuclei theoretically and compares the results with the experimentally observed R1 orientation dependency.


| INTRODUCTION
Modeling and reconstruction of quantitative magnetic resonance (MR) parameters are powerful tools for extracting physical and structural properties of neural tissue in vivo. [1][2][3] Several studies have demonstrated that the extracted tissue features are directly linked to pathological processes such as edema formation, inflammation, or demyelination and that they can serve as biomarkers to assess the course of neurological diseases or to study brain development. [4][5][6] A technique commonly employed to extract tissue features is based on the reconstruction of the MR signal evolution employing suitable biophysical models. The accuracy of this technique, however, is constrained by the limited knowledge of the tissue-specific MR properties and their relative influence on the accessible MR signal on a macroscopic scale. A possible way to address this issue is based on a modeling of brain microstructures on an atomistic scale and employing a molecular dynamics (MD) simulation to theoretically determine the relevant MR properties. We recently introduced such a framework for a trajectory-based determination of the dipolar induced R 1 relaxation rate of hydrogen nuclei within the myelin sheath. 7 In a previous study, we also showed that the R 1 rate in white matter (WM) structures is influenced by the orientation of nerve bundles relative to the main magnetic field. 8 Specifically, a monotonically increasing R 1 relaxation rate from 1.277 Hz (T 1 = 783 ms) for nerve bundles oriented parallel to 1.311 Hz (T 1 = 763 ms) for nerve bundles oriented perpendicular relative to the B 0 -field was observed at 3 T in WM regions with a fractional anisotropy (FA) >0. 4. If one compares the observed T 1 difference of ≈ 20 ms with the contrast variance of approximately 20 ms to 80 ms in identical brain regions, it is obvious that identifying the source of this phenomenon is important for a more comprehensive understanding of the R 1 heterogeneity in WM regions. [9][10][11] Moreover, given the rapid decrease of dipolar interactions with increasing distance, the R 1 orientation dependency likely originates from hydrogen nuclei located within the myelin sheath. Investigating and understanding this effect on a molecular basis may thus enable the extraction of tissue characteristics for investigating myelin integrity. Approximately 40% of these nuclei exist in an aqueous phase trapped between lipid membranes termed myelin water (MW). In contrast to interstitial and axonal water, MW predominantly diffuses close to the lipid membranes, which influences its dynamical properties in several ways. The grained and complex membrane surface, for example, decreases the mobility of the individual water molecules and thus enhances their longitudinal relaxation rate. 12,13 In addition, the hydrophilic character of the phospholipids binds the hydrogen atoms of water molecules to the membrane surface. This effect decreases the configurational entropy in the interfacial water and restricts molecular tumbling along the membrane normal as already demonstrated in scattering and nuclear magnetic resonance experiments. 14,15 Intramolecular induced R 1 relaxation in water, however, is governed by dipole-dipole (DD) couplings between the two hydrogen nuclei and is thus directly related to the water dynamics relative to the B 0 -field.
The working hypothesis of the current work was, therefore, that the structural polarization of water molecules leads to an anisotropic DD interaction, and thus, to an orientation-dependent R 1 rate in MW. Because myelin water molecules make up about 10%-15% of the total white matter water content, this effect could be strong enough to alter the detectable MR signal and contribute to the observed R 1 orientation dependency. Moreover, the structural configuration of hydrogen nuclei within the lipid membrane of the myelin sheath may play an important role. These nuclei reside in a spatially and temporally asymmetric spin environment that also causes an anisotropic DD interaction with the neighboring protons. This results in an orientation-dependent relaxation, as already demonstrated in previous nuclear magnetic resonance studies on model lipids. [16][17][18] Even though the hydrogen nuclei within the membrane are invisible in conventional MRI, protons of the membrane surface are in direct contact to MW and influence the detectable signal indirectly 19 through magnetization transfer (MT) effects such as cross-relaxation or chemical exchange. An accurate estimation of the relevance of the structural orientation of the hydrogen nuclei in both MW and the membrane surface on the MR observable R 1 orientation dependency, however, requires an explicit analysis of the molecular dynamics within the myelin sheath, combined with a realistic signal weighing scheme.
The aim of the current work was, therefore, to apply the MD-framework of myelin presented in Ref. 7 to investigate the relevance of the effects described previously on the experimentally observed R 1 orientation dependency in human white matter. 8 For this purpose, we first carried out molecular dynamics simulations of a lipid-water system representing a small fraction of the membrane structure inside the myelin sheath. The trajectories obtained from the simulation were employed to determine the DD induced R 1 relaxation rate of the relevant 1 H protons at different orientations of the membrane structure relative to the B 0 -field. Based on the calculated relaxation rates, a three-pool model comprising MW, free water (FW), and the hydrogen nuclei of the membrane surface representing the solid myelin (SM) pool was employed to simulate the effective R 1 rate in WM voxels at different fiber orientations. The results obtained clearly demonstrate that the calculated R 1 rate is significantly higher in both SM and MW for perpendicular oriented fibers. More important, the theoretically determined R 1 rate matches the experimental results reported in Ref. 8

| System preparation and molecular dynamics simulations
To simulate the dynamics of hydrogen nuclei within the myelin sheath, two MD simulations were performed using the GROMACS package 5.1.2 20 in combination with the CHARMM36 all-atom force field. 21 The molecular system was prepared according to a realistic lipid/water composition as described in Ref. 7 In short, the myelin sheath was represented by two lipid layers solvated in 2276 TIP4p-FB water molecules. 22 23 After a 50 ns thermodynamic equilibration to body temperature (310.15 K) and normal pressure, two production runs with identical settings but different simulation times, T MD , were performed. The shorter first run was conducted for T S MD = 50 ns with a sampling frequency of ΔT −1 = 0.25 −1 ps −1 and provided the trajectory set to determine the R 1 relaxation rate of hydrogen nuclei located in MW. The second run was performed for T L MD = 250 ns and a sampling frequency of ΔT −1 = 2 −1 ps −1 . The trajectory set of the longer run was used to determine the R 1 rate of motionally restricted hydrogen nuclei located within the membrane surface. After the production runs, the trajectories were recomposed to the initial molecular topology to avoid broken structures caused by the periodic boundary conditions in the simulation using the subroutines provides by the GROMACS software package. Finally, the hydrogen density profile of the full SM pool was calculated from the long-term trajectory set in order to label those hydrogen protons within the membrane that are in physical contact with water molecules. The density profile is shown in Figure 1A and the analysis reveals that approximately 32% (blue shaded area in Figure 1A) of the membrane-associated hydrogen nuclei were able to effectively interact through short-range interactions with MW, which hereinafter is referred to as the effective interacting solid myelin fraction, SM = 0.32.

| Computation of the R 1 orientation dependence
Based on the trajectory set of the water molecules from the short-term (T S MD = 50 ns) simulation, the longitudinal relaxation rate of the MW, R 1,MW , was determined as a function of the fiber orientation as follows: First, a hydrogen nucleus was randomly selected and the relative positions in relation to the remaining hydrogen nuclei in the MW were computed F I G U R E 1 Hydrogen density profile of different lipid species, the total membrane (Memb.) and water molecules (A). The blue shaded area represents the 1 H protons of the membrane surface that are in direct contact with water molecules, which accounts for approximately 32% of the total hydrogen fraction of the SM pool. Three-pool model for WM with labeled parameters used in the main text (B) for each time step. Then, the trajectories were transformed in polar coordinates and the R 1 relaxation rate was calculated from the spectral densities by explicit integration of the autocorrelation function in the interval = 0, … ,0.3T S MD according to the method described in Ref. 7 To sample the R 1 rate of the selected water molecule at various segment locations around the axon, the trajectory set was successively rotated in 16 steps counterclockwise by Δ = 22.5 • around the z-axis as illustrated in Figure 2A. At each rotation step, the relative positions to the other hydrogen nuclei were redetermined, from which the R 1 rate was calculated for each rotation angle, . Then, to determine R 1 for different fiber orientations, the trajectory set was rotated by an angle around the y-axis of the laboratory frame ( Figure 2A) and the procedure described previously was repeated by again rotating trajectories in Δ = 22.5 • steps around the z'-axis of the new coordinate system. The rotation around the y-axis (angle ) was performed 12 times in steps of Δ = 8.1°, resulting in a maximum angle of max = 90°, which corresponds to a perpendicular orientation of axons relative to B 0 . The procedure described previously was performed for 26% (N = 1200) of all hydrogen nuclei in MW, resulting in relaxation rates, R 1,MW ( , ), for 192 combinations of and for each proton. From this data set, the -dependent average relaxation rate, R N 1,MW ( , ) = ⟨R 1,MW ( , )⟩ N , and the effective relaxation rate, R 1,MW ( ) = ⟨⟨R 1,MW ( , )⟩ ⟩ N , were calculated, where <...> denotes either the average or the median of the N = 1200 hydrogen nuclei and the rotation angle , respectively. In addition, to exclude outliers and demonstrate the robustness of the effect, the effective relaxation rate was also calculated using the 5%-trimmed mean value. Finally, the orientation-dependent relaxation shifts, were computed from the effective relaxation rates.
The determination of the relaxation rate of hydrogen nuclei located in the surface area of the SM,R 1,SM , was completely analogous to the approach described for MW. Here, a single hydrogen nucleus from the membrane surface was randomly selected from the long-term (T L MD = 250 ns) trajectory set for which the autocorrelation function and the relaxation rate were calculated as described previously. This step was repeated for 1289 randomly selected hydrogen nuclei using the same combinations of and as for MW. Finally, the -dependent average (median) relaxation rate, R N 1,SM ( , ), the effective relaxation rate, R 1,SM ( ), and the orientation-dependent relaxation shifts, ΔR 1,SM ( ), were computed in complete analogy to protons located in MW.

| Three-pool model and R 1 ( ) fitting
To determine the MR observable relaxation rate, the threepool model illustrated in Figure 1B was employed. This model is based on the four-pool model introduced by Harrison et al, 24 which has already been used in previous studies to simulate the MR signal evolution in WM regions, considering MT effects between adjacent tissue compartments. [25][26][27] For the original model, however, the authors suggested that the relaxation behavior in WM structures can be described with four different tissue compartments: MW, the interstitial/axonal water termed as FW, immobile macromolecules and nonaqueous myelin that is replaced in F I G U R E 2 Schematic illustration of the laboratory frame (black), the rotated frame (red), and a snapshot of the last frame from the long-term MD simulation (A). Flowchart of the optimization procedure described in the Methods section (B) the current study by the effectively interacting hydrogen fraction in SM. However, because macromolecules in the immobile macromolecular pool are predominantly associated with neurons and glia cells located outside the myelin sheath, it can be assumed that (1) their structural orientation is on average isotropic when compared to the lipids within the myelin sheath and that (2) they interact exclusively with water molecules in the corresponding space outside of myelin. For the current study, we therefore assumed that macromolecules only lead to an orientation-independent R 1 offset in such a way that the fourth pool can effectively be included in the FW pool. In the case of a linear magnetization transfer rate between the pools, the time evolution of the longitudinal magnetization can then be written as set of coupled Bloch equations given by where M i (t, ) are the time-and orientation-dependent pool magnetizations, M Eq i the corresponding equilibrium magnetizations, k ij the directional magnetization exchange rates from the i-th to the j-th pool and R offset 1,i the orientation-independent offsets of the pool-specific relaxations rates. The factor in Equations 1 and 2 scales the magnitude of the orientation-dependent shift and was introduced to take into account the imperfect alignment of fiber trajectories within the WM voxels. The physical meaning and adequacy of this shape factor is discussed in detail in the Discussion section.
The three-pool model given by Equations (1-3) has 11 degrees of freedom (three-pool magnetizations (M i ) and the respective R 1 offsets (R offset 1,i ), four exchange rates (k ij ), and the shape factor ). The number of free parameters can be reduced by assuming that the net magnetization transfer rates between two pools as well as the individual compartment sizes are time independent. Using these assumptions, the magnetization transfer rates are mutually related and satisfy 28 where p i are the probabilities to find a hydrogen nuclei in pool i, and T MW cr (T D Cr ) is the cross-relaxation time describing the magnetization exchange process between two pools. In addition, because the total water content in a WM voxel is time independent, the probability of finding a hydrogen nuclei in the FW pool can be expressed by p FW = 1 − p MW . A similar relationship can also be used for p SM when the sizes of the SM and the MW pool are linearly correlated. Here, a fixed ratio of 0.4 between the density of water hydrogen protons and lipid hydrogen protons is assumed for myelin, consistent with the value employed in the MD simulation. 29 The probability of finding a hydrogen nuclei in the surface area of the membrane can then be written as 32 is the effective interacting solid myelin fraction as described previously. The assumption of strictly correlated pool sizes for the myelin sheath, together with Equations 4 and 5, reduces the degrees of freedom from 11 to 7. An iterative optimization procedure was performed (see Figure 2B Because the net magnetization in WM comprises components from different molecular environments that are in mutual exchange, the decay process is more complex and should be described by a multicomponent model. 30 In our experimental study, however, the decay of the MR signal was reconstructed assuming a monoexponential decay. Therefore, Equation 6 is adequate to compare the theoretical predictions obtained in the current work with the experimental data.
As can be seen from Equation 6, the optimization was performed by minimizing the mean squared residual between the simulated and the experimentally observed R 1 values at different angles using the gradient descent method. In order to obtain realistic tissue characteristics, the pool parameters were constrained within the intervals shown in Table 1, whereby the exchange rate between SM and MW was adjusted to the surface hydrogen fraction of 32% to take into account the reduced pool size of the SM pool. The optimization procedure was carried out using both the average (Case I) and the median (Case II) of the orientationdependent relaxation shifts, ΔR 1,SM ( ) . The latter was chosen to exclude the influence of outliers. A potential bias might be introduced by restricting the analysis to surface hydrogen protons of SM without any extension of the three-pool-model. This effectively neglects the magnetization transfer between head and tail groups due to spin diffusion along the alkyl chain. The whole analysis was thus repeated by including all 1 H protons within the SM compartment. In accordance with previous simulation studies, it was assumed that water molecules interact equally with all protons of the SM pool, which corresponds to an effectively interacting hydrogen fraction of SM = 1. R 1,SM ( ) was then determined for 2845 protons randomly selected from both tail and head groups. Based on the calculated R 1 rates, the average (Case III) and the median (Case IV) R 1 orientation dependencies were calculated and the optimization procedure was repeated using identical starting conditions and fitting constraints as described previously. Figure 3A shows the -dependent average relaxation rate of MW as a function of segment position ( ) and fiber orientation ( ). As is evident from this figure, the relaxation rate is nearly constant for < 10 • (R N 1,MW ( < 10 • , ) ≈ 0.91 Hz) and oscillates periodically in between 0.91 Hz and 1.02 Hz for > 10 • . Figure 3B depicts the effective R 1 rate of MW averaged over the different segment positions ( ). The average effective relaxation rate increases with sinusoidal shape from R 1,MW ( = 0 • ) = 0.91 Hz to R 1,MW ( = 90 • ) = 0.96 Hz with increasing , which corresponds to a total increase of 0.05 Hz. Comparable results, but with a reduced rate at = 0 • , were also found for the 5%-trimmed mean and the median relaxation rate of the individual protons, respectively. The relaxation rate of hydrogen nuclei located on the membrane surface of SM shows an angular dependence similar to that observed for MW ( Figure 4). Here, the -dependent average relaxation rate approaches 3.4 Hz for small and oscillates between 3.4 Hz and 4.1 Hz as a function of segment position at = 90 • . The effective relaxation rate increases from R 1,SM ( ) = 3.4 Hz for parallel fibers ( = 0 • ) to approximately 3.76 Hz for fibers orientated perpendicular ( = 90 • ) to the B 0 -field, as can be seen in Figure 4B. For the 5%-trimmed mean, a slightly larger increase from 3.15 Hz to 3.61 Hz is observed whereas the median relaxation rate increases from 2.6 Hz to 3.25 Hz. The analysis in which all protons of the SM pool were included yielded comparable results, apart from an overall reduction of R 1,SM (≈ 0.8 Hz) and a slightly reduced orientation-dependent offset (by ≈ 0.2 Hz). Figure 5 shows the predictions of the three-pool model for all four cases (Case I: mean, SM = 0.32; Case II: median, SM = 0.32; Case III: mean, SM = 1; and Case IV: median, SM = 1) along with the data measured in the experimental study. 8 The corresponding fit parameters are reported in Table 1. An increase of the relaxation rate with is predicted regardless of the model employed. All four cases are consistent with the experimental data within the error bars, although cases II-IV fit the data with a larger r 2 , as is evident from Figure 5.

| DISCUSSION
In this study, we investigated the orientation dependence of the dipole-dipole induced R 1 relaxation within myelin and its surrounding water compartment. For this purpose, we carried out atomistic MD simulations of a myelin-alike environment and calculated the dipolar induced R 1 relaxation rate of myelin water and the proton compartment in solid myelin that is in direct contact with myelin water for different orientations of the trajectory sets relative to the main magnetic field. As is evident from Figures 3 and 4, R 1 rates of both solid myelin and myelin water clearly depend on the location within the axon ( ) and the nerve fiber orientation ( ) relative to the B 0 -field. This can be attributed to the anisotropic molecular tumbling within myelin water and the heterogeneous spin environment within the solid myelin compartment. Due to the hydrophilicity of the lipid heads, water molecules form an ordered hydration network near the membrane surfaces. In this network, the reorientation of the individual water molecules is restricted, which leads to an asymmetric residual coupling between the hydrogen protons and causes the R 1 orientation dependency shown in Figure 3. The R 1 orientation dependence of myelin water can indeed be attributed exclusively to intramolecular dipoledipole-induced relaxation (see Supporting Information Figure  S1 and Supporting Information Figure S2). A similar mechanism also causes the orientation dependence of the R 1 rate of the hydrogen nuclei within the membrane surface. The anisotropic coupling leads to an asymmetrical interaction Hamiltonian that causes the orientation-dependent R 1 relaxation. [16][17][18] However, if

F I G U R E 4
-dependent average relaxation rate of hydrogen nuclei located in the membrane surface area of SM as a function of segment position ( ) and fiber orientation ( ) (A) and effective relaxation rates averaged over (B) one compares the different predicted R 1 orientation dependencies shown in Figure 4, it is obvious that the results are sensitive to outliers. This also affects the fitting of the experimental data (see Figure 5 and subsequent discussion).
Based on the calculated pool-specific R 1 rates, a threepool model was employed to determine the MR-detectable signal from which the experimentally accessible R 1 rate was predicted for different orientations.
Fitting the model to the experimentally observed data from a previous study resulted in a good agreement for all four cases (r 2 = 0.851 to 0.998).
However, the question of which of the pools (SM or MW) predominantly governs the observable R 1 orientation dependencies remains open. This question can be answered if one assumes a fast magnetisation transfer (ie, 1/T D cr ≫ R 1,MW ) between free water and myelin water, neglecting the direct magnetization transfer with the solid myelin pool. The observable R 1 rate is then given by the sum of the pool-specific R 1 rates weighted by the respective hydrogen fraction. This fastexchange model, however, would yield an orientationdependent shift that is approximately one order smaller than the shift of the myelin water alone, given the small (approx. 10%-15%) fraction of myelin water. This would clearly contradict the experimental data. Therefore, the orientation dependence of R 1 in WM should be attributed to the asymmetrical spin environment within the myelin sheath and its interaction with the myelin water compartment.
Concerning the functional relationship between R 1 and fiber orientation, it is obvious that the model-predicted relaxation rate clearly increases monotonically with increasing angle. This observation, however, is in contrast to the findings reported in a recent MRI study by Knight et al. 31 Here, the authors reported a decreasing R 1 rate for 0 • < < 55 • followed by an increase for > 55 • . This contradicts the theoretical predictions of the current work as well as the experimental results reported in Ref. 8

| Predicted model parameters
As already discussed in Ref. 7, a quantitative comparison of the predicted and the calculated pool-specific R 1 rates with experimental findings at 3 T is rather complicated. The main challenges in the experimental quantification of R 1 in myelin are the accurate determination of the individual compartment sizes from the raw signal and the reconstruction of the intrinsic relaxation rates considering magnetization transfer between the pools. Nonetheless and as stated in Ref. 7, relaxation rates of R 1,SM ≈ 3-4 Hz for SM and R 1,MW ≈ 0.9 Hz to 1.1 Hz for MW at 3 T are in accordance with experimental observations on model lipids. 7 However, the relaxation rates of R offset 1,MW = 0.51 Hz (0.57 Hz) and R offset 1,SM = 2.51 Hz (2.57 Hz) predicted for cases I (III) are slightly too small. Moreover, when comparing the other pool characteristics predicted for those cases, it becomes clear that most parameters were adjusted to the predefined fitting constraints. Although the parameters are physically reasonable, results obtained appear to be unrealistic in terms of a realistic tissue composition, which indicates that the orientation-dependent offset calculated from the averaged R 1 rates might be too small to provide a satisfactory explanation for the experimental data.
In contrast, the predicted model parameters for cases II and IV are within reasonable ranges and represent typical nuclear magnetic resonance characteristics of WM structures. Concerning the relaxation rates of the free water, it is also obvious that R offset 1,FW is too large when compared with the R 1 relaxation rate of bulk water at 3 T (1.1 Hz vs. 0.3 Hz). However, as already mentioned in the Methods section, the FW pool includes the immobile macromolecules located in water compartments outside the myelin sheath. Even though these macromolecules are invisible in routine MRI, they interact through MT with the surrounding water and increase the effective R 1 . The application of a three--instead of a four--pool model is also important for the comparison of MT exchange rates reported by Kalantari et al. 27 Here, the authors investigated the influence of different R 1 scenarios on the predicted cross-relaxation times and show that the pool-specific R 1 rate is a crucial parameter affecting the predicted exchange rates. This has to be kept in mind when comparing our results with the findings of other studies. Nonetheless, the results of T D Cr ≈ 0.5 s and T MW Cr = 0.1 s to 0.14 s are in accordance with the cross-relaxation times of T D Cr = 0.5 s to 1.0 s and T MW Cr = 0.08 s to 0.227 s reported in 25,26 using fourpool models. Comparing cases I and II, it seems likely that outliers have a significant effect on the results obtained. This is not unexpected as only a relatively small proton sample (1200 in MW and 2845 in SM) was simulated due to constraints with the available hardware components. Given the fact that the interaction Hamiltonian scales with 1∕r 6 , where r is the distance between two protons, the average value is expected to be overestimated by a small fraction of nuclei that are in very close contact, consistent with the results shown in Figures  3 and 4. Therefore, the median value should be less biased and the corresponding results are likely to be more reliable. This is especially true for head-group protons that are relatively mobile. Inclusion of the full SM compartment resulted in a reduced difference between mean and median relaxation rates as can be seen from Figure 5. This issue should be further addressed using larger sample sizes.

| Interpretation of the shape factor and limitations
To adapt the parameters of the three-pool model to the experimental data, the shape factor α was introduced in Equations 1 and 2. The adequacy of this parameter can be explained by the fact that the fiber pathways in a single voxel are in general not purely parallel, which decreases the experimentally measured R 1 orientation shift. For example, in voxels with isotropic internal fiber orientations (FA = 0), the measured R 1 rate should be independent of orientation. In contrast, for a voxel with exclusively parallel fibers (FA = 1), or for a single nerve fiber with a comparable lipid-water ratio as simulated above, the measurable orientation-dependent R 1 shift should be maximal. In our experimental study, 8 the R 1 values were determined from voxels with FA >0.4 (FA = 0.59 ± 0.016), which corresponds to a weighted intermediate between the two extreme cases. In order to address this issue, the shape factor was introduced to scale the magnitude of ΔR 1,SM ( ) and ΔR 1,MW ( ) for the experimental case where 0.4 < FA < 1.
However, as can be seen in Table 1, the predicted values of α are larger than the average experimental FA. Even though this seems to be an apparent contradiction, two issues have to be considered when comparing to the average FA of 0.59 reported in Ref. 8: First, the FA in the experimental study was computed using the standard single-tensor model 32 provided by the FSL software package. 33 Even though this model is widely employed to investigate the diffusivity in brain tissue, it significantly underestimates the FA in crossing fiber regions. [34][35][36] Second, it should be noted that the two quantities are not inevitably correlated as they describe different physical/geometrical characteristics. In contrast to the shape factor, which indicates the orientational dispersion of the nerve fibers, FA is derived from the diffusion tensor and describes the isotropy of the molecular motion. The diffusion tensor, however, also includes the molecular motion orthogonal to the main fiber orientations, and therefore, depends on tissue characteristics such as the axonal diameter, fiber packing density, and degree of myelination. A large axonal diameter or a reduced fiber packing density, however, increases the diffusivity orthogonal to the fiber directions, resulting in FA <1 even in voxels with purely parallel fibers. In contrast, in this case the measurable orientation dependency should be comparable to the calculated results, and the R 1 fitting would yield a shape factor of = 1. Nonetheless, even though these points might partly explain the discrepancies between FA and the shape factor, the mismatch between both quantities is a limitation of the current study that has to be considered.
Apart from this point, the current study is subject to some limitations that are important to mention. One point to note is that a three-pool model was applied to simulate the MR signal in WM. Even though this model allows a proper reconstruction of the experimental data, it greatly simplifies the tissue composition in WM. A more realistic model, for example, should separate the FW pool into an intracellular and extracellular compartment and include at least one supplementary pool for the immobile macromolecules outside the myelin sheath. Moreover, the SM pool in the current work comprises only hydrogen nuclei within the lipid membrane that are in close contact with water molecules. Although this representation seems to be reasonable because MT is mainly governed by short-range interactions, it must be considered that hydrogen nuclei located in the tail groups interact with protons in the head groups by spin diffusion along the chain and could thus affect the magnetization transfer to MW. However, inclusion of the full solid myelin pool into the exchange model (Equations 1-3) did not significantly alter the results obtained, although the fitting parameters varied slightly (see Table 1). This variation represents the typical systematic error that can be attributed to uncertainties of the MT dynamics between MW and SM.
Furthermore, it should be mentioned that the calculated R 1 rates are only due to DD interactions between hydrogen nuclei from the same compartment, for example, MW or SM. A full and comprehensive analysis, however, should consider that the hydrogen nuclei in one compartment also interact via through-space DD interactions with the hydrogen nuclei in adjacent compartments. Because the pools are spatially separated and DD interactions are directly related to the relative position of the nuclei, it can be assumed that such intercompartmental interactions are also orientation-dependent. As we have shown in Ref. 7, however, these effects are smaller by approximately one order of magnitude than the intermolecular and intramolecular induced R 1 relaxation within a single compartment, and therefore, negligible for the current work.
Another effect that was not explicitly considered in the current work is the contribution of susceptibility differences between water and myelin. Several studies investigating the orientational dependence of R * 2 have shown that susceptibility differences between water and solid myelin affect the coherence of spin ensembles depending on the fiberto-field angle. For example, the work of Knight et al suggests that diffusion of 1 H in an inhomogeneous magnetic field in the vicinity of cylindrical structures results in an orientation-dependent nonlinear phase evolution and phase decoherence. 37 However, it was concluded that this effect is not responsible for a reduction of the transverse relaxation time. Apart from different assumptions made about fiber geometry, fiber distribution, or the compartment being analyzed, results obtained in studies on R 2 or R * 2 cannot be directly transferred to R 1 . This is due to the fact that R 1 is governed by fluctuations of the transverse field around the Larmor frequency whereas R 2 is sensitive to longitudinal field fluctuations at ω = 0 Hz. However, a recent molecular dynamics study has shown that diffusion-mediated R 1 relaxation in myelin water can indeed be neglected as compared to the dipolar induced mechanism. 38

| CONCLUSION
The aim of this work was to investigate whether dipolar induced R 1 relaxation within the myelin sheath contributes to the orientation-dependent R 1 contrast in WM regions. For this purpose, atomistic MD simulations were carried out and the dipole-dipole induced R 1 relaxation rate was determined for different fiber orientations relative to the B 0 -field. The results obtained clearly demonstrate that the R 1 rate in both myelin water and solid myelin depends on the fiber orientation relative to B 0 . The model employed here indeed allows the reconstruction of the experimental data of R 1 using realistic tissue characteristics and two different exchange scenarios between solid and aqueous myelin pools. Even though the applied model is, as discussed previously, subject to some limitations, the results strongly suggest that the orientation-dependent R 1 contrast in WM originates primarily from the anisotropic molecular environment within the myelin sheath and can be explained by dipolar induced relaxation. In addition, the outcomes show that the R 1 heterogeneity in WM is not determined just by the ratio of bounded and free hydrogen nuclei but also by the structural order and orientation relative to the external field. The characterization of this relationship may be useful for future studies investigating WM structures and could enable the more precise extraction of tissue features, that is, by employing more accurate T 1 mapping schemes such as the one presented in Ref. 39