Multi-compartment relaxometry and diffusion informed myelin water imaging – Promises and challenges of new gradient echo myelin water imaging methods

Myelin water fraction (MWF) mapping based on data fitting of a 3-pool exponential model of multi-echo gradient echo (mGRE) data using MRI shows great promises for in vivo myelin quantification. However, this multi-exponential fitting is ill-conditioned because of the similar relaxation times and frequency shifts of the various compartments. Additionally, the bound water residing in the myelin sheath of white matter is expected to have a faster longitudinal magnetisation recovery than that of the free water in the axonal and extracellular space. When the Ernst angle is used to achieve maximum SNR and improve fitting, this will introduce a T1-weighting effect to the derived MWF. In this study, we first demonstrate that diffusion-weighted imaging can be used to infer the compartmental signal properties using an analytical fibre model to achieve a robust MWF estimation. Second, we show that by incorporating a variable flip angle scheme to the mGRE acquisition with a multi-compartment relaxometry model, not only the MWF is corrected from the T1 dependency but also the fitting procedure becomes less ill-conditioned and more SNR efficient. Finally, we demonstrate these two approaches can be combined allowing higher spatial resolution MWF maps than what has been reported to date with robust MWF estimation on a small cohort.


Introduction
Myelin water imaging (MWI) is a promising MRI tool for indirect quantification of myelin content, which can find applications in studying brain plasticity (Lebel and Deoni, 2018) and neurodegenerative diseases such as multiple sclerosis (Kolind et al., 2012;Laule et al., 2004;Tozer et al., 2005;Vargas et al., 2015;Vavasour et al., 1998). It was first implemented using a multi-echo spin echo (MSE) sequence (MacKay et al., 1994), from which the multi-component nature of the T 2 distribution was analysed. Myelin water (MW), representing the water molecules residing in-between the intra-phospholipid layers of the myelin sheath, is regarded as the short yet detectable T 2 components (T 2 < 25 ms), while the long T 2 components represent the free water molecules in the intra-axonal and extra-axonal space. Good correlations between the myelin water fraction (MWF) estimated with MSE and the myelin concentration were found in histological studies (Laule et al., 2016(Laule et al., , 2008), yet the high-energy deposition and long scan times when attempting whole-brain coverage have prevented its wider dissemination. Image acquisition developments over the last decade using gradient and spin echo (GRASE) acquisitions have allowed circumventing the whole volume coverage and long acquisition time limitations (Alonso-Ortiz et al., 2015;Lee et al., 2020;Prasloski et al., 2012), although at the cost of very anisotropic voxels (typically~5 mm slice thickness).
Multi-compartment modelling using multi-echo gradient echo (mGRE) imaging provides an alternative method for MWI (Du et al., 2007;Gelderen et al., 2011;Hwang et al., 2010;Nam et al., 2015b) that has low RF power deposition, high scan time efficiency and isotropic resolution . Good correlation has been found between the spin echo and the gradient echo methods when doing a region of interest (ROI) analysis that included various grey matter (GM) and white matter (WM) regions (Alonso-Ortiz et al., 2017). In GRE MWI, WM signal is described as a 3-pool system, with the signal originating from the water in the myelin sheath, axon and extra-axonal space respectively. Similar to the MSE approach, MW has the shortest apparent transverse relaxation time T 2 *(¼1/R 2 *) among the three compartments, and additionally each compartment has a distinct frequency shift as predicted by the hollow cylinder fibre model (HCM) (Sati et al., 2013;Wharton and Bowtell, 2012). Fitting such a 3-pool model is challenging for various reasons. The difference in frequency shifts between the compartments is small when compared to the signal decay rates. MW signal not only decays quickly but also accounts for less than 20% of the total signal (Lee et al., 2016), making data acquired at short echo times with high SNR crucial for the measurement. Additionally, intra-axonal water (IW) and extra-axonal water (EW) have relatively slow and similar R 2 *s, and their frequency difference is relatively small compared to MW.
The semi-solid nature of the myelin sheath also impacts the (apparent) longitudinal relaxation rate R 1 (¼1/T 1 ) of MW (Gelderen et al., 2016;Gelderen and Duyn, 2019), making it significantly faster than that of the free water. This phenomenon was clearly illustrated in a double inversion experiment aimed at suppressing signals with slow longitudinal recovery R 1 (with T 1 range of 750-2000 ms) (Oh et al., 2013). In that context, it showed that the remaining fast longitudinal recovering signal had similar T 2 * distribution to that in MW (Oh et al., 2013). The Ernst angle and short repetition time (TR) are often used to achieve high SNR and scan time efficiency for GRE experiments. However, given the discrepancy of the signal recovery rates between various compartments, the measured compartmental signal intensities can have different steady-state saturation, making MWF subject to a T 1 -weighting effect (Lee et al., 2017a), as illustrated in Fig. 1.
Recent advances in diffusion-weighted imaging (DWI) acquisition strategies, such as simultaneous multi-slice imaging which can accelerate the acquisition of multiple b-shell schemes, and WM signal modelling allow estimations of not only diffusion properties (such as apparent diffusion coefficient and fibre orientation) but also multi-compartment information such as the volume fraction of different water pools (Kaden et al., 2016;Zhang et al., 2012). WM fibre orientation affects the frequency shifts observed in the 3 different compartments modelled in MWI (Sati et al., 2013;Bowtell, 2013, 2012). In the past, we used the richness of the WM contrast in GRE signal to improve DWI fibre tractography (Chan et al., 2018) but to limited gains. Now, noting that in MWI, the signal ratio between intra-and extra-axonal spaces is one of the least robust fitting results (Lee et al., 2016) and we propose that the information obtained from DWI, although not being specific to myelin, has the potential to improve the GRE MWI results.
In this study, we tested two pathways to improve model-based GRE MWI. First, we investigated the feasibility of integrating DWI information to the MWI signal model thus improving the robustness of MWI. This approach can be particularly advantageous if DWI data is already present in the imaging protocol. In this case, the gradient echo add-on can be particularly short compared to the standard MWI. As a second pathway, by using a variable flip angle (VFA) acquisition strategy that changes the weighting of the various compartments, the fitting was made more robust while additionally correcting the T 1 -dependency in the MWF estimation and providing various metrics that reflect the micro-environment of MW.

Model 1: Standard myelin water imaging
Standard GRE MWI considers 3 water pools (MW, IW and EW) in WM. Each water pool has a distinct signal amplitude A, a transverse relaxation decay rate R * 2 , and a frequency shift ω (Gelderen et al., 2011;Nam et al., 2015b): where ω b is a background frequency offset and ϕ is an initial phase that comes from B 1 þ phase offset (Kim et al., 2014;Schweser et al., 2011).
Subsequently, by fitting the complex-valued mGRE data to Eq. (1), myelin concentration can be indirectly quantified as the MWF, which is a ratio of the MW signal to the total water signal.

Model 2: Diffusion informed myelin water imaging (DIMWI)
Since myelinated axons are surrounded by multiple phospholipid layers composed of highly organised radially-oriented molecules with anisotropic magnetic susceptibility, the contrast observed within WM is fibre orientation dependent (Wharton and Bowtell, 2012;Yablonskiy and Sukstanskii, 2015). Such observation can be explained by the HCM proposed by Bowtell, 2013, 2012). By approximating a myelinated axon as an infinite long hollow cylinder, the average frequency offsets induced in the myelin and axonal spaces can be characterised by the isotropic and anisotropic magnetic susceptibility of the myelin sheath (χ I and χ A ), and are given by: Due to the fast R 1 of MW, the change of the MW signal is less sharp than that of the free water compartment, resulting in a shift of the Ernst angle of the total signal to a higher value. (b) The logarithm of the simulated multi-echo GRE signal from the 3-pool MWI model as a function of echo time and flip angle when the R 1 of the fast MW R 2 * compartment assumed to be faster than that of free water. Dashed lines represent the single exponential fit using only the later echoes when the signal of the MW is (nearly) fully relaxed. The increasing discrepancy between the dotted and solid lines with the increase of flip angle is caused by the T 1 -weighting effect on the apparent MWF in the signal.
where γ is the proton gyromagnetic ratio; g is the ratio between inner and outer axonal radii, which can be obtained from the relationship g ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jung et al., 2018;Stikov et al., 2015) and ρ MW is the relative MW density; c ¼ 1=4 À 3g 2 lnð1 =gÞ=ð2ð1 À g 2 ÞÞ; θ is the angle between the fibre and the main magnetic field (B 0 ) directions and can be obtained from DWI; E describes the contribution of chemical exchange within the myelin sheath to the frequency offset (Wharton and Bowtell, 2013). Therefore, if the myelin sheath properties (i.e. water density, magnetic susceptibility and exchange) can be assumed to be relatively constant throughout the brain then, given the fibre orientation from DWI, the frequency shifts of IW and MW can be computed directly from the standard MWI model parameters while the EW shift is set to be zero, as suggested both analytically and numerically by the HCM (Sati et al., 2013;Wharton and Bowtell, 2012).
HCM simulations in previous studies (Lee et al., 2016;Sati et al., 2013) have shown a narrow frequency dispersion inside the axon (full width of half maximum of the spectral line ≪ 0.05 ppm), which has led us to assume the R 2 * of IW can be approximated by its R 2 . The signal decay in EW can be modelled as the sum of the R 2 of IW and an extra dephasing effect, D E , caused by the local inhomogeneous field in the extra-axonal space at time t (Wharton and Bowtell, 2013): Àg 2 Þ is the effective susceptibility of the nerve fibre in the HCM, α ¼ 3=ðjχ D jγB 0 sin 2 θÞ is the transition time of dephasing from the quadratic regime to the linear regime and FVF is the fibre volume fraction and can be derived from FVF ¼ A IW =ðA IW þ A EW g 2 Þ.
Recent development in DWI acquisition strategies and signal modelling can provide more detailed microstructural information. Utilising a multi-shell acquisition scheme, techniques such as neurite orientation dispersion and density imaging (NODDI) (Zhang et al., 2012) and spherical mean technique (SMT) (Kaden et al., 2016) can infer microstructural information such as fibre dispersion and volume fractions between various water compartments. In particular, the volume fraction of intracellular (intra-axonal) water v ic , which is a volume ratio between intra-cellular and extracellular space, is directly related to MWI. Incorporating the DWI information to the standard MWI model, the multi-echo GRE signal, S DIMWI , with a given fibre orientation (θ) can be written as: While standard GRE MWI can be seen as a multi-compartment transverse relaxometry fitting problem, it does not cover the whole range of relaxation properties and mechanisms of these compartments. With a VFA acquisition scheme (Deoni et al., 2003), the T 1 -weighted steady-state of a spoiled GRE signal intensity of individual compartment N (N2{MW, IW, EW}) can be modelled as: where M 0;N is the thermal equilibrium magnetisation of compartment N, which is directly proportional to its proton density, and α is the flip angle of the GRE sequence. This equation can be incorporated in the standard MWI model of Eq. (1) to replace the compartmental signal amplitudes A MW , A IW and A EW . In this way, the proton density of each compartment and their T 1 weightings can be decoupled by changing the acquisition parameters. Consequently, true MWF can be computed without any T 1weighting effect with an additional longitudinal relaxation map of the various compartments.
Preliminary analysis and simulations (data not shown) suggested that imperfect gradient and RF spoiling can have large effects on the measured signal when large flip angles are used. Furthermore, intercompartmental magnetisation exchange cannot be ignored in the steady-state calculation as has been extensively reported in the literature Ou and Gochberg, 2008;Zhang et al., 2015). We also found that the R 1 of the IW cannot be distinguished from that of the EW on the GRE data. To account for these three observations, we considered a two-pool exchange system consisting of a myelin sheath pool (sum of myelin macromolecules and water molecules) which has a fast R 1 (R 1;Myelin ), and a free water pool (IEW¼IW þ EW) which has a slow R 1 (R 1;IEW ). The concept of magnetisation exchange between the myelin sheath and the free water protons is included in the Bloch-McConnell equations governing the longitudinal magnetisation M z (Spencer and Fishbein, 2000). In the absence of RF pulses, they can be written as: where k IEWM and k MIEW are the directional exchange rates between the free water pool and the myelin sheath pool. In thermal equilibrium, the exchange rates are constrained by k IEWM M 0;IEW ¼ k MIEW M 0;Myelin . By solving Eq. (7) and considering the RF pulse with flip angle α, applied at every TR, the steady-state signal of the myelin sheath and free water (A Myelin and A IEW ) right after the excitation pulse is given by: A Myelin ðα; TRÞ A IEW ðα; TRÞ ¼ sin αðI À e ΛL TR cos αÞ where I is an identity matrix. Eq. (8) can be computed with the consideration of incomplete RF spoiling using the full EPG-X formalisms described in (Malik et al., 2017). To be able to incorporate the EPG-X framework for MWI, the following assumptions were made in the proposed model: (1) M 0;Myelin is approximated by the M 0;Myelin ¼ M 0;MW =ρ MW , (2) M 0;IEW ¼ M 0;IW þ M 0;EW , (3) R 2 of the free water compartment is approximated by IW R 2 * based on the observation of narrow IW frequency spectrum in HCM simulations (Lee et al., 2016;Sati et al., 2013) and is a free parameter in data fitting, (4) R 2 of the myelin sheath is approximated by MW R 2 *, its fast relaxation (when compared to TR) minimises its impact on the EPG-X simulated signal and is also a free parameter to be estimated, and (5) the frequency offset between free water and myelin sheath is approximated by the frequency offset between the IW and MW based on the zero frequency offset of EW and the frequency difference between IW and EW is small (Lee et al., 2016;Sati et al., 2013;Wharton and Bowtell, 2012 The multi-compartment relaxometry approach provides a framework to model the steady state of the multi-compartment signals, governed by their T 1 s and exchange processes, but is independent of the signal decay across echoes. The framework can, therefore, be expanded to include the microstructure information derived from multi-shell DWI data when it is available, through the incorporation of Eqs. (8) and (9) to Eqs. (2)-(5) from the DIMWI model, achieving further reduction on fitting parameters.
In the study conducted by Nam et al. (2015b), the frequency shift of each water pool was combined with the background frequency offset, resulting in 10 parameters to be estimated in the standard MWI (standard). For the diffusion informed model (DIMWI), the frequency shifts of IW and MW, the signal decay for inhomogeneous field dephasing on EW and the volume fraction of IW are derived from DWI, so that the number of fitting parameters is reduced to 6. For the multi-compartment relaxometry model (MCR), one R 1 (the myelin R 1 was subsequently fixed in the following processing; see Methods section) and one exchange constant have to be estimated simultaneously with the 8 compartmental properties (3x proton density, 3x R 2 * and 2x frequency shift) and the B 1 þ phase offset, in addition to an individual background frequency component for the data of each flip angle to account for the small variation of head position or scanner drifts. This results in 10 tissue properties, plus 1þ number of acquisition parameters. With the diffusion-derived microstructure properties, the number of fitting parameters of the MCR-DIMWI models can be reduced by 4. The number of fitted parameters, model assumptions together with the subsequent processing steps for each model are further summarised in Fig. 2.

Material and methods
MWF is the key parameter in MWI, thus the robustness of estimating this parameter from the proposed signal models was studied using numerical simulations, including the consideration of the model assumptions deviating from the ground truth. MWF maps generated from various methods were also evaluated using in vivo imaging data in a small cohort.

Simulation 1: Search space behaviour of the MCR model
Incorporating inter-compartmental exchange in multi-compartment signal modelling, which is also considered in the MCR model, can lead to inaccurate parameter estimation due to the non-convex nature of the fitting residual (West et al., 2019). To thoroughly understand how the exchange term can affect the cost function of the MCR model, a noise-free signal was simulated using the MCR model (Eqs. (1), (8) and (9)) based on the tissue properties and acquisition scheme shown in Tables 1 and 2. A three-dimensional search grid was created comprising the R 1 of the free (intra-axonal and extra-axonal) water (T 1 range from 800 to 1500 ms, with a spacing of 25 ms), R 1 of the myelin sheath (T 1 range from 100 to 500 ms, with a spacing of 10 ms) and the magnetisation exchange rate from the free water to the myelin sheath (range from 0 to 6 s À1 , with a spacing of 0.5). Data fitting was subsequently conducted using the noise-free signal as an input to the MCR model to obtain the minimum cost of each candidate in the search space, by fixing the T 1 s and exchange rate to the values in the search grid and allowing the rest of the parameters in the model to be estimated in data fitting. The cost function was then visualised as a normalised root-mean-square residual η, similar to the approach in (West et al., 2019): where N is the number of the data points, S j ðξÞ is the signal in the jth data for a parameter set ξ in the search space, b ξ is the set of parameters used to simulate the signal, and σ is the notional standard deviation of the noise and was set to the value such that the mean SNR at TE ¼ 0 for all the flip angles in the acquisition scheme is 200, regarded as a high SNR condition. When η is equal or smaller than 1, it means that the residual of the solution cannot be differentiated from the noise if there was noise present in the data (the simulation was noise-free). If a model has a unique solution in the search space, the area with η 1 will be confined in a narrow range of parameters in the grid that is close to the solution.

Monte Carlo simulations
3.1.2.1. Simulation 2: Impact of fitting tolerance in MWF. Monte Carlo simulation is a powerful way to estimate the accuracy and precision of the various models. To avoid data overfitting, an early stopping approach was implemented in a previous study via setting a relatively large fitting tolerance value (Nam et al., 2015b), assuming the initial starting parameter set is reasonably close to the actual solution. To study the fitting condition of the DIMWI model, Monte-Carlo simulations were conducted to simulate the 3-pool WM signal using the HCM formalism (Wharton and Bowtell, 2013) with a clinically feasible mGRE acquisition scheme summarised in Tables 1 and 2 Random Gaussian noise was added to 4000 replicas of the data to simulate SNR(TE ¼ 0) of 150 and fitted with the standard and DIMWI models to study the bias and precision of MWF estimation using the lsqnonlin function in Matlab (MathsWork, Natick, USA) at three different tolerance levels (both StepTolerance and FunctionTolerance options) corresponding to 10 À5 , 10 À6 and 10 À7 , respectively.  (Stikov et al., 2015). b . c (Wharton and Bowtell, 2012). d (Nam et al., 2015b). e (Du et al., 2014). f . 3.1.2.2. Simulation 3: Model performance to reflect a range of MWF. In practice, it is important to distinguish the differences in MWF with the change of WM microstructure. By using the extreme FVF and g-ratio values in 6 different WM fibre tracts previously reported , mGRE signal was simulated based on the HCM formalism with T 1 -weighting (MCR-DIMWI) corresponding to 4 MWF conditions: For standard and DIMWI fitting, the signal was simulated from an Ernst angle acquisition scheme illustrated in Table 2, with random Gaussian noise added to the 4000 replicas of the data to simulate SNR(TE ¼ 0) of 150. The same noise level was used in MCR and MCR-DIMWI fitting with signal simulated based on the VFA scheme in Table 2. Mean and standard deviation of the fitted MWF were derived to evaluate the ability of the signal models for distinguishing the MWF of different WM microstructures. The mean MWF of all conditions were also fitted to a linear function of y ¼ mx þ c to determine the linearity of the estimated MWF.
3.1.2.3. Simulation 4: Realistic consideration of the fixed myelin sheath properties in the proposed models. To be able to integrate HCM without simply adding new fitting parameters in DIMWI, the isotropic and anisotropic magnetic susceptibility of the myelin sheath, the relative MW density and the contribution of exchange to phase were assumed to be constant which are unlikely to hold throughout the whole brain (see Theory). Therefore, it is important to understand the MWF estimation bias when these quantities deviate from the assumptions. The simulations were repeated for SNR(TE ¼ 0) of 150 while varying the values of the parameters in the simulated signal, assuming 30% underestimation and overestimation of the fixed values in our model assumptions. A scenario of EW signal possessing an independent decay rate instead of decaying with the sum of IW R 2 * and the dephasing effect D E was also considered. Multi-echo GRE signal was simulated by replacing the decay terms (R 2 * IW and D E ) of EW in Eq. [4] by an R 2 * EW , while the exact formalism of Eq.
[4] was used in the DIMWI model fitting in the Monte Carlo simulations. Three conditions were tested respectively corresponding to EW T 2 * of 34 ms, 48 ms and 62 ms.
Given the search space behaviour of the MCR model, we decided to fix the T 1 of the myelin sheath to 234 ms (Du et al., 2014) to reduce the number of possible solutions. Monte Carlo simulations were subsequently conducted to investigate the model performance when the R 1 of myelin in the simulated signal deviates from the fixed value in the MCR model, utilising the same approach as in the aforementioned DIMWI simulations.

Data processing
The MP2RAGE combined contrast image with the background being suppressed using a mask based on thresholding the image acquired at the second inversion time was processed with the anatomical processing pipeline of the FMRIB Software Library (FSL, https://fsl.fmrib.ox.ac.uk/f sl/fslwiki/fsl_anat) (Jenkinson et al., 2012), comprising bias-field correction, linear and nonlinear registration to MNI152 space, brain extraction, tissue-type segmentation and subcortical structure segmentation. The resulting tissue masks and T 1 map were subsequently registered to the GRE space.
Susceptibility-induced distortion and eddy current-induced distortion were first corrected from the DWI data using topup (Andersson et al., 2003) and eddy (Andersson and Sotiropoulos, 2016) of FSL. Accordingly, a ball-and-stick model (Hernandez et al., 2013) was used to compute 3 fibre orientations per voxel, a DTI model was used to compute the fractional anisotropy and the volume fraction of IW was estimated by the multi-compartment DWI based on SMT (Kaden et al., 2016). All results were registered to the GRE space for data fitting of the standard and the proposed MWI models.
Gibbs-ringing artefact on GRE data was first mitigated using a local subvoxel-shifts method (Kellner et al., 2016). Motion between scans was corrected by coregistration of all scans to the middle GRE acquisition using a rigid body transform. Complex-valued based fitting is sensitive to phase errors (Lee et al., 2017b). Eddy current-induced phase artefacts, which can be prevalent in the early echoes of the mGRE phase data, were corrected by removing the first-order polynomial components from the phase difference between the phase of the average of successive odd echoes and the even echoes in between (Li et al., 2015). Because of the non-linear phase behaviour expected in WM, the fitting was performed on the GM mask obtained from the MP2RAGE results. MWI data fitting was performed in a voxel-wise manner using software developed in-house. Our preliminary analysis indicated that the noise level in the GRE data increases with echo time and the fitting residual of various models is substantially smaller than the noise (see Supplementary material Fig. S1). This observation, possibly resulting from physiological noise, supports the use of protocols with relatively short echo times when compared to the measured T 2 * of the IW and EW. Similar to (Nam et al., 2015b), the cost function of the least-square fitting was weighted by the amplitude of the signal. To further mitigate the overfitting effects, an optimal tolerance value was determined for each model as described in detail on the Supplementary material (Figs. S2-5). The 7 measurements of the 1.8 mm isotropic resolution GRE data with the flip angle of 20 were averaged as proposed in (Nam et al., 2015b) and fitted to the standard MWI model. For the DIMWI model, only one GRE measurement of the flip angle of 20 was used with the DWI results. The MCR-based models were used to fit the complex-valued data from all 7 flip angles as described in the Theory section.

Data analysis
Histogram analysis was performed to study the reproducibility of the MWF estimated by various models in both WM and GM in the small cohort of 6 subjects based on 1.8 mm isotropic mGRE data, using the tissue masks obtained from the combination of the segmentation results from FAST and FIRST of FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki). The masks were further eroded by 1 voxel in all directions to reduce partial volume effect at tissue boundaries. The MWF histograms were averaged across subjects on the 1.8 mm isotropic data for all MWI methods. In addition, linear correlations were performed between the 4 MWI methods based on the 1.8 mm isotropic mGRE data to quantitatively evaluate if the different methods are indeed measuring similar distributions of myelin across the brain.
To be able to compare the various parameters derived from our methods to those found in literature, 5 fibre bundles were studied, including forceps major (FMj), forceps minor (FMn), corticospinal tract (CST, both hemispheres), cingulum (CG, both hemispheres) and inferior longitudinal fasciculus (ILF, both hemispheres), on the 1.8 mm isotropic data acquired on the 6 subjects for the 4 MWI models. ROIs of these fibre bundles were obtained from each subject using probabilistic tractography tool PROBTRACKX of FSL on the DWI data as previously reported (Gil et al., 2016), in turns, the derived ROI masks were registered to the GRE space. For each ROI and subject, the mean and standard deviation of each parameter were computed and the mean of these quantities across subjects were further evaluated.

Simulations
4.1.1. Simulation 1: Search space behaviour of the MCR model Fig. 3 shows the normalised root-mean-square residual between the simulated signal using the parameters described in Table 1 and different parameter solutions in the 3D search grid of the T 1 of myelin, T 1 of free water and exchange rate when the mean notional SNR is 200 at TE ¼ 0. The points in colour indicate the possible solutions at this SNR level. A wide range of T 1 and exchange rate combinations can minimise residual to a level that is significantly smaller than the target noise level even though the notional SNR is high: ranging from 925 ms to 1500 ms for T 1 of the free water (ground truth ¼ 1050 ms); from 100 ms to 380 ms for T 1 of the myelin (ground truth ¼ 234 ms) and from 1.5 s À1 to 5 s À1 for the exchange rate (ground truth ¼ 2 s À1 ), suggesting that the ground truth parameter set is difficult to obtain. Consequently, the estimated R 1 s and exchange rate are heavily affected by the initial starting values in data fitting.

Simulation 2: Impact of fitting tolerance in MWF
The definition of the fitting tolerance values can have a strong impact on the variability of the fitted parameters. Fig. 4 shows that when changing the fitting tolerance to smaller values, no observable change in MWF bias (deviation of the median from the ground truth) can be found from the standard model. On the other hand, the increase of interquartile range (IQR) denotes the ill-conditioning problem of the model. With the DIMWI model, MWF bias can be lowered by using a smaller fitting tolerance when the fibre is parallel to the B 0 , yet this introduces more variability (wider IQR) in the result. The MWF of perpendicular fibre obtained from DIMWI has the lowest bias and variability, demonstrating that in this case the ill-conditioned nature of the problem was overcome and the results are very consistent in all settings.

Simulation 3: Monte Carlo simulations on a range of MWF estimation
The ability of a signal model to distinguish different MWF levels is of uttermost importance in MWI. The mean MWF derived from Monte Carlo simulations suggests that all four models (standard, DIMWI, MCR and MCR-DIMWI) produced reasonably accurate estimation in a range of MWF with the slopes of the linear fit being close to one (Fig. 5). The relatively greater discrepancy between true MWF and the results of the standard model and DIMWI is expected due to the T 1 weighting effect. Significant reduction of standard deviation can be observed with DIMWI when the WM fibre is perpendicular to B 0 (see changes from Fig. 5a-d) and with the MCR-based models (see changes from Fig. 5a-c).

Simulation 4: Monte Carlo simulations on biases introduced by the fixed parameters in the signal model
In reality, the actual myelin sheath properties can deviate from the assumptions we made for the HCM in DIMWI. When the assumptions of the DIMWI model do not agree with the simulated values, DIMWI can still potentially reduce MWF bias and improve precision for fibres that are perpendicular to B 0 (Fig. 6a-e). Minor improvements in bias reduction and comparable IQR are observed with DIMWI compared to the standard result for parallel fibre, which is when the HCM derived frequency shifts and signal decay rates of the IW and EW compartments become indistinguishable.
MWF obtained from the MCR model gives a more stable measurement when the myelin R 1 is fixed regardless of fibre orientation compared to the standard model, despite the values in the simulated signal disagreeing with the fixed parameter in the model (Fig. 6f). This improvement comes at the cost of introducing systematic bias in the estimation when the actual myelin R 1 deviated from the fixed value in the MCR model.

1.8 mm isotropic resolution GRE data
MWI results were first compared across various fitting models using a GRE protocol similar to those reported in the literature in terms of echo times and resolution. The 1.8 mm isotropic resolution MWF maps obtained from DIMWI using only one GRE measurement demonstrated similar quality as the result of averaging 7 measurements using the standard model (Fig. 7). Despite the addition of the diffusion data, the overall acquisition time was still reduced. Although the DIMWI result is affected by the image artefact close to the frontal region in the GRE data, the resulting maps are generally less noisy in WM, especially for the fibre bundles perpendicular to the B 0 such as optic radiation and the splenium of the corpus callosum (green arrows). A significant improvement in MWF map quality can be obtained from the MCR model (obtained in the same total acquisition time as the standard MWF maps). The MCR results are less affected by image/fitting noise and are more robust to the image artefact in the data (red arrows). For example, the globus pallidus, a subcortical GM structure with high iron deposition (fast R 2 * decay), usually appears as myelin rich using the standard model (see Fig. 5 and (Nam et al., 2015b)) yet in MCR it appears dark (blue arrows) suggesting a lower myelin concentration as most of the remaining GM. This can also be seen in supplemental Fig. S6, where a coronal slice crossing the substantia nigra shows that standard MWI overestimates the local myelin water, while this overestimation gets corrected by all proposed methods. Combining the DIMWI and the MCR models can further improve the MWF map quality, especially in perpendicular fibre bundles (yellow arrows).
In multiple parameter fitting problems, it is important to also evaluate some of the remaining by-products of the fitting process. It is interesting to note that MWF maps are one of the most robust metrics derived by all the different approaches. Additionally, all the proposed methods produced robust IW R 2 * maps, although with faster decay rates compared to the standard model. The impact of small imaging artefacts on the GRE measurements was more pronounced on both the MW R 2 * map (Fig. 8a) and the frequency maps ( Fig. 8c and d) when using both the standard model and the DIMWI. The most affected regions were located in regions close to the frontal sinuses (highlighted with red arrows) where the frequency maps across the various acquisitions are most likely to vary as a result of small movements. It is important to note that the MW R 2 * maps from the diffusion informed models were close to the fitting boundaries in GM (will be discussed later). The MW and IW frequency shifts ( Fig. 8c and d) derived from the HCM are very similar between DIMWI and MCR-DIMWI. The standard and MCR MW frequency shifts are also similar to the HCM predictions but it is not the case for the IW frequency shift. Similar quality of free water R 1 (Fig. 8f) and exchange rate (Fig. 8g) can also be observed between the two MCR-based methods. Clearly, the exchange rate map derived from the MCR-based models has a varying bias from the anterior to the posterior part of the brain shown in Fig. 8g. This is likely related to result from the normalised RMSE behaviour shown in Fig. 3c. The derived R 1 maps have strong GM/WM contrast but are not as high SNR as a typical VFA R 1 map when the exchange is neglected. Interestingly, there is more contrast within WM than on standard single compartment R 1 maps at 3T (Gil et al., 2016;Marques et al., 2010) where WM bundles are not usually distinguishable. It is possible to appreciate many of these features also in Fig. S6.

1.4 mm isotropic resolution GRE data
Encouraged by the improved robustness to noise amplification of the proposed MWI methods, we evaluated their robustness in the context of higher resolution images. Increasing the resolution from 1.8 mm isotropic to 1.4 mm isotropic reduced the SNR of the GRE data by approximately a factor of 2 and, unsurprisingly, the MWF derived from Fig. 4. Box plots (central line: median; height: IQR; whiskers: inclusion of the most extreme data that are not considered outliers, approximately AE 2.7 times the standard deviation) of MWF estimated with the standard and DIMWI models at 3 tolerance levels (10 -5 , 10 -6 and 10 -7 ) and in 2 fibre orientations (0 and 90 ) with respect to B 0 . The IQR becomes wider for smaller tolerance values due to the ill-conditioned standard model fitting, and it does not guarantee a reduction of MWF bias. Estimation bias (deviation of the median from the true MWF) can be reduced by setting a smaller tolerance value for DIMWI, with the expense of estimation stability when the fibre is parallel to the B 0 . Accurate results with high precision (narrower IQR) can be achieved when the fibre is perpendicular to the B 0 from DIMWI for all tolerance levels.
the standard model is extremely noisy (Fig. 9, first column). Incorporating the diffusion information can significantly improve the MWF quality, leading to more distinguishable fibre tracts in WM (Fig. 9, second column). The MWF maps from the MCR model possess similar high quality as the 1.8 mm isotropic lower resolution results shown in Fig. 7  (Fig. 9, third column), and using the DIMWI formalism can further reduce the noise in the maps, resulting in a more homogeneous MWF within the same fibre bundles (Fig. 9, fourth column) and suggesting that the resolution could be further enhanced in future studies.

Histogram analysis
To study the reproducibility of the MWF result of each model across the six subjects, MWF histograms under the WM and GM masks were evaluated and are displayed in Fig. 10a (based on the 1.8 mm isotropic resolution data). The subject-specific histograms are slightly different from each other, and inter-subject MWF variation is more observable in WM for all models. At the group level, both the standard and DIMWI model generated broad MWF spectra in WM (Fig. 10b), while the MCRbased models show more distinguishable and narrower peaks in both WM and GM. This can be seen both at the subject (Fig. 10a) and group levels (Fig. 10b). The width of the histogram of the MCR-DIMWI model is clearly narrower than that of the MCR model (Fig. 10b), supporting the observation in Fig. 7 that MCR-DIMWI produced the least noisy MWF maps among all models. Finally, the smaller standard deviation (shaded region) in the group level suggests that the results are highly reproducible.
The correlation analysis between different MWI models shows the similarity of the MWF contrast observed within tissues. Note that we have separated the analysis of WM only and the sum of WM and GM. The highest correlation was found between MCR and MCR-DIMWI in both ROIs (R ¼ 0.82 in WM þ GM and R ¼ 0.56 in WM) among all methods, while moderate correlations were also found in WM between DIMWI and  . Three scenarios were considered for each parameter, including the simulated values are 30% less than (blue), exactly as (red) and 30% more than (yellow) the fixed value. The T 2 * EW simulations were slightly different from the others whereas the EW signal was simulated with a T 2 * decay instead of using the HCM formalism. Ground truth MWF is indicated by the horizontal dashed lines, and the true MWF changes with the simulated value in (d) resulting in individual ground truth for each scenario. Bias of the estimation is assessed as the deviation of the median to the ground truth and precision is compared from the width of IQR. MCR-DIMWI (R ¼ 0.46), and between standard method and DIMWI (R ¼ 0.38). Fig. 11a shows the barplots of mean MWF for the 5 different ROIs evaluated with the standard and the proposed methods. Good agreement was found on the CST, CG and ILF between the left and right hemispheres and, therefore, the statistics of these fibre bundles were computed using the masks from both hemispheres. The error bars demonstrate clearly the reduced noise of the MCR derived metrics. All methods evaluated here suggest that the highest MWF is found on FMj and ILF and the lowest for the CG. Not surprisingly, when evaluating the R 1 of the free compartment, a similar pattern is observable between both MCR methods (Fig. 11b). Here, the CG shows the lowest R 1 , whereas the ILF has the highest R 1 , and the two are distinguishable at the subject level. The ROI analysis on the R 2 * of MW (Fig. 11c) demonstrates that this is a measurement with high noise where no reliable differences can be seen between fibres. Finally, it should be noted that the R 2 * of IW benefits from the addition of the diffusion information, particularly for the MCR method (note the reduced standard deviation of MCR-DIMWI in Fig. 11d). Yet the differences remain quite small with the increased R 2 * of the FMn being attributed to the artefacts arising from B 0 inhomogeneities that can induce a faster decay of the signal.

Discussion
We have demonstrated that GRE-based MWI can be improved by either using the microstructural information derived from DWI and the Fig. 7. MWF maps obtained from (left to right) the standard model (with the average of 7 mGRE repetitions at flip angle 20 ), the DIMWI model (with 1 mGRE and 1 DWI measurements), the MCR model (with 7 mGRE of different flip angles) and the MCR model incorporated with diffusion information. The 4 rows represent 4 slices at different positions from the 3D acquisition of one subject. The right column is the volume fraction of IW estimated from DWI. Fibre bundles such as optic radiation and the splenium of the corpus callosum that are perpendicular to the B 0 are smoother and more distinguishable in the MWF map of DIMWI compared to those in the standard results (green arrows), which is also the case in MCR-DIMWI when compared to MCR MWF maps (yellow arrows). Unrealistic MWF in the globus pallidus showed in the standard method does not appear in the proposed models (blue arrows). Comparing to the other methods, DIMWI is more susceptible to the image artefacts especially close to the frontal region affected by the strong B 0 inhomogeneity resulting in close to zero MWF in a part of the forceps minor (red arrows). MCRbased results show less ringing artefact when compared to the single flip angle counterparts.
HCM formalism or modelling the T 1 -weighting effect with a VFA scheme. Firstly, by reducing the number of fitting parameters in the standard 3pool signal model using the HCM formalisms and the compartmental volume fraction derived from DWI, it is possible to obtain a more robust MWF estimation, especially for those fibre tracts that are not exactly parallel to the B 0 (See Fig. 4). Secondly, with the MCR model and VFA strategy, we introduced more variation in the data which allows not only to correct the T 1 weighting effect but also to have a more robust and reproducible MWF by effectively improving the fitting condition of the MWF estimation. Furthermore, we demonstrated that the two models can be combined to further improve the quality of the MWF result.
Overfitting is one of the major issues when using model-based processing in GRE MWI, which was illustrated in the Monte Carlo simulations with various tolerances shown in Fig. 4. Given the relatively large number of estimates in the 3-pool model of GRE MWI and, in particular, the similarity between the IW and EW compartments slow signal decay rates and frequency shifts, it is difficult to have a stable separation of these two compartments with a short echo time GRE acquisition, which can degrade the accuracy and precision of the MWF estimation. Theoretically, acquiring data at longer echo times can improve the fitting condition . However, propagation of phase errors due to breathing, resulting in a signal intensity-dependent noise amplitude as shown in Supplementary Fig. S1, limits the usefulness of the later echoes when performing in vivo imaging without dynamic shimming. Ways to overcome this limitation could be based on using field cameras (Duerst et al., 2014) or data-driven consistency estimation (Meineke and Nielsen, 2018).
Another way to improve the robustness of the obtained results is to reduce the degrees of freedom of the fitting process. This was our motivation to introduce WM microstructural information from other MRI techniques such as DWI. The multi-compartment models of DWI provide crucial information about the signal ratio of the free water compartments (IW/EW), which is a less robust parameter in MWI (Lee et al., 2016). Together with the HCM formalism, DIMWI reduces the number of fitting parameters from 10 to 6, and the advantages can be seen from the lower bias and narrower IQR in the simulations in Figs. 4 and 6, and less noisy in vivo MWF maps in Figs. 7 and 9. The VFA strategy with the MCR model paves an alternative direction to alleviate the ill-conditioning problem of MWI. By capitalizing on the observation that different water compartments have different T 1 relaxations, the VFA acquisition introduces wider signal diversity and, therefore, narrows down the number of possible solutions in the search space in data fitting. This is also supported by the fact that the MWF maps obtained from the MCR methods were very robust to the fitting tolerance when compared to the standard methodology (see supplementary Figs S2-S5). The advantages were demonstrated by: the MWF maps with excellent noise suppression on the MCR-based models (Figs. 7 and 9); the notable contrast between GM and WM; and realistic MWF in the globus pallidus (similar to other GM structures rather than myelin rich). These are further supported by the observation of MWF histograms in both WM and GM with the MCR-based methods being significantly narrower (Fig. 10), benefiting from a more stable data fitting. The VFA acquisition, from which the compartmental T 1 steady-state saturation can be modelled, provides MWF measurements that should be independent of sequence parameters such as flip angle and TR (although there might be some residual MT effects (Teixeira et al., 2019a(Teixeira et al., , 2019b). 2D multi-slice acquisition with long TR can also reduce the effect of T 1 weighting, as shown by , but will be necessarily less efficient from an SNR perspective .
To be able to reduce the number of fitting parameters without adding new unknowns in the DIMWI model, we fixed some of the myelin sheath properties in the HCM across the whole brain. Such choices, if incorrect, will adversely affect the results of other parameters in the model. For example, in GM, the MW R 2 * was close to the fitting boundaries when using DIMWI (Fig. 8). We attribute this observation to a more random arrangement of the GM microstructure which contradicts the assumptions of the HCM. We have tried to mitigate this inconsistency by considering multiple fibre orientations per pixel (3 orientations weighted by the corresponding anisotropic volume fraction) and including this in the signal model. The relatively low MW concentration in GM implied only a minimal impact on the IW R 2 * maps, which were nevertheless of higher quality compared to the standard model. Another noticeable difference between the models with and without diffusion information is the IW frequency map, where a more positive shift was predicted by the Fig. 8. Parameter maps estimated from all models on the 1.8 mm isotropic data in 1 subject. R 2 * maps of (a) MW and (b) IW, frequency maps of (c) MW and (d) IW, (e) the colour-coded fibre orientation map of the largest anisotropy estimated from DWI, (f) R 1 of free water and (g) exchange rate maps. The strong B 0 inhomogeneity in the frontal regions resulting in errors in some of the fitting parameters is known to be one of the limitations of GRE-based approach (Alonso-Ortiz et al., 2017Nam et al., 2015a) (red arrows). The IW R 2 * of the iron-rich globus pallidus is clearly faster in the proposed methods (blue arrows).
HCM. One possible explanation is that the frequency was overdetermined due to the high degree of freedom in the standard signal model, which can be seen as the IW frequency shift became less negative when the MCR model was deployed. Moreover, the IW frequency shift derived from the HCM is dependent on the anisotropic magnetic susceptibility of the myelin sheath which is fixed to À0.1 ppm (Wharton and Bowtell, 2012) in the DIMWI model. Therefore, any deviation of the myelin sheath susceptibility from the DIMWI assumption can result in systematic differences in the frequency maps as suggested in Fig. 6. The geometry of myelinated axons can further complicate the compartmental frequency shifts in contrast to the HCM that assumed the axons are circular (Xu et al., 2017). Interestingly, using the mGRE signal simulated with more realistic axonal geometry (Hedouin et al., 2019;Xu et al., 2017) can potentially serve as an alternative to substitute the analytical HCM formalism in DIMWI, which might also be able to explain the differences in the IW frequency maps between the current models.
Measuring and fixing some of the multi-compartment longitudinal relaxation and exchange rates was one of the main challenges when using the MCR models due to the lack of consensus in the literature regarding their true values. Detailed modelling of R 1 with exchange effect has been proposed using up to a four-pool model (Barta et al., 2015;Levesque and Pike, 2009) or a realistic myelin sheath model with the consideration of the exchange effect between the water molecules and the macromolecule protons in different compartments (Gelderen and Duyn, 2019). Such models require rich data from an inversion recovery sequence at multiple inversion times and with various inversion pulses which are difficult to obtain when doing full brain coverage and aiming at clinically acceptable scan times. Furthermore, these studies have been typically demonstrated Fig. 9. High-resolution MWF maps derived from 1.4 mm isotropic resolution GRE data on 2 subjects. The standard results are noisier than the corresponding maps in Fig. 7 due to lower SNR in the GRE data. With DIMWI, WM fibre bundles are more distinguishable and the MCR model results are superior to the models using a single flip angle across the whole brain. Further noise reduction can be seen in the MCR-DIMWI results. through signal averaging in large ROIs to achieve high SNR and, even with that increase in SNR, some of the model parameters were fixed to be able to robustly fit the remaining compartment parameters. This led us to choose the two-pool model to account for the multi-component nature of T 1 when describing the steady-state signal which, having fewer variables, is less SNR demanding. Recently, West et al. demonstrated that fitting a two-pool steady-state model is robust and reproducible when neglecting magnetisation exchange between the two proton pools (West et al., 2019). When the chemical exchange between bound and free protons is considered, the MWF estimation in mcDESPOT is inherently inaccurate (West et al., 2019). We observed that a two-pool model without exchange does not fit our data satisfactorily, resulting in a systematic signal overestimation when high flip angles are used (data not shown). To account for the magnetisation exchange observed in the GRE data without deteriorating the fitting condition, we considered an exchange system between 2 proton pools (free water and myelin sheath) that can be related to the biophysical environment, but even with such a simplified system, the search space displayed in Fig. 3 indicated that a wide range of parameters can explain the same signal behaviour. Fixing the myelin R 1 improved fitting stability at the expense of introducing a systematic bias when the actual value deviated from our assumption (see Fig. 6f). For example, if the fixed R 1 of myelin is higher than the actual value, a faster exchange will be needed to correctly fit the data and a higher MWF is assumed. Alternatively, if both relaxation and exchange rates were left as free parameters, the ill-posed nature of the model (see Fig. 4) would result in no significant improvements on the MWF estimation in respect to the standard model (see supplemental Fig. 7). Thus, more work is still needed to study accurately the multi-compartment R 1 relaxation when using a VFA approach. On the other hand, the magnetisation transfer (MT) effect is not modelled explicitly in the current MCR model although it is accepted to be one of the major sources of T 1 contrast in WM (Gelderen et al., 2016). These effects are partially explained by the magnetisation exchange effect between myelin and free water, and a relatively fast R 1 in the free water compartments compared to pure water R 1 . Nevertheless, balancing the RF power across acquisitions will still be important to control the MT effect (Teixeira et al., 2019a). Another alternative would be to fit the MCR model to multiple inversion time acquisitions with alternating inversion pulses (Gelderen et al., 2016;Gelderen and Duyn, 2019), rather than the VFA currently used. In such acquisitions, the change of the interference patterns between myelin and free water is greater than those observed in VFA because of the different polarity of those signals across inversion times.
The ROI analysis in 5 major WM fibre bundles allowed us to compare our results with other findings in the literature. This analysis showed that the largest MWF among the 5 studied fibres was found on the FMj and ILF. This finding is consistent with the observations of R 1 values measured by the MP2RAGE sequence studying the same subset of fibres, see Fig. 6 of Gil et al. (2016). Strong correlations between histological quantification of myelin and a single compartment R 1 were previously found in several studies (Mottershead et al., 2003;Stüber et al., 2014) focusing both on GM and WM, so this finding might not come as a surprise. Liu et al. have created an MWF atlas and reported mean values in similar regions (Liu et al., 2019). While both methods agree that the FMj (including the splenium) has higher MWF than the FMn (including the genu), the corticospinal tracts (referred in the atlas from (Liu et al., 2019) as the posterior internal capsule) showed very different results, appearing as the highest MWF on MSE methods and low on GRE methods. Interestingly, this is a pattern that was observed on our R 1 plots (Fig. 11b,   Fig. 10. (a) Individual and (b) group-averaged MWF histograms in GM (blue) and WM (red). All subjects have similar yet different histograms for each method. The group histograms of both the standard and DIMWI models have a board spectrum in both tissues. The histograms of MCR are well-characterised and there is a clear MWF difference between WM and GM. The WM histogram of MCR-DIMWI is even narrower and has a lower standard deviation (shaded region) across subjects, because of the further noise reduction achieved. (c) MWF correlation matrix in only WM (red) and both GM and WM (blue) between all methods. The number indicates the corresponding correlation coefficient. A relatively moderate-strong correlation was often found between methods utilising similar principles (e.g. MCR vs MCR-DIMWI: R ¼ 0.56 in WM and R ¼ 0.82 in WM þ GM; and DIMWI vs MCR-DIMWI: R ¼ 0.46 in WM).
although it lacks statistical significance). In summary, we expected some discrepancies between the evaluated methods shown in Fig. 11 , given the correlation results in Fig. 10c. While the weaker correlation in Fig. 10c could be related to the method robustness to noise amplification, the results in Fig. 11 point to more fundamental differences. When comparing our results to the literature findings using the MSE method (Liu et al., 2019), the inclusion criteria and extent of each WM mask can also play an important role, for example, our defined CST mask extends towards the brainstem. One should nevertheless note that the strong correlation (R 2 ¼ 0.84) reported in (Alonso-Ortiz et al., 2017) in MWF between the MSE and GRE approaches was obtained when pooling GM and WM ROIs together. When reproducing the analysis of (Alonso-Ortiz et al., 2017) using only ROIs of WM, the R 2 drops to 0.44. Such a reduction in correlation between MSE MWI and GRE MWI is in line with the differences we observe between our WM MWF and the work from Liu et al. (2019).
Lastly, this study has some limitations. The long processing time (~6 s per voxel) of the MCR model is a limitation if the technique is used with a larger sample size. The main hurdle is the computational cost of recomputing the extended phase graph in the steady state, despite our optimisation of the EPG-X implementation for data fitting. Another limitation has to do with the number of parameters to be fitted. Because of the current sequence implementation, each acquisition is acquired minutes apart, resulting in a need of an increased number of fitting parameters: separate frequency and phase offset for each T 1 -weighted data point (8 parameters in the case of our 7 flip angles used). If the acquisitions with the various T 1 weightings were interleaved, similar to what is done in a multi-echo MP2RAGE (Metere et al., 2017), 2 parameters would be sufficient to describe the phase and frequency offset for all acquisitions. Another methodological approach would be replacing the data fitting approach by a deep learning-based parameter estimation (Hedouin et al., 2019;Lee et al., 2019) which can potentially reduce the processing time significantly. To ensure we have sufficient signal variation for the MCR model, we acquired GRE data with 7 different flip angles, though proved to be sufficient, the scan time was inevitably long. Using highly-accelerated techniques such as echo planar time-resolved imaging (EPTI) with model-based subspace reconstruction (Dong et al., 2020;Wang et al., 2019) has potential to achieve significant scan time reduction without making a compromise on the number of flip angles being used for MCR.

Conclusion
Improved GRE MWI is achieved by incorporating the microstructural information derived from DWI together with analytical GRE signal behaviour characterised by the HCM, and multi-compartment relaxometry with VFA acquisition scheme. DIMWI is an interesting add-on to DWI to study WM microstructure. A relatively short acquisition (3.5 min added to a 15-min DWI scan) can provide information on myelin integrity that is not directly accessible in DWI. The MCR method significantly improves the robustness of the GRE-based MWI, supporting high quality, high resolution MWF maps to be acquired, which can be combined with DIMWI to achieve further improvements. Measuring compartmental R 1 and R 2 * could potentially reveal properties of the myelin environment that are more specific to structural changes than MWF and will be the focus of future studies.

Funding
This work is part of the research programme with project number FOM-N-31/16PR1056 supported by the Netherlands Organisation for Scientific Research (NWO).