Front. Neurosci. Frontiers in Neuroscience Front. Neurosci. 1662-453X Frontiers Media S.A. 10.3389/fnins.2023.1133086 Neuroscience Original Research Fiber-orientation independent component of R2* obtained from single-orientation MRI measurements in simulations and a post-mortem human optic chiasm Fritz Francisco J. 1 * Mordhorst Laurin 1 Ashtarayeh Mohammad 1 Periquito Joao 2 Pohlmann Andreas 2 Morawski Markus 3 4 Jaeger Carsten 3 4 Niendorf Thoralf 2 Pine Kerrin J. 4 Callaghan Martina F. 5 Weiskopf Nikolaus 4 6 Mohammadi Siawoosh 1 4 7 * 1Department of Systems Neurosciences, University Medical Center Hamburg-Eppendorf, Hamburg, Germany 2Berlin Ultrahigh Field Facility (B.U.F.F.), Max-Delbrueck-Center for Molecular Medicine in the Helmholtz Association, Berlin, Germany 3Paul Flechsig Institute – Center for Neuropathology and Brain Research, University of Leipzig, Leipzig, Germany 4Department of Neurophysics, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany 5Wellcome Centre for Human Neuroimaging, UCL Queen Square Institute of Neurology, University College London, London, United Kingdom 6Felix Bloch Institute for Solid State Physics, Faculty of Physics and Earth Sciences, Leipzig University, Leipzig, Germany 7Max Planck Research Group MR Physics, Max Planck Institute for Human Development, Berlin, Germany

Edited by: Helene Ratiney, CREATIS, Unité CNRS UMR 5220 – INSERM U1294 – Université Lyon 1 – INSA Lyon, France

Reviewed by: José P. Marques, Radboud University, Netherlands; Christoph Birkl, Innsbruck Medical University, Austria

 *Correspondence: Francisco J. Fritz, f.lagosfritz@uke.de Siawoosh Mohammadi, mohammadi@mpib-berlin.mpg.de 25 08 2023 2023 17 1133086 28 12 2022 04 08 2023 Copyright © 2023 Fritz, Mordhorst, Ashtarayeh, Periquito, Pohlmann, Morawski, Jaeger, Niendorf, Pine, Callaghan, Weiskopf and Mohammadi. 2023 Fritz, Mordhorst, Ashtarayeh, Periquito, Pohlmann, Morawski, Jaeger, Niendorf, Pine, Callaghan, Weiskopf and Mohammadi

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The effective transverse relaxation rate (R2*) is sensitive to the microstructure of the human brain like the g-ratio which characterises the relative myelination of axons. However, the fibre-orientation dependence of R2* degrades its reproducibility and any microstructural derivative measure. To estimate its orientation-independent part (R2,iso*) from single multi-echo gradient-recalled-echo (meGRE) measurements at arbitrary orientations, a second-order polynomial in time model (hereafter M2) can be used. Its linear time-dependent parameter, β1, can be biophysically related to R2,iso* when neglecting the myelin water (MW) signal in the hollow cylinder fibre model (HCFM). Here, we examined the performance of M2 using experimental and simulated data with variable g-ratio and fibre dispersion. We found that the fitted β1 can estimate R2,iso* using meGRE with long maximum-echo time (TEmax ≈ 54 ms), but not accurately captures its microscopic dependence on the g-ratio (error 84%). We proposed a new heuristic expression for β1 that reduced the error to 12% for ex vivo compartmental R2 values. Using the new expression, we could estimate an MW fraction of 0.14 for fibres with negligible dispersion in a fixed human optic chiasm for the ex vivo compartmental R2 values but not for the in vivo values. M2 and the HCFM-based simulations failed to explain the measured R2*-orientation-dependence around the magic angle for a typical in vivo meGRE protocol (with TEmax ≈ 18 ms). In conclusion, further validation and the development of movement-robust in vivo meGRE protocols with TEmax ≈ 54 ms are required before M2 can be used to estimate R2,iso* in subjects.

 effective transverse relaxation rate biophysical model R2* orientation-independent R2* myelin water fraction g-ratio fibre dispersion multi-echo gradient recalled echo section-at-acceptance Brain Imaging Methods Introduction

The effective transverse relaxation rate (R2* = 1/T2*) is a nuclear magnetic resonance (NMR) relaxation property (Tofts, 2004) that enables non-invasive characterisation of the microstructure of the human brain (MacKay et al., 2006; Does, 2018; Weiskopf et al., 2021). The microstructural sensitivity of R2* makes it particularly interesting for neuroscience and clinical research studies (Langkammer et al., 2010; Draganski et al., 2011; Callaghan et al., 2014; Kirilina et al., 2020). This is because R2* is sensitive not only to free and myelin water pools in the brain (MacKay et al., 2006; Dula et al., 2010; Weiskopf et al., 2021) but also to microscopic perturbations in the main magnetic field ( B 0 → ) (Chavhan et al., 2009). These microscopic perturbations are caused by the different magnetic susceptibilities of biological structures (Duyn and Schenck, 2017) like the diamagnetic myelin sheath (Kucharczyk et al., 1994; Duyn, 2014; Lee et al., 2017; Alonso-Ortiz et al., 2018) and paramagnetic iron deposits in glial cells (Ordidge et al., 1994; Li et al., 2009; Yao et al., 2009). Moreover, it has been shown that R2* is also strongly dependent on the angular orientation of the white matter fibre tracts relative to B 0 → (Lee et al., 2011, 2012) confounding the mapping of R2* to the underlying microstructure. The impact of this confounding factor can be attenuated by decomposing the angular orientation dependence of R2* into an isotropic, i.e., angular-independent component (R2,iso*), and an angular-dependent component using either complex gradient-recalled echo (GRE) acquisitions at several angular orientations (Oh et al., 2013; Wharton and Bowtell, 2013; Rudko et al., 2014) or hybrid diffusion weighted imaging (DWI) and GRE acquisitions with reduced numbers of angular orientations (Gil et al., 2016). However, both methods are impractical for clinical research due to the constrained and inconvenient positioning of the patient’s head in the radiofrequency receiver coil needed to achieve the required distinct angular orientations.

A practical approach to estimate R2,iso* was recently proposed by Papazoglou et al. (2019). They showed that R2,iso* can be estimated from the magnitude signal of a single multi-echo GRE (meGRE) measurement using a second-order model in time hereafter denoted as M2. The model was derived from a two-pool system based on the hollow cylinder fibre model (HCFM) (Wharton and Bowtell, 2012, 2013). In M2, the linear component in time (β1) is a proxy for R2,iso* and the orientation-dependent part is regressed out by the second-order term in time (β2). Although M2 is just an approximation of the original HFCM multi-compartment model and thus less accurate, it is, to our knowledge, the only way of estimating R2,iso* from magnitude-only meGRE data with a single orientation of the head.

Another advantage of the M2 is its relation to the HCFM model, allowing for direct translation of the M2-proxy for R2,iso* (i.e., the β1 parameter) into microscopic tissue properties. However, a drawback of this model is the assumption in M2 that the signal contribution of the myelin water can be neglected, limiting the microscopic interpretability of the estimated β1 parameter. For example, the M2-based prediction of β1 depends only on the transverse relaxation rate of the free water molecules of the non-myelinated compartments (R2N) and thus is independent of any changes associated with the myelin water signal or the myelin water fraction (MWF). This model’s prediction could contradict experimental observations reporting that R2* (and presumably R2,iso*) is linearly dependent on MWF (see Lee et al., 2017; Kirilina et al., 2020; Milotta et al., 2023). Moreover, M2 assumes that axonal fibres are perfectly aligned or even described by one representative axon. However, most of the fibre bundles in the human brain possess a diverse range of topographies, i.e., show fanning and bending, or mildly to acute crossing (e.g., Schmahmann et al., 2007, 2009; Jeurissen et al., 2019) and different levels of relative myelination (e.g., Mohammadi et al., 2015). Besides that, the performance of M2 in estimating R2,iso* via β1 has only been tested with data acquired at very long maximum echo times up to ≈ 54 ms (Papazoglou et al., 2019). Such a long maximum echo time is unusual for in vivo meGRE measurements with whole-brain coverage (Weiskopf et al., 2013; Ziegler et al., 2019), because it increases the total scan time as well as the propensity for bulk and physiological motion.

This work explores the potential and pitfalls of using M2 to estimate R2,iso* via β1, from a single-orientation meGRE, while varying biological fibre properties and maximum echo times. To this end, we use simulated (hereafter in silico) data and ex vivo MRI. The in silico data are simulated using a three-pool system based on the HCFM to generate realistic meGRE datasets from an ensemble of myelinated axons, for which the ground truth biophysical parameters (i.e., g-ratio, fibre dispersion and angular orientation) are known and can be varied. The ex vivo dataset combines high-resolution DWI and multi-orientation meGRE imaging of a human optic chiasm to generate gold-standard datasets where the fibre orientation and dispersion can be estimated. Both datasets are used to perform the following analyses: First, we assess the performance of M2 to estimate R2,iso* via β1 for varying g-ratio values and fibre dispersions. Second, we assess the microstructural interpretability of β1. To this end, we test the model-prediction of M2 that β1 is independent of MWF by evaluating the deviation between the biophysically predicted β1 by M2 and the fitted β1 using the in silico data. Additionally, we perform the same comparison to the fitted β1 as above using a novel heuristic expression that incorporates the MWF dependence into the predicted β1. Third, we use the heuristic expression for β1 to calculate MWF from the β1 of the ex vivo data for two sets of compartmental R2 values, i.e., in vivo and ex vivo configurations. And fourth, we assess the performance of M2 to estimate R2,iso* via β1 for different echo times ranges.

 Background Overview of the hollow cylinder fibre model and the approximated log-quadratic model

The HCFM (Wharton and Bowtell, 2012, 2013) proposes an analytical approximation describing the dependence of the GRE signal on the angular orientation ( θ μ → ) defined as the angle between the main magnetic field B 0 → and the hollow-cylinder fibre ( μ → ). This approximation establishes that the total MR signal comes from water molecules in an infinitely long and perfectly aligned hollow cylinder affected by the diamagnetic myelin sheath (Liu, 2010). The diamagnetic myelin sheath magnetically perturbs the water molecules in three distinct compartments: (1) the intra-axonal (SA), (2) myelin (SM) and (3) extra-cellular (SE) compartments (details in Supplementary material, section 2). When the signal of the water molecules in the myelin compartment is neglected (i.e., at long echo times: TE > T2 of myelin, T2M), the signal magnitude of the HCFM can be approximated by a log-quadratic model (M2) in time (Papazoglou et al., 2019):

M2 : ln ( | S N ( t , θ μ → ) | ) ≈ β 0 − β 1 t − β 2 ( θ μ → ) t 2 ,

where β 0 , 1 , 2 are the model parameters. In this model, the slope β1 is considered as a proxy for R2,iso* because it does not possess any θ μ → dependence, whereas β2 contains all the θ μ → dependent information of R2* (detailed derivation can be found in section 4, Supplementary Equations S17b,c).

Classically, R2* is estimated by the slope (α1) of the log-linear model (M1) (Elster, 1993):

M 1 : ln ( | S ( t ) | ) ≈ α 0 − α 1 ( θ μ → ) t

where α1 is a function of R2,iso* and the θ μ → dependent components of R2* (e.g., see Lee et al., 2011, 2012).

In this model, the offset parameter α 0 captures a large portion of the remaining information like contrast parameters, e.g., magnetisation transfer and longitudinal relaxation rate; and experimental parameters, e.g., transmit field. For M2, we assume that β 0 behaves in an identical fashion, i.e., this parameter captures all the remaining information as α 0 , even though this assumption is not explicitly shown in the HCFM [see Discussion in Wharton and Bowtell (2013)].

 Myelin independence of <italic>β</italic><sub>1</sub> parameter as predicted by the log-quadratic model (M2)

The slope β1 of M2, which is a proxy for R2,iso*, is derived from the HCFM by assuming a two-pool system in the slow-exchange regime: a fast decaying water pool consisting of the myelin water with a relaxation rate R 2 M and a non-myelin water pool with a relaxation rate R 2 N = R 2 A = R 2 E . In this work, we assumed that this non-myelin water pool consisted of the intra and extra cellular water, based on the findings and simplifications of Dula et al. (2010) and Wharton and Bowtell (2013). The only source of dephasing in the HCFM is caused by the magnetic properties of the hollow-cylinder fibre. All potential other perturbers are ignored (e.g., non-local effects of susceptibility inhomogeneities due to cavities, vessels, iron molecules, and diffusion) as well as other anisotropic magnetic properties, e.g., the magnetisation transfer effects (Pampel et al., 2015), influencing transverse relaxation rate (Knight et al., 2017; Birkl et al., 2021; Tax et al., 2021) or longitudinal relaxation rate (Labadie et al., 2014; Schyboll et al., 2019; Chan and Marques, 2020; Kleban et al., 2021). In white matter, this simplification seems to be reasonable since the HCFM describes the orientation dependence of the meGRE signal to a great extend (Wharton and Bowtell, 2012). Consequently, in the approximation of M2 (Eq. 1), the predicted β 1 parameter (hereafter β 1 , n m , where nm represents ‘no myelin contribution’) is given by the transverse relaxation rate of the non-myelin water pool ( R 2 N ):

β 1 ≈ β 1 , n m = R 2 N .

We hypothesise here that for realistic tissue composition where the myelin compartment cannot be neglected (i.e., g-ratio equal to or smaller than 0.8), Eq. 3 is invalid. This hypothesis is supported by previous observations showing that R2* (and presumably R2,iso*) depends on the myelin water fraction, MWF (e.g., Lee et al., 2017; Weber et al., 2020; Milotta et al., 2023).

Here, we propose an alternative heuristic biophysical expression of the predicted β1 parameter (hereafter β 1 , m , where m denotes ‘with myelin contribution’):

β 1 , m = ρ N R 2 N ( 1 − V M ) + V M ρ M R 2 M ( ρ N ( 1 − V M ) + V M ρ M ) ,

under the assumption that the proton densities of the non-myelinated compartments are equal (i.e., ρA = ρE ≡ ρN) and the volume of the non-myelinated compartment is defined as one minus the myelin compartment’s volume (= 1 – VM). The heuristic expression in Eq. 4 can be analytically derived under the condition that TE < T2 of the myelin compartment based on the HCFM (details can be found in section 4, Supplementary Equation S18).

In this case, β 1 , m can be re-written as a function of the myelin water fraction (MWF, Supplementary Equation S19a, section 4), R2N and R2M:

β 1 ≈ β 1 , m = ( 1 − MWF ) R 2 N + MWF · R 2 M .

Consequently, the MWF can be calculated by re-ordering Eq. 5 as a function of the R2’s values:

M W F = β 1 , m − R 2 N R 2 M − R 2 N .

Based on our hypothesis, we expect that the heuristic expression for β 1 , m can better describe the fitted β1 when varying the g-ratio, and thus is a better proxy of R2,iso*.

 Materials and methods

This section explains the approaches used for data acquisition, data analysis and for comparing the results obtained from the ex vivo data and the findings derived from the in silico data.

<italic>Ex-vivo</italic>: optic chiasm Sample and data acquisition

A human optic chiasm (OC) from a patient without any diagnosed neurological disorder was measured (male, 59 years, multi-organ failure, 48 h postmortem interval, ~80 days of fixation in phosphate buffered saline (PBS) pH 7.4 with 0.1% sodium acide NaN3 containing 3% paraformaldehyde +1% glutaraldehyde) with prior informed consent (Ethical approval #205/17-ek). Two MR techniques were used: multi-echo GRE (meGRE) and diffusion-weighted MRI (dMRI). The meGRE data used here have been in parts used already in Papazoglou et al. (2019).

All meGRE acquisitions were performed on a 7 T Siemens Magnetom MRI scanner (Siemens Healthcare GmbH, Erlangen, Germany) using a custom 2-channel transmit/receive circularly polarised (CP) coil with a diameter of 60 mm. The OC sample was fixed within an acrylic sphere of 60 mm diameter filled with agarose (1.5% Biozym Plaque low melting Agarose, Merck, Germany) dissolved in PBS (pH 7.4 + 0.1% sodium acide) and scanned at sixteen orientations (covering a solid angle, with azimuthal and elevation angles from 0° to 90°, Figure 1A) using the 3D meGRE MRI (hereafter: GRE dataset). For each angular meGRE measurement (Figure 1B), sixteen echoes were acquired at equally spaced echo times (TE) ranging from 3.4 to 53.5 ms (increment 3.34 ms) with a repetition time (TR) of 100 ms, a field of view (FoV) of (39.00 mm)3, a matrix size of 1123, resulting in an isotropic voxel resolution of (0.35 mm)3, non-selective RF excitation with a flip angle of 23° and a gradient readout bandwidth of 343 Hz/px.

Acquisition of the multi-angular multi-echo gradient recalled echo (meGRE) ex vivo data. (A) An illustration of the different angular measurements performed on the optic chiasm (OC) specimen. The red dots show the position of the optical tracts (see inset) for the different measurements. The different coordinates (spatial, x-y-z and anatomical, anterior-head-right, A-H-R) are shown (adapted illustration from Papazoglou et al., 2019). (B) Illustration of the first echo meGRE image acquired at the first and last angular measurement. The 3D view shows the specimen position to the main magnetic field B 0 → and the position of the optical tract (red dot). The yellow line shows the same coronal slice image.

The multi-shell dMRI data (hereafter: dMRI dataset), suitable for NODDI analysis, were acquired, three months later, on a 9.4 T small animal MR system (Bruker Biospec 94/20; Bruker Biospin, Ettlingen, Germany) using a 2-channel receiver cryogenically cooled quadrature transceiver surface RF coil (Bruker Biospin, Ettlingen, Germany) and a gradient system with Gmax = 700 mT/m per gradient axis. This dataset was acquired with a slice-selective (2D) pulsed-gradient spin-echo (PGSE) technique, consisting of four diffusion-weighting shells (number of directions) of b = 1,000 s/mm2 (60), 4,000 s/mm2 (60), 8,000 s/mm2 (60) and 12,000 s/mm2 (60) with 35 non-diffusion-weighted volumes (~ 0 s/mm2). The fixed diffusion parameters were diffusion time Δ = 13 ms, diffusion gradient duration δ = 6 ms. The remaining sequence parameters were TE = 27 ms, TR = 30 s (to acquire all the slices), FoV = 20.75 × 16.00 × 12.50 mm3, matrix size = 83 × 64 × 50, isotropic voxel resolution = (0.25 mm)3, slice selective pulses with flip angles of 90° (excitation) and 180° (refocusing) and receiver bandwidth of 9,411 Hz/px.

Note that we used different MR systems for dMRI and meGRE measurements, because it was intended to use the optimal system for the respective measurement. The dMRI dataset was acquired on a Bruker Biospec with cryo-coil to take advantage of the scanner’s high gradient strength, allowing for acquisition of high-resolution images at optimal b-values for ex vivo tissue [up to 12,000 s/mm2 (Roebroeck et al., 2018)]. Moreover, we used a TR of 30 s for dMRI measurements to ensure full magnetisation recovery and to reduce possible biases for diffusion analysis (section 3.1.2). The meGRE data was acquired on a human 7 T Siemens Magnetom MRI scanner because an optimised meGRE sequence was available on this system, including a self-built ex vivo sample coil.

Since the tissue sample was already fixed and immersed in PBS and sodium acide solution to preserve the tissue quality (Minassian and Huang, 1979), the time between the acquisitions (less than 2 weeks) did not affect the tissue quality.

 Dispersion and mean fibre orientation estimation from dMRI dataset

To incorporate the voxel-wise information regarding the angular orientation of the fibres to B 0 → and fibre’s dispersion, the dMRI datasets were corrected using a simple rigid-body registration to remove a potential drift of the sample during measurements. The dMRI data were analysed with two diffusion models: Neurite Orientation Dispersion and Density Imaging (NODDI) (Zhang et al., 2012) and Diffusion Tensor Imaging (DTI) (Basser et al., 1994). The NODDI toolbox was adjusted for ex vivo analysis (Wang et al., 2019) and used all the diffusion shells. The main neurite (hereafter fibre) orientation ( μ → ), a measure of the fibre dispersion (κ), and fibre density (volume fraction of the intracellular compartment, ICVF) maps were estimated from this analysis. The DTI model used the first two diffusion shells (b-values: 1000 s/mm2 and 4,000 s/mm2) and was used only for estimating the fractional anisotropy (FA) map, which in turn was used only for diffusion-to-GRE coregistration (section 3.1.3). Note that eddy currents were small in this dMRI data because the data were acquired with a small FoV in the gradient-coil centre using the standard Bruker gradient coil. Moreover, the cryo-coil provided sufficiently high SNR values for unbiased dMRI model parameters (the mean SNR across the specimen was approximately 57 for the b-value = 0 s/mm2 images). The SNR was calculated by dividing the MR signal by the standard deviation of the background voxels of its corresponding image (Kellman and McVeigh, 2005).

 Coregistration of the GRE angular measurements and dMRI results

To establish a voxel-to-voxel relationship between the meGRE signal at different angular orientations and the properties estimated from dMRI, i.e., κ, μ → and ICVF, we coregistered the angular meGRE measurements and the dMRI measurement. To this end, we estimated two sets of transformation matrices: first, transformation matrices that coregister the i-th angular measurements in GRE space, T G R E : i , 1 (with i = 2… 16, see Figure 2A); and second, a transformation matrix that coregisters from GRE space to dMRI space, T D i f f , G R E (see Figure 2B). The coordinate system of GRE space was defined by the first meGRE angular measurement. This reference was chosen due to the alignment of the optical nerves to B 0 → and following the procedure adopted in a previous study (Papazoglou et al., 2019). The meGRE coregistration and estimation of T G R E : i , 1 were performed using the 3D Slicer software1 (Fedorov et al., 2012), while the GRE-to-diffusion transformation was performed using the coregistration module in SPM 12.2

Coregistration of the ex vivo GRE and dMRI measurements. (A) A transformation matrix (TGRE) is obtained by coregistering all other multi-echo gradient-recall-echo (meGRE) datasets (I2.0.16) to the first measurement (I1, TGRE: i,1). This transformation matrices not only align, voxel-wise, the images of the meGRE datasets (I’2.0.16) to the first dataset, but also adjusts the directions of the main magnetic field ( B 0 → ) per angular measurement to preserves their relative orientation with respect to the first meGRE dataset. (B) A transformation matrix (TDiff,GRE) is obtained by coregistering the diffusion MRI (dMRI) image to the first angular GRE measurement. This transformation will allow the coregistration of the NODDI analysis results to the GRE data.

 Voxel-wise estimation of the angular orientation, <inline-formula> <mml:math id="M36"> <mml:mrow> <mml:msub> <mml:mi>θ</mml:mi> <mml:mover accent="true"> <mml:mi>μ</mml:mi> <mml:mo>→</mml:mo> </mml:mover> </mml:msub> </mml:mrow> </mml:math> </inline-formula>, between fibres and <inline-formula> <mml:math id="M37"> <mml:mrow> <mml:mover accent="true"> <mml:mrow> <mml:msub> <mml:mi>B</mml:mi> <mml:mn>0</mml:mn> </mml:msub> </mml:mrow> <mml:mo stretchy="true">→</mml:mo> </mml:mover> </mml:mrow> </mml:math> </inline-formula>

The angular orientation θ μ → between fibres and B 0 → for each meGRE angular measurement was calculated in dMRI space and the resulting θ μ → maps mapped onto GRE space. For that, the arccosine of the inner product between B 0 → ( θ i ) and μ → , i.e., θ μ → = a r c c o s B 0 → θ i • μ → was computed (Figure 3C). In this computation, B 0 → ( θ i ) is the resulting B 0 → after the transformation from the i-th meGRE angular measurement to the first meGRE angular measurement ( T G R E : i , 1 ), and the transformation from GRE to dMRI space ( T D i f f , G R E − 1 ) (Figure 3A). The main fibre direction was obtained by the μ → map from the NODDI analysis (Figure 3B).

Estimation of the voxel-wise angular θ μ → map. This estimation needed the B0 direction per angular GRE measurement ( B 0 → ( θ i ) ) in diffusion space and the main fibre direction. (A) The B 0 → ( θ i ) was estimated by applying to B 0 → , first, the transformation matrix between GRE volumes (TGRE:i,1) and later from GRE-to-diffusion (T−1Diff,GRE). (B) The main fibre direction ( μ → ) was acquired by analysing the dMRI data with the NODDI model. (C) Then, the voxel-wise θ μ → per angular measurement was computed by the arccosine of the scalar product between the projected B 0 → ( θ i ) and the main diffusion direction ( μ → ), θ μ → = a c o s B 0 → θ i • μ → . This sketch shows the steps for the last GRE angular measurement.

Note that θ μ → was computed in dMRI space instead of GRE space to avoid undersampling and interpolation caused by transforming the dMRI-based μ → to GRE space. These sources of error do not occur by transforming B 0 → to dMRI space, i.e., computing T D i f f , G R E − 1 · T G R E : i , 1 · B 0 → , for each GRE angular measurement, since it is a global rather than a per-voxel measure. Finally, the θ μ → maps together with the ICVF and κ maps (not shown in Figure 3) were transformed using T D i f f , G R E . Exemplary θ μ → maps in GRE space are shown in Supplementary Figure S1 (first row).

 Masking and pooling the <italic>ex vivo</italic> data

Before analysis, the ex vivo data required further pre-processing to remove outliers and to ensure a robust assessment of the effect of fibre dispersion and θ μ → on R2*. For that, the ex vivo data were masked using the coregistered ICVF map and later pooled across the sixteen coregistered meGRE angular measurements.

In this process, all voxels in the OC with an ICVF >0.8 were selected and pooled across all the meGRE angular measurements, hereafter referred to as cumulated data. The ICVF threshold was used because the extra-axonal space in the ex vivo specimen is reduced (e.g., Stikov et al., 2011). The application of this threshold reduced the number of voxels in the OC by 7.2% (~ 600 over 8,744 voxels). By pooling the data, the resulting cumulated data dependent on TE but also on θ μ → from 0° to 90°, and on fibre dispersion assessed by κ from 0 to 6.

 Simulated R<sub>2</sub>* signal decay from the HCFM

Multi-echo GRE signal decay was simulated as ground truth (hereafter, in silico data) to assess the impact on M2 of variable fibre orientation, dispersion and myelination (i.e., g-ratio). For that, we estimated averaged MR signals calculated from an ensemble of 1,500 hollow cylinders. The cylinders were evenly distributed on a sphere with defined spherical coordinates: an azimuthal angle φ rotating counter-clockwise from 0° to 359° starting aligned with the +X axis, and elevation angle θ rotating from 0° (+Z) to 180° (−Z). The signal contribution per hollow cylinder was modelled using the frequency-averaged signal equations from the HCFM for all the compartments including the myelin compartment (Supplementary Equations S1, S3, section 2).

In the simulation framework, three assumptions were made. First, the B 0 → was fixed and oriented parallel to +Z (Figure 4A). Second, the approximated piece-wise function of DE in the SE signal was replaced by its analytical solution (Supplementary Equations S2b, S4, S7, S8, section 3; respectively), because a discontinuity in this piece-wise function was observed at the so-called critical time (Yablonskiy and Haacke, 1994; Wharton and Bowtell, 2013). See section 3 in Supplementary material for a detailed discussion. And third, we ignored the near-field effects between cylinders, therefore the total signal is the sum of all the complex signals from each cylinder as defined in Supplementary Equations S1–S3.

Schematics of the simulated in silico data: (A) Simulation: 1500 hollow cylinders, each of them defined by the vector x i → , were distributed evenly on a sphere (see the blue dots). A mean orientation μ → of the cylinders is defined, with the external magnetic field ( B 0 → ) oriented parallel to the Z-axis. The signal contribution per cylinder was modelled using the Hollow Cylinder Fibre Model (HCFM) with the intra-axonal (SA), extra-axonal (SE) and myelin (SM) compartments (inset). (B) Addition of cylinder’s dispersion: the dispersion effect was added by weighting the signal coming from the cylinders by the parameter κ from the Watson distribution and μ → (Eq. 8b). The parameter κ is limited from κ = 0 for isotropically dispersed to κ = infinity to fully parallel fibres. Here, μ → is parallel to B 0 → .

To incorporate the effect of fibre dispersion in the in silico data, the ensemble-average signal was calculated by weighting Sc with the Watson distribution (W) (Sra and Karp, 2013 and Eq. 8b). This weight from the Watson distribution was calculated using the position of each simulated cylinder, x i → , and a mean fibre orientation μ → , both defined with spherical coordinates (φ, θ) and ( ϕ μ → , θ μ → ), respectively (Figure 4A). For simplification, μ → was restricted to an azimuthal angle of zero ( ϕ μ → = 0°). Then, the analytical expression of the ensemble-average signal, SW, is defined as follows:

S W ( κ , t , θ μ → ) = ∑ i ( S C ( t , θ i ) [ W ( κ , θ μ → , ϕ i , θ i ) ] ) ∑ i W ( κ , θ μ → , ϕ i , θ i ) , where W ( κ , θ μ → , θ i , ϕ i ) = C 1 ( 1 2 , 3 2 , κ ) − 1 e κ ( μ → ( θ μ → ) · x i → ( ϕ i , θ i ) ) 2 .

In Eq. 7b, C 1 () is the confluent hypergeometric function, which is the normalisation factor of the Watson distribution, and the exponent holds the norm of the inner product between each individual i-th cylinder x i → and μ → . The level of dispersion was modulated by the parameter κ (Zhang et al., 2012; Sra and Karp, 2013) as shown in Figure 4B for a few cases. It is important to note that the notation θ μ → for the elevation angle of μ → used here is equal to the one used to describe the fibre’s angular orientation in the ex vivo data (section 3.1.4). This is intentional since they stand for the same concept in both datasets. This simulation approach was used in previous conference publications [Fritz et al. (2020, 2021)].

With the ensemble averaged signal equation (Eq. 7a), a meGRE signal decay can be computed based on the relevant parameters listed in the following Table 1.

We tested the validity of M2 and its microstructural derivatives like MWF (Eq 6) for varying microstructural parameter settings. This included two different sets of compartmental R2 values (i.e., R2N and R2M) (Dula et al., 2010; Wharton and Bowtell, 2013). The sets represent two possible extremes of compartmental R2 values at 7 T: (1) ex vivo compartmental R2 values measured from an ex vivo rat spinal cord with similar tissue preparation procedure as the optic chiasm in this work [i.e., fixed with 4% PFA and hydrated in PBS (Dula et al., 2010)], and (2) in vivo compartmental R2 values obtained from in vivo human measurements (Wharton and Bowtell, 2013). The ex vivo compartmental R2 values are reported in the main manuscript whereas the in vivo compartmental R2 values are reported in the Supplementary material, Section 7.

Finally, each simulated meGRE signal decay was replicated 5,000 times with an additive Gaussian complex noise (Gudbjartsson and Patz, 1995) to approximate the SNR of the experimental ex vivo data (see Supplementary material, section 5). The experimental SNR was calculated by dividing the MR signal acquired at the first echo by the standard deviation of the background voxels of its corresponding image (Kellman and McVeigh, 2005), resulting in a mean SNR across the selected voxels of the OC of 112.

This simulation framework is publicly and freely available in Github.3

 Data analysis Data fitting and binning

The ex vivo data (section 3.1) and in silico data (each of 5,000 replicas per simulated meGRE signal decay, section 3.2) were analysed with the log-linear and log-quadratic models, M1 (Eq. 2) and M2 (Eq. 1), respectively. In both models, the α’s (α0 in arbitrary units, α1 in units of 1/s) from M1, and β’s (β0 in arbitrary units, β1 in units of 1/s and β2 in units of 1/s2) from M2, hereafter referred to as the α-parameters and β-parameters, were estimated. To fit the data, ordinary Least Square (OLS) optimization was used for both models in custom-made Matlab code. Three fittings were performed using three different meGRE subsets, that varied by their maximum TE (TEmax) values: TEmax = 54 ms (all 16 time points), TEmax = 36 ms (first 10 points) and TEmax = 18 ms (first 5 time points). The first meGRE subset with TEmax of 54 ms replicated the meGRE protocols of the ex vivo studies, while the meGRE subset with TEmax of 18 ms could be considered as a typical meGRE protocol for in vivo studies [at least with regards to the sample size and TE range used in the multi-parametric mapping protocol (Weiskopf et al., 2013)]. The meGRE subset with TEmax of 36 ms was chosen as an intermediate subset between both protocols.

To compare the α- and β-parameters between datasets as a function of fibre dispersion (κ) and θ μ → , the fitted parameters were binned and averaged for the ex vivo cumulated data (section 3.1.5) and for the in silico data. The binning on the fitted parameters was performed to ensure: (1) a reduced effect size bias in the ex vivo cumulated data, given the unequal number of voxels at specific θ μ → and κ (Figure 5); and (2) a better comparison between in silico and ex vivo data.

Preparation of the ex vivo data for analysis. (A) The cumulated ex vivo data were distributed first as a function of κ parameter, to ensure similar fibre dispersion. Heuristically it was divided in highly dispersed (κ < 1), mildly dispersed (1 ≤ κ < 2.5) and negligibly dispersed (κ ≥ 2.5) fibres. Coincidentally, this division enclosed specific areas in the OC (red, green and blue ROIs). (B) After division, the cumulated data were binned irregularly as a function of the estimated voxel-wise angular orientation ( θ μ → ) per κ range (orange bars), to avoid a possible effect size bias caused by its non-uniform distribution (blue bars). The first angular irregular bin or angular offset θ 0 was obtained and showed to be κ range dependent (Supplementary Table S1, section 5).

In the binning process, both datasets were distributed first as a function of κ, and later as a function of θ μ → . The first distribution was performed to ensure a similar degree of fibre dispersion as observed in Figure 4B and in the work of Fritz et al. (2020). For that, three different fibre dispersion ranges were defined as a function of κ: κ < 1 for the highly dispersed fibres, 1 ≤ κ < 2.5 for the mildly dispersed fibres, and κ ≥ 2.5 for the negligibly dispersed fibres. Coincidentally, these fibre dispersion ranges depicted specific areas in the OC (Figure 5A). However, the in silico data required two extra averages on the fitted parameters to bin it as a function of the different fibre dispersion ranges: first, across the 5,000 replicas and, second, across the κ values within each fibre-dispersion range. The average across κ was performed in such a way that it resembled the frequency distribution of κ observed in the ex vivo cumulated data (for more detail, see Supplementary material, section 5). After separating the fitted parameters per fibre dispersion range for both datasets, the data were irregularly binned as a function of θ μ → bins per defined κ range.

The irregular θ μ → bins were introduced to avoid a bias due to the uneven distribution of voxels with azimuthal orientations across the 16 angular measurements (Figure 5B, blue bars). To determine the irregular θ μ → bin sizes, a cumulated θ μ → distribution of voxels was estimated and divided into 20 equally populated bins (Figure 5B, orange bars). The mean of the first angular irregular bin was defined as the angular offset θ 0 . The range of θ μ → values contained in each irregular bin and the θ 0 values are shown in Supplementary Table S1, section 5.

After binning, the average and standard deviation (sd) of the α- and β-parameters between datasets was calculated per irregular θ μ → bin in the ex vivo cumulated data. For the in silico data, the average and sd of the same parameters were obtained by weighting the distribution of θ μ → in each bin in a similar way to that seen in the irregular bins in the ex vivo cumulated data (for more detail, see Supplementary material, section 5).

 Quantitative analysis

Four different analyses were performed in order to study: (1) the effect of g-ratio and fibre dispersion, via κ, on the estimated angular-independent β1 parameter in M2, (2) the microstructural interpretability of β1 via the deviation between fitted β1 and its predicted counterparts from M2 (β1,nm, Eq. 3) and from the heuristic expression (β1,m, Eq. 4), (3) the possibility of calculating the MWF (Eq. 5) from the fitted β1 using the heuristic expression β1,m, and (4) the effect of TE, via the different meGRE subsets, on the performance of M2. The last analysis was divided into two parts, testing: (A) its capability to reduce the orientation dependence in β1 (and thus be a valid proxy for R2,iso*), and (B) if M2 can be better explained by the different meGRE subsets than M1. Using the simulation framework, the validity of M2 and its derived microstructural parameters were tested based on analyses 1, 2 and 4. While both datasets were used for the first and fourth analyses, only the in silico data were used for the second analysis while ex vivo data were only used for the third analysis.

First analysis: ability of M2 to obtain the angular-independent β<sub>1</sub> parameter for varying g-ratio and fibre dispersion values

For the first analysis, the ability of M2 to estimate an orientation-independent effective transverse relaxation rate, R2,iso*, via the β1 parameter was assessed. Since R2,iso* by definition is the angular independent part of R2 * and according to the HCFM should be given by β1 parameter at θ μ → = 0 ≡ θ 0 , we assessed the residual θ μ → dependence of the β1 parameter with respect to θ 0 and compared it with its counterpart for α1, i.e., the proxy for the θ μ → dependent R2*.

For this, we first calculated the θ μ → dependence of each parameter with respect to θ 0 using the normalised-root-mean-squared deviation (nRMSD, in %):

n R M S D γ κ j = ∑ l = 1 N − 1 γ κ j , θ 0 − γ κ j , θ l 2 N − 1 γ κ j , θ 0 · 100 % ,

where θ 0 varied slightly for each κ range (sub-index j) but was close to zero (see Supplementary Table S1) with 𝛾 ∈ {𝛼1, 𝛽1}.

To compare the nRMSD of each parameter, we calculated the difference between them, ΔnRMSD, as:

Δ n R M S D ( κ j ) = n R M S D ( β 1 ( κ j ) ) − n R M S D ( α 1 ( κ j ) )

in percentage-points (%-points). If the ΔnRMSD is positive or higher than 0%-points, this implies that the θ μ → dependency of β1 is similar or higher, in magnitude, to α1. The latter says therefore that M2 failed in estimating an angular-independent parameter from R2*. A negative ΔnRMSD in turn implies that the θ μ → independence of β 1 ( κ j ) has been reduced. A perfect orientation independence is achieved if n R M S D ( β 1 ( κ j ) ) = 0 and, consequently, Δ n R M S D ( κ j ) = − n R M S D ( α 1 ( κ j ) ) .

 Second analysis: assessment of the microstructural interpretability of β<sub>1</sub>

For the second analysis, the microstructural interpretation of β1 was quantitatively assessed by comparing the relative difference (ε) between estimated β1 at the angular orientation θ μ → for the fitted in silico data ( β 1 ( θ μ → ) ) and the predicted β1 ( β 1 , p ) using M2 or the heuristic expression:

ϵ ( θ μ → , κ j ) p = ( 1 − β 1 ( θ μ → , κ j ) β 1 , p ) · 100 % ,

where β 1 , p ∈ { β 1 , n m , β 1 , m } were defined in Eqs 5, 6, respectively. Additionally, the mean ϵ ( θ μ → , κ j ) p across angles was calculated as ϵ κ j p ≡ 1 N ∑ l = 1 N ϵ θ l , κ j p .

 Third analysis: myelin water fraction and g-ratio estimation from ex vivo data using the heuristic expression of R<sub>2,iso</sub>* via β<sub>1,m</sub>

For the third analysis, the MWF was estimated from the fitted β1 in ex vivo data using the analytical expression for β1,m (Eq. 6). For that, we used the two sets of R2 values for the non-myelinated (R2N) and myelinated (R2M) compartments reported in Table 1. Only the ex vivo R2 values were reported in this section, while the in vivo R2 values were reported in Supplementary material, section 7.2.3.

Microstructural parameters used to generate the in silico data.

Parameter Value Reference
Anisotropic and isotropic susceptibilities (χA and χI) −0.1 ppm Wharton and Bowtell (2013)
Exchange (E) 0.02 ppm Wharton and Bowtell (2013)
Proton density intra-and extra- axonal compartments (ρA and ρE)* 5,000 a. u. Wharton and Bowtell (2013)
Larmor frequency at 7 T (ω0) 1.873 ∙ 106 rad/ms –
Fibre volume fraction (FVF) 0.5 n. u. Wharton and Bowtell (2013)
Proton density myelin compartment (ρM)* 3,500 a. u. Wharton and Bowtell (2013)
R2 intra-and extra- axonal compartments (R2A = R2E = R2N) 18.53 s−1 (ex vivo)27.8 s−1 (in vivo) Dula et al. (2010) Wharton and Bowtell (2013)
R2 myelin compartment (R2M) 75.41 s−1 (ex vivo)125 s−1 (in vivo) Dula et al. (2010) Wharton and Bowtell (2013)
Angular orientation ( θ μ → ) 2°:2°:90° –
Index of fibre dispersion (κ) 0.001:0.1:6.0 –
g-ratio 0.66, 0.73, 0.8 Emmenegger et al. (2021) only for 0.66 and Wharton and Bowtell (2013) for 0.8.**
Time (i.e., echo time) 3.25:3.25:53.5 ms –

The last four parameters are the parameter or simulation space. *Proton densities were scaled by a factor of 5,000 but they kept the same proton density proportion between the non-myelinated and myelinated compartments (1:0.7). ** The mean g-ratio value of 0.73 was arbitrarily defined.

 Fourth analysis: the effect of echo time ranges on the performance of M2

In the fourth analysis, the performance of M2 was tested in two sub-analyses when using the meGRE datasets with different TE ranges (see section 3.3.1).

First sub-analysis: assessing the residual <inline-formula> <mml:math id="M123"> <mml:mrow> <mml:msub> <mml:mi>θ</mml:mi> <mml:mover accent="true"> <mml:mi>μ</mml:mi> <mml:mo>→</mml:mo> </mml:mover> </mml:msub> </mml:mrow> </mml:math> </inline-formula> dependence in β<sub>1</sub> for meGRE subsets with different maximum echo times

For the first sub-analysis, the orientation dependence of β1 was assessed for the different meGRE subsets from the ex vivo dataset and the in silico data for variable g-ratio. For that, α1 and β1 from M1 and M2 were compared once again as in the first analysis and the ΔnRMSD was calculated to assess the residual θ μ → dependence of β1 in comparison to the θ μ → dependence of α1.

 Second sub-analysis: assessing if M2 is better explained by the data using meGRE subsets with different maximum echo times

For the second sub-analysis, the weighted-corrected Akaike Information Criterion (wAICc, Eq. 12) was introduced [more details can be found in Supplementary Equation S30, section 6 and Burnham et al. (2011)]. According to Burnham et al. (2011), the wAICc can be used to assess whether a given model (here M2) is better explained [or “supported” as introduced in Burnham et al. (2011)] by the data than a set of other models (here M1). In this work, we used the AICc (i.e., Akaike Information Criterion, AIC, with a correction for small sample sizes) instead of the AIC or the Bayesian Information Criterion (BIC) to better account for the small sample size in comparison to the number of model’s parameters. Note that to use the AIC the ratio between the sample size and the number of parameters (n/k) should be above 40 (Burnham and Anderson, 2002) and this condition was not always fulfilled in our data.

The wAICc for M2 is defined by:

wAICc = 1 1 + exp − 0.5 Δ AIC ,

where ΔAICc in Eq. 12 is the difference of the AICc for models M1 and M2:

Δ AICc = AICc ( M 1 ) − AICc ( M 2 )

The AICc and wAICc were estimated per voxel (for ex vivo) and replica (for in silico) from the previous analysis using the sum-of-squares error (SSE) from the fitting of each model (see Supplementary Equation S32, section 6) and for each of the three meGRE subsets. Note that the AICc and wAICc were estimated only ex vivo and in silico data with negligible fibre dispersion (κ ≥ 2.5). Then, the averaged wAICc as well as its standard deviation (sd) were calculated. In this work, we interpreted the range of possible wAICc values in a more conservative manner. Hereby, we mainly focused on the case AICc(M1) > AICc(M2) (Eq. 12), where the resulting wAICc (Eq. 11) is greater than 0.5: a wAICc >0.73 implies that M2 is better explained by the meGRE data than M1, and a wAICc between 0.5 and 0.73 implies that M2 and M1 are ambiguously explained by the data but M2 is still preferred. For the case of AICc(M1) ≤ AICc(M2), where wAICc ≤0.5, M2 was not explained by the data as compared to M1. More details regarding the calculations as well as the threshold of 0.73 can be found in the Supplementary material, section 6. Note that we were only reporting the average wAICc, thus the wAICc for some voxels (for ex vivo) or replicas (for in silico) might belong to a different range than the average wAICc, which can be observed by the estimated sd wAICc.

In the following sections, the dependence of the parameters under study, i.e., nRMSD( α ( κ j ) ), nRMSD( β ( κ j ) ), ΔnRMSD ( κ j ), α1( θ μ → , κ j ), β1( θ μ → , κ j ) (Eqs 8, 9), ϵ ( θ μ → , κ j ) p (Eq. 10) and ϵ ( κ j ) p , to θ μ → and κ were simplified for readability purposes. Therefore, these parameters will be hereafter nRMSD(α1), nRMSD(β1), ΔnRMSD, α1, β1, ϵ p and ϵ p , respectively.

 Results First analysis: ability of M2 to obtain the angular-independent <italic>β</italic><sub>1</sub> parameter for varying g-ratio and fibre dispersion values

Figure 6 shows the performance of M2 when estimating R2,iso* via β1 for variable g-ratio and fibre dispersion. To visualise this, we compared the θ μ → dependence of α1 from M1 to the residual θ μ → dependence of β1 from M2 (Figures 6A,B). Both θ μ → dependencies were quantified in Figure 6C using their respective nRMSD (Eq. 8). The results are from the analysis performed on the ex vivo and in silico data. The in silico data was generated using the ex vivo compartmental R2 values (the corresponding results for the in vivo compartmental R2 values are presented in Supplementary Figure S4).

Orientation dependence of linear model parameters (α1 and β1) for varying g-ratio and fibre dispersion values. (A,B) Depicted is the α1 parameter of M1 (proxy for R2*) and β1 parameter of M2 (proxy for the isotropic part of R2*) as a function of the angle between the main magnetic field and the fibre orientation ( θ μ → ) for different fibre dispersion and g-ratio values. The different columns depict different dispersion regimes: highly dispersed (κ < 1, first column), mildly dispersed (1 ≤ κ < 2.5, second column) and negligibly dispersed (κ ≥ 2.5, third column) fibres. Note that the smallest angle ( θ 0 ) varied across dispersion regimes: 17.3° (κ < 1), 20.4° (1 ≤ κ < 2.5) and 22.9° (2.5 ≤ κ). This was caused by the irregular binning (see section 3.1.4) (C) Depicted is the normalised root-mean-squared deviation (nRMSD, Eq. 11 in %) of the α1 parameter of M1 (proxy for R2*) and β1 parameter of M2 (proxy for the isotropic part of R2*) for different fibre dispersion and g-ratio values. Across the entire figure, the distinct colours (blue and green curves and bars) distinguish between in silico data with variable g-ratios (increasing blue hue with increasing g-ratio) and ex vivo data (olive curve).

The ability of M2 to reduce the θ μ → dependency of β1 varied with g-ratio and fibre dispersion. The θ μ → dependency of α1 (and residual θ μ → dependency of β1) was also strongly influenced by g-ratio and fibre dispersion: smaller g-ratio values and reduced fibre dispersion increased the θ μ → dependency of α1 and (the residual θ μ → dependency) of β1 (Figures 6A,B, respectively).

The fibre dispersion affected the performance of M2 the same between in silico and ex vivo datasets (Figure 6C). In both datasets, the improvement is largest for negligible dispersion (starting from ΔnRMSD = −12.0%-points for the in silico data with a g-ratio of 0.8 and ΔnRMSD = −37.4%-points for the ex vivo data). For the ex vivo data, the nRMSD(β1) was the lowest for the negligibly dispersed fibres (nRMSD(β1): 1.3% at κ ≥ 2.5). For the in silico data, the nRMSD(β1) was the lowest for the highly dispersed fibres and for a g-ratio of 0.73 (nRMSD(β1): 0.1%), and it increased with decreasing fibre dispersion (nRMSD(β1) up to 2.7%). For the g-ratios of 0.66 and 0.8, the nRMSD(β1) was higher but still below 12%.

The θ μ → dependence of α1 on fibre dispersion was the same between in silico and ex vivo datasets (Figure 6C, top): the lower the dispersion the higher the nRMSD(α1). The θ μ → dependence of α1 increased as the g-ratio decreased.

 Second analysis: assessment of the microstructural interpretability of <italic>β</italic><sub>1</sub>

Figures 7A,B report the angular-orientation ( θ μ → ) dependent relative differences ( ϵ n m and ϵ m , Eq. 10) between the fitted β1 from the in silico data and its predicted counterparts using M2 (Eq. 3) and the heuristic expression (Eq. 4). Figure 7C shows the mean and standard deviation of ϵ n m and ϵ m across angles for ex vivo compartmental R2 values (the corresponding results for the in vivo R2 values are presented in Supplementary Figure S5).

Assessment of the microstructural interpretability of β1 by the deviation between fitted and biophysically predicted β1. (A–B) The relative difference (ε, Equation 10) was calculated between the fitted β1 to the in silico data and two biophysically-modelled expressions for β1 based on the HCFM. The two expressions for β1 values were calculated from the original expression for M2, β1,nm (Equation 3, resulting in εnm, A) and the heuristic expression, β1,m (Equation 4, resulting in εm, B). This was calculated per g-ratio and fibre dispersion. (C) The corresponding mean, <ε>, and standard deviation, sd(ε), of the relative differences across the angular orientations ( θ μ → ) were estimated. The hue intensity coding represents increasing g-ratio value for both error estimations.

ϵ n m was large, between −100% and − 40%, and varied strongly with g-ratio and fibre dispersion. ϵ n m showed the largest θ μ → dependence where the largest deviation was observed (i.e., for the g-ratio of 0.66 and the lowest fibre dispersion, Figure 7A). ϵ m was always smaller than ϵ n m and showed a smaller θ μ → dependence across all the studied fibre dispersions and g-ratios. It varied between −20 and 20% and had the largest values and variation for the smallest g-ratio and negligibly fibre dispersion. For the average across angles, we found that negligibly dispersed fibres showed the smallest ϵ n m and ϵ m per g-ratio.

The mean across angles for ϵ n m , ϵ n m , was up to −85% whereas the mean across angles for ϵ m , ϵ m , was only up to −12% (Figure 7C). On average across all g-ratios and fibre dispersion arrangements, ϵ n m was approximately 8 to 9 times larger than ϵ m . Both relative mean differences became more negative with increasing g-ratio and decreasing fibre dispersion. The ϵ m for the negligibly dispersed fibres at g-ratio 0.66 was close to −2% but accompanied by a large standard deviation across θ μ → due to the strong θ μ → -dependency of the corresponding fitted β 1 parameters. For both ϵ n m and ϵ m , the variability (Figure 7C) across different θ μ → values, s d ( ϵ n m ) and s d ( ϵ m ) respectively, was highest when the fibre dispersion and g-ratio were lowest.

 Third analysis: myelin water fraction estimation from <italic>ex vivo</italic> data using the heuristic expression of R<sub>2,iso</sub>* via <italic>β</italic><sub>1,m</sub>

Figure 8 reports the MWF estimated from the ex vivo data by inverting the heuristic expression for β1,m (Eq. 6), using the ex vivo compartmental R2 values (the corresponding results for the in vivo R2 values are presented in Supplementary Figure S6). Figure 8A shows the estimated MWF as a function of θ μ → while Figure 8B shows the median and standard deviation (sd) of the estimated MWF across θ μ → .

Dependence of the MWF estimation on angular orientation for three different fibre dispersion ranges in ex vivo data. (A) The MWF was estimated by using the heuristic analytical expression of β1 (β1,m, Eq. 4) and the fitted β1 for the ex vivo data using the compartmental R2 values from Dula et al. (2010) (hues of green) in Table 1. This calculation was performed per angle ( θ μ → ) and for the three different fibre dispersion ranges: highly dispersed, mildly dispersed and negligibly dispersed. The increasing green hue represents decreasing fibre dispersion. (B) The corresponding median and standard deviation (sd) were estimated across θ μ → per fibre dispersion range.

The estimated MWF was larger with decreasing fibre dispersion (Figure 8A). Moreover, there was a trend towards larger estimated MWF for larger θ μ → . Across θ μ → , the estimated median ex vivo MWF was 0.14 for fibres with negligible dispersion but moved towards to even lower and unrealistically small values (MWF: 0.069) for dispersed fibres (Figure 8B). The standard deviation across MWF was similar for different fibre dispersions, ranging from 0.0068 to 0.0104.

 Fourth analysis: the effect of echo time ranges on the performance of M2

In this section, two sub-analyses were performed for in silico data at variable g-ratio and ex vivo data, both with negligibly dispersed fibres (i.e., κ ≥ 2.5), using the three meGRE subsets with different maximum echo time (TEmax) for ex vivo compartmental R2 values (the corresponding results for the in vivo R2 values are presented in Supplementary Figures S7, S8). In the first sub-analysis, its result is depicted similarly as in Figure 6, but for different TEmax and κ ≥ 2.5. In the second sub-analysis, it was assessed whether M2 was better explained by the different meGRE subsets than M1 using the average wAICc of M2 (Eq. 11).

First sub-analysis: assessing the residual <inline-formula> <mml:math id="M186"> <mml:mrow> <mml:msub> <mml:mi>θ</mml:mi> <mml:mover accent="true"> <mml:mi>μ</mml:mi> <mml:mo>→</mml:mo> </mml:mover> </mml:msub> </mml:mrow> </mml:math> </inline-formula> dependence in β<sub>1</sub> for meGRE subsets with different maximum echo times

Using the meGRE subsets with smaller TEmax (36 ms and 18 ms), M2 was less effective across all g-ratios (Figures 9A,B, second and third column). For some microstructural parameter settings, even an increased θ μ → dependence was observed for β1 compared to α1: nRMSD(β1) went up by 5.6%-points at 36 ms (in silico, g-ratio: 0.8) and by 14.1%-points at 18 ms (ex vivo). Moreover, for the meGRE subset with the smallest TEmax (18 ms), an atypical θ μ → dependence of β1 (and α1) was found in the ex vivo data: β1 (and α1) decreased with increasing θ μ → up to approximately 55° (magic angle, dashed magenta lines in Figures 9A,B) and then slightly increased again. The θ μ → dependence up to the magic angle was not observed in the in silico data at any investigated meGRE subset. Moreover, the θ μ → dependence of α1 in the ex vivo data decreased when meGRE subsets with decreasing TEmax were used. This trend was mostly also observable in the in silico data (Figure 9A). Note that we investigated the orientation dependence of α1 and β1 also for mildly and highly dispersed fibres but did not find new trends in those datasets (data not shown).

Effect of the maximal echo time, i.e., meGRE subsets with different maximum echo times, on the θ μ → dependency of α1 and β1. (A,B) Angular orientation ( θ μ → ) dependence of α1 in M1 and β1 in M2 for the three meGRE subsets with varying maximum TE (TEmax: 54 ms, 36 ms and 18 ms). Two datasets are compared: ex vivo (green curve) and in silico (blue curve) data at variable g-ratios. Only datasets of the negligibly dispersed fibres (κ ≥ 2.5) are presented. The magenta vertical lines in some of the subplots indicates the magic angle ( θ μ → = 55°). (C) Depicted is the normalised root-mean-squared deviation (nRMSD, Eq. 8 in %) of the α1 parameter of M1 (proxy for R2*) and β1 parameter of M2 (proxy for the isotropic part of R2*) shown in (A) and (B), respectively.

 Second sub-analysis: assessing if M2 is better explained by the data using meGRE subsets with different maximum echo times

The average wAICc showed different trends across the different meGRE subsets with varying TEmax for both datasets. For the ex vivo data, the average wAICc decreased when meGRE subsets with smaller TEmax were used. Using the meGRE subsets with the largest and intermediate TEmax (54 and 36 ms), the average wAICc indicated that M2 was better explained than M1 by the data with wAICc values in the ranges of wAICc >0.73 (TEmax = 54 ms) and 0.73 > wAICc >0.5 (TEmax = 36 ms), respectively. Interestingly, for the in silico data, the average wAICc decreased as a function of g-ratio for the meGRE subset with the largest TEmax, from wAICc: 0.71 to 0.44; but increased with increasing g-ratio for the meGRE subset with intermediate TEmax, from wAICc: 0.31 to 0.59. However, none of the highest wAICc overpassed the threshold of 0.73. Note that the large standard deviation of the reported wAICc per dataset indicates that the results are only valid on average whereas the wAICc for single voxels (ex vivo data) or replicas (in silico data) can be outside the reported ranges (see Figure 10).

Assessing if model M2 is better explained by the meGRE signal decay than M1, quantified by the averaged wAICc for M2 (Eq. 11). This quantification was done per meGRE subsets with different maximum echo time (TEmax) for the in silico data at variable g-ratios (increased blue hue in bars, higher g-ratio) with R2 values from Dula et al. (2010); and ex vivo data (green bar) for negligibly dispersed fibres (κ > 2.5). The magenta and orange lines mark the following ranges: over the magenta line (wAICc = 0.73), M2 is better explained by the data; between the magenta and orange (wAICc = 0.5) lines, there is a preference for M2 but it is ambiguous whether M2 is better explained than M1 by the data; and bellow the orange line, M2 is not better explained by the data than M1.

 Discussion

This work quantitatively evaluated the performance of the log-quadratic model (M2) for estimating the orientation-independent part of R2* (R2,iso*) via its linear parameter, β1, using a single-orientation multi-echo GRE (meGRE) measurement in simulations and in a human optic chiasm. We found that M2 can estimate R2,iso* via β1 when using meGRE with long maximum echo time (TEmax ≈ 54 ms) for all investigated fibre dispersion and g-ratios. Our simulation results show that the proposed heuristic expression for β1 better explained the fitted β1 for ex vivo compartmental R2 values than the M2-based prediction. Using this heuristic model, we estimated realistic MWF values from β1 fitted to the ex vivo data. However, we found that its validity depends on the choice of compartmental R2-values and we found that the heuristic model cannot be used for tissue with dispersed fibres. We created an openly available simulation framework to test the validity of the heuristic expression for different microstructural arrangements. We found that M2 cannot reduce the orientation dependence of β1, and therefore cannot be used as a proxy of R2,iso* when the meGRE subsets with shorter maximum echo times were used (TEmax ≈ 36 ms or 18 ms). For the meGRE subset with the shortest TEmax of 18 ms, we found that the orientation-dependence of the classical R2* showed the highest deviation between ex vivo and in silico data for angles below the magic angle (55°), indicating that, at short echo times, the mechanism for the orientation-dependence of R2* is not captured by our HCFM-based simulation.

Ability of M2 to estimate the angular independent <italic>β</italic><sub>1</sub> for varying g-ratio and fibre dispersion values

Our results show that M2 has the potential to estimate R2,iso* from a single-orientation meGRE via β1 for the ex vivo data of an optic chiasm tissue sample and the in silico data. We found that the performance of M2, assessed by the residual θ μ → dependence of β1, varied for different g-ratios and fibre dispersions (Figure 5 and Supplementary Figure S4). For the ex vivo compartmental R2 values (Figure 5), the residual θ μ → dependence of β1 was always less than 12% even if the θ μ → dependence of the original R2* (using the α1 parameter of M1) was up to 50%. For the in vivo compartmental R2 values (Supplementary Figure S4), the residual θ μ → dependence of β1 was always less than 20%. The comparison of the performance of M2 for different compartmental R2 values indicates that the performance of M2 might vary for tissue with different microstructural tissue properties such as the compartmental R2 values or the fibre volume fraction.

 Assessment of the microstructural interpretability of <italic>β</italic><sub>1</sub>

As hypothesised in the introduction, the fitted β1 parameter is an unsuitable proxy for estimating microscopic tissue parameters via the dependency of M2 on the biophysical HCFM (Eq. 3). Using the ex vivo compartmental R2 values to generate the in silico data, we obtained an error of up to −70% for the fibres with negligible dispersion (Figure 7C) between the fitted β1 and the β1 predicted using the biophysical relation in M2 (Eq. 3). With the proposed heuristic expression for β1 (Eq. 4), the relative error was reduced by a factor of about 10 and more for fibres with negligible dispersion (e.g., from −65% to −6% for a g-ratio of 0.73, Figure 7C), indicating that this expression is better suited for the biophysical interpretation of β1 than the M2-based expression. However, we also found that the heuristic expression is not valid for all microstructure parameters, e.g., for in vivo compartmental R2 values the error switched signed, e.g., it changed from −35 to 20% for a g-ratio of 0.73 and negligible fibre dispersion (Supplementary Figure S6). This shows that the validity of the new heuristic expression for β1 as a sum of the relaxation rates of the myelin and non-myelin water pools weighted by their signal fractions is constrained to a specific range of relaxation rate values.

In this manuscript, we provide a simulation framework that allows to test whether for a given set of microscopic parameters the validity of the heuristic expression is given.

Note that neither the proposed heuristic correction nor the previous M2-based expression account for the effect of fibre dispersion which might explain why the accuracy of the predictions decreased with increasing fibre dispersion (Figure 7). While the influence of fibre dispersion has been successfully incorporated into M2 in another study (Fritz et al., 2020), it remains an open task for future studies to also do this for the heuristic expression of β1.

 Myelin water fraction estimation from <italic>ex vivo</italic> data using the heuristic expression for <italic>β</italic><sub>1</sub>

Under the condition that M2 estimates an orientation-independent β1 and that the heuristic expression of β1 provides a valid biophysical interpretation, the myelin water fraction (MWF) can be estimated from the fitted β1 (Eq. 6). When using the ex vivo compartmental R2 values, we found a median (across orientation) MWF of 0.14 for fibres with negligible dispersion (Figure 8B), which is congruent with the mean value reported in white matter of 0.10 (Uddin et al., 2019). In the Supplementary materials section 7.2.3, we exemplified what happens if the MWF is calculated for a set of microscopic parameters for which the heuristic expression is invalid. We found that the resulting MWF is negative and thus implausible. As such, the estimation of the MWF through β1 seems a less effective method than existing MWF estimation approaches but might still be useful to estimate the MWF if magnitude-only meGRE data with a single head orientation are available.

 The effect of echo time on the performance of M2

Our findings revealed that the ability of M2 to estimate β1 was reduced for meGRE subsets with shorter maximum echo time (TEmax). This was evidenced by: (i) an increased residual orientation dependence of β1, and (ii) M2 not being better explained by the meGRE data than M1. The performance of M2 decreased when the maximum TE (TEmax) also decreased. This was not only observed for meGRE subsets with TEmax values typically used for in vivo studies (i.e., TEmax = 18 ms), but also at the intermediate TEmax (= 36 ms). Note that these observations could also be driven by the reduced time points of the meGRE subsets at shorter TEmax: while the meGRE subset at TEmax = 54 ms contained 16 time points, the meGRE subset at TEmax = 18 ms only contained five time points. A limited sample size or number of time points, however, is an unsolved challenge for in vivo application of M2 because typical in vivo meGRE protocols, specifically MPM protocols, use short TEmax (~ 18 ms) and few echo times only (~ 6–8 echoes). Therefore, future studies should aim at increasing the TEmax and/or the time points. This will require highly accelerated acquisitions [e.g., like in Han et al. (2014) for spin echo sequences or Kim et al. (2019) for 3D-GRE sequences] and the correction of motion artefacts (Magerkurth et al., 2011), B0 fluctuations due to breathing (e.g., Vannesjo et al., 2015) and susceptibility artefacts (e.g., Port and Pomper, 2000), which are particularly strong at later echo times.

Interestingly, the biggest discrepancy between in silico and ex vivo results for β1 was seen for the meGRE subset with the shortest TEmax value at θ μ → smaller than the magic angle (55°, Figure 9B). This is because β1 and α1 of the measured ex vivo data showed an atypical θ μ → dependence in this θ μ → range: they decreased as a function of increasing θ μ → up to the magic angle. A similar observation was also made by Bartels et al. (2022) for the orientation dependence of R2. They suggested that a mechanism that could explain a reduction in R2 at the magic angle would be the Magic Angle Effect in highly structured molecules like myelin sheaths (see Bydder et al., 2007). Since, in our experiment, this phenomenon would be superimposed on the orientation dependence of R2*, it may be particularly evident when the latter effect is negligible, i.e., at low θ μ → . Note that our finding was observed only for one tissue sample. Thus, further testing on different tissue samples is necessary to verify the generalisability of our finding.

 Considerations

Our results indicate that the ability of M2 to estimate the orientation-independent component of R2* varies with echo time and strongly depends on microstructural parameters. As the space of parameters in the simulations are large, not all possible combinations could be investigated here. In future studies, we will test the performance of M2 in scenarios that map directly to in vivo meGRE experiments as opposed to the ex vivo case that was the focus of this study.

M2 can separate the orientation dependence of R2* leaving an orientation-independent parameter β1, but at the same time this estimated β1 cannot be predicted accurately based on the current analytical derivation of M2 (Supplementary Equations S15, S16, section 4). Future studies should aim to find a better derivation of M2 from the HCFM that does not neglect the contribution of the myelin water as well as incorporating other sources of dephasing, e.g., due to diffusion and near-field interactions. In fact, an analytical derivation without neglecting the contribution of the myelin compartment was performed in this manuscript (Supplementary material, section 4). However, this derivation is mathematically valid only for meGRE subsets with a maximal TE smaller than the T2 of the myelin compartment. Thus, the derived expression (Supplementary Equations S13, S14, section 4) does not hold for our simulated datasets because TEmax > T2 myelin for all meGRE subsets. This might also explain why the heuristic expression does not work for the in vivo compartmental R2 values, for which the T2 myelin is smaller. Nevertheless, it can be used to motivate our heuristic expression for β1 (Eq. 4) because it is the same expression as in Supplementary Equation S15b. This derivation might also be relevant for studies that are performed at lower magnetic fields, e.g., at 3 T, where the condition TEmax > T2 myelin could be fulfilled because the R2 from the myelinated and non-myelinated compartments are different (e.g., shorter) from the ones used in our current simulation.

Our simulations did not always show the same trend as the ex vivo data (e.g., Figures 6, 9) and were occasionally quantitatively different. This could be related to simplifications that were employed in our simulations and/or the underlying simplifications of the HCFM. The most important simplifications in our simulations were: First, the assumption that the R2 was the same for both intra- and extracellular compartments. Although, these R2 have been found to be different (e.g., Beaulieu et al., 1998; Assaf and Cohen, 2000; Does and Gore, 2000; Does, 2018; Veraart et al., 2018; Tax et al., 2021), we expect the differences not to play a substantial role at the short TEs that were used here [e.g., TEmax: 54 ms < T2 of the extra-axonal compartment ≈ 58 ms in Tax et al. (2021)]. Second, we assumed that the signal coming from multiple dispersed hollow cylinders is a superposition of the complex signal of multiple single hollow cylinders at different orientations, neglecting the near-field interaction of the cylinders. As compared to previous studies where near-field interaction was more faithfully described in two dimensions (Xu et al., 2018; Hédouin et al., 2021), our simulation framework allowed for better control over the fibre dispersion in three dimensions via the Watson distribution parameter κ. The most important simplifications of the HCFM are: (1) neglecting the orientation dependence of R2 with respect to the external magnetic field (Knight et al., 2017; Birkl et al., 2021; Tax et al., 2021) and (2) the different longitudinal magnetisation of the compartments which affects the longitudinal relaxation rate (R1) (see, e.g., Labadie et al., 2014; Shin et al., 2019; van Gelderen and Duyn, 2019; Chan and Marques 2020; Kleban et al., 2021). While the anisotropic part of R2 is three times smaller than the anisotropic part of R2* at 3 T (Gil et al., 2016) and could explain residual orientation dependence of β1, other assumptions requires further study, for example removing the R1 dependence in the estimated R2* (Milotta et al., 2023). Nevertheless, even with all the simplifications, the HCFM-based in silico data described the θ μ → dependence of α1 and β1 similarly to the ex vivo data across all dispersion regimes when using the long maximal TE protocol.

The ex vivo data require further discussion. First, we investigated only one human optic chiasm tissue sample with relatively long postmortem interval of 48 h, which could explain parts of the differences that we found when comparing with the in silico dataset. Second, the coregistration of the diffusion and meGRE datasets (see section 3.1.4) might lead to image interpolation artefacts affecting the κ and θ μ → estimates. Moreover, coregistration between meGRE images at different orientations could lead to additional blurring of the data. However, these coregistration steps are necessary to ensure maximal correspondence between the same voxels across maps. We expect that the additional coregistration-related blurring will only slightly reduce the variability when binning the data (e.g., the standard deviation along R2* in Figure 6). Third, the Watson dispersion from the NODDI model cannot describe all existing fibre arrangements in the brain accurately, e.g., the crossing fibre arrangement. However, in the optic chiasm specimen crossing-fibre arrangements were only found in a few regions, e.g., at the crossing of the optical tract and optic nerve. Therefore, the contribution of such crossing-fibre voxels with estimated κ values in the range of highly to mildly dispersed fibres will be averaged-out with the single-fibre orientation voxels with similar κ values during the irregular binning pre-processing (section 3.3.1). However, this could result in an increasing standard deviation in the estimated α-parameters in the log-linear model and β-parameters in the log-quadratic model.

 Conclusion

We showed that our recently introduced biophysical log-quadratic model (M2) of the multi-echo gradient-recall echo (meGRE) signal can estimate the fibre-angular-orientation independent part of R2* (R2,iso*) for varying g-ratio values and fibre dispersions. Thus, the estimated linear time-dependent parameter of M2, β1, provides an attractive alternative for estimating R2,iso* to standard methods that require multiple acquisitions with distinct positioning of the sample in the head-coil. We also showed that β1 can be used to estimate the myelin water fraction (MWF) for ex vivo compartmental R2 values using a newly proposed heuristic expression relating β1 to microstructural tissue parameters including the myelin water signal. We provide a freely available simulation framework to test the validity of the heuristic expression for varying sets of microstructural parameters. We found that the heuristic expression cannot be used for in vivo compartmental R2 values.

Importantly, we found that an angular-independent β1 (and thus R2,iso*) cannot be estimated with the log-quadratic model for meGRE measurements with maximum shorter echo times, that are typically used for whole-brain in vivo meGRE experiments. Therefore, it indicates that we need to develop new meGRE protocols with longer echo times that remain time efficient and motion robust. This could be achieved by using highly accelerated acquisitions with a higher data sampling for shorter echo times. Finally, at echo time ranges of about 18 ms, an unexpected R2* orientation-dependence was found in the ex vivo dataset at angles below the magic angle: a decrease of R2* for increasing angles. However, more testing is required to confirm that our finding can be generalised to other brain regions and specimens since our results are based on thorough measurements of one human optic chiasm tissue sample.

 Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Author contributions

FF: conceptualization, MRI data analysis, in-silico data analysis and manuscript’s writer. LM and MC: manuscript review. MA: MRI data pre-processing and acquisition. JP and AP: MRI data acquisition and protocol design. MM: Ex-vivo specimen preparation and containment, manuscript review. CJ: Ex-vivo specimen preparation and containment. TN and NW: resources and manuscript review. KP: MRI data acquisition and protocol design, manuscript review. SM: conceptualization, funding acquisition, co-writer of the manuscript and manuscript review, supervision. All authors contributed to the article and approved the submitted version.

 Funding

This work was supported by the German Research Foundation (DFG Priority Program 2041 “Computational Connectomics,” [AL 1156/2-1;GE 2967/1-1; MO 2397/5-1; MO 2249/3-1; MO 2397/5-2], by the Emmy Noether Stipend: MO 2397/4-1, MO 2397/4-2) and by the BMBF (01EW1711A and B) in the framework of ERA-NET NEURON and the Forschungszentrums Medizintechnik Hamburg (fmthh; grant 01fmthh2017). The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013) / ERC grant agreement n° 616,905. MFC is supported by the MRC and Spinal Research Charity through the ERA-NET Neuron joint call (MR/R000050/1). The Wellcome Centre for Human Neuroimaging is supported by core funding from the Wellcome [203,147/Z/16/Z]. The Max Planck Institute for Human Cognitive and Brain Sciences has an institutional research agreement with Siemens Healthcare. NW holds a patent on acquisition of MRI data during spoiler gradients (US 10,401,453 B2). NW was a speaker at an event organized by Siemens Healthcare and was reimbursed for the travel expenses.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.1133086/full#supplementary-material

 References Alonso-Ortiz E. Levesque I. R. Pike G. B. (2018). Impact of magnetic susceptibility anisotropy at 3 T and 7 T on T2*-based myelin water fraction imaging. NeuroImage 182, 370–378. doi: 10.1016/j.neuroimage.2017.09.040, PMID: 28958882 Assaf Y. Cohen Y. (2000). Assignment of the water slow-diffusing component in the central nervous system using q-space diffusion MRS: implications for fiber tract imaging. Magn. Reson. Med. 43, 191–199. doi: 10.1002/(SICI)1522-2594(200002)43:2<191::AID-MRM5>3.0.CO;2-B, PMID: 10680682 Bartels L. M. Doucette J. Birkl C. Zhang Y. Weber A. M. Rauscher A. (2022). Orientation dependence of R2 relaxation in the newborn brain. NeuroImage 264:119702. doi: 10.1016/j.neuroimage.2022.119702, PMID: 36272671 Basser P. J. Mattiello J. LeBihan D. (1994). MR diffusion tensor spectroscopy and imaging. Biophys. J. 66, 259–267. doi: 10.1016/S0006-3495(94)80775-1 Beaulieu C. Fenrich F. R. Allen P. S. (1998). Multicomponent water proton transverse relaxation and T2-discriminated water diffusion in myelinated and nonmyelinated nerve. Magn. Reson. Imaging 16, 1201–1210. doi: 10.1016/s0730-725x(98)00151-9, PMID: 9858277 Birkl C. Doucette J. Fan M. Hernández-Torres E. Rauscher A. (2021). Myelin water imaging depends on white matter fiber orientation in the human brain. Magn. Reson. Med. 85, 2221–2231. doi: 10.1002/mrm.28543, PMID: 33017486 Burnham K. P. Anderson D. R. (2002). Model selection and multimodel inference: a practical information-theoretic approach. Springer New York, NY: Springer Verlag. Burnham K. P. Anderson D. R. Huyvaert K. P. (2011). AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behav. Ecol. Sociobiol. 65, 23–35. doi: 10.1007/s00265-010-1029-6 Bydder M. Rahal A. Fullerton G. D. Bydder G. M. (2007). The magic angle effect: a source of artifact, determinant of image contrast, and technique for imaging. J. Magn. Reson. Imaging 25, 290–300. doi: 10.1002/jmri.20850, PMID: 17260400 Callaghan M. F. Freund P. Draganski B. Anderson E. Cappelletti M. Chowdhury R. . (2014). Widespread age-related differences in the human brain microstructure revealed by quantitative magnetic resonance imaging. Neurobiol. Aging 35, 1862–1872. doi: 10.1016/j.neurobiolaging.2014.02.008, PMID: 24656835 Chan K.-S. Marques J. P. (2020). Multi-compartment relaxometry and diffusion informed myelin water imaging – promises and challenges of new gradient echo myelin water imaging methods. NeuroImage 221:117159. doi: 10.1016/j.neuroimage.2020.117159, PMID: 32663644 Chavhan G. B. Babyn P. S. Thomas B. Shroff M. M. Haacke E. M. (2009). Principles, techniques, and applications of T2*-based MR imaging and its special applications. Radiographics 29, 1433–1449. doi: 10.1148/rg.295095034, PMID: 19755604 Does M. D. (2018). Inferring brain tissue composition and microstructure via MR relaxometry. NeuroImage 182, 136–148. doi: 10.1016/j.neuroimage.2017.12.087, PMID: 29305163 Does M. D. Gore J. C. (2000). Compartmental study of diffusion and relaxation measured in vivo in normal and ischemic rat brain and trigeminal nerve. Magn. Reson. Med. 43, 837–844. doi: 10.1002/1522-2594(200006)43:6<837::aid-mrm9>3.0.co;2-o, PMID: 10861878 Draganski B. Ashburner J. Hutton C. Kherif F. Frackowiak R. S. J. Helms G. . (2011). Regional specificity of MRI contrast parameter changes in normal ageing revealed by voxel-based quantification (VBQ). NeuroImage 55, 1423–1434. doi: 10.1016/j.neuroimage.2011.01.052, PMID: 21277375 Dula A. N. Gochberg D. F. Valentine H. L. Valentine W. M. Does M. D. (2010). Multiexponential T2, magnetization transfer, and quantitative histology in white matter tracts of rat spinal cord. Magn. Reson. Med. 63, 902–909. doi: 10.1002/mrm.22267, PMID: 20373391 Duyn J. H. (2014). Frequency shifts in the myelin water compartment. Magn. Reson. Med. 71, 1953–1955. doi: 10.1002/mrm.24983, PMID: 24700549 Duyn J. H. Schenck J. (2017). Contributions to magnetic susceptibility of brain tissue. NMR Biomed. 30:e3546. doi: 10.1002/nbm.3546, PMID: 27240118 Elster A. D. (1993). Gradient-echo MR imaging: techniques and acronyms. Radiology 186, 1–8. doi: 10.1148/radiology.186.1.8416546, PMID: 8416546 Emmenegger T. M. David G. Ashtarayeh M. Fritz F. J. Ellerbrock I. Helms G. . (2021). The influence of radio-frequency transmit field inhomogeneities on the accuracy of G-ratio weighted imaging. Front. Neurosci. 15:770. doi: 10.3389/fnins.2021.674719, PMID: 34290579 Fedorov A. Beichel R. Kalpathy-Cramer J. Finet J. Fillion-Robin J.-C. Pujol S. . (2012). 3D slicer as an image computing platform for the quantitative imaging network. Magn. Reson. Imaging 30, 1323–1341. doi: 10.1016/j.mri.2012.05.001, PMID: 22770690 Fritz F. J. Ashtarayeh M. Periquito J. Pohlmann A. Morawski M. Jaeger C. ., (2021). Effects of fibre dispersion and myelin content on R2*: Simulations and post-mortem experiments. In: Proceedings of the 2021 ISMRM & SMRT Virtual Conference & Exhibition. Presented at the International Society of Magnetic Resonance Imaging (ISMRM). Fritz F. J. Edwards L. J. Streubel T. Pine K. J. Weiskopf N. Mohammadi S., (2020). Simulating the effect of axonal dispersion and noise on anisotropic R2* relaxometry in white matter. In: Proceedings of the 2020 ISMRM & SMRT Virtual Conference & Exhibition. Presented at the International Society of Magnetic Resonance Imaging (ISMRM). Gil R. Khabipova D. Zwiers M. Hilbert T. Kober T. Marques J. P. (2016). An in vivo study of the orientation-dependent and independent components of transverse relaxation rates in white matter. NMR Biomed. 29, 1780–1790. doi: 10.1002/nbm.3616, PMID: 27809376 Gudbjartsson H. Patz S. (1995). The rician distribution of noisy mri data. Magn. Reson. Med. 34, 910–914. doi: 10.1002/mrm.1910340618, PMID: 8598820 Han S. H. Cho F. H. Song Y. K. Paulsen J. Song Y. Q. Kim Y. R. . (2014). Ultrafast 3D spin-echo acquisition improves gadolinium-enhanced MRI signal contrast enhancement. Sci. Rep. 4:5061. doi: 10.1038/srep05061, PMID: 24863102 Hédouin R. Metere R. Chan K.-S. Licht C. Mollink J. van Walsum A.-M. C. . (2021). Decoding the microstructural properties of white matter using realistic models. NeuroImage 237:118138. doi: 10.1016/j.neuroimage.2021.118138, PMID: 33964461 Jeurissen B. Descoteaux M. Mori S. Leemans A. (2019). Diffusion MRI fiber tractography of the brain. NMR Biomed. 32:e3785. doi: 10.1002/nbm.3785 Kellman P. McVeigh E. R. (2005). Image reconstruction in SNR units: a general method for SNR measurement. Magn. Reson. Med. 54, 1439–1447. doi: 10.1002/mrm.20713, PMID: 16261576 Kim T. H. Bilgic B. Polak D. Setsompop K. Haldar J. P. (2019). Wave-LORAKS: combining wave encoding with structured low-rank matrix modeling for more highly accelerated 3D imaging. Magn. Reson. Med. 81, 1620–1633. doi: 10.1002/mrm.27511, PMID: 30252157 Kirilina E. Helbling S. Morawski M. Pine K. Reimann K. Jankuhn S. . (2020). Superficial white matter imaging: contrast mechanisms and whole-brain in vivo mapping. Sci. Adv. 6:eaaz9281. doi: 10.1126/sciadv.aaz9281, PMID: 33028535 Kleban E. Gowland P. Bowtell R. (2021). Probing the myelin water compartment with a saturation-recovery, multi-echo gradient-recalled echo sequence. Magn. Reson. Med. 86, 167–181. doi: 10.1002/mrm.28695, PMID: 33576521 Knight M. J. Dillon S. Jarutyte L. Kauppinen R. A. (2017). Magnetic resonance relaxation anisotropy: physical principles and uses in microstructure imaging. Biophys. J. 112, 1517–1528. doi: 10.1016/j.bpj.2017.02.026, PMID: 28402893 Kucharczyk W. Macdonald P. M. Stanisz G. J. Henkelman R. M. (1994). Relaxivity and magnetization transfer of white matter lipids at MR imaging: importance of cerebrosides and pH. Radiology 192, 521–529. doi: 10.1148/radiology.192.2.8029426, PMID: 8029426 Labadie C. Lee J.-H. Rooney W. D. Jarchow S. Aubert-Frécon M. Springer C. S. Jr. . (2014). Myelin water mapping by spatially regularized longitudinal relaxographic imaging at high magnetic fields. Magn. Reson. Med. 71, 375–387. doi: 10.1002/mrm.24670, PMID: 23468414 Langkammer C. Krebs N. Goessler W. Scheurer E. Ebner F. Yen K. . (2010). Quantitative MR imaging of brain iron: a postmortem validation study. Radiology 257, 455–462. doi: 10.1148/radiol.10100495, PMID: 20843991 Lee J. Shmueli K. Kang B.-T. Yao B. Fukunaga M. van Gelderen P. . (2012). The contribution of myelin to magnetic susceptibility-weighted contrasts in high-field MRI of the brain. NeuroImage 59, 3967–3975. doi: 10.1016/j.neuroimage.2011.10.076, PMID: 22056461 Lee J. van Gelderen P. Kuo L.-W. Merkle H. Silva A. C. Duyn J. H. (2011). T2*-based fiber orientation mapping. NeuroImage 57, 225–234. doi: 10.1016/j.neuroimage.2011.04.026, PMID: 21549203 Lee J. Shin H.-G. Jung W. Nam Y. Oh S.-H. Lee J. (2017). An R2* model of white matter for fiber orientation and myelin concentration. NeuroImage 162, 269–275. doi: 10.1016/j.neuroimage.2017.08.050, PMID: 28823825 Li T.-Q. Yao B. van Gelderen P. Merkle H. Dodd S. Talagala L. . (2009). Characterization of T2* heterogeneity in human brain white matter. Magn. Reson. Med. 62, 1652–1657. doi: 10.1002/mrm.22156, PMID: 19859939 Liu C. (2010). Susceptibility tensor imaging. Magn. Reson. Med. 63, 1471–1477. doi: 10.1002/mrm.22482, PMID: 20512849 MacKay A. Laule C. Vavasour I. Bjarnason T. Kolind S. Mädler B. (2006). Insights into brain microstructure from the T2 distribution. Magn. Reson. Imaging 24, 515–525. doi: 10.1016/j.mri.2005.12.037, PMID: 16677958 Magerkurth J. Volz S. Wagner M. Jurcoane A. Anti S. Seiler A. . (2011). Quantitative T*2-mapping based on multi-slice multiple gradient echo flash imaging: retrospective correction for subject motion effects. Magn. Reson. Med. 66, 989–997. doi: 10.1002/mrm.22878, PMID: 21381108 Milotta G. Corbin N. Lambert C. Lutti A. Mohammadi S. Callaghan M. F. (2023). Mitigating the impact of flip angle and orientation dependence in single compartment R2* estimates via 2-pool modeling. Magn. Reson. Med. 89, 128–143. doi: 10.1002/mrm.29428, PMID: 36161672 Minassian H. Huang S. (1979). Effect of sodium azide on the ultrastructural preservation of tissues. J. Microsc. 117, 243–253. doi: 10.1111/j.1365-2818.1979.tb01180.x, PMID: 399299 Mohammadi S. Carey D. Dick F. Diedrichsen J. Sereno M. I. Reisert M. . (2015). Whole-brain in-vivo measurements of the axonal G-ratio in a group of 37 healthy volunteers. Front. Neurosci. 9:441. doi: 10.3389/fnins.2015.00441, PMID: 26640427 Oh S.-H. Kim Y.-B. Cho Z.-H. Lee J. (2013). Origin of B0 orientation dependent R2* (=1/T2*) in white matter. NeuroImage 73, 71–79. doi: 10.1016/j.neuroimage.2013.01.051, PMID: 23376494 Ordidge R. J. Gorell J. M. Deniau J. C. Knight R. A. Helpern J. A. (1994). Relative assessment of brain iron levels using MRI at 3 tesla. Magn. Reson. Mater. Phys. Biol. Med. 2, 449–450. doi: 10.1007/BF01705295 Pampel A. Müller D. K. Anwander A. Marschner H. Möller H. E. (2015). Orientation dependence of magnetization transfer parameters in human white matter. NeuroImage 114, 136–146. doi: 10.1016/j.neuroimage.2015.03.068, PMID: 25862261 Papazoglou S. Streubel T. Ashtarayeh M. Pine K. J. Edwards L. J. Brammerloh M. . (2019). Biophysically motivated efficient estimation of the spatially isotropic component from a single gradient-recalled echo measurement. Magn. Reson. Med. 82, 1804–1811. doi: 10.1002/mrm.27863, PMID: 31293007 Port J. D. Pomper M. G. (2000). Quantification and minimization of magnetic susceptibility artifacts on GRE images. J. Comput. Assist. Tomogr. 24, 958–964. doi: 10.1097/00004728-200011000-00024, PMID: 11105718 Roebroeck A. Miller K. L. Aggarwal M. (2018). Ex vivo diffusion MRI of the human brain: technical challenges and recent advances. NMR Biomed. 32:e3941. doi: 10.1002/nbm.3941 Rudko D. A. Klassen L. M. de Chickera S. N. Gati J. S. Dekaban G. A. Menon R. S. (2014). Origins of <em>R</em><sub class="stack">2</sub><sup class="stack">∗</sup> orientation dependence in gray and white matter. Proc. Natl. Acad. Sci. 111, E159–E167. doi: 10.1073/pnas.1306516111, PMID: 24374633 Schmahmann J. D. Pandya D. N. Wang R. Dai G. D’Arceuil H. E. de Crespigny A. J. . (2007). Association fibre pathways of the brain: parallel observations from diffusion spectrum imaging and autoradiography. Brain 130, 630–653. doi: 10.1093/brain/awl359, PMID: 17293361 Schmahmann J. D. Schmahmann J. Pandya D., (2009). Fiber pathways of the brain. Oxford University Press, USA. Schyboll F. Jaekel U. Petruccione F. Neeb H. (2019). Fibre-orientation dependent R1(=1/T1) relaxation in the brain: the role of susceptibility induced spin-lattice relaxation in the myelin water compartment. J. Magn. Reson. 300, 135–141. doi: 10.1016/j.jmr.2019.01.013, PMID: 30743171 Shin H.-G. Oh S.-H. Fukunaga M. Nam Y. Lee D. Jung W. . (2019). Advances in gradient echo myelin water imaging at 3T and 7T. NeuroImage 188, 835–844. doi: 10.1016/j.neuroimage.2018.11.040, PMID: 30476624 Sra S. Karp D. (2013). The multivariate Watson distribution: maximum-likelihood estimation and other aspects. J. Multivar. Anal. 114, 256–269. doi: 10.1016/j.jmva.2012.08.010 Stikov N. Perry L. M. Mezer A. Rykhlevskaia E. Wandell B. A. Pauly J. M. . (2011). Bound pool fractions complement diffusion measures to describe white matter micro and macrostructure. NeuroImage 54, 1112–1121. doi: 10.1016/j.neuroimage.2010.08.068 Tax C. M. W. Kleban E. Chamberland M. Baraković M. Rudrapatna U. Jones D. K. (2021). Measuring compartmental T2-orientational dependence in human brain white matter using a tiltable RF coil and diffusion-T2 correlation MRI. NeuroImage 236:117967. doi: 10.1016/j.neuroimage.2021.117967, PMID: 33845062 Tofts P. S. (2004). “Concepts: measurement and MR” in Quantitative MRI of the brain. ed. Tofts, P. (Wiley, New Jersey, U.S.: John Wiley & Sons, Ltd), 1–15. Uddin M. N. Figley T. D. Solar K. G. Shatil A. S. Figley C. R. (2019). Comparisons between multi-component myelin water fraction, T1w/T2w ratio, and diffusion tensor imaging measures in healthy human brain structures. Sci. Rep. 9:2500. doi: 10.1038/s41598-019-39199-x, PMID: 30792440 van Gelderen P. Duyn J. H. (2019). White matter intercompartmental water exchange rates determined from detailed modeling of the myelin sheath. Magn. Reson. Med. 81, 628–638. doi: 10.1002/mrm.27398, PMID: 30230605 Vannesjo S. J. Wilm B. J. Duerst Y. Gross S. Brunner D. O. Dietrich B. E. . (2015). Retrospective correction of physiological field fluctuations in high-field brain MRI using concurrent field monitoring. Magn. Reson. Med. 73, 1833–1843. doi: 10.1002/mrm.25303, PMID: 24903278 Veraart J. Novikov D. S. Fieremans E. (2018). TE dependent diffusion imaging (TEdDI) distinguishes between compartmental T2 relaxation times. NeuroImage 182, 360–369. doi: 10.1016/j.neuroimage.2017.09.030, PMID: 28935239 Wang N. Zhang J. Cofer G. Qi Y. Anderson R. J. White L. E. . (2019). Neurite orientation dispersion and density imaging of mouse brain microstructure. Brain Struct. Funct. 224, 1797–1813. doi: 10.1007/s00429-019-01877-x, PMID: 31006072 Weber A. M. Zhang Y. Kames C. Rauscher A. (2020). Myelin water imaging and R2* mapping in neonates: investigating R2* dependence on myelin and fibre orientation in whole brain white matter. NMR Biomed. 33:e4222. doi: 10.1002/nbm.4222, PMID: 31846134 Weiskopf N. Edwards L. J. Helms G. Mohammadi S. Kirilina E. (2021). Quantitative magnetic resonance imaging of brain anatomy and in vivo histology. Nat. Rev. Phys. 3, 570–588. doi: 10.1038/s42254-021-00326-1 Weiskopf N. Suckling J. Williams G. Correia M. M. Inkster B. Tait R. . (2013). Quantitative multi-parameter mapping of R1, PD(*), MT, and R2(*) at 3T: a multi-center validation. Front. Neurosci. 7:95. doi: 10.3389/fnins.2013.00095 Wharton S. Bowtell R. (2013). Gradient Echo based Fiber orientation mapping using R2* and frequency difference measurements. NeuroImage 83, 1011–1023. doi: 10.1016/j.neuroimage.2013.07.054, PMID: 23906549 Wharton S. Bowtell R. (2012). Fiber orientation-dependent white matter contrast in gradient echo MRI. Proc. Natl. Acad. Sci. 109, 18559–18564. doi: 10.1073/pnas.1211075109, PMID: 23091011 Xu T. Foxley S. Kleinnijenhuis M. Chen W. C. Miller K. L. (2018). The effect of realistic geometries on the susceptibility-weighted MR signal in white matter. Magn. Reson. Med. 79, 489–500. doi: 10.1002/mrm.26689, PMID: 28394030 Yablonskiy D. A. Haacke E. M. (1994). Theory of NMR signal behavior in magnetically inhomogeneous tissues: the static dephasing regime. Magn. Reson. Med. 32, 749–763. doi: 10.1002/mrm.1910320610, PMID: 7869897 Yao B. Li T. -Q. Gelderen P. Van Shmueli K. de Zwart J. A. Duyn J. H., (2009). Susceptibility contrast in high field MRI of human brain as a function of tissue iron content. NeuroImage 44, 1259–1266. doi: 10.1016/j.neuroimage.2008.10.029, PMID: 19027861 Zhang H. Schneider T. Wheeler-Kingshott C. A. Alexander D. C. (2012). NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61, 1000–1016. doi: 10.1016/j.neuroimage.2012.03.072, PMID: 22484410 Ziegler G. Hauser T. U. Moutoussis M. Bullmore E. T. Goodyer I. M. Fonagy P. . (2019). Compulsivity and impulsivity traits linked to attenuated developmental frontostriatal myelination trajectories. Nat. Neurosci. 22, 992–999. doi: 10.1038/s41593-019-0394-3, PMID: 31086316 Glossary
Acronyms
Biophysical terms and model parameters
AWF (Intra-) Axonal water fraction
EWF Extra-axonal water fraction
FVF Fibre volume fraction
HCFM Hollow cylinder fibre model
ICVF Intra-cellular volume fraction (from NODDI)
MWF Myelin water fraction (Eq. 6)
gratio g-ratio
Magnetic resonance imaging and sequence acronyms
dMRI Diffusion-weighted Magnetic Resonance Imaging
DWI Diffusion-weighting Imaging
GRE Gradient-recalled echo
meGRE Multi-echo gradient-recalled echo
OC Optic chiasm
R2* Effective transverse relaxation rate
R2,iso* Orientation independent or isotropic part of R2*
TE Echo time
TEmax Maximal echo time
Hollow cylinder fibre model parameters
SA Signal of the intra-axonal compartment
SE Signal of the extra-axonal compartment
SM Signal of the myelin compartment
SN Sum of the signals of the non-myelinated (SA and SE) compartments
SC Sum of all the signal compartments (SA, SE and SM)
R2A Transverse relaxation rate of the intra-axonal compartment
R2E Transverse relaxation rate of the extra-axonal compartment
R2N Transverse relaxation rate of the non-myelinated compartments
R2M Transverse relaxation rate of the myelin compartment
ρA Proton density of the intra-axonal compartment
ρE Proton density of the extra-axonal compartment
ρN Proton density of the non-myelinated compartments
ρN Proton density of the myelin compartment
VM Volume fraction of the myelin compartment
Symbols
In silico and ex vivo data descriptors
θ μ → Angular orientation of the mean fibre bundle
θ 0 First angular orientation or angular offset
B 0 → Main magnetic field
κ Coefficient of dispersion (from Watson Distribution and NODDI)
μ → Vector of the mean fibre bundle
x → Vector of the individual cylinder in the simulated in silico data
T Diff , GRE Transformation matrix from dMRI to GRE images
T GRE : i , 1 Transformation matrix from GRE images at the i-th angular orientation measurement to the first angular orientation measurement
Model parameters and analysis metrics
α0 Intercept parameter of M1
α1 Slope or linear parameter of M1
β0 Intercept of M2
β1 Slope of linear parameter of M2
β1,nm β1 ground-truth value without myelin signal contribution (Eq. 3)
β1,m β1 ground-truth value with myelin signal contribution (Eq. 4)
β2 Quadrature or second order parameter of M2
εm Relative difference between fitted β1 and predicted β1,nm
εnm Relative difference between fitted β1 and predicted β1,m
AIC Akaike Information Criterion
AICc Akaike Information Criterion corrected
ΔAICc Difference of Akaike Information Criteria (Equation 13)
wAICc Weighted Akaike Information Criterion corrected (Eq. 12)
nRMSD Normalised root-mean-squared deviation (Eq. 9)
ΔRMSD Normalised root-mean-squared deviation difference (Eq. 10)
M1 Log-linear model (Eq. 2)
M2 Log-quadratic model (Eq. 1)

1http://www.slicer.org

2http://www.fil.ion.ucl.ac.uk/spm

3https://github.com/quantitative-mri-and-in-vivo-histology/r2s_iso_estimation