Promise and pitfalls of g-ratio estimation with MRI

TL;DR

This paper reviews the development, methodology, and challenges of MRI-based g-ratio imaging, highlighting its potential for advancing understanding and diagnosis of neurological health and disease.

Contribution

It provides a comprehensive overview of the g-ratio MRI imaging technique, including its methodology, limitations, and future prospects.

Findings

MRI g-ratio imaging has potential for in vivo brain studies

Methodological challenges limit current accuracy

Future applications could enhance neurological diagnostics

Abstract

The fiber g-ratio is the ratio of the inner to the outer diameter of the myelin sheath of a myelinated axon. It has a limited dynamic range in healthy white matter, as it is optimized for speed of signal conduction, cellular energetics, and spatial constraints. In vivo imaging of the g-ratio in health and disease would greatly increase our knowledge of the nervous system and our ability to diagnose, monitor, and treat disease. MRI based g-ratio imaging was first conceived in 2011, and expanded to be feasible in full brain with preliminary results in 2013. This manuscript reviews the growing g-ratio imaging literature and speculates on future applications. It details the methodology for imaging the g-ratio with MRI, and describes the known pitfalls and challenges in doing so.

Equations

$$g=\sqrt{\frac{1}{1+MVF/AVF}}.$$

$$AVF=(1-MVF)(1-v_{iso})v_{ic},$$

$$g=\sqrt{1-MVF/FVF}.$$

$$MVF=cF+b,$$

Full text

††journal: NeuroImage

1 Introduction

The g-ratio is an explicit quantitative measure of the relative myelin thickness of a myelinated axon, given by the ratio of the inner to the outer diameter of the myelin sheath. Both axon diameter and myelin thickness contribute to neuronal conduction velocity, and given the spatial constraints of the nervous system and cellular energetics, an optimal g-ratio of roughly 0.6-0.8 arises (Rushton, 1951; Waxman, 1975; Chomiak & Hu, 2009). Spatial constraints are more stringent in the central nervous system (CNS), leading to higher g-ratios than in peripheral nerve (Chomiak & Hu, 2009). Study of the g-ratio in vivo is interesting in the context of healthy development, aging, learning, and disease progression and treatment. In demyelinating diseases such as multiple sclerosis (MS), g-ratio changes and axon loss occur, and the g-ratio changes can then partially recover during the remyelination phase (Albert et al., 2007). The possibility that the g-ratio is dependent on gender during development, driven by testosterone differences, has recently been proposed (Paus & Toro, 2009) and investigated (Pesaresi et al., 2015; Perrin et al., 2009). Possible clinical ramifications of a non-optimal g-ratio include “disconnection” syndromes such as schizophrenia (Paus & Toro, 2009), in which g-ratio differences have been reported (Uranova et al., 2001; Du et al., 2014).

The g-ratio is expected to vary slightly in healthy neuronal tissue. The relationship between axon size and myelin sheath thickness is close to, but not exactly, linear (Berthold et al., 1983), with the nonlinearity more pronounced for larger axon size (Hildebrand & Hahn, 1978), where the g-ratio is higher (Graf von Keyserlingk & Schramm, 1984). During development, axon growth outpaces myelination, resulting in a decreasing g-ratio as myelination catches up (Schröder et al., 1988). There is relatively little literature on the spatial variation of the g-ratio in healthy tissue. Values in the range 0.72-0.81 have been reported in the CNS of small animals (mouse, rat, guinea pig, rabbit) (Benninger et al., 2006; Duval et al., 2016a). Other primary pathology and disorders may lead to an abnormal g-ratio. These include leukodystrophies and axonal changes, such as axonal swelling in ischemia.

There are many outstanding questions in demyelinating disease that could be best answered by imaging the g-ratio in vivo. For example, in MS, disease progression is still the topic of active research. Most histopathological data are from patients at the latest stages of the disease. Therapies designed to promote remyelination are being developed to augment immunomodulatory and immunosuppressive treatments. Since remyelination of chronically demyelinated axons would be neuroprotective, this may help slow progression in MS. Detailed longitudinal study of the extent of remyelination can therefore aid in choosing avenues for therapy. While techniques exist for measurement of the g-ratio ex vivo (Graf von Keyserlingk & Schramm, 1984), measurement of the g-ratio in vivo is an area of active research.

Many MR imaging contrasts are sensitive to changes in the g-ratio, in that they are sensitive to changes in, e.g., the total myelin content or total fiber content. The purpose of the g-ratio imaging formulation is to decouple the fiber density from the g-ratio, such that a more complete and specific picture of the microstructural detail can be achieved. In recent work (Stikov et al., 2011, 2015b), it has been shown that the combination of an MRI marker that is sensitive to the myelin volume fraction (MVF) and an MRI marker that is sensitive to the intra-axonal volume fraction or axon volume fraction (AVF) is sufficient to compute a g-ratio for each voxel, i.e., an aggregate g-ratio, without explicit estimation of axon diameter and myelin sheath thickness. The aggregate g-ratio is a function of the ratio of the MVF to the AVF. g-Ratio imaging does not involve the acquisition of a novel contrast, but is a specific computation using the parameters (MVF and AVF) extracted from existing contrasts. The challenge then becomes how to estimate the MVF and the AVF precisely and accurately with MRI. The g-ratio imaging framework, coupled with independent microstructural measures such as axon diameter (Assaf et al., 2008; Zhang et al., 2011), comprises the field of in vivo histology of white matter. The ultimate goal is to describe microstructure in detail on a scale much finer than an imaging voxel, as a distribution or single aggregate value for each the voxel.

2 Methodology

2.1 The g-ratio formulation

As previously defined, the g-ratio is the ratio of the inner to the outer diameter of the myelin sheath of a myelinated axon (see Fig. 1). It has been shown in recent work (Stikov et al., 2011, 2015b) that the aggregate g-ratio can be expressed as a function of the myelin volume fraction and the axon volume fraction, and hence can be estimated without explicit measurement of these diameters:

[TABLE]

This formulation applies to any imaging modality (e.g., electron microscopy (EM) and scanning electron microscopy (SEM), where the MVF and AVF can be measured after segmentation of the image - see Fig.1), but it is of particular interest to be able to estimate the g-ratio in vivo. MRI provides us with several different contrast mechanisms for estimation of these volume fractions, and given MVFMRI and AVFMRI, we can estimate gMRI. Hereon, we sometimes refer to gMRI as “the g-ratio”, but note that it is derived from MRI images with certain contrasts sensitive but not equivalent to the MVF and AVF. We forego the subscript MRI on the MVF and AVF acronyms, but it should be clear from context when these are MR estimates. Estimation of these quantities is discussed in the next sections.

2.2 Axon volume fraction

Diffusion MRI is particularly well suited to aid in the estimation of the axon volume fraction. As we will see later, it cannot estimate the absolute AVF, because there is very little signal from the myelin in a typical diffusion MRI experiment. However, the relative AVF can be estimated, because diffusion MRI is sensitive to the displacement distribution of water molecules moving randomly with thermal energy, and this displacement distribution is affected by the cellular structure present in the tissue. As the molecules impinge on the cellular membranes, organelles, and cytoskeleton, the displacement distribution takes on a unique shape depending on the environment. Intra-axonal diffusion is said to be restricted, resembling free Gaussian diffusion at short diffusion times, but departing markedly from Gaussianity at longer times, where the displacement distribution is limited by the pore shape. There will be a sharp drop in the probability of displacement beyond the cell radius. Extra-axonal diffusion is said to be hindered, resembling free Gaussian diffusion, but with a smaller variance due to impingement of motion.

Many diffusion models exist for explicit estimation of the relative cellular compartment sizes. These include neurite orientation density and dispersion imaging (NODDI) (Zhang et al., 2012), the composite hindered and restricted model of diffusion (CHARMED) (Assaf & Basser, 2005), diffusion basis spectrum imaging (DBSI) (Wang et al., 2011), restriction spectrum imaging (RSI) (White et al., 2013), white matter tract integrity (WMTI) from diffusion kurtosis imaging (DKI) (Fieremans et al., 2011), temporal diffusion spectroscopy (Xu et al., 2014), microscopic anisotropy obtainable from multiple pulsed field gradient MRI (Shemesh et al., 2012; Zhou et al., 2013; Avram et al., 2013), the spherical mean technique (Kaden et al., 2016), the distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND) (Scherrer et al., 2016), and many others (Rokem et al., 2015; Stanisz et al., 1997). It is also possible to perform NODDI with relaxed constraints (NODDIDA) (Jelescu et al., 2016a), and to do this calculation analytically (LEMONADE) (Novikov et al., 2016).

Another approach, termed the apparent fiber density (AFD), uses high diffusion weighting to virtually eliminate the hindered diffusion signal, leaving only intra-axonal water (Raffelt et al., 2012). It has been used to estimate the relative axon volume fraction of different fiber populations in a voxel. A modification, termed the tensor fiber density (TFD), can be performed with lower diffusion weighting (Reisert et al., 2013).

The simplest diffusion MRI models do not differentiate between the tissue compartments. For instance, the diffusion tensor (Basser et al., 1994) models the entire displacement distribution as an anisotropic Gaussian function. The parameters defining this function will change if the intra-axonal volume fraction changes, but to what extent is it practical to extract meaningful quantitative compartment volume fractions from the tensor? Recently, a framework called NODDI-DTI has been developed, in which the proximity of DTI-based parameters to the computed NODDI parameters is assessed (Edwards et al., 2016), given certain assumptions. The FA and the mean diffusivity (MD) are highly correlated in straight, parallel fiber bundles, and will change with changing AVF, leading to estimates of the relative intra-axonal volume. However, this formulation is probably an oversimplification of the microstructural situation, and more detailed modeling is a better choice to ensure specificity to white matter fibers.

All of the diffusion models above are potential candidates for use in g-ratio imaging. They have strengths and weaknessess that will not be detailed here. The original full-brain gMRI demonstration (Stikov et al., 2015b) employed the NODDI model of diffusion. It was chosen because of its suitability in the presence of complex subvoxel fiber geometry, including fiber divergence, which may occur to a significant scale in almost all imaging voxels (Ghosh et al., 2016), and its suitability on clinical scanners with relatively low gradient strength. Having a fast implementation of the model fitting with numerical stability is important for large studies, hence, the convex-optimized AMICO implementation is beneficial (Daducci et al., 2015).

While diffusion MRI is a modality of choice for imaging microstructure, it can only measure the displacement distribution of water molecules that are visible in a diffusion MRI experiment. This limits us to water that is visible at an echo time (TE) on the order of 50-100 ms, and therefore excludes water that is trapped between the myelin bilayers, which has a T2 on the order of 10 ms. Hence, the estimates provided by these models are of the intra-axonal volume fraction of the diffusion visible volume. Myelin does not figure in the models. Given, e.g., the NODDI model outputs, a complementary myelin imaging technique must be used to estimate the absolute axon volume fraction. The AVF is given by

[TABLE]

where $v_{iso}$ and $v_{ic}$ are the isotropic and restricted volume fractions from the NODDI model, and the MVF is obtained from one of many possible myelin mapping techniques, examples of which are discussed below.

Diffusion contrast may not be our only window onto the axon volume fraction. Recent work has shown that it is possible to disambiguate the myelin, intra-axonal, and extra-axonal water compartments using complex gradient echo (GRE) images (Sati et al., 2013; Wu et al., 2016). The myelin water is separable from the combined intra- and extra-axonal water using multicomponent T${}_{2}^{*}$ reconstruction, providing a myelin marker (see below). However, incorporation of the phase of the GRE images potentially allows us to separate all three compartments based on frequency shifts. Challenges include the fact that the frequency shift is dependent on the orientation of the axon to the main magnetic field B0. When the axon is oriented perpendicular to B0, the myelin water will experience a positive frequency shift, the intra-axonal water a negative frequency shift, and the extra-axonal water will not experience a frequency shift.

Note that the AVF as defined by these diffusion MRI models is specific to white matter. While it makes sense to define an axon volume fraction in grey matter, the models in general cannot distinguish between axons and dendrites. The NODDI model’s $v_{ic}$ parameter, for example, is “neurite density”, i.e., the density of all cellular processes that can be assumed to have infinitely restricted diffusion in their transverse plane. Hence, gMRI from such MRI data is undefined in grey matter.

2.3 Fiber volume fraction

The fiber volume fraction (FVF), or fiber density, is the sum of the AVF and the MVF. Can diffusion MRI, or any other MRI contrast mechanism, measure the total fiber volume fraction itself? Clearly, GRE images have potential, as discussed above. Is diffusion imaging sensitive to the FVF, as opposed to the AVF? While myelin water is virtually invisible in diffusion MRI, diffusion MRI is not insensitive to myelin. First, the ratio of the intra- to extra-axonal diffusion MRI visible water in a voxel will change as the myelin volume fraction in that voxel changes. In early work on the fiber g-ratio, it was shown that assuming a simple white matter model of straight, parallel cylinders, the fractional anisotropy (FA) of the diffusion tensor, which weighs both intra-axonal and extra-axonal anisotropy, is proportional to the total fiber volume fraction, with a quadratic relationship (Stikov et al., 2011). The NODDI parameter $v_{ic}$ also changes with demyelination, even if all axons remain intact. Second, diffusion acquisitions are heavily T2 weighted, and T2 is myelin-sensitive. The total diffusion weighted signal thus decreases as myelin content increases. However, to robustly quantify myelin volume fraction, it is necessary to add a second contrast mechanism, even if it is additional T2 weighted images, to the scanning protocol. This is discussed below.

Despite the nomenclature, as noted above, even the Apparent Fiber Density and Tensor Fiber Density are in fact relative axon densities. They would provide a relative FVF only if the g-ratio is constant. In a recent study of the g-ratio (Mohammadi et al., 2015), the TFD was equated with the FVF, not the AVF, for input to the g-ratio formula. The g-ratio is a function of the ratio of the MVF to the AVF (Eq. 1), or alternately, the ratio of the MVF to the FVF:

[TABLE]

This means that conclusions reached about the variation of the g-ratio found by equating the TFD with the FVF will be robust in the above case. Absolute g-ratios in the above case were calibrated to have a mean of 0.7 in healthy white matter.

We note that if diffusion MRI were capable of estimating the absolute FVF (or AVF) as well as the ratio of the intra-axonal to extra-fiber water, the g-ratio could immediately be estimated from these two quantities, without further myelin imaging. This has yet to be done robustly, and it is therefore preferable to use a more robust independent myelin marker.

2.4 Myelin volume fraction

There are many different contrasts and computed parameters that are sensitive to myelin (Laule et al., 2007). The possible sources of signal from the myelin compartment are the ultra-short T2 protons in the macromolecules of the myelin sheath itself (T${}_{2}\sim 10\mu s$) and the short T2 water protons present between the phospholipid bilayers (T${}_{2}\sim 10ms$, see Fig. 1). Most MRI contrast mechanisms are sensitive to myelin content, but few are specific, for reasons that are detailed in section 3.3. The myelin phospholipid bilayers create local Larmor frequency variations for water protons in their vicinity due to diamagnetic susceptibility effects. This results in myelin content modulated transverse relaxation times T${}_{2}^{*}$ (Hwang et al., 2010) and T2, and longitudinal relaxation time T1. It has been shown that macromolecular content is the dominant source of variance in T1 in the brain (Rooney et al., 2007). The local Larmor field shift (fL) and the susceptibility itself ($\textstyle\chi$) can be computed as well (Liu et al., 2015). Ultra-short TE (UTE) imaging can be used to image the protons tightly bound to macromolecules (Wilhelm et al., 2012; Du et al., 2014). An alternate approach to isolating the myelin compartment is magnetization transfer (MT) imaging, where the ultra-short T2 macromolecular proton pool size can be estimated by transfer of magnetization to the observable water pool.

MT based parameters sensitive to macromolecular protons include the magnetization transfer ratio (MTR) (Wolff et al., 1991), the MT saturation index (MTsat) (Helms et al., 2008, 2010), the macromolecular pool size (F) from quantitative magnetization transfer (Sled & Pike, 2001; Yarnykh, 2002; Ramani et al., 2002), single-point two-pool modeling (Yarnykh, 2012), and inhomogeneous MT (Varma et al., 2015).

Alternately, the myelin water can be imaged with quantitative multicomponent T2 (MacKay et al., 1994) or T${}_{2}^{*}$ (Du et al., 2007; Alonso-Ortiz et al., 2016) relaxation, which yields the myelin water fraction (MWF) surrogate for myelin density. Variants include gradient and spin echo (GRASE) MWF imaging (Does & Gore, 2000; Prasloski et al., 2012), linear combination myelin imaging (Jones et al., 2004; Vidarsson et al., 2005), T2 prepared MWF imaging (Oh et al., 2006), multi-component driven equilibrium single point estimation of T2 (mcDESPOT) (Deoni et al., 2008), direct visualization of the short transverse relaxation time component via an inversion recovery preparation to reduce long T1 signal (ViSTa) (Oh et al., 2013), and a fast adiabatic T2-prep and spiral readout approach (FAST-T2) (Nguyen et al., 2016). Other alternate approaches exploiting myelin-modulated relaxation times include combined contrast imaging (T1W/T2W) (Glasser & Van Essen, 2011) or independent component analysis (Mangeat et al., 2015). Proton density is also sensitive to macromolecular content, and the proton-density based macromolecular tissue volume (MTV) (Mezer et al., 2013) has been used as a quantitative myelin marker.

While these MRI measures have been shown to correlate highly with myelin content (Thiessen et al., 2013; Schmierer et al., 2007; Gareau et al., 2000; Laule et al., 2006; Mottershead et al., 2003), they have not been incorporated in a specific tissue model in a manner similar to the diffusion signal, and hence some calibration is needed. This is still a topic of research. Caveats of improper calibration of the MVF are discussed in section 3.3.

Above, we have discussed imaging techniques for both diffusion-visible microstructure and myelin. Any multi-modal modal imaging protocol with contrasts such as these, sensitive to the axon and myelin volume fractions, is sensitive to the g-ratio (e.g., (Molina-Romero et al., 2016; Nossin-Manor et al., 2015; De Santis et al., 2016a; Bells et al., 2011)). The purpose of the explicit g-ratio formulation is to create a measure that is specific to the g-ratio. It provides us with a metric that is independent of the fiber density, which none of the MRI contrasts sensitive to the MVF and AVF separately accomplish. It is interesting to ask whether we could use a technique such as deep learning (Bengio, 2009) to estimate the g-ratio, skipping explicit modeling completely.

2.5 MRI acquisition protocol: an example

Given the large number of MRI contrasts that could potentially be used in the g-ratio imaging framework, we cannot prescribe an exact protocol for g-ratio imaging here. However, in the following sections, we illustrate several important points about g-ratio imaging using experimental data acquired at our site. The following describes the acquisition protocol we have used.

We acquired data from healthy volunteers and from multiple sclerosis patients. These data were acquired on a Siemens 3T Trio MRI scanner with a 32 channel head coil. A T1W structural MPRAGE volume with 1 mm isotropic voxel size was acquired for all subjects. For diffusion imaging, the voxel size was 2 mm isotropic. For most experiments, the NODDI diffusion protocol consisted of 7 b=0 s/mm2, 30 b=700 s/mm2, and 64 b=2000 s/mm2 images, 3x slice acceleration, 2x GRAPPA acceleration, all acquired twice with AP-PA phase encode reversal. For the other experiments, as detailed below when they are introduced, the slice acceleration and phase encode reversal were not employed. For a dataset optimized for diffusion tensor reconstruction, a dataset with 99 diffusion encoding directions at b=1000 s/mm2 and 9 b=0 s/mm2 images was acquired.

For magnetization transfer images, we also used 2 mm iso-tropic voxels to match the diffusion imaging voxel size. For MTR, one 3D non-selective PD-weighted RF-spoiled gradient echo (SPGR) scan was acquired with TR=30 ms and excitation flip angle $\alpha=5^{\circ}$, and one MT-weighted scan was acquired with the same parameters and an MT pulse with 2.2 kHz frequency offset and 540*∘* MT pulse flip angle. For MTsat computation, these same MT-on and MT-off scans were used, with one additional T1-weighted scan with TR=11 ms and excitation flip angle $\alpha=15^{\circ}$. For qMT computation, 10-point logarithmic sampling of the z-spectrum from 0.433-17.235 kHz frequency offset was acquired, with two MT pulse flip angles for each point, 426*∘* and 142*∘, and excitation flip angle $\alpha=4.5^{\circ}$. The qMT acquisition was accelerated with 2x GRAPPA acceleration. Additional scans for correction of the maps included B1 field mapping using the double angle technique (Boudreau et al., 2017; Stollberger & Wach, 1996), with 60∘* and 120*∘* flip angles, B0 field mapping using the two-point phase difference technique, with TE1/TE2 = 4.0/8.48 ms, and T1 mapping using the variable flip angle technique (Fram et al., 1987), with flip angles 3*∘* and 20*∘. Additional T2-FLAIR and PDW* images were acquired for the MS subjects to aid in lesion segmentation.

3 The pitfalls: outstanding technical challenges

In this section, we discuss pitfalls and outstanding issues in g-ratio imaging. These include the challenge of combining multiple different contrast mechanisms, limitations in diffusion modeling, specificity and calibration of myelin markers, and the limitation of estimating a single g-ratio metric in the presence of subvoxel heterogeneity. Experimental results are included in these sections to illustrate these problems. In the experiments, we focus on the NODDI and DTI models and MT-based myelin imaging because they were used in the early g-ratio imaging publications (Stikov et al., 2011, 2015b), but many of these concerns apply to any chosen AVF and MVF markers.

3.1 Confounds in multi-modal image acquisition

The computation of the g-ratio metric includes several pre-processing steps, including distortion and field inhomogeneity correction, that deserve further discussion. The MT-based contrasts are acquired with spin-warp acquisition trains, and the diffusion-based contrasts are acquired with single-shot EPI. When any acquisition details are changed, the distortions in the images change, and co-registration of voxels for voxelwise quantitative computations becomes more difficult.

The blip-up blip-down phase encode strategy (section 2.5) allows for precise correction of susceptibility-induced distortion in the diffusion images (Andersson et al., 2003). Lack of correction for this distortion leads to visible bands of artifactually high gMRI near tissue-CSF interfaces (see, e.g, (Campbell et al., 2014; Cercignani et al., 2016a)). This was illustrated by Mohammadi et al. (Mohammadi et al., 2015) (see Fig. 2). Uncorrected diffusion MRI data leads to g-ratios in the vicinity of unity at the edge of the genu of the corpus callosum, caused by voxels where the AVF is artifactually high (containing little or no CSF), and the MVF low (because the correctly localized voxels actually contain CSF). The white matter - CSF boundary is a region of obvious misregistration, but much of the frontal lobe suffers from susceptibility induced distortion, and would therefore have incorrect g-ratios.

Multi-modal imaging protocols are a powerful tool for investigation of microstructure. We have thus far discussed combining multiple images with partially orthogonal contrasts in order to estimate the g-ratio. However, problems such as the above misregistration issue arise. Can a single acquisition train provide multiple contrasts? One such approach was described recently for simultaneous mapping of myelin content and diffusion parameters (De Santis et al., 2016a). It consists of an inversion-recovery preparation before a diffusion weighted sequence, allowing for fitting of a model that incorporates both T1 (a myelin marker (Stüber et al., 2014)) and axonal attributes. This approach is conceptually extensible to other myelin-sensitive preparations or modifications of a diffusion weighted sequence, such as quantitative T2 estimation (Kim et al., 2016) or MT preparation (Gupta et al., 2003).

Can the g-ratio be estimated using a single contrast mechanism? This could also offer inherent co-registration, as well as possibly increasing the acquisition speed. As discussed in section 2.2, analysis of complex GRE images may lead us to a technique for estimating the MVF and AVF volume fractions from one set of images.

Making a multi-modal imaging protocol short enough for the study of patient populations and use in the clinic is a considerable challenge. In section 3.3, we investigate the MT saturation (Helms et al., 2008) as a more time efficient replacement for qMT. There exist other short MT-based approaches, such as single-point two-pool modeling (Yarnykh, 2012) and inhomogeneous MT (Varma et al., 2015). Another approach could be to use compressed sensing (Lustig et al., 2007) for MT-based acquisitions (McLean et al., 2017). Recent advances in myelin water fraction imaging (Nguyen et al., 2016) may make MWF estimation more efficient, and GRE-based myelin water fraction approaches (Alonso-Ortiz et al., 2016) may also offer a faster approach for estimating the MWF, with the possibility, as mentioned above, to eliminate the diffusion imaging part of the protocol.

Diffusion imaging has benefited from many acceleration approaches in recent years, including parallel imaging, which can also be used in the myelin mapping protocols, slice multiplexing (Setsompop et al., 2012), and hardware advances such as the CONNECTOM gradient system.

3.2 Diffusion modeling

3.2.1 Model parameters

The diffusion MRI post-processing techniques described in section 2.2 give a range of outputs. Some are physical quantities (such as the diffusion displacement distribution; kurtosis), while some are parameters of detailed biological models (such as the intra-axonal volume fraction). Models are valuable, but the user has to be aware of the assumptions made.

The parameter space in existing models ranges from three free parameters in NODDI to six (Assaf & Basser, 2005), twenty three (Jespersen et al., 2010), and thirty one (Novikov et al., 2016; Jelescu et al., 2015) in other models, with the maximum dependent on acquisition details. Recent analysis hypothesizes that the lower number of free parameters in, e.g., NODDI and CHARMED, may be matched to the level of complexity possible on current clinical systems (Ferizi et al., 2015), while high gradient strength, high b-values, and more b-shells may be necessary for more complex models (Jelescu et al., 2015), and would make them more optimal. This is a general problem with multi-exponential models when diffusion weighting is weak (Kiselev & Il’yasov, 2007). On standard MR systems, relaxing the constraints on fixed parameters has been shown to lead to degeneracy of solutions (Jelescu et al., 2016a). Regularization approaches such as the spherical mean technique (SMT) (Kaden et al., 2016) have been employed in an attempt to make the problem less ill-posed.

One of the fixed parameters in the NODDI model is the parallel diffusivity in the intra- and extra-axonal space, both set to the same fixed value. Other models explicitly model these as unequal; for instance, WMTI assumes that the intra-axonal diffusivity is less than or equal to the extra-axonal diffusivity. The actual values are unknown, however simulations have shown that the assumption of equal parallel diffusivities leads to a 34-53% overestimation of the intra-axonal compartment size if the diffusivities are in fact unequal (Jelescu et al., 2015), with the intra-axonal diffusivity either greater than or less than the extra-axonal diffusivity. If the fixed diffusivities are incorrect, the intra-axonal compartment size estimated by the NODDI model will be non-zero even if there is no anisotropy (Lampinen et al., 2017).

Independent of whether the intra- and extra-axonal parallel diffusivities are equal, another source of this bias is the tortuosity model (Szafer et al., 1995) employed by many models, including NODDI, DIAMOND, and the SMT. This model computes the perpendicular extra-axonal diffusivity as a function of the diffusion-visible intra-axonal volume fraction of the non-CSF tissue (vic in the NODDI model). This tortuosity estimate is bound to be inaccurate because the tortuosity is expected to vary as the absolute fiber volume fraction of the non-CSF tissue, not the diffusion-visible fiber volume fraction (Jelescu et al., 2015). These two quantities are very different, as the myelin and axon volume fractions are almost equal in healthy tissue (Stikov et al., 2015b). The MVF could be explicitly included in the equation, and would be expected to result in a myelin volume dependent reduction in $v_{ic}$. However, in healthy tissue, where the FVF should scale roughly as the vic parameter, the tortuosity model of Szafer et al. does not appear to hold when applied to experimental data with independent estimates of the parallel and perpendicular extra-axonal diffusivities (Novikov et al., 2016). The model is probably not correct for varying packing density on the sub-voxel scale (Novikov et al., 2011); it has been shown to depend on the packing arrangement and breaks down for tight axon packing (Fieremans et al., 2008; Novikov & Fieremans, 2012). This might explain the discrepancy between model and experiment mentioned above, because the geometry of axonal packing can vary considerably for a given average volume fraction.

Another fixed parameter in the NODDI model is the T2 relaxation time of all tissue, assumed to be the same, even in CSF. This leads to an overestimation of viso, which can be corrected (Bouyagoub et al., 2016) given T2 estimates from, e.g., a T2 mapping technique such as mcDESPOT (Deoni et al., 2008).

3.2.2 Geometry in diffusion imaging

Diffusion MRI is exquisitely sensitive to fiber geometry. The fractional anisotropy may be more sensitive to geometry than to any microstructural feature (Hutchinson et al., 2016). Hence, microstructural models must be careful to take geometry (crossing, splaying, curving, microscopic packing configuration) into account. A typical diffusion imaging voxel is roughly 8 mm3, while the axons probed by microstructural models are on the order of one micron.

Early work on the fiber g-ratio investigated the human corpus callosum, where it was assumed that the white matter fibers were effectively straight and parallel (Stikov et al., 2011). The use of this model that assumes straight, parallel fibers suffers from several problems. First, the regions of the brain where this model can be expected to hold at all are very limited, as there are crossing or splaying fibers in up to 95% of diffusion MRI voxels in parenchyma (Jeurissen et al., 2010; Behrens et al., 2007; Ghosh et al., 2016), and curvature is almost ubiquitous at standard imaging resolution. Even the axons of callosal fibers are not straight and parallel, with splay up to 18*∘* (Ronen et al., 2014; Mollink et al., 2016). Second, the model assumes a relatively uniform, if random, packing of axons on the scale of the MRI voxel. Due to the nonlinear nature of the FA, it will depend strongly on the packing geometry. If two voxels, one with densely packed axons and one with sparsely packed axons, are combined into one, the FA for that voxel will be less than the average of the two original voxels, whereas the fiber density will be the average of the fiber densities. Third, FA is in practice acquisition and b-value dependent.

3.2.2.1 Experiments: AVF from NODDI in the presence of crossing fibers

The NODDI model that has been used in several g-ratio imaging studies to date assumes there is a single fiber population with potential splay or curvature, but does not explicitly model crossings. To what extent does the fiber dispersion model of NODDI handle crossing fiber bundles? We have employed the diffusion MRI simulator dSim (Sveinsson & Dougherty, 2011) to investigate this question (Campbell et al., 2014). We simulated realistic axonal packing (Aboitiz et al., 1992) in voxels with straight, parallel fibers and with two equal size bundles of straight fibers crossing at 90*∘* (See Fig.3). Fiber volume fractions were set equal for both configurations and were varied from 0.3 to 0.7. g-Ratios were varied from 0.7 to 0.9. The diffusion weighted signal was generated, and the NODDI model parameters computed using the NODDI Matlab toolbox (NODDI Matlab Toolbox, 2013). The FVF was computed from the NODDI parameters using the known MVF. The computed FVF was $3.8\pm 0.3\%$ lower in the crossing fiber case for the NODDI-based FVF. This demonstrates that the NODDI model, while not explicitly designed for crossing fibers, gives acceptable results in this case, and can be used for full-brain g-ratio estimation at standard voxel size with significant subvoxel fiber crossing, with only a small decrease in the estimated FVF due to partial volume averaging of fiber orientations.

3.2.2.2 Experiments: Comparison of DTI and NODDI for FVF estimation

NODDI works optimally with diffusion MRI measurements made on at least two shells in q-space, i.e., two different nonzero b-values, although recent work has proposed solutions for single shell data, at least where certain assumptions can be made about the tissue, or where high b-values are used (Grussu et al., 2014; Magnollay et al., 2014). In contrast, the diffusion tensor can be robustly fitted and the fiber volume fraction inferred (Stikov et al., 2011) (see section 2.3) using a much more sparsely sampled, single shell dataset. Many research programs have large databases of single-shell diffusion data, often with limited angular sampling of q-space as well. It is therefore of interest to explore to what extent such data, using the diffusion tensor model, can be used in investigation of the g-ratio.

In the simulations described above, we also computed the diffusion tensor using in-house software (mincdiffusion, 2013). The FVF (FVFDTI) was computed from the FA using the quadratic relationship determined from previous simulations (Stikov et al., 2011). As expected, the FA is not a predictor of FVF in the presence of crossing fibers: the computed DTI-based FVF was $58.7\pm 3.2\%$ lower in the crossing fiber case compared to the parallel fiber case.

To compare NODDI and DTI in vivo, diffusion and qMT data were acquired as described in section 2.5 for one healthy volunteer, without slice acceleration or phase encode reversal. The qMT data were processed with in-house software and the NODDI parameters as described above. Additionally, the diffusion tensor was calculated using the b=1000 s/mm2 diffusion shell. The AVF, MVF, and gMRI were computed voxelwise from the diffusion and qMT data as described in section 2.1: the NODDI-based FVF (FVFNODDI) is the sum of the MVF and the AVF computed from Eq. 2, and FVFDTI was calculated from the fractional anisotropy of the diffusion tensor using the quadratic relationship (Stikov et al., 2011). The corpus callosum was skeletonized on the FA image and a voxel-wise correlation between the FVF computed from DTI and from NODDI and qMT was performed for these voxels. The coefficient of proportionality between F and MVF was determined from previous EM histological analysis (Stikov et al., 2015b).

Fig. 4 shows the FVF computed using both NODDI and the FA from DTI in the skeleton of the healthy human corpus callosum. The Pearson correlation coefficient between the FVF measured using the two techniques was r=0.79, with $\mbox{FVF}_{\mbox{NODDI}}=1.15*\mbox{FVF}_{\mbox{DTI}}+0.00$. This indicates a slight discrepancy between the FVF using NODDI compared to DTI, and a reasonably high correlation between techniques on the skeleton.

Possible explanations for the higher estimates using NODDI appear in section 3.2.1, although without ground truth is is difficult to say which approach is more accurate. Additionally, because the FA does not explicitly model compartments, it is subject to partial volume effects. While partial volume averaging with CSF will decrease the FA, the FA-based quadratic FVF model appears to break down in this case. This effect could possibly be reduced by applying the free water elimination technique (Pasternak et al., 2009) to obtain the correct FVF for the non-CSF compartment and then scaling to reflect the partial volume averaging with CSF afterward. To conclude, while the FA is generally a poor indicator of FVF, it may be a reasonable surrogate in certain special cases when data are limited.

It is interesting to consider how useful imaging a cross-section of a white matter fascicle may be, regardless of the model used. If the g-ratio can be assumed to be constant along an axon, measurement of a cross-section is useful. However, in many pathological situations, such as Wallerian degeneration, it is of interest to study the entire length of the axon.

3.2.3 Restricted extracellular diffusion

Most existing diffusion models assume that extra-axonal diffusion is Gaussian, hindered by the structures present, but not restricted. However, observation of tightly packed axons in microscopy (e.g., Fig.1) indicates that the intra- and extra-axonal spaces may not be as distinguishable as the models assume. It is unclear to what extent the extra-axonal diffusion is non-Gaussian. If axons are packed tightly together, is extra-axonal diffusion non-Gaussian? It is not clear whether the water mobility through the tight passageways between fibers is distinguishable from the restricted diffusion within spaces surrounded by contiguous myelin. If signal from the extra-axonal space is erroneously attributed to the intra-axonal space, the model output will be incorrect.

Some models make no attempt to distinguish between intra- and extra-cellular restricted diffusion, meaning the pore size estimates may reflect a mixture of the two (Ning et al., 2016). Time-dependent (i.e., non-Gaussian) diffusion has recently been observed in the extra-axonal space using long (Fieremans et al., 2016; De Santis et al., 2016b) and short (Xu et al., 2014) diffusion times. This may be due to axon varicosity, axonal beading, or variation in axonal packing (Fieremans et al., 2016). In order for the intra- and extra-axonal compartment fractions to be estimated, the diffusion MRI signal from these compartments must be distinguishable with the experimental paradigm and modeling used.

Diffusion modeling is an active field, and advances in the near future will hopefully improve precision and accuracy of AVF estimates using diffusion MRI. Histological validation may aid in understanding the strengths and limitations of these estimates. At present, the limitations of these models propagate to the g-ratio, as do the limitations of MVF estimates, which are discussed below.

3.3 MVF calibration

In this section, we investigate the effect of imperfect calibration of MRI markers of the MVF. As an example, we present detailed experimental analysis of the magnetization transfer ratio (MTR), which is a commonly used myelin marker that is implemented on clinical scanners and would be a logical first choice for g-ratio imaging in the clinic. Computing a precise and accurate MVF from the MTR is, however, challenging, as is detailed below. We also investigate the performance of MTsat (Helms et al., 2008) as a potential time-efficient replacement for qMT in the estimation of the MVF. Care must be taken in the use of any myelin marker in quantitative g-ratio computation.

How do we make a quantitative estimate of the MVF from myelin sensitive MRI markers? Linear correlations have been shown between the individual myelin sensitive metrics (such as F (Thiessen et al., 2013; Schmierer et al., 2007), MTR (Gareau et al., 2000; Schmierer et al., 2004), R1 (Mottershead et al., 2003), MWF (Laule et al., 2006), and MTV (Duval et al., 2016b)) with the MVF from histology. Given the linear correlations that have been established, a logical first approximation is to assume a linear relationship between the chosen myelin-sensitive metric and the MVF. Then, using the macromolecular pool size F as an example, the relationship is

[TABLE]

with $c$ and $b$ constants. While a non-zero value for $b$ has been indicated by some studies (Thiessen et al., 2013; West et al., 2017), this may be an artifact due to the inherent bias in linear regression. The assumption of a linear relationship hinges on the assumption that non-myelin macromolecular content scales linearly with myelin content, but this relationship can break down in disease, or even in healthy tissue. In disease, astrocyte scarring, glial cell processes, and inflammatory cell swelling could all modulate the relationship between a marker of macromolecular content and the MVF. If the myelin and non-myelin macromolecular content do scale linearly, as is assumed here, a theoretical prior that $b=0$ is reasonable. However, $b$ could reflect a fixed population of macromolecules that are uncorrelated with myelin content. It would most likely be negative, i.e., F could be positive with MVF=0.

There is evidence that even if a simple scaling relationship exists between F and MVF, it is dependent on acquisition and post-processing details. For instance, a recent study calibrated F at two different sites, and found a different scaling factor for each (Cercignani et al., 2016a). These scaling factors in turn differ from those obtained from other investigations (Stikov et al., 2015b; Dula et al., 2010; Thiessen et al., 2013). Hence, careful calibration for each study must be performed. Several studies have calibrated scaling factors based on a given expected g-ratio in healthy white matter (Cercignani et al., 2016a; Mohammadi et al., 2015). However, the g-ratio in healthy white matter is not precisely known.

None of the myelin-sensitive MRI markers is 100% specific to myelin, and most are sensitive to myelin in a slightly different way. Magnetization transfer contrast is specific to macromolecules, and more specific to lipids than to proteins (Kucharczyk et al., 1994). Macromolecules in the axon membrane itself, in neurofilaments within the axons, and in glial cell bodies, will contribute to the MT signal, with myelin constituting only 50% of the macromolecular content in healthy white matter (Bjarnason et al., 2005). Additionally, MT-based metrics such as the magnetization transfer ratio will have residual contrast from other mechanisms. We expect the MTR contrast to vary linearly with macromolecular content, but also with T1 (Vavasour et al., 2011). T1 has the opposite sensitivity to myelin than does the MT effect (Mottershead et al., 2003), meaning that these effects work against each other, reducing the dynamic range and power of MTR as a marker of myelin. Furthermore, T1 is sensitive to iron and calcium content, intercompartmental exchange, and diffusion, and hence sensitive to axon size (Harkins et al., 2016) and axon count (Schmierer et al., 2008). This means the relationship between MTR and MVF may not be monotonic, and is certainly nonlinear. This nonlinearity is evident in published plots of MTR vs. F, e.g., that shown by Levesque et al. (Levesque et al., 2005), and the lack of dynamic range of MTR is also evident (Levesque et al., 2005; Garcia et al., 2012). The MTsat technique aims to remove the T1 dependence in MTR. Both MTR and MTsat depend on the offset frequency used in the acquisition. ihMT shows promise as a more myelin-specific MT marker due to its sensitivity to specific molecules in myelin that broaden the z-spectrum asymmetrically, although it has recently been shown that asymmetric broadening is not essential to generate a non-zero ihMT signal (Manning et al., 2016), and the technique suffers from low signal. qMT is the most comprehensive of the MT-based myelin markers, although its use is impeded by long acquisition times, and its parameters appear to be sensitive to the specific model and fitting algorithm.

Proton-density based techniques (Mezer et al., 2013) will, like MT, be sensitive to all macromolecules, with a different weighting on these macromolecules compared to the lipid dominated MT signal. Relaxation-based myelin markers are also not 100% specific to myelin. The confounds with using T1 directly were mentioned above, and the dependence on iron and calcium concentration, intercompartmental exchange and diffusion will also affect T2. T${}_{2}^{*}$ is also sensitive to iron concentration, as well as fiber orientation (Cohen-Adad, 2014). Isolating the short T2 or short T${}_{2}^{*}$ compartment enhances specificity to myelin, but MWF estimates vary nonlinearly with myelin content as the sheath thins and exchange and diffusion properties are modulated (Levesque & Pike, 2009; West et al., 2014; Harkins et al., 2012). Variants may suffer from reduced accuracy or precision, for example, the mcDESPOT technique has been shown to overestimate the MWF (Bouhrara et al., 2016) and lack precision (Lankford & Does, 2013). Combining T1 and T2 in various ways (Glasser & Van Essen, 2011; Mangeat et al., 2015) may increase specificity, although this approach relies on myelin being the dominant source of contrast. In the UTE technique, it is as yet unclear how to map the signal directly to myelin content. In addition to these confounds, most of these myelin markers have recently been shown to have orientation dependence. These include T${}_{2}^{*}$, $\textstyle\chi$ (Rudko et al., 2014; Liu et al., 2015), and T2 of the macromolecular pool (Pampel et al., 2015).

While these myelin imaging techniques are certainly powerful tools in the study of healthy and diseased brain, can they be used reliably in the g-ratio imaging framework? As an illustration of the effects of miscalibration of myelin markers, consider the following scenario (Campbell et al., 2016). We investigate three MT-based myelin markers: MTR, MTsat, and macromolecular pool size F. We assume a simple linear scaling between our MRI marker and the MVF. As described in previous work (Stikov et al., 2015b), we calibrate F using combined in vivo MRI acquisition ex vivo electron microscopy in the macaque. We then calibrate MTR and MTsat to match the mean F-based MVF in white matter. We subsequently compute gMRI, using the NODDI model of diffusion and the MVF derived from the myelin markers (Eq.s 2,1).

3.3.1 MVF calibration: Experimental Methods

MTR, MTsat, qMT, and NODDI data were acquired for five healthy volunteers and one MS patient, as described in section 2.5. Additionally, for two of the healthy subjects, we acquired one MT-on image with an offset frequency of 1.2 kHz, which is the standard MTR sequence used at our site. For the MS patient, the MTsat images were computed from the qMT MT-off and MT-on (MT pulse offset 2.732 kHz, flip angle 142*∘) images and one additional T1W* image with TE=3.3 ms, TR=15 ms, and excitation flip angle $\alpha=20^{\circ}$. The diffusion images were preprocessed using FSL (Smith et al., 2004), and the NODDI parameters were computed using the NODDI matlab toolbox (NODDI Matlab Toolbox, 2013). The qMT computation of the macromolecular pool size F was performed using in-house software (Sled & Pike, 2001; Cabana et al., 2015), including B0 and B1 correction. MTsat was computed according to Helms et al. (Helms et al., 2008, 2010), using the 2.2 kHz offset frequency data, and also the 1.2 kHz offset frequency data where available. MTR was computed for both offset frequencies where available. A semi-empirical B1 correction was made (Weiskopf et al., 2013) to correct for higher order B1 effects. Binary segmentation of white and grey matter was performed using an in-house pipeline, using the MPRAGE image only. Lesion segmentation for the MS subject was performed with in-house software.

The combined MRI/histology dataset (Stikov et al., 2015b) was used to scale each myelin marker (MTR, MTsat, and F) to give the MVF, with the assumption of a linear relationship (Eq. 4) with intercept b=0. Correlations between the three myelin markers were computed in brain parenchyma. Percent differences were computed between healthy white matter and healthy grey matter for each of the three myelin markers. The AVF was computed using Eq. 2, and g-ratios were computed in the MS and healthy brains using Eq. 1. Average g-ratios were computed in healthy white matter, normal appearing white matter (NAWM), and MS lesions.

A theoretical computation was also performed, varying the mapping of an arbitrary myelin metric to MVF using Eq. 4. We separately varied the slope (c) and the intercept (b) for a range of fiber volume fraction values and mapped the computed g-ratio as a function of FVF. When varying the slope, the intercept was fixed at the origin.

3.3.2 MVF calibration: Results

For the 2.2 kHz offset frequency (Helms et al., 2008), the average correlation of MTR with F was r=0.62 (p$<$0.001), and of MTsat with F was r=0.80 (p$<$0.001) in parenchyma. The relationship between MTsat and F is more linear than the relationship between MTR and F (p$<$0.0001). Fig. 5 shows plots of MTR versus F (left) and MTsat versus F (right) in parenchyma for one subject at 2.2 kHz offset frequency. Of note, the plot of MTR versus F has a distinctive nonlinear shape, similar to that seen in the literature (Levesque et al., 2005). When T1 effects are reduced using MTsat, the linearity and dynamic range increase. For the two subjects in which the lower, 1.2 kHz offset frequency was also used to compute both MTR and MTsat, the average correlation of MTR with F was 0.55 and of MTsat with F was 0.73.

In healthy brain, the percent difference between white and grey matter was 15.02% for MTR, 40.08% for MTsat, and 45.86% for F. The narrower dynamic range of the MVF derived from MTR can also be seen in Fig. 6, where grey matter has markedly higher values. If this simple scaling to obtain the MVF is used in the g-ratio formula, the g-ratio in healthy white matter is relatively constant. However, when lesions exist, the contrast using the different MVF markers is very different. In the MS patient, the mean g-ratio in normal appearing white matter (NAWM) was 0.76 for all three MVF markers. In MS lesions, the mean g-ratio was 0.65, 0.80, and 0.80, for MTR, MTsat, and F, respectively. Fig. 7 shows the spatial distribution of g-ratios in the MS patient for the three MVF markers.

Fig. 8 shows the theoretical effect of having an improper slope (top) or intercept (bottom) in the relationship between an arbitrary myelin marker and the MVF, in the case where the (theoretical) relationship is in fact linear. The plots show that the computed g-ratio becomes fiber density dependent, in addition to being incorrect.

3.3.3 MVF calibration: Discussion

The MTR is a commonly used myelin marker, however, due to T1 sensitivity, it lacks dynamic range. This results in unrealistic g-ratios in MS lesions that are lower than in NAWM. T1 is one of the possible MR-based myelin markers, but in the context of the MTR experiment, it confounds the contrast, because the MT effect dominates but is diminished by the T1 contrast, which works against it. MTsat correlates more highly with F, which is obtained from an explicit qMT model designed to isolate the macromolecular tissue content. It is important to note, however, that this correlation may be driven to some extent by the different B1 sensitivities of the techniques as the MTR was not corrected for B1 induced variability (Volz et al., 2010; Yarnykh & Khodanovich, 2015). Independent of this demonstration of the potential of MTsat for myelin mapping, researchers have found that MTsat may be more sensitive to tissue damage than MTR in multiple sclerosis, with higher correlation with disability metrics (Lema et al., 2017). MTsat has recently been used by other groups in g-ratio imaging of healthy adults (Mohammadi et al., 2015).

If the MVF is miscalibrated in this g-ratio imaging formulation, there will be a residual dependence on fiber volume fraction in our formulation. This reduces the power of the g-ratio metric, which ideally is completely decoupled from the fiber density. Independent of specificity of the myelin marker, if the myelin calibration is inaccurate, this residual dependence on fiber volume fraction occurs. It is clear that the g-ratio metric we will compute is g-ratio weighted, and the better the calibration, the more weighted to the g-ratio it will be. Until quantitative myelin mapping is accurate, the g-ratio metric will not be specific to the g-ratio.

There is evidence that fiber density drops precipitously in some MS lesions (Stikov et al., 2015b). In Fig. 7, this was the case, and we can see that MTR does not drop enough, making MS lesions appear to have a lowered g-ratio instead of a higher g-ratio as expected. Inspecting the bottom (red) curve in Fig. 8, we see that even if there is a linear relationship between the myelin marker of choice and the MVF, miscalibration leads to an apparent g-ratio metric that is elevated in regions of lower fiber density, and significantly lower in regions of healthy fiber density. This occurs when in fact all of the fibers have the same g-ratio, and could easily be interpreted as hypomyelination in an MS subject or population.

As noted above, the problem of miscalibration of the myelin marker exists independent of the specificity of the myelin marker. However, specificity is itself a major confound, as previously detailed. In the case of demyelinating disease, one must also consider that, with many myelin markers, all myelin will affect the MR signal, even if it is not part of an intact fiber. Research indicates that in MS, there is acute demyelination followed by a period of clearance of myelin debris, followed by effective remyelination. During clearance, remyelination can occur, but this myelin is of poor quality (Lampron et al., 2015). On the scale of an MRI voxel, there can be myelin debris, poor remyelination, and higher quality remyelination. The extent to which myelin debris affects the myelin volume estimates may depend on the myelin mapping technique chosen. It is also not clear how well the estimates of MVF, and also AVF, behave at very low fiber density.

One possible solution for MVF calibration is to calibrate the g-ratio to a known value in certain regions of interest (Mohammadi et al., 2015; Cercignani et al., 2016a), as mentioned above. However, care must be taken that this step is not adjusting for differences in the diffusion part of the pipeline (e.g., different implementations of the diffusion model (Cercignani et al., 2016a)), and therefore still leaving a fiber density dependence. Additionally, the correct value in these regions of interest must be known. Calibration based on expected MVF would remove this sensitivity, but is subject to error due to partial volume averaging of white matter with other tissue. If the relationship between the myelin-sensitive metric and the MVF is not a simple scaling, such calibration will fail. Particular care needs to be taken when studying disease.

If the assumed relationship between the myelin marker and the MVF is incorrect, the computed g-ratio will be incorrect. Is it possible to compute a g-ratio that is correct to within a scaling factor, and not sensitive to the fiber density? This would require that the AVF or FVF be estimated independent of the MVF. Simple models such as the diffusion tensor, apparent fiber density (Raffelt et al., 2012), and tensor fiber density (Reisert et al., 2013), are indicators of fiber or axon density, but detailed modeling is most likely superior. Consideration of contrasts other than diffusion MRI, such as gradient-echo based approaches (Sati et al., 2013), might also help with this problem. The g-ratio is a function of the ratio of the MVF to the AVF, and a technique that measures this ratio directly would be optimal. However, due to the extremely short T${}_{2}^{*}$ of myelin, GRE based estimates would be of the myelin and axon water fraction, and hence would still need to be calibrated using the volumetric occupancy of water in these tissues.

In summary, both specificity and accuracy are important for both AVF and MVF estimation. For both AVF and MVF estimation, more sophisticated models may be required. For example, we have thus far ignored cell membranes. The axon membrane should technically be included in the AVF, and its volume is up to 4% of the AVF (Sepehrband et al., 2015), but it would most likely be included in the MVF using MT-based MVF estimation. It is also important to keep in mind the different sensitivities of different MRI markers (e.g., MWF from relaxometry vs. MT-based parameters), because these will not be expected to give identical gMRI metrics. The determination of which MRI markers are most optimal for g-ratio imaging is still the topic of active research (Ellerbrock & Mohammadi, 2017).

3.4 g-Ratio distribution

The g-ratio imaging paradigm extracts a single g-ratio metric per voxel. At typical imaging resolution feasible for the constituent MR images, a voxel contains hundreds of thousands of axons. As with the axon diameter, the g-ratio really applies to an individual axon, and takes on a broad distribution of values in tissue (Graf von Keyserlingk & Schramm, 1984) (see Fig. 9, which shows the g-ratio distribution in the macaque corpus callosum, measured using electron microscopy). The range of myelination includes some unmyelinated axons within healthy white matter. The g-ratio distribution may broaden and become bi-modal in disease. Even within a single axon with intact myelin, the g-ratio may vary due to organelle swelling. Fiber bundles that cross within one voxel may have different g-ratio distributions. In development, some fibers within one fiber bundle will fully develop, while others will be pruned, resulting in an interim bimodal g-ratio distribution within the fascicle. The current MRI-based g-ratio framework will not be able to distinguish these cases, as it reports only an intermediate g-ratio value. It is robust to crossing fibers, in that it will report the same intermediate g-ratio value whether the separate bundles cross or lie parallel to each other. The broad g-ratio distribution is in part a resolution problem, but the g-ratio is expected to be heterogeneous on a scale smaller than we can hope to resolve with MRI.

The aggregate g-ratio we compute in the case of a distribution of values is not precisely fiber- or axon-area weighted, but is close to axon area weighted within a reasonable range of values (West et al., 2016b). Larger axons will have a greater weight in the aggregate g-ratio metric we measure. Simply put, the aggregate g-ratio is the g-ratio one would measure if all axons had the same g-ratio.

In the case of an ambiguous g-ratio distribution, what techniques can we use to infer what situation is occurring? In multiple sclerosis, for example, two possible scenarios probably occur frequently. One is patchy demyelination, on a scale much smaller than a voxel and smaller than the diffusion distance, and the other is more extensively and uniformly distributed thin myelin. These two scenarios could give rise to equal AVF, MVF, and aggregate g-ratio measurements. One possible way to differentiate these cases could be to look more closely at parameters available to us from diffusion models. It has been shown that the extra-axonal perpendicular diffusivity is relatively unchanged by patchy demyelination in a demyelinating mouse model (Jelescu et al., 2016b), because diffusing molecules encounter normal hindrance to motion on most of their trajectory, whereas the axon water fraction is sensitive to this patchy demyelination. Hence, the discrepancy between these two measures can be taken as a measure of patchy demyelination. Alternatively, one can scan subjects longitudinally and infer disease progression. From the ambiguous timepoint described above, the axons in the patches that are demyelinated may die, leaving a decreased AVF and MVF, and a return to a near-healthy g-ratio. In the case of uniformly thin myelin, the remyelination may continue, leaving a near-healthy AVF, MVF, and g-ratio. Note that the g-ratio metric still does not distinguish these pathologically distinct cases. There are two unknowns - the fiber density and the g-ratio (or, alternately, the MVF and the AVF), and one must consider both to have a full picture of the tissue. Looking at the time courses, one can hypothesize what the g-ratio distribution was at the first timepoint.

It would be technically challenging to measure the g-ratio distribution in vivo. Even with an estimate of a distribution of diffusion properties, and an estimate of the distribution of a myelin-sensitive metric, the g-ratio distribution is ill-defined. However, several recent acquisition strategies may help us get closer to this aim. One approach is to take advantage of the distinguishable diffusion signal between different fiber orientations. In the IR-prepared diffusion acquisition described above (De Santis et al., 2016a), the model specifies multiple fiber populations with distinct orientations, each with its own T1 value. This means the diffusion properties, including the restricted pool fraction (a marker of intra-axonal signal from the CHARMED model), are paired with a corresponding T1 for each fiber orientation. Hence, a g-ratio metric could be computed for each fiber orientation. This could be of benefit in, e.g., microstructure informed white matter fiber tractography (e.g., (Girard et al., 2015)) of fiber populations with distinct g-ratios. “Jumping” from one fiber population to another is very common in tractography (Campbell & Pike, 2013; Descoteaux et al., 2016; Maier-Hein et al., 2016), and constraining tractography to pathways with consistent microstructural features could help reduce false positives in regions of closely intermingling tract systems.

It may be possible, conceptually, to estimate the g-ratio distribution via a 2D diffusion-relaxation spectroscopic approach. While extremely acquisition intensive, 2D spectroscopy of T2 and the diffusion coefficient (Callaghan et al., 2003) has been demonstrated recently as a probe of microstructure (Kim et al., 2016). The acquisition involves making all diffusion measurements at different echo times. If a distribution of a myelin volume sensitive metric (here, T2) can be estimated simultaneously in 2D with a distribution of a diffusion-based metric sensitive to the axon volume, it may be possible to infer the distribution of g-ratios.

This has been an incomplete but useful list of pitfalls. Now, we will consider the promise of imaging the aggregate g-ratio weighted metric, despite its pitfalls. g-Ratio imaging is being explored in many different contexts, described below.

4 The promise: g-ratio imaging studies

The promise of g-ratio imaging is its potential to provide us with valuable in vivo estimates of relative myelination. In the last few years, studies showing the potential of this framework have begun to emerge.

4.1 Healthy white matter

Fig. 10 shows an image of gMRI in healthy white matter using our qMT and NODDI g-ratio protocol (see section 2.5). With our protocol, in healthy subjects, the gMRI map is relatively flat, with a mean gMRI of 0.76 (SD=0.05). Other groups have explored these and other MVF and AVF sensitive contrasts for g-ratio mapping in healthy white matter. These include a study of the effects of age and gender in a population of subjects aged 20 to 76 using qMT and NODDI (Cercignani et al., 2016b), studies of healthy adults using MTsat and the TFD (Mohammadi et al., 2015) and MTV and DTI (Berman et al., 2017), and a study of healthy subjects using the ViSTa myelin water imaging technique and NODDI (Jung et al., 2016).

A variation of the g-ratio with age appears to be detectable with this methodology (Cercignani et al., 2016b). A variation with gender has not been seen, and if it exists in adolescence (Paus & Toro, 2009), a study designed for sufficient statistical power at a precise age will be required to detect it. In addition to exploring the effect of age and gender, spatial variability of the g-ratio has been investigated. An elevated g-ratio at the splenium of the corpus callosum has been seen (Stikov et al., 2015b; Mohammadi et al., 2015). The splenium has been reported to contain axons of very large diameter (LaMantia & Rakic, 1990), and these would be expected, due to the nonlinearity of the g-ratio (Hildebrand & Hahn, 1978), to have relatively thinner myelin sheaths. Electron microscopy in the macaque (Stikov et al., 2015b) (see Fig. 11) confirms this; the “super-axons” dominate the aggregate g-ratio measure, which was seen to be elevated in the splenium using both EM measurements and MRI of the same tissue (Stikov et al., 2015b).

g-Ratio imaging has also been performed in the healthy human spinal cord (Duval et al., 2015), where there are considerable technical challenges, such as motion, susceptibility, and the need for significantly higher resolution than we have described for cerebral applications. Duval et al. acquired g-ratio data at 0.8 mm x 0.8 mm inplane voxel size. This study used the CHARMED model of diffusion, more accessible on scanners with high gradient strength, on a CONNECTOM scanner. It used the MTV myelin marker. Of interest, the g-ratio was not found to vary significantly across white matter tracts in the spinal cord, while the diffusion metrics (restricted fraction, diffusivity of the hindered compartment, and axon diameter) and the MTV metric did vary across tracts. This is expected, as heterogeneity in packing and axon diameter is expected to be greater than heterogeneity of the g-ratio, and the g-ratio is also robust to partial voluming effects.

Multiple groups have studied the g-ratio in vivo in the developing brain (Dean et al., 2016; Melbourne et al., 2016). Axon growth outpaces myelination during development, and therefore a decreasing g-ratio is expected as myelination reaches maturity, as was seen in these studies.

4.2 Multiple sclerosis

Imaging the g-ratio in vivo in multiple sclerosis has been explored by several groups (Stikov et al., 2015b; Cercignani et al., 2015; Hori et al., 2016) and is of interest for several reasons. It can possibly help assess disease evolution, and can help monitor response to treatment. It has the potential to aid in the development of new therapies for remyelination. It can also help us understand which therapies might be more fruitful avenues of research. While currently available therapies are immunomodulatory or immunosuppressive in nature, several remyelinating therapies are in clinical trials. Preventing MS-related demyelination (via immunomodulatory therapies) is always preferred at earlier disease stages. Once demyelination has occurred, remyelinating therapies may help protect demyelinated (but not transected) axons from delayed degeneration due to the loss of trophic factors no longer received from myelin. Given that remyelination can occur only when demyelinated but still viable axons are present, interpretation of MRI markers of remyelination would be improved by including a marker of axonal integrity. Thus, the g-ratio framework would be useful in the evaluation of remyelinating therapies. In particular, dynamic changes of g-ratio over time could be measured within new and chronic lesions (as detected on conventional MRI) and the (temporally aligned) timecourses compared between treatment arms.

Despite the promise of imaging the g-ratio in vivo with MS, it is important to remember the pitfall of specificity and miscalibration of the myelin metric when interpreting g-ratio estimates in MS. In our own preliminary experience imaging MS, the relationship between the g-ratio in lesions and in NAWM appears to be complex. In section 3.3, we detailed our observations in one MS patient for whom we acquired a full qMT protocol. In a total of four subjects for whom we have acquired MTsat data, we observed heterogeneity in the gMRI values within and across subjects. The average ($\pm$ standard deviation) g-ratio in lesions compared to NAWM was not always higher than that in NAWM, with (gMRI(lesions),gMRI(NAWM))=(0.80$\pm 0.07$,0.76$\pm$0.05), (0.72$\pm$0.05,0.77$\pm$0.03),

(0.74$\pm$0.06,0.75$\pm$0.05), (0.75$\pm$0.08,0.78$\pm$0.04) for the four subjects. On average in these patients, gMRI was 0.75 in lesions and 0.76 in NAWM. This variability of gMRI in lesions compared to NAWM is consistent with preliminary findings from other groups using MT contrast for MVF (Cercignani et al., 2015), and could indicate that variable levels of non-myelin macromolecular content may confound the g-ratio metric. The g-ratio itself is not expected to be lower than in NAWM, but may appear so because other macromolecular content is confounding the MT measurements.

Further studies of MS are ongoing, including pediatric populations, optic neuritis, and studies investigating whether gadolinium enhancing lesions have a distinct g-ratio.

4.3 Other potential applications

g-Ratio imaging has potential to aid in the understanding and treatment of multiple other diseases. White matter abnormalities may underlie many developmental disorders. These include Pelizaeus Merzbacher disease and Sturge-Weber syndrome (Hori, 2016). An increased apparent g-ratio could result from axonal changes that occur with intact myelin (for example, axonal swelling due to infarction). g-Ratio differences have been seen in schizophrenia using electron microscopy (Uranova et al., 2001), and researchers hope to be able to study such changes in vivo in schizophrenia and other psychiatric disorders. Another potential application of g-ratio imaging is bridging the gap between microstructure and large-scale functional measures such as conduction delays. Adding a g-ratio measure to the human connectome could result in a framework for evaluating delays. Finally, g-ratio imaging offers the possibility to quantify plastic changes in myelin thickness due to learning and adaptation to injury or disease. MRI has been used to measure such changes, with the hypothesis that myelin thickness is changing (Reid et al., 2016), but quantitative g-ratio measurement could help describe the neural changes in plasticity in more detail.

Of note, when examining the g-ratio compared to conduction velocity or functional metrics, linear correlation would not be expected to be a good statistic. The g-ratio has an optimal value, hence, both decreases and increases from this optimum would be expected to decrease function. In the future, assuming sufficient reproducibility, establishment of an age-dependent atlas of normal g-ratio values could help determine whether an individual lies within the normal range.

4.4 Ex vivo g-ratio imaging

One major application of ex vivo g-ratio imaging is validation of the technique for use in vivo, which should help elucidate the true promise of g-ratio imaging. Ex vivo validation has been performed, in order to investigate the g-ratio explicitly (Stikov et al., 2015b, a; West et al., 2016a), or one or both of the individual metrics used to compute it (Duval et al., 2016b; Jelescu et al., 2016b; Jespersen et al., 2010; West et al., 2017, 2014; Schmierer et al., 2008; Scherrer et al., 2016; Sepehrband et al., 2015; Wood et al., 2016). These studies compare in vivo or ex vivo MRI metrics to electron microscopy, optical microscopy, myelin staining, immunohistochemistry, and coherent anti-Stokes Raman spectroscopy (CARS). While no microscopy technique is perfect, microscopy provides a reasonable validation for imaging techniques, taking into account the possibility for tissue shrinkage and distortion, limitations in contrast and resolution, and segmentation techniques (Zaimi et al., 2016).

Interpretation of findings of demyelinating models should take into account the particularities of the demyelinating challenge. Jelescu et al. have shown that the extra-axonal diffusivity perpendicular to axons correlates with the g-ratio in a cuprizone demyelinating model in mice (Jelescu et al., 2016b). This is probably driven by a fiber volume fraction decrease, because little axon loss would be expected in this model. In other words, the extra-axonal diffusivity is not specific to the g-ratio per se, but to the fiber volume fraction (i.e., size of the extracellular space), but these two quantities correlate highly in this particular case. Similarly, West et al. have shown a correlation between the discrepancy between F and MWF and the g-ratio in a knockout model in mice (West et al., 2014). This is probably a correlation with absolute myelin thickness, via exchange effects, as opposed to the g-ratio per se.

Ex vivo g-ratio imaging could also be of use for the study of neural tissue independent of validation of in vivo imaging. MRI is far more suited to imaging large samples of tissue than are other ex vivo imaging techniques, such as electron microscopy.

5 Conclusion

Computing a g-ratio metric is an effective way to interpret any combination of myelin-weighted and axon/fiber-weighted MR data. In this article, we have discussed the considerable promise of g-ratio imaging to help us understand disease, develop therapies, and monitor disease progression. Additionally, we have shown how the g-ratio framework can provide a window onto the study of normal brain variability, development, aging, plasticity, and functional dynamics. We have also explored the pitfalls of g-ratio imaging, which include MR artifacts, lack of specificity, low spatial resolution, and long acquisition times. Keeping these confounds in mind, it is clear that what we are currently measuring is an aggregate g-ratio weighted metric that is strongly dependent on the MRI markers used to compute it. The framework described in this article provides information on two quantities: the fiber density and the g-ratio (equivalently, the myelin and axon volume fractions), and attempts to decouple these two quantities to the best of the ability of our current imaging technology. To improve upon this description of the microstructure, one needs to fully decouple the g-ratio from the fiber density, and to provide precise and accurate measures of the myelin and the axon volume fractions. This task is at the frontier of microstructural MRI, showing the way for the future of multi-modal brain imaging.

Acknowledgements

The authors would like to thank Tomas Paus, Robert Dougherty, Eva Alonso-Ortiz, J.F. Cabana, Christine Tardif, Jessica Dubois, Dmitry Novikov, Ofer Pasternak, Atef Badji, Robert Brown, Masaaki Hori, and David Rudko for their insights and contributions to this work. This work was supported by grants from Campus Alberta Innovates, the Canadian Institutes for Health Research (GBP, FDN-143290, and JCA, FDN-143263), the Natural Science and Engineering Research Council of Canada (NS, 2016-06774, and JCA, 435897-2013), the Montreal Heart Institute Foundation, the Fonds de Recherche du Québec - Santé (JCA, 28826), the Quebec BioImaging Network (NS, 8436-0501), the Canada Research Chair in Quantitative Magnetic Resonance Imaging (JCA), and the Fonds de Recherche du Québec - Nature et Technologies (JCA, 2015-PR-182754).