Deciphering the fibre-orientation independent component of R2* (R2,iso*) in the human brain with a single multi-echo gradient-recalled-echo measurement under varying microstructural conditions

The effective transverse relaxation rate (R2*) is sensitive to the microstructure of the human brain, e.g. the g-ratio characterising the relative myelination of axons. However, R2* depends on the orientation of the fibres relative to the main magnetic field degrading its reproducibility and that of any microstructural derivative measure. To decipher its orientation-independent part (R2,iso*), a second-order polynomial in time (M2) can be applied to single multi-echo gradient-recalled-echo (meGRE) measurements at arbitrary orientation. The linear-time dependent parameter, β1, of M2 can be biophysically related to R2,iso* when neglecting the signal from the myelin water (MW) in the hollow cylinder fibre model (HCFM). Here, we examined the effectiveness of M2 using experimental and simulated data with variable g-ratio and fibre dispersion. We showed that the fitted β1 effectively estimates R2,iso*when using meGRE with long maximum echo time (TEmax ≈ 54 ms) but its microscopic dependence on the g-ratio was not accurately captured. This error was reduced to less than 12% when accounting for the MW contribution in a newly introduced biophysical expression for β1. We further used this new expression to estimate the MW fraction (0.14) and g-ratio (0.79) in a human optic chiasm. However, the proposed method failed to estimate R2,iso* for a typical in-vivo meGRE protocol (TEmax ≈ 18 ms). At this TEmax and around the magic angle, the HCFM-based simulations failed to explain the R2*-orientation-dependence. In conclusion, estimation of R2,iso* with M2 in vivo requires meGRE protocols with very long TEmax ≈ 54 ms.

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 (R2A = R2E = R2N) R2M Transverse relaxation rate of the myelin compartment B. Symbols a. In silico and ex vivo data

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 (Does, 2018;MacKay et al., 2006;Weiskopf et al., 2021). The microstructural sensitivity of R2* makes it particularly interesting for neuroscience and clinical research studies Draganski et al., 2011;Kirilina et al., 2020;Langkammer et al., 2010). This is because R2* is sensitive not only to free and myelin water pools in the brain (Dula et al., 2010a;MacKay et al., 2006;Weiskopf et al., 2021) but also to microscopic perturbations in the main magnetic field ( 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 (Alonso-Ortiz et al., 2018;Duyn, 2014;Kucharczyk et al., 1994;Lee et al., 2017) and paramagnetic iron deposits in glial cells Ordidge et al., 1994;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 0 ⃗⃗⃗⃗ (Lee et al., 2011(Lee et al., ,2012 confounding the mapping of R2* to the underlying microstructure. The angular orientation dependence of R2* can be decomposed into an isotropic, i.e. angular-independent, component (R2,iso*) and an angular-dependent component using either gradient-recalled echo (GRE) acquisitions at several angular orientations (Oh et al., 2013;Rudko et al., 2014;Wharton and Bowtell, 2013) or hybrid diffusion weighted imaging (DWI) and GRE acquisitions with reduced numbers of distinct 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 diverse range of topographies, i.e. show fanning and bending, or mildly to acute crossing, e.g. (Jeurissen et al., 2019;Schmahmann et al., 2009Schmahmann et al., , 2007 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 incorporating very long maximum echo times of ≈ 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 singleorientation meGRE, while varying biological fibre properties and maximum echo times. Moreover, it tests the counter-intuitive hypothesis, based on M2, that the estimated β1 is independent of the MWF.
To this end, we use simulated (hereafter in silico) data and ex vivo MRI. The in silico data were simulated using 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 multiorientation meGRE imaging of a human optic chiasm to generate gold-standard datasets where the fibre orientation and dispersion are known. Both datasets are used to perform the following analyses: First, we assess the capability 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 hypothesis that β1 is independent of MWF by evaluating the deviation between fitted β1 using the in silico data and the biophysically-predicted β1 by M2. Additionally, we perform the same comparison as above using a novel heuristic expression that incorporates the MWF dependence into the predicted β1. Third, we demonstrate that the heuristic expression for β1 can be used to calculate MWF and the g-ratio from the β1 of the ex vivo data. And, fourth, we assess the capability of M2 to estimate R2,iso* via β1 for shorter maximal echo times more typical of in vivo meGRE applications.

Background
2.1 Overview of the hollow cylinder fibre model and the approximated log -quadratic model.
The hollow cylinder fibre model (HCFM, Bowtell, 2013, 2012)) proposes an analytical approximation describing the dependence of the GRE signal on the angular orientation ( ⃗ ⃗ ) defined as the angle between the main magnetic field 0 ⃗⃗⃗⃗ and the hollow-cylinder fibre ( ). This approximation establishes that the total MR signal comes from water molecules in an infinitely long 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 Appendix, section 9.1). When the signal of the water molecules in the myelin compartment is neglected, the signal magnitude of the HFCM can be approximated by a log-quadratic model (M2) in time (Papazoglou et al., 2019): , where 0,1,2 are the model parameters and SN is the non-myelin signal (i.e., SN = SA + SE). In this model, the slope β1 is considered as a proxy for R2,iso* because it does not possess any ⃗ ⃗ dependence (Equations A16b), whereas β2 contains all the ⃗ ⃗ dependent information of R2* (Equations A16c, details in Appendix, section 9.3).
In contrast to M2, the slope (α1) in the classic log-linear model (M1, (Elster, 1993)) of R2* is a function of R2,iso* and the ⃗ ⃗ dependent components of R2* (e.g. see (Lee et al., 2012b(Lee et al., , 2011): 2.2. Myelin dependence of β1 parameter as predicted by the log -quadratic model (M2) The slope β1 of M2, which is considered to be a proxy for R2,iso*, is derived from the HCFM of a two-pool system in the slow-exchange regime: a fast decaying water pool consisting of the myelin water with a relaxation rate 2 and a slow decaying water pool with a relaxation rate 2 consisting on the intra and extra cellular water. In this model, the only source of dephasing is caused by the hollow-cylinder fibre and all potential other perturbers are ignored (e.g. non-local effects of susceptibility inhomogeneities due to cavities, vessels, and iron molecules). Consequently in the approximation of M2 (Equation A16b, section 9.3), the predicted 1 parameter (hereafter 1, ) is given by the transverse relaxation rate of the non-myelin water pool ( 2 ): We hypothesise here that for realistic tissue composition (i.e. g-ratio equal to or smaller than 0.8), where the myelin compartment cannot be neglected, Equation 3 is invalid. This hypothesis is supported by previous observations showing that R2,iso* depends on the myelin water fraction, MWF (e.g. (Lee et al., 2017;Weber et al., 2020)).
Here, we propose an alternative heuristic biophysical expression of the predicted β1 parameter (hereafter 1, ), also based on the HCFM (Equation A17) but without assuming the myelin compartment to be negligible. In this case, the expected dependence of R2,iso* on variation in the MWF remains: where 2 is the relaxation rate of the myelin water pool. It follows from the heuristic model that the fitted β1 is a weighted sum of the relaxation rates of the two pools, and that Equation 3 ( 1, ) is a special case of Equation 4 ( 1, ) when the MWF is equal to 0.
Based on our hypothesis, we expect that the heuristic expression for 1, can better describe the fitted β1 when varying the g-ratio, and thus is a better proxy of R2,iso*.
Finally, we describe the dependence of MWF on the g-ratio by redefining the compartmental volumes: intra-axonal VA, extra-axonal VE and myelin VM (Equation 18a), as a function of the gratio and fibre volume fraction (FVF) as: VA = FVF • g 2 ratio, VE = 1 -FVF, and VM = FVF • (1 -g 2 ratio). If the proton densities of the non-myelinated compartments are equal (ρA = ρE = ρN), then the MWF can be rewritten as: Therefore, the g-ratio could be estimated from the MWF if the proton densities and FVF were known.

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.
All R2*-weighted GRE 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) and scanned at sixteen orientations (covering a solid angle, with azimuthal and elevation angles from 0° to 90°, Figure 1A) using the 3D multi-echo GRE (meGRE) MRI (hereafter: GRE dataset). For each angular meGRE measurement, 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 112 3 , 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. The acquired MR data are the same as reported in (Papazoglou et al., 2019).

Dispersion and mean fibre orientation estimation from dMRI dataset
To incorporate the voxel-wise information regarding the angular orientation of the fibres to 0 ⃗⃗⃗⃗ and fibre's dispersion, the dMRI datasets 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/mm 2 and 4000 s/mm 2 ) and only the fractional anisotropy (FA) map was estimated, because this map was used only for diffusion-to-GRE coregistration (section 3.1.3).

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 angular measurements in GRE space (see Figure 1A); and second, a transformation matrix that coregister from GRE space to dMRI space (see Figure 1B). The coordinate system of GRE space was defined by the first meGRE angular measurement.
To estimate the transformation matrices that coregister the angular meGRE measurements to the first angular meGRE measurement, a manual coregistration was performed and refined later with an automatic coregistration. This pair of coregistrations resulted in the transformation matrix : ,1 (i = 2 … 16). When aligning the meGRE volumes, the respective 0 ⃗⃗⃗⃗ directions had to be adjusted accordingly. This was done by aligning the 0 ⃗⃗⃗⃗ direction of the -th meGRE using the respective transformation matrix : ,1 (see insets). Both coregistrations were performed using the 3D Slicer software (http://www.slicer.org and (Fedorov et al., 2012)).
To estimate the transformation matrix from GRE to dMRI space, the FA map from the DTI analysis (section 3.1.2) was coregistered to the first meGRE angular measurement ( Figure 1B). This transformation matrix, , , was applied to coregister the NODDI results, i.e., κ, and ICVF, to the GRE space. However, this transformation matrix was also used to align the sixteen new 0 ⃗⃗⃗⃗ Same coronal slice at 1st and last GRE angular measurements (1st echo) (1) 3.1.4. Voxel-wise estimation of the angular orientation, ⃗ ⃗ , between fibres and 0 ⃗⃗⃗⃗ : The angular orientation ⃗ ⃗ between fibres and 0 ⃗⃗⃗⃗ for each meGRE angular measurement was calculated in dMRI space by computing the arccosine of the inner product between 0 ⃗⃗⃗⃗ ( ) and , i.e., ⃗ ⃗ = ( 0 ⃗⃗⃗⃗ ( ) • ) ( Figure 3C). In this computation, 0 ⃗⃗⃗⃗ ( ) is the resulting 0 ⃗⃗⃗⃗ after the transformation from the i-th meGRE angular measurement to the first meGRE angular measurement ( : ,1 ), and the transformation from GRE to dMRI space ( , −1 ) ( Figure 3A). The main fibre direction was obtained by the map from the NODDI analysis ( Figure 3B).
Note that ⃗ ⃗ was computed in dMRI space instead of GRE space to avoid undersampling and interpolation because of transforming to GRE space. These sources of error do not occur by transforming 0 ⃗⃗⃗⃗ to dMRI space, i.e. computing , −1 ⋅ : ,1 ⋅ 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 , . Exemplary ⃗ ⃗ maps in GRE space are shown in Supplementary Material, Figure S1 (first row).

Masking and pooling the ex vivo 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 a 7.2% (~ 600 over 8744 voxels). By pooling the data, the resulting cumulated data has the signal decays as a function of TE but also of ⃗ ⃗ , and fibre dispersion assessed by κ.
3.2. Simulated R2* signal decay from the HCFM R2* 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, an averaged MR signal was calculated from an ensemble of 1500 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 at +X axis, and elevation angle θ rotating from 0° (+Z) to 180° (-Z). The signal contribution per hollow cylinder was modelled with the hollow cylinder fibre model (HCFM) for all the compartments, SC (Equation A1).
In this work, two considerations were taken. First, the 0 ⃗⃗⃗⃗ was fixed and oriented parallel to +Z this piece-wise function was observed in the so-called critical time (Wharton and Bowtell, 2013;Yablonskiy and Haacke, 1994). See section 9.2 for a detailed discussion.
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 Equation 6b). This weight from the Watson distribution was calculated using the position of each simulated cylinder, ⃗⃗⃗ , and a mean fibre orientation , both defined with spherical coordinates (φ, θ) and ( ⃗ ⃗ , ⃗ ⃗ ), respectively. 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: In Equation 9b, 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 ⃗⃗⃗ and . The level of dispersion was modulated by the parameter κ (Sra and Karp, 2013;Zhang et al., 2012) 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 for both datasets. This simulation approach was used in previous conference publications (Fritz et al., 2020) and .

from the Watson distribution and (Equation 9b
). The parameter κ is limited from κ = 0 for isotropically dispersed to κ = infinity to fully parallel fibres. Here, is parallel to 0 ⃗⃗⃗⃗ .
With the ensemble averaged signal equation (Equation 6a), R2*-weighted signal decay can be created. In this work, the R2*-weighted signal decay was dependent on ⃗ ⃗ , κ and g-ratio, as a function of time. The values used are reported in Appendix (Table A3). The remaining fixed parameters required by the ensemble-averaged equation were obtained from (Dula et al., 2010a;Wharton and Bowtell, 2013) and are listed in section 9.4 in Appendix (Table A1 and A2).
Finally, each simulated R2*-weighted signal decay was replicated 5000 times with an additive Gaussian complex noise (Gudbjartsson and Patz, 1995) to approximate the SNR of the experimental ex vivo data (see section 9.4). 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 (section 3.1.4).

Data fitting and binning
The ex vivo data (section 3.1) and in silico data (each of 5000 replicas per simulated R2*weighted signal decay, section 3.2) were analysed with the log-linear and log-quadratic models, M1 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 the in silico data. However, the in silico data required two extra averages on the fitted parameters: first, across the 5000 replicas and, second, across the κ values used for simulation. 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 section 9.4).
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).
After separating the fitted parameters per fibre dispersion range for both datasets, the data was irregularly binned per ⃗ ⃗ bin per defined κ range. This was performed to avoid any bias due to effect size, since a non-uniform distribution of voxels was found in the ex vivo cumulated data as a function of ⃗ ⃗ ( Figure 5B, blue bars). To estimate the irregular ⃗ ⃗ bin, 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 Table A.4 in section 9.4.

Figure 5: Preparation of the ex vivo data for analysis. (A) The cumulated ex vivo data was 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 (Table A4,  After binning, the average and standard deviation (sd) was calculated per irregular ⃗ ⃗ bin in the ex vivo cumulated data. For the in silico data, the average and sd 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 section 9.4).

Quantitative analysis
Four different analyses were performed for both datasets in order to study: (1) the effect of gratio and fibre dispersion, via κ, on the estimated angular-independent β1 using M2, (2) the microstructural interpretability of β1 via the deviation between fitted β1 and its predicted counterparts from M2 (β1,nm, Equation 3) and from the heuristic expression (β1,m, Equation 4), (3) the possibility of calculating the MWF and g-ratio from the fitted β1 using the heuristic expression β1,m (Equation 4), and (4) the effect of TE on the performance of M2 in estimating R2,iso* from β1. The first two analyses were aimed to test whether β1 can be used as a proxy of R2,iso*.
For the first analysis, the capability 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 normalized-root-mean-squared deviation (nRMSD): where 0 varied slightly for each κ range (sub-index j) but was close to zero (see Table A.4 in section 9.4).
To compare the nRMSD of each parameter, we calculated the difference between them, ΔnRMSD, as: 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, or disentangling the ⃗ ⃗ dependency from R2*.
The above analyses were performed per g-ratio across κ range using the same values as for in silico data (Tables A1, A2 and A3, section 9.4). For the third analysis, β1,m (Equation 4) was rearranged to estimate MWF from the fitted β1 in ex vivo data. For that, the R2 values of the non-myelinated (R2N) and myelinated (R2M) compartments are reported in Table A1. After estimating the MWF, the g-ratio values were also estimated by rearranging Equation 5. For that, the fibre volume fraction, FVF, and proton density values, ρN and ρM, required for this calculation are reported in Table A1.
For the last analysis, the effect of TE on the capability of M2 to estimate R2,iso* via β1 was assessed by comparing the ex vivo dataset and the in silico data with similar g-ratio, obtained from the previous analysis. For that, α1 and β1 from M1 and M2 were compared once again as in the first analysis. However, now the models were fitted to meGRE datasets with different longest echo times (TEmax): 54 ms, 36 ms and 18 ms. Again, the ΔnRMSD was calculated to assess the residual ⃗ ⃗ dependence of β1 in comparison to the ⃗ ⃗ dependence of α1.

First analysis:
Capability of M2 to obtain the angular-independent β1 parameter for varying g-ratio and fibre dispersion values Figure 6 and 7 show the capability of M2 to estimate R2,iso* via β1 for variable g-ratio and fibre dispersion. To visualise this, the ⃗ ⃗ dependency of α1 from M1 was compared to the residual ⃗ ⃗ dependency β1 from M2 ( Figure 6A and 6B). Both ⃗ ⃗ dependencies were quantified in Figure 7A and 7B using their respective nRMSD (Equation 7a) and the ΔnRMSD in Figure 7C (Equation 7b). The results are from the analysis performed on the in silico and ex vivo data.
The capability of M2 to reduce the ⃗ ⃗ dependency of β1 varied with g-ratio and fibre dispersion, the ⃗ ⃗ 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 ( Figure 6A and 6B, respectively).  The fibre dispersion affected the performance of M2 differently between in silico and ex vivo datasets. For the ex vivo data, the nRMSD(β1) was the lowest for the negligibly dispersed fibres A B C (nRMSD(β1): 1.3% at κ ≥ 2.5) but less so for the highly dispersed fibres (nRMSD(β1) : 4.1% at κ < 1). 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% to 2.7% with decreasing fibre dispersion). The nRMSD(β1) was higher, but still below 12%, for g-ratios of 0.66 and 0.8. The ⃗ ⃗ dependency of α1 on fibre dispersion was the same between in silico and ex vivo datasets: the lower the dispersion the higher the nRMSD(α1). The ⃗ ⃗ dependency of α1 increased the lower the g-ratio was. When comparing the residual ⃗ ⃗ dependency of β1 with the ⃗ ⃗ dependency of α1, the improvement is large for negligible dispersion (from ΔnRMSD = -11.9%-points to ΔnRMSD = -37.4%-points) for both datasets.  Figure 8A and 8B report the angular-orientation ( ⃗ ⃗ ) dependent relative differences ( and ,Equation 8) between the fitted β1 from the in silico data and its predicted counterparts using M2

(Equation 3) and the heuristic expression (Equation 4
). Figure 8C shows the mean and standard deviation of and across angles.

A B C
was large, between -100% and -40%, and varied strongly with g-ratio and fibre dispersion. Even more, showed an ⃗ ⃗ dependence where the largest deviation was observed for the smallest g-ratio (0.66) and the lowest fibre dispersion ( Figure 7A). By contrast, was smaller, between -20% and 20%, and showed a smaller ⃗ ⃗ dependence, which was largest for the smallest g-ratio and lowest fibre dispersion. On average, we found that negligibly dispersed fibres showed the smallest and .
The mean across angles for , ⟨ ⟩, was smaller than 85% whereas the mean across angles for , ⟨ ⟩,was smaller than 12% ( Figure 8C). On average across all g-ratios and fibre dispersion arrangements, ⟨ ⟩ was approximately 8 to 9 times larger than ⟨ ⟩. Both relative mean differences decreased with increasing g-ratio and decreasing fibre dispersion for almost all ⟨ ⟩ and ⟨ ⟩. ⟨ ⟩ for the negligibly dispersed fibres at g-ratio 0.66 was close to 2% but accompanied by a large standard deviation across ⃗ ⃗ , indicating strong ⃗ ⃗ -dependency of the corresponding fitted 1 parameters. For both and , the variability ( Figure 8C) across different ⃗ ⃗ values, ( ) and ( ) respectively, was highest when the fibre dispersion and g-ratio were lowest.

Third analysis: Myelin water fraction (MWF) and g-ratio estimation from ex vivo data using
the heuristic expression of R2,iso* via β1,m Figure 9 reports the MWF estimated from the ex vivo data by inverting the heuristic expression for β1,m (Equation 4). Figure 9A shows the estimated MWF as a function of ⃗ ⃗ while Figure 9B shows the median and standard deviation (sd) of the estimated MWF across ⃗ ⃗ .

A B C
The estimated MWF was larger with decreasing fibre dispersion ( Figure 9A). Moreover, there was a trend towards larger estimated MWF for larger ⃗ ⃗ . On average across ⃗ ⃗ (Figure 9B), the estimated median ex vivo MWF decreased by 98% from highly dispersed fibres (MWF: 0.0028) and by 50.8% from mildly dispersed fibres (MWF: 0.069) in comparison to negligibly dispersed fibres (MWF: 0.14). The standard deviation across MWF was similar for different fibre dispersions, ranging from 0.0068 to 0.0104.  Table A2. This calculation was performed per angle ( ⃗ ⃗ ) and for the three different fibre dispersion ranges: highly dispersed (green olive), mildly dispersed (cyan) and negligibly dispersed (magenta). (B) The corresponding median and standard deviation (sd) were estimated across ⃗ ⃗ per fibre dispersion range.
With the estimated median MWF, a median g-ratio can be estimated using Equation 5. The estimated median g-ratios were 0.996, 0.895 and 0.785 for highly, mildly and negligibly dispersed fibres, respectively. Figures 10A and 10B show the ⃗ ⃗ dependency of α1 and β1 as a function of TEmax for the ex vivo data and the in silico data, for the negligibly dispersed fibres (i.e., κ ≥ 2.5). Figures 11A and 11B show the corresponding nRMSD (Equation 7a) for both parameters at different TEmax, while Figure 11C shows the difference between both nRMSD (ΔnRMSD, Equation 7b). Note that only negligibly dispersed fibres and the in silico data at g-ratio of 0.8 are studied here, because it is known from the results in Figure 8 that those possess the smallest relative difference and sd to the model predictions.

Fourth analysis: the effect of echo time on the performance of M2
At the largest TEmax ex vivo and in silico data showed the same trend ( Figure 10, first column). M2 could greatly reduce the ⃗ ⃗ dependency of β1 when compared to the ⃗ ⃗ dependency of α1 ( Figure 11A-B): nRMSD(α1) of 15.7% (in silico) and 37.9% (ex vivo) was reduced to 3.8% (in silico) and 1.3% (ex vivo). At smaller TEmax (36 ms and smaller), M2 was less effective (Figures 10, second and third column). Even

A B
an increased ⃗ ⃗ dependency was observed for β1 when compared to α1 ( Figure 11C): ΔnRMSD = 5.3%points at 36 ms (in silico) and ΔnRMSD = 14.1%-points at 18 ms (ex vivo). Moreover for 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 blue lines in Figure 10A and 10B) and then slightly increased again. The ⃗ ⃗ dependency up to the magic angle was not observed in the in silico data at any investigated TEmax. Moreover, the ⃗ ⃗ dependency of α1 in the ex vivo data decreased with decreasing TEmax also this trend was not observable in the in silico data ( Figure 10A).

Discussion
This work quantitatively explored the efficiency of the log-quadratic model (M2) in deciphering the orientation-independent part of R2* (R2,iso*) via its linear parameter, β1, from a single-orientation multi-echo GRE (meGRE) while varying microstructural fibre properties, i.e. fibre dispersion and gratios. Our findings demonstrated that M2 was effective in estimating R2,iso* via β1 when using meGRE with long maximum echo time (TEmax ≈ 54 ms) for all investigated microscopic arrangements in both simulations and ex vivo measurements. Moreover, we confirmed our hypothesis that the fitted β1 of the meGRE signal from tissue with different g-ratios cannot be predicted using the biophysical relation in M2, which is derived from the hollow cylinder fibre model (HCFM) and neglects the contribution of the myelin water compartment. We proposed a heuristic expression that predicts β1 to be a weighted sum of the relaxation rates of the myelin and non-myelin water pools and showed that this expression can better describe the data. Using the heuristic expression we estimated the MWF and the g-ratio from the fitted β1 of the ex vivo data achieving plausible results. Lastly, we found that M2 was not capable of estimating R2,iso* correctly when using shorter maximum echo times (TEmax < 36 ms) that are typical of whole-brain in vivo protocols. We made another unexpected observation at the shortest investigated TEmax (18 ms): Here, the orientation-dependency 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 this orientation-dependency of R2* is not captured by the HCFMbased simulation used here.
Capability of M2 to estimate the angular independent β1 for varying g-ratio and fibre dispersion values M2 has the potential to estimate R2,iso* from a single-orientation meGRE via β1. By assessing the residual ⃗ ⃗ dependency of β1 we found that M2 was effective although its capability varied for different g-ratios and fibre dispersions ( Figure 6). Nevertheless, the residual ⃗ ⃗ dependency of β1 was always less than 12% even if the ⃗ ⃗ dependency of the original R2* (via the α1 parameter of M1) was up to 50% ( Figure 7A). The residual ⃗ ⃗ dependence of β1 was smallest at a g-ratio of 0.73 ( Figure 7B).
The highest performance of M2 was found for negligibly dispersed fibres at the lowest g-ratio (0.66, Figure 7C), where an original ⃗ ⃗ dependence of R2* (i.e. via the α1 parameter of M1) of almost 50% was reduced to less than 12% in the ⃗ ⃗ dependency of β1. Note that the g-ratio value at which the performance of M2 is maximal might also depend on the compartmental R2 values. In the simulations, we used the R2 values from (Dula et al., 2010b) (Table A2). It is possible that using different R2 values would result in a different g-ratio values for which the model's performance was maximal.

Assessment of the microstructural interpretability of β1
M2 is derived from the biophysical HCFM Bowtell, 2013, 2012) and thus the fitted β1 parameter can be related to microscopic tissue parameters (Equation 16b, section 9.3). We found an error of up to 85% ( Figure 8C) between the β1 obtained by fitting M2 to the in silico data and the predicted β1 parameter using the biophysical relations in M2 (Equation 3). This confirms our hypothesis that neglecting the myelin contribution in the derivation of M2 results in an invalid biophysical expression for β1. We showed that the proposed heuristic expression for β1 is better suited for biophysical interpretation than the original version from M2, resulting in a relative error that was less than 12% ( Figure 8C). The newly found heuristic expression for β1 implies that in the slow-exchange regime, the linear time-dependent component of the logarithm of the meGRE signal can be described as a sum of the relaxation rates of the myelin and non-myelin water pools weighted by their signal fractions (Equation 4). This can be understood when considering that the β1 parameter in M2 captures the contribution of the linear component of the logarithmic-signal decay, which is the equivalent to an effective mono-exponential decay. The effective relaxation rate of a mono-exponential decay, however, can be expressed as the sum of compartmental relaxation rates weighted by their corresponding signal fractions as is well-known from the fast exchange regime. In other words, we showed here that when approximating the logarithm of the signal of a two-pool model in the slowexchange regime by a second order polynomial in time, the first-order term in time captures the fastexchange regime behaviour of the signal decay whereas the second-order term in time accounts for deviations from the linear (i.e. mono-exponential) signal decay.
Note that the proposed heuristic correction does not account for the effect of fibre dispersion which might explain why the accuracy of the prediction was reduced with increasing fibre dispersion. 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 do this as well for the heuristic expression of β1.

Myelin water fraction and g-ratio estimations from ex vivo data using the heuristic expression of β1
Assuming an effective performance of M2 in estimating the angular-independent β1 and using the heuristic expression for its biophysical interpretation, MWF and g-ratio can be estimated from the fitted β1 of the ex vivo data (Equation 5). For the negligibly dispersed fibres we found a median (across orientation) MWF of 0.14 ( Figure 9B), which is congruent with the mean value reported in white matter of 0.10 (Uddin et al., 2019). By using the FVF and proton density values from the in silico data (Table  A1), we found a median g-ratio of 0.79. The estimated g-ratio value is higher than typical MRI-based gratio values reported for the in vivo brain, which ranges between 0.65 and 0.70 (Berman et al., 2018;Emmenegger et al., 2021;Stikov et al., 2015) but is close to the value used in the work by (Wharton and Bowtell, 2013). The reason for the dissimilarity between predicted g-ratio and its counterpart from literature might be related to the additional assumptions that were made to estimate the g-ratio: while the estimation of MWF only requires knowledge of the compartmental R2 values, the g-ratio estimation requires additional knowledge of fibre volume fraction (FVF) and proton density values (Equation 6).
In this study, we used the same parameters as reported in (Wharton and Bowtell, 2013) for the FVF, and the proton densities, whereas the compartmental R2 values were based on (Dula et al., 2010b). Particularly, the value employed for the FVF (0.5, Table A1, section 9.4) is considerably lower than the values reported in literature, e.g. 0.75 in Stikov et al., 2011, which might explain the similarity between our reported g-ratio value and the one from (Wharton and Bowtell, 2013). Furthermore, we found that increasing the amount of fibre dispersion leads to a substantial underestimation of the MWF and an overestimation of the g-ratio. This might be of practical importance when using this method for MWF or g-ratio estimation.

The effect of echo time on the performance of M2
We found that M2 is less effective in estimating R2,iso* when the maximum TE is reduced to values more typically used for in vivo studies (i.e., maximal TE of 18 ms) but also at intermittent echo time ranges (36 ms, Figure 11). This observation is at first glance in contradiction with the validity range of M2 because the mathematical approximations, especially the approximation of the dephasing component of the extra-cellular compartment DE (section 9.2, Figure A1), are valid only in the regime TE ≤ 36 ms (see section 9.3). This apparent contradiction can be resolved by acknowledging the contribution of the myelin compartment signal which is neglected in M2 but non-negligible in the short TE regime. For in vivo application of M2, new in vivo meGRE protocols need to be developed that allow for a longer maximal TE (e.g. ≥ 54 ms), where sufficient data points with dominant extra-cellular signal compartment are collected. However, estimating robust in vivo R2,iso* maps at these large TEmax values requires the correction of motion artefacts (Magerkurth et al., 2011), which is an interesting future project for itself.
Interestingly, the biggest discrepancy between in silico and ex vivo results for β1 was seen for the smallest maximal TE value at ⃗ ⃗ smaller than the magic angle (55°, Figure 10B). This is because β1 and α1 of the measured ex vivo data showed in this ⃗ ⃗ range an atypical ⃗ ⃗ dependence: they decreased as a function of increasing ⃗ ⃗ up to the magic angle. One might suspect that this atypical ⃗ ⃗ dependence of β1 has been artificially introduced by the higher-order model M2. However, the fact that it was also found for α1 makes it more likely that we have found a new orientation dependence of R2* that cannot be explained by the HCFM. A mechanism that could explain a reduction in R2* at the magic angle would be shortening of R2 due to the Magic Angle Effect in highly structured molecules like myelin sheaths (see (Bydder et al., 2007)). Since this phenomenon would be superimposed on the orientation dependence of R2* being investigated here, it may be particularly evident when the latter effect is negligible, i.e. at low ⃗ ⃗ .

Considerations
Two of our findings might appear contradictory, at first glance: On the one hand, M2 can effectively capture the orientation independent component of R2* (i.e. estimate R2,iso*) at long maximal TEs (~ 54 ms), indicating that the myelin-water pool can be neglected in this TE range. On the other hand, M2 cannot predict the fitted β1 whereas a heuristic expression that incorporates the contribution of the myelin water in β1 substantially improves the prediction. This contradiction can be resolved when considering the time-dependency of the orientation dependent and independent parameters of M2, i.e. β2 and β1 respectively (Equation 1). The β2 parameter scales with the square of time and thus captures the logarithm of the signal-decay at the higher-TE time-points where the contribution of the myelin water is negligible, whereas β1 scales linearly with time and thus captures the logarithm of the signal-decay at the smaller-TE time-points where the contribution of the myelin water cannot be neglected. Consequently, M2 can effectively separate-out the orientation dependency of R2* but at the same time might fail to predict β1 accurately. Future studies should aim to find a better derivation of M2 from the HCFM that does not neglect the contribution of the myelin water. From the perspective of interpretation, our heuristic expression of β1 might be particularly helpful because we found that it is a good proxy for the fitted β1.
To generate the in silico data, we employed simplifications to the HCFM when extended to multiple cylinders contained in an MRI volume with varying degree of dispersion. We assumed that the signal coming from multiple dispersed hollow cylinders is a super-position of the complex signal of multiple single hollow cylinders with a specific orientation to an arbitrary main orientation of the fibres. As a result, the near-field interaction of the cylinders was neglected. Moreover, the dephasing due to the myelin compartment was also assumed to be negligible (i.e., DM ≈ 0). Nevertheless, the in silico data described the ⃗ ⃗ dependence of α1 and β1 as in the ex vivo data. This is seen across all dispersion regimes when using the long maximal TE protocol. As compared to previous studies where the dephasing process was more faithfully described in two dimensions (Hédouin et al., 2021;Xu et al., 2018), our model allowed for better control over the fibre dispersion in three dimensions via the Watson distribution parameter κ. Future work should investigate whether the validity of the in silico data could be improved by combining the approach of (Hédouin et al., 2021) in a three-dimensional simulation environment where the degree of fibre dispersion can be changed as well.
The ex vivo data possess two issues warranting discussion: (i) the coregistration of the diffusion results from the NODDI model into the GRE space and (ii) the use of the Watson dispersion from the NODDI model as a descriptor of the different fibre arrangements in the brain. Regarding the coregistration of the diffusion results from NODDI (see section 3.1.4), image interpolation could introduce a bias on the κ and, especially, ⃗ ⃗ maps. The latter can be more affected in areas with strong angular gradients (e.g. a 0°-90° between two adjacent voxels), resulting in local over-or underestimation of the ⃗ ⃗ values. One solution is to ensure the acquisition of both MRI techniques (R2*w GRE and dMRI) occur with the specimen in the same MR system, with identical positioning, field-ofview and image resolution. For our study, this was not possible because we used a preclinical MR system to acquire the high-resolution diffusion MRI data whereas the meGRE data were acquired on a human 7 T system. Regarding the use of Watson dispersion from the NODDI model, this distribution cannot describe all existing fibre arrangements in the brain accurately, e.g., the crossing fibre arrangement. In the optic chiasm specimen such 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 of the multi-echo gradient-recall echo (meGRE) signal can effectively estimate the fibre-angular-orientation independent part of R2* (R2,iso*) for varying g-ratio values and fibre dispersions. Thus, it provides an attractive alternative to standard methods for deciphering the orientation-dependence of R2* that requires multiple acquisition with distinct positioning of the sample in the head-coil. Doing so would provide a more robust marker for neuroscientific studies in a broadly accessible manner. We also showed that the estimated linear time-dependent parameter of M2, β1, can be used to estimate the myelin water fraction (MWF) and g-ratio using a newly proposed heuristic expression relating β1 to microstructural tissue parameters including the myelin water signal. Importantly, we found that the proxy of R2,iso* cannot be estimated effectively with the log-quadratic model at lower echo time ranges (i.e. at maximal echo times smaller than 36 ms) that are typically used for whole-brain in vivo meGRE experiments. To make M2 usable for in vivo applications, future studies need to develop new meGRE protocols with longer TEmax (54 ms) that remain time efficient and motion robust. Finally, at echo time ranges smaller than 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. Our HCFM based simulations were not able to model this angular dependence, which points towards a distinct mechanism in white matter that cannot be explained by the HCFM.  Yablonskiy and Haacke, (1994) proposed the analytical expression for the magnetic dephasing of a medium due to the presence of cylindrical dephasors (defined as cylinders with a different magnetic susceptibility than the medium) oriented with an angle ⃗ ⃗ to 0 ⃗⃗⃗⃗ defined as:

Acknowledgment
Where Vc is the cylinder's volume (equal to the fibre volume fraction, FVF), J0 is the zerothorder Bessel's function of the First Kind, u is the variable of integration and ωE is the frequency offset in the extracellular space. The latter is defined as: Where χE is the mean susceptibility of the myelin sheet, defined as (χI + 0.25 χA)(1 -g 2 ratio). In their work, Equation A3 was approximated for two-time scales divided by the so-called critical time (α in (Wharton and Bowtell, 2013)), defined as: For times shorter than the critical time, the dephasing function is approximated by a quadratic function, while for times longer than the critical time this function becomes linear. The corresponding analytical expressions (Yablonskiy and Haacke, 1994) are: Where and are expressions having all the parameters that are not time dependent, including sin 4 ( ⃗ ⃗ ) and sin 2 ( ⃗ ⃗ ), respectively. This simplified expression, especially the quadratic approximation, is used later (section 9.3). However, this piecewise approximation has a discontinuity at this critical time, as observed in Figure A1. To avoid this discontinuity when DE overpasses the critical time for the in silico data, we used an analytical solution to Equation A4. This solution was performed in Mathematica 12 (Wolfram Research, Inc., Champaign, IL (2020)), giving the following expression: In where J1 is the first-order Bessel's function of the First-Kind, and H0 and H1 are the zeroth and first-order Struve functions ( (Struve, 1882) and (Aarts and Janssen, 2016)), respectively. The offset frequency in the extracellular space (ωE) is dependent on the mean angular orientation and the g-ratio, as defined in (Yablonskiy and Haacke, 1994) and (Wharton and Bowtell, 2012).

, exp(-DE) in Equation A2b) was evaluated in function of time (in ms) and at three different angular orientations (0°, 60° and 90°). Two expressions for the DE function (Equation A4) were used: the analytical solution given in Equation A8
(Integrated DE, blue curve) and the piece-wise approximation proposed in the work of Yablonskiy et al. 1994

in Equation A7
(approximated DE, orange curve). Both functions were evaluated using the simulation values (section 9.4, Tables A1 to A3).

Analytical interpretation of the log -quadratic model (M2) and approximation with myelin compartment added
The log-quadratic model (M2) is derived from the signal equation from the HCFM with neglected myelin-water signal (SM, Equation A4). This signal is neglected due to short T2 and small volume of this compartment (see Wharton and Bowtell, 2013). The magnitude of the remaining signal of the nonmyelin compartments (SN) is defined as: , where and are the real and imaginary components of SN and SN is defined as follows: Evaluating Equation A10 with Equations A2a-b resulted in: | | ≈ √ 2 2 −2 2 + 2 2 −2 2 −2 + 2 −( 2 + 2 ) ( ) Using the natural logarithm function (ln(x)) of the above equation results in: (| |) ≈ 1 2 ( 2 2 −2 2 + 2 2 −2 2 −2 This expression can be linearised if the three functions related to time, i.e. the transverse relaxation rates (e.g. − 2 ), the frequency offset of the intra-axonal compartment (cos(ωAt)) and dephasing of the extra-axonal compartment (DE), are sufficiently small to be approximated using the 1 st and 2 nd order of the Taylor expansion, respectively, as follows: If these conditions are fulfilled, the logarithm function of Equation A12 can be approximated by a 2 nd order Taylor expansion in time, resulting in: If the quadratic approximation for DE is used (DE = t 2 , Equation A7), this expression can be summarized as: In the scenario where R2A is equal to R2E, the analytical expression for β1 becomes 1, ) turned out to explain better the in silico fitted 1 than the 1, (Equation 3 and A16b, see Figure 7), the validity range of the second-order approximation of the entire signal SC with the added myelin compartment is highly restrictive as a function of time and cannot be used for the experimental parameters used here (data not shown).
Equation A17 can be re-written as a function of the myelin water fraction (MWF), axonal water fraction (AWF) and extra-axonal water fraction (EWF), defined as: . Since the sum of the water fractions are equal to 1, β1 becomes: Where Equation 4 (or A19b) is obtained if we assume in Equation A19a that the relaxation rate in the intra-and extra-cellular water is the same: 2 = 2 = 2 .
9.4. In silico data setup: simulation parameters, SNR and anisotropic binning The in silico MR data was simulated using the HCFM for each hollow cylinder. The fixed parameters were obtained from (Wharton and Bowtell, 2013) (Wharton and Bowtell, 2013) in section 3.2. *Proton densities were scaled by a factor of 5000 but they kept the same proton density proportion between the nonmyelinated and myelinated compartments (1:0.7).
Other fixed parameters were obtained from (Dula et al., 2010b) and they are listed as follows:   (Emmenegger et al., 2021) and (Wharton and Bowtell, 2013), respectively. The mean value of 0.73 was arbitrarily defined.
To make the in silico as-similar-as possible to the ex vivo data, noise was added to the signal decay of the in silico data, in such a way that the in silico SNR is like the SNR seen in the ex vivo GRE data. For that, the ex vivo SNR was calculated by dividing the signal of the white matter region of the OC and the standard deviation of the background in (image) magnitude space. No noise correlation correction was performed in this calculation given the coil having 2 receiver channels. As a result, an average SNR value of 112 was obtained for this region (section 3.2), and with this SNR value, a complex random Gaussian noise was added to the in silico data as follows: Where N(0,σ) is the Normal distribution with mean 0 and the standard deviation defined by: Where the magnitude signal is divided by the desired SNR at time 0 (|Ssilico(t = 0)|).
With the noise added, the magnitude of the in silico MR signal at SNR = 112 was obtained: To compare the in silico data analysis across the 5000 signal decays per simulated g-ratio, sampled κ and ⃗ ⃗ to the irregularly binned ex vivo data analysis (section 3.3.1 and Figure 5B), the αparameters and β-parameters from the in silico data required three consecutive averaging-steps: (1) an averaging across the 5000 samples, resulting in the sampled-averaged ̂( : 0,1), ̂( : 0,1,2) and their standard deviations ( ), ( ) per sampled κ value and ⃗ ⃗ .
(2) A weighted averaging across κ values per each ⃗ ⃗ irregular bin of the ex vivo data in each κ range. For that, it was obtained the distribution of the κ values from the voxels contained in each of the 20 defined ⃗ ⃗ irregular bins. The ⃗ ⃗ range per bin and κ range is given in Table A.4. Then, all the obtained distributions were averaged per κ range ( Figure A2 from A to C) to remove possible influence of the irregular ⃗ ⃗ bins on κ. The standard deviation from this average was calculated, normalised and used later (referred as the sd( ( )) in Equation A28). Next, a probability distribution, ( ), was fitted accordingly ( Figure A2 from D to F) and the weighted averaging on ̂ (the same procedure is performed for ̂) was calculated as follows: Where the expression for ( ) was heuristically chosen and varied per each fibre dispersion (κ range): a Beta distribution for the highly dispersed fibres (κ < 1 range, Figure A2-D), defined as: Is the Gamma function. The coefficients a and b estimated for this range were 3.145 and 1.234, respectively. Given the clear half-shaped normal distribution, a Half-Normal distribution for the mildly dispersed fibres (1 ≤ κ < 2.5 range, Figure A2-E) was used, defined as: (1 ≤ < 2.5) = √2 √ ( −( − 1) 2 2 2 ) (A26) . The coefficients μ and σ were 0 and 0.4498, respectively. And given the fast decay of the values at the beginning of the distribution, an Exponential distribution for the highly aligned fibres (2.5 ≤ κ range, Figure A2-F) was used, defined as: (2.5 ≤ ) = (− ( − 2.5)) (A27) . The coefficient λ was 0.2241. The standard deviation of ⟨ ⟩ was also estimated by error-propagating the ( ) weighted by ( ) and its standard deviation sd( ( )), as follows: (⟨ ⟩ ) = √ ∑ ( ( ( )) ( ) ∑ ( ) ) 2 + (̂( ) ( ( ))) 2 (A28) While the first squared term requires the normalisation factor (∑ ( )) because the weights ( ) are not normalised, the second is not needed since sd( ( )) is already normalised. Finally, the ⟨ ⟩ and (⟨ ⟩ ), and the ⟨ ⟩ and (⟨ ⟩ ) (as in Equation A28) were averaged and errorpropagated, respectively, as a function of the ⃗ ⃗ values for each irregular bin. Figure A 2: Assembling the in silico data across the simulated κ ranges and angular ( ⃗ ⃗ ) anisotropic bins. To make the in silico data comparable to the ex vivo data, the frequency of voxels as a function of κ was obtained per defined ⃗ ⃗ irregular bin in the ex vivo data (Figure 3). This was performed for the 20 ⃗ ⃗ bins (AnisoBin X, with X the corresponding bin from 1 to 20, see Table A.4) and per fibre dispersion (κ range): highly dispersed (κ < 1, A), mildly dispersed (1 ≤ κ < 2.5, B) and negligibly dispersed (κ ≥ 2.5, C) fibres. The mean and standard deviation across histograms were obtained (error bars). The means were normalised with respect to the cumulated value (i.e., sum of all the mean values) and fitted with a continuous function ( ( )) per κ range, previously normalised: a beta distribution for κ < 1 (D), half-normal distribution for 1 ≤ κ < 2.5 (E) and exponential distribution for κ ≥ 2.5 (F). The standard deviation was also normalised by the cumulated value per κ range and used as the standard deviation of the continuous distributions (sd( ( ))).