Orientationally invariant indices of axon diameter and density from diffusion MRI

NeuroImage 52 (2010) 1374–1389

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g

Orientationally invariant indices of axon diameter and density from diffusion MRI Daniel C. Alexander a,⁎, Penny L. Hubbard b, Matt G. Hall a, Elizabeth A. Moore c, Maurice Ptito d, Geoff J.M. Parker b, Tim B. Dyrby e a

Centre for Medical Image Computing, Department of Computer Science, University College London, Gower Street, London WC1E 6BT, UK Imaging Science and Biomedical Engineering, School of Cancer and Imaging Sciences, University of Manchester, Manchester, M13 9PT, UK Philips Healthcare, Best, Eindhoven, Netherlands d School of Optometry, University of Montreal, Montreal, Canada e Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark b c

a r t i c l e

i n f o

Article history: Received 3 December 2009 Revised 14 May 2010 Accepted 16 May 2010 Available online 23 May 2010 Keywords: Diffusion MRI Active imaging Axon diameter Fibre density White matter Microstructure

a b s t r a c t This paper proposes and tests a technique for imaging orientationally invariant indices of axon diameter and density in white matter using diffusion magnetic resonance imaging. Such indices potentially provide more specific markers of white matter microstructure than standard indices from diffusion tensor imaging. Orientational invariance allows for combination with tractography and presents new opportunities for mapping brain connectivity and quantifying disease processes. The technique uses a four-compartment tissue model combined with an optimized multishell high-angular-resolution pulsed-gradient-spin-echo acquisition. We test the method in simulation, on fixed monkey brains using a preclinical scanner and on live human brains using a clinical 3 T scanner. The human data take about one hour to acquire. The simulation experiments show that both monkey and human protocols distinguish distributions of axon diameters that occur naturally in white matter. We compare the axon diameter index with the mean axon diameter weighted by axon volume. The index differs from this mean and is protocol dependent, but correlation is good for the monkey protocol and weaker, but discernible, for the human protocol where greater diffusivity and lower gradient strength limit sensitivity to only the largest axons. Maps of axon diameter and density indices from the monkey and human data in the corpus callosum and corticospinal tract reflect known trends from histology. The results show orientationally invariant sensitivity to natural axon diameter distributions for the first time with both specialist and clinical hardware. This demonstration motivates further refinement, validation, and evaluation of the precise nature of the indices and the influence of potential confounds. © 2010 Elsevier Inc. All rights reserved.

Introduction White matter is the cabling that supports communication throughout the brain. It consists of bundles of axons packed to densities often over 105 mm−2 (Waxman et al., 1995; Aboitiz et al., 1992a). Most axons have a diameter between 0.2 and 20 μm and lengths vary from millimeters to over a meter (Waxman et al., 1995). Axon diameter determines conduction velocity (Ritchie, 1982). Largediameter axons, for example, in the corticospinal tract (CST) or midbody of the corpus callosum, transfer information quickly and support rapid communication for synchronous processing of sensorimotor stimuli. Small-diameter axons, for example, in the genu, pack more densely allowing fibre bundles to carry more diverse information between higher processing areas that favour quantity of information over transmission speed (Ptito, 2003; Aboitiz et al.,

⁎ Corresponding author. E-mail address: [email protected] (D.C. Alexander). 1053-8119/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2010.05.043

1992a; Lamantia and Rakic, 1990). Thus, axon diameter and density provide information about the role and performance of white matter pathways. Diffusion magnetic resonance imaging (MRI) offers unique insight into live tissue microstructure through its sensitivity to displacements of water molecules over millisecond timescales. Semipermeable barriers within tissue affect the dispersion pattern of water and, consequently, diffusion MRI measurements are sensitive to the geometry and organization of the barriers. Diffusion tensor imaging (DTI) (Basser et al., 1994) models the distribution of particle displacements with a Gaussian distribution. The apparent diffusion tensor (DT) provides statistics of mean diffusivity and diffusion anisotropy (Basser and Pierpaoli, 1996), which provide some insight into tissue microstructure and have become popular as markers of white matter integrity. However, a limitation of these markers is that they do not relate directly to features of tissue microstructure and are sensitive to a variety of different effects simultaneously. For example, the size and packing density of cells, the permeability of cell walls and membranes, and the distribution of orientations of anisotropic cells all

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

affect the mean diffusivity and anisotropy. Hence, a change in these statistics is difficult to associate with specific changes in tissue microstructure. More sophisticated approaches model the geometry and material properties of tissue microstructure to predict the dispersion pattern of water within and thus the diffusion MR signal. Modelling approaches can provide sensible estimates of more specific microstructural parameters. Stanisz et al. (1997) model bovine optic nerve tissue as a three-compartment system with one population of water inside elongated ellipsoidal axons, another inside spherical glial cells, and a third in the extracellular space. Each compartment has its own dimensions, volume fraction, membrane permeability, and internal diffusivity and relaxivity. Stanisz et al. (1997) fit the threecompartment model to 800 nuclear magnetic resonance (NMR) measurements of a fixed sample using a stimulated echo sequence (Merboldt et al., 1991) with 3 ms pulses and varying diffusion time and gradient strength. Estimates of the cell sizes and densities agree with manual estimates from microscope images of the same samples. Permeability, diffusivity, and relaxivity estimates are all sensible but unvalidated. The simpler CHARMED model (Assaf et al., 2004, 2008; Assaf and Basser, 2005) has only two compartments: impermeable, parallel cylindrical axons in a homogeneous extracellular space. Each compartment has its own diffusivity, but relaxation times are the same in both. The diameter of the cylindrical axons has a twoparameter gamma distribution. The original CHARMED framework (Assaf et al., 2004; Assaf and Basser, 2005) provides maps of volume fraction between the intra- and extracellular compartments. The technique uses pulsed-gradient spin-echo (PGSE) (Stejskal and Tanner, 1965) in high-angular-resolution diffusion imaging (HARDI) (Jones et al., 1999; Tuch et al., 2002) with several b-values (multishell HARDI) to allow model fitting for arbitrary fibre orientation. Assaf and Basser (2005) demonstrate the method in live humans using a clinical MR scanner. However, the technique does not attempt to estimate the axon diameter. The axon diameter parameters of the model are fixed during fitting, and the acquisition protocol uses only a single diffusion time, so it is not designed to support such an estimate. The volume fraction estimates may correlate with axon density, but a true axon density estimate requires knowledge of the axon diameters. Later work (Assaf et al., 2008) uses the CHARMED model to estimate the axon diameter distribution, but abandons the HARDI acquisition in favour of a fixed gradient direction to allow multiple combinations of diffusion time and gradient strength for sensitivity to axon diameter. They perform a similar NMR experiment to Stanisz et al. (1997) and estimate distributions of axon diameters in bovine optic and sciatic nerve samples that agree with distributions measured by hand on microscope images of the same samples. They also combine the method with imaging and show that segmentation of a pig spinal cord image from clustering the gamma distribution parameters shows consistency with a segmentation using various histological stains. Barazany et al. (2009) test the method in the corpus callosum of a live rat. They add an isotropic–diffusion compartment to account for partial volume with cerebrospinal fluid (CSF) for in vivo imaging. Results show spatial variation in the estimated axon diameter distribution that reflects subsequent histological evaluation. Several limitations prevent direct translation of the existing axon diameter estimation techniques (Stanisz et al., 1997; Assaf et al., 2008; Barazany et al., 2009) above to whole-brain mapping in live human subjects: 1. The imaging protocols require orders of magnitude higher magnetic field gradient strength than current human MRI scanners can provide. The NMR experiments in Stanisz et al. (1997) and Assaf et al. (2008) use maximum gradient strength |G|max around 1 T m−1 and the imaging experiments in Assaf et al. (2008) and

1375

Barazany et al. (2009) use 300 mT m−1, whereas current human systems provide |G|max between 40 and 80 mT m−1. 2. The acquisition time is too long for live human volunteers or patients to tolerate. The MRI experiments in Assaf et al. (2008) require around a day of imaging time. Even the in vivo experiment in Barazany et al. (2009) requires 2 hours of imaging. Human volunteers tolerate little more than 1 hour. 3. The methods require prior knowledge of the fibre orientation, because the gradient direction must be perpendicular to the fibre in all the measurements. This limits parameter estimation to specific targeted structures, since the orientation of white matter fibres varies widely over the brain. The simulation study in Alexander (2008) tests feasibility of relaxing these limitations using clinically feasible HARDI combined with a simplified version of the CHARMED model that has a singleaxon diameter rather than a gamma distribution. Unlike previous multishell HARDI acquisitions, all three variables of the PGSE sequence (gradient strength |G|, pulse width δ, and separation Δ) vary between different shells, which provides the sensitivity to axon diameter. The simulation study uses an optimization algorithm to identify the combination of HARDI shells that maximize sensitivity to the model parameters. The experiment design optimization considers only protocols that divide the total number of measurements into M shells of N measurements. The algorithm optimizes the M combinations of PGSE settings with the directions held fixed. The objective function is the normalized sum of Cramer–Rao lower bounds on the model parameters averaged over a range of a priori settings. The optimization also identifies the best tradeoff between N and M with the total number of measurements NM fixed. Orientational invariance proves feasible, because M is small, typically 3 or 4, in contrast to previous protocols (Stanisz et al., 1997; Assaf et al., 2008; Barazany et al., 2009), which have several hundred distinct combinations. Simulations with N = 30, M = 4 and maximum gradient strength |G|max = 70 mT m− 1 suggest feasibility of recovering axon diameter and density without knowledge of the fibre orientation. More precisely, they show that, with protocols optimized for typical constraints of a human scanner, posterior distributions on single axon-diameters have low variance when the true diameter is between 10 and 40 μm. However, posterior distributions for true diameters between 2 and 4 μm are all similar and close to uniform over the range 0 to 5 μm. The result suggests that axon diameters of about 5 μm or less are indistinguishable from one another, although they are identifiable as small compared to diameters greater than 5 μm. Diameters over 5 μm (up to around 40 μm) can be identified more precisely. Sensitivity extends to lower axon diameters as gradient strength increases and to higher diameters as T2 increases (allowing longer diffusion times). As diffusivity decreases, the window of sensitivity shifts to lower diameters. The aim of this article is to test the ideas in Alexander (2008) for mapping orientationally invariant indices of axon diameter and density in biological tissue. An important difference between real white matter and the simple model underlying the earlier simulation experiments is that the tissue contains a distribution of axon diameters. The main contribution here is an evaluation of the method in the presence of distributed axon diameter, which includes preliminary results from brain data. The method differs from Alexander (2008) in a few ways: (i) slightly improved a priori parameter settings for the experiment design optimization; (ii) an additional CSF compartment, as in Barazany et al. (2009), and an isotropically restricting compartment, similar to Stanisz et al. (1997), in the model; and (iii) grid search and gradient descent stages in the fitting procedure to initialize the Markov Chain Monte Carlo (MCMC). Key differences of the method with Assaf et al. (2004) and Assaf and Basser (2005) are (i) simplification of the two-compartment CHARMED model and inclusion of the CSF and isotropically restricting

1376

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

compartments; (ii) estimation of the axon diameter; (iii) optimized, rather than ad hoc, experiment design that explores all possible PGSE combinations; (iv) the Bayesian parameter estimation. Key differences with Stanisz et al. (1997), Assaf et al. (2008), and Barazany et al. (2009) are (i) orientational invariance, (ii) the model adaptations, (iii) the optimized experiment design, (iv) the Bayesian estimation, and (v) rather than a model of the full axon diameter distribution, as in Assaf et al. (2008) and Barazany et al. (2009), we estimate a single summary statistic, which we call the axon diameter index; Stanisz et al. (1997) also estimates a single axon-diameter. The methods section details the experiment design, data acquisition, tissue model, preprocessing, and parameter estimation procedures. Experiments then evaluate the method in simulation and show parameter maps from scan–rescan experiments on human subjects and fixed monkey brains. Finally, we conclude and discuss limitations, applications, and areas for further work. Methods This section describes the samples and imaging protocols for fixed and live brain imaging. It then specifies the tissue model, details the modelfitting procedure, and defines the indices of axon diameter and density. Tissue samples Test data sets for in vivo imaging come from two healthy human volunteers. Subject 1 is a female, 28 years old, and weighs 54 kg; subject 2 is a male, 51 years old, and weighs 105 kg. The subjects were scanned with informed consent and the approval of the local research ethics committee. Monkey data sets come from two normal appearing young-adult female vervet monkey (Chlorocebus aethiops) brains with matched age of 32 months. Perfusion fixated postmortem brains obtained from the Montreal Monkey Brain Bank were prepared for scanning, as described by Dyrby et al. (2010), in Copenhagen, Denmark. The monkeys were handled and cared for on the island of St. Kitts according to a protocol approved by the local ethics committee (The Caribbean Primate Center of St. Kitts). Imaging protocols Live human brain imaging uses a clinical 3 T Philips Achieva system with |G|max = 60 mT m− 1. The Philips scanner can acquire around 360 diffusion-weighted images in 1 hour, which we set as the acquisition time limit for in vivo imaging. The experiment design follows the procedure in Alexander (2008) with the following a priori model parameter settings: intracellular volume fraction f1 = 0.7, free diffusivity d∥ = 1.7 × 10− 9 m2 s− 1, axon diameter a = 10, 20, and 40 μm, and perpendicular extracellular diffusivity of d⊥ = 1.2 × 10− 9 m2 s− 1. The optimized protocol divides the 360 measurements into M = 4 HARDI shells each with N = 90 gradient directions. Fig. 1 shows the four combinations of PGSE settings in the human protocol together with example images for each. The protocol also includes four b = 0 images. It uses sagittal echo-planar imaging (EPI) with cardiac gating, in-plane image matrix 128 × 128 and 1.8 × 1.8-mm2 voxels. Images have ten slices, centred on the midsagittal plane, with thickness 3.9 mm and no gap. The echo time TE = 122 ms, and the repetition time TR = 6 s, depending on heart rate. Echo and repetition times are the same for all measurements. Fixed monkey brain imaging uses a 4.7-T Varian experimental system with maximum gradient strength |G|max =140 mT m− 1. The optimization is similar to the human protocol optimization, but uses lower a priori diffusivity (d∥ = 0.45 × 10− 9 m2 s− 1 and d⊥ = 0.25× 10− 9 m2 s− 1) for fixed tissue. Higher gradient strength and lower diffusivity shift the sensitivity range to lower diameters and we adjust the a priori settings to a =4 and 10 μm.

Fig. 1. Individual pulse sequences in the human protocol and example images for each. The gradient direction is inferior–superior for each image. The pulse sequence variables are the diffusion gradient vector G, the pulse length δ and the pulse separation Δ. The gradient strength |G| in the title of each sequence diagram refers to the strength during the pulse. The diffusion-weighting factor b = (Δ − δ/3)(γδ|G|)2 for the PGSE sequence, where γ is the gyromagnetic ratio. The echo time is the same for each measurement.

Fig. 2 details the monkey protocol, which also has M=4 and N=90. The protocol includes 12 b=0 images. It uses a standard 2D spin-echo with a single-line readout, in-plane image matrix 256×128, 0.4×0.4 mm2 voxels, and three 0.4-mm sagittal slices, centered on the mid-sagittal plane, with an interslice gap of 2.0 mm; TE=60 ms and TR=5 s. Data sets We repeat the acquisition on each human subject in a separate session on a separate day. Each session takes around an hour depending on heart rate. Motion and eddy-current distortion correction among the human images uses FLIRT (Jenkinson et al., 2002) with careful visual monitoring to avoid errors. White matter signal to noise ratio (SNR) at b = 0 is about 20. We repeat the acquisition twice on each monkey brain in the same session and the total acquisition time for each brain is 150 h. The monkey brain data require no alignment (Dyrby et al., 2010). White matter SNR at b = 0 is also around 20.

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

1377

Fig. 2. As Fig. 1, but for the monkey protocol.

Tissue model The white matter model has four tissue compartments or populations of water molecules that each provides a separate normalized MR signal S1, …, S4: 1. Signal S1 comes from intra-axonal water trapped inside parallel cylinders with equal diameter, a. 2. Signal S2 comes from extra-axonal water adjacent to, but outside, the cylinders. 3. Signal S3 comes from CSF, where the cylinders do not affect diffusion. 4. Signal S4 comes from stationary water trapped in glial cells and other small compartments or bound to membranes and other subcellular structures. No exchange between the populations of water molecules occurs. The full model for the diffusion MRI signal S⋆ is then ⋆

S =

⋆ S0

4

∑ fi S i ;

i=1

ð1Þ

where S⋆0 is the MR signal with no diffusion weighting, fi is the proportion of water molecules in population i, 0 ≤ fi ≤ 1 and ∑4i= 1fi = 1.

The model for S1 is the Gaussian phase distribution approximation (Murday and Cotts, 1968) of the signal from particles trapped inside a cylinder (Van Gelderen et al., 1994). Intra-axonal diffusion is unhindered along the axis of the cylinders. The model for S2 is the DT model (Basser et al., 1994), i.e., anisotropic Gaussian distributed displacements, with diffusion coefficient d∥ in the direction of the axons and d⊥ in all perpendicular directions. The parallel diffusivity, d∥, is the same as the intrinsic diffusivity inside the cylinders in the model for S1, following (Alexander, 2008). A simple tortuosity model (Szafer et al., 1995) sets d⊥ = d∥ ð1−νÞ;

ð2Þ

where ν = f1 = ðf1 + f2 Þ

ð3Þ

is the intra-axonal volume fraction within the white matter, i.e., ignoring the CSF and stationary compartments, which do not contribute to S2. The model for S3 is isotropic Gaussian displacements, as in Barazany et al. (2009), with diffusion coefficient dI. The signal from the stationary population remains unattenuated by diffusion weighting so S4 = 1. The full set of parameters of the model is fi, i = 1, …, 4, d∥, dI, a, the axis of the cylinders n, and S⋆0.

1378

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

Previous diffusion studies on fixed tissue (Stanisz et al., 1997; Panagiotaki et al., 2009) motivate the inclusion of the “stationary” compartment. Both studies demonstrate clear departures from Gaussian displacements along the fibre direction. Stanisz et al. (1997) suggest that the departure comes from water trapped within glial cells and Panagiotaki et al. (2009) suggest it may also come from axons with different orientation to the main bundle. Interestingly, however, our in vivo data show the effect much less than our postmortem data. This suggests that other mechanisms underlie the departure. Those mechanisms may include restriction within glial cells and nonparallel axons, but changes in the chemical environment during fixation and loss of active function after death may also contribute. Here, we choose the simplest possible model that captures the effect in the data. We assume the effect is negligible in vivo and fix f4 = 0 in the model for live tissue, which reduces the number of parameters and stabilizes the fitting procedure. Model fitting Model fitting uses the Rician MCMC procedure in Alexander (2008), with the minor adaptations specified below, after a new initial grid search and gradient descent to identify a good starting point. The initial grid search finds a starting point for the gradient descent by picking the combination of parameters that maximize the likelihood of the data, using a Rician noise model, over a regular grid of combinations of physically plausible settings. To reduce the number of combinations, we set n to the principal direction of the best fit DT and fix d∥ and dI to 1.7 and 3.0 × 10−9 m2 s− 1, respectively, for the in-vivo human, and to 0.6 and 2.0 × 10− 9 m2 s− 1 for the fixed monkey brain. The gradient descent algorithm keeps d∥ and dI fixed, but allows n and all the other parameters to vary. Gradient descent uses a Gauss– Newton method, with the same objective function as the grid search, to refine the maximum likelihood parameter estimates. The gradient descent provides a starting point for the MCMC, which keeps all parameters fixed except a, f1, and f2. We collect 40 samples at intervals of 200 iterations after a burn in of 2000 iterations. Final parameter estimates come from the mean over the MCMC samples. We distinguish parameter estimates from the parameters themselves with a tick, so a′ is the estimate of a from the mean of the MCMC samples. Scalar statistics Our index of axon diameter is a′. The model assumes a single a, although tissue contains a distribution of axon diameters. The fitted a′ therefore is the single diameter that produces attenuated signals that best match the average attenuated signals over the distribution. The mean axon diameter E(a) = ∫p(a)ada, where p is the true distribution of axon diameters. However, large diameter axons contain more water and so contribute more to the MR signal than small diameter axons. Cylinder volume is proportional to a2 and so the contribution to the signal of each axon is approximately proportional to the square of its diameter (Packer and Rees, 1972). Thus, we might expect a′ to correlate with 3

α=

∫pðaÞa da ∫pðaÞa2 da

;

ð4Þ

where the normalization in the denominator is simply the intraaxonal volume. Our index of axon density is ρ=

4ν′ ; πa′2

using πa′2/4 as the cross-sectional area of the “average” axon.

ð5Þ

Preprocessing The white matter model assumes all cylinders are parallel so is inappropriate in regions of complex fibre architecture such as crossings or significant splaying, fanning or bending. We limit investigation to voxels in which we expect homogeneous fibre orientation (straight parallel fibres). Thus, we first compute the linearity and planarity (Westin et al., 1999) of the best-fit DT and process only voxels with linearity greater than 0.6 and planarity less than 0.2. We further exclude voxels where S⋆0 is over twice the mean S⋆0 in white matter to avoid large partial volume effects with CSF. The processed region includes most of the central corpus callosum, and the parts of the CST that are not adjacent to or intersecting other fibres. Before fitting, a mean filter smoothes each diffusion-weighted and b = 0 image independently. The basic smoothing kernel is 5 × 5 pixels in plane, but it excludes voxels outside the high-linearity–lowplanarity mask. It also excludes voxels in which the absolute dot product of the principal direction of the best-fit DT with that of the central pixel is less than 0.9 to avoid contributions from fibres with different orientations. Experiments and results This section tests the acquisition protocols and fitting procedures above, first in simulation and then on the human and monkey data sets. Simulations The simulation experiments test the extent to which a′ and ρ distinguish the distributions of axon diameters and densities we expect to see in white matter. Data synthesis The diffusion simulation system in Hall and Alexander (2009), which is freely available in the Camino toolkit (Cook et al., 2006; http://www.camino.org.uk), generates synthetic diffusion MRI data by simulating random walks of molecules within restricting geometries. Here we use packed cylinders with distributions of diameters consistent with anatomical measurements of axon diameter distributions, as in Hall and Alexander (2009). The experiments use a total of 44 virtual white matter environments (substrates) consisting of packed nonabutting impermeable cylinders. Aboitiz et al. (1992a) and Lamantia and Rakic (1990) together provide 11 histograms of axon diameters from different regions of white matter. Shrinkage during histological preparation reduces the axon diameters causing the histograms to exclude the largest naturally occurring diameters. Doubling the diameters in each histogram produces a further 11 histograms and ensures that the whole set covers the realistic range of diameters. The distribution of cylinder diameters in each substrate approximates one of those 22 histograms and each histogram provides two substrates with different packing densities. The first is the maximum packing density we can achieve using one run of the cylinder-packing algorithm in Hall and Alexander (2009) and the other has intracylinder volume fraction (analogous to f1) approximately 0.1 less than that maximum. The simulations include no isotropic or stationary compartments (by analogy, f3 = f4 = 0 so f1 = ν). Each substrate consists of 100 parallel nonabutting cylinders with axes along the z-axis. The substrates are periodic and tile perfectly. We generate synthetic signals for both the human and monkey protocols from each substrate. Each simulation uses 50,000 walkers and 20,000 time steps, which ensures repeatability of synthetic signals well within the noise levels we study. Free diffusivity is

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

1.7 × 10−9 m2 s− 1 for the in vivo simulation and 0.6 × 10− 9 m2 s− 1 for the fixed-tissue simulation. Qualitative assessment Fig. 3 compares signals from substrates with different axon diameter distributions, derived from Aboitiz et al. (1992a), to demonstrate qualitatively that the protocols discriminate the substrates. The substrate cross sections in the top row show larger axons in the middle of the corpus callosum than at either end, which is consistent with the histograms in Aboitiz et al. (1992a). The signal plots from each substrate show differences, which confirm that the signals discriminate these different distributions of cylinder diameters. Differences are clearest when the gradient direction is perpendicular to the fibre direction, i.e., to the left of each plot where |Gz|/|G| = 0. The monkey protocol shows greater perpendicular

1379

signal attenuation in the midbody substrate than the genu or splenium, because the larger axons in the midbody accommodate larger internal displacements. Signals from the genu and splenium substrates appear similar, although the high b signal still shows contrast reflecting the slightly higher α in the genu. Signal variation in the human protocol is weaker, but discernible. The perpendicular signals show little difference between the genu and splenium substrates. However, the three highest b signals all show greater attenuation in the mid-body substrate than the genu and splenium. Parameter estimates Figs. 4 and 5 demonstrate that a′ and ρ discriminate substrates with different diameter distributions and densities. Fig. 4 plots α for all 44 substrates against a′ for both protocols over 10 trials with independent Rician noise added to make the signal to noise ratio with

Fig. 3. Top: cross sections of virtual white matter environments for the genu (left), midbody (middle), and splenium (right) of the human corpus callosum approximating the data in Aboitiz et al. (1992a) after doubling. The substrates are periodic and each image shows a patch with sides of length 100 μm. The smaller squares outline the 100 cylinder basic substrate. The fraction of the volume within the cylinders is the same for each substrate. Middle: plots of normalized synthetic signals from the monkey protocol in Fig. 2 against absolute dot product between the gradient direction and the cylinder axis orientation; signals from perpendicular gradient directions are towards the left and from parallel directions to the right. The colours distinguish the different HARDI shells. Bottom: as middle but for the human protocol in Fig. 1. In both sets of signal plots, two curves overlap and appear as one, because two of the four combinations of PGSE settings are close to identical.

1380

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

no diffusion weighting 20. For completeness, Fig. 4 also plots E(a) against a′. Fig. 5 shows a similar plot for ν′ against the true intracylinder volume fraction, ν, for each substrate. Fig. 4 reveals that, while a′ is not a perfect predictor for α, it does correlate. The correlation is strong for the monkey protocol for substrates with α greater than 2 or 3 μm, although a′ underestimates high α and overestimates low α. The variance of a′ among the 10 trials is higher with the human protocol. However, correlation with α is still noticeable, particularly for the mean over the 10 trials with α greater than 3 or 4 μm. With the human protocol, a′ consistently overestimates α. For both protocols, the estimate a′ is more robust to noise for higher α, as the earlier simulations in Alexander (2008) suggest. As expected, the true mean axon diameter E(a) correlates less well with a′ than α. Fig. 5 shows that ν′ is a good predictor of ν for both protocols. The precision and accuracy of the density index, ρ, therefore depend largely on those of a′. The human protocol produces slightly better estimates of ν than the monkey protocol, which occasionally underestimates ν. Departures of a′ from α reflect differences in the sensitivity of the measurements to different axon diameters. The low |G|max on the human scanner means that only the largest diameters attenuate the signal. Although the precise mechanism requires further study, this effect seems to skew a′ towards higher values than α, probably because α retains contributions from lower axon diameters. The monkey protocol has sensitivity to smaller axon diameters, but lacks sensitivity to very large axon diameters, which causes a′ to underestimate large α and also causes the underestimation of ν. Longer diffusion times and greater diffusivity reveal restriction in larger pores better allowing the human protocol to avoid underestimation of ν. We emphasize that α is merely a useful statistic for comparison rather than the ideal axon

diameter index. The desirable feature of a′ is that it discriminates realistic axon diameter distributions. Partial volume effects Fig. 6 demonstrates the effect of nonzero f3 on a′ and ν′. We simulate partial volume by averaging the synthetic signal from the substrate and S3 = exp(−bdI) with weightings f3 2 {0.1, 0.2, …, 0.9}. For each substrate, we compare estimates a′ and ν′ when f3 N 0 with those when f3 = 0 to reveal any bias that partial volume with CSF may introduce. The top row in Fig. 6 shows that f3′ from the monkey protocol is a good estimate of f3; f3′ from the human protocol correlates with but consistently underestimates f3. The plots in the second row show that the average change in a′ when 0 b f3 ≤ 0.7 compared to f3 = 0 is close to zero for both protocols, which suggests that moderate nonzero f3 does not bias a′. The bottom row shows that, for the monkey protocol, low to moderate f3 has little effect on ν′, but for f3 N 0.5, ν′ becomes biased downward. However, even low nonzero f3 affects ν′ from the human protocol. We conclude that a′ is not biased by moderate CSF contamination, although f3 N 0 may increase the variance of a′ slightly. Thus, a′ and ρ are robust at low f3 with the monkey protocol, but ρ is sensitive to f3 with the human protocol. The human protocol is less able to distinguish the effects of CSF contamination from changes in ν than the monkey protocol. A similar experiment with nonzero f4 (not shown) reveals similar robustness of a′ and ν′ from the monkey protocol to signal from stationary particles. Human data This section shows results from the in vivo human data.

Fig. 4. Top: plots of a′ fitted to synthetic data against α calculated from the histograms of cylinder diameters in the substrates that produce the data. Bottom: a′ against E(a). The plots on the left are for the monkey protocol and on the right for the human protocol. Red circles show estimates from 10 independent trials with noisy data for each substrate; blue crosses show the mean estimates over the 10 noise trials.

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

1381

Fig. 5. As Fig. 4, but plots ν′ against the true intracylinder volume fraction ν for each substrate.

Parameter maps Fig. 7 shows maps of a′ and ρ over midsagittal slices of each acquisition of each human brain. Fig. 8 shows maps of a′, ρ, ν′, and f3′

Fig. 6. Plots showing the effect of nonzero f3 on the fitted model parameters. The left column is for the monkey protocol and the right for the human protocol. The top row plots the estimated isotropic volume fraction f′3 against the true f3. The second row plots the difference Δa′ between the estimate a′ with nonzero f3 and the estimate with f3 = 0, against f3. The third row plots the difference Δν′ between ν′ with and without f3 = 0, against f3.

Fig. 7. Maps of a′ (left) and ρ (right) in the midsagittal slices of the human data sets. The top two rows show results from the two scans of subject 1. The bottom two rows show the two scans from subject 2.

1382

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

Fig. 8. Full parameter maps for a′ (left), ρ (left middle), ν′ (right middle) and f′3 (right) in scan 1 of human subject 1. We show only slices 2 to 9 of the 10 slices, since the two extreme slices contain registration artifacts.

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

over the whole volume for scan 1 of human subject 1. Fig. 9 (left) compares the trends in a′ across regional subsections of the midsagittal corpus callosum for each scan with estimates of α in each region from the histology data in Fig. 4 of Aboitiz et al. (1992a). The midsagittal corpus callosum shows consistent anterior– posterior trends of low–high–low a′ and high–low–high ρ. The pattern persists in slices surrounding the midsagittal slice, but changes in more lateral slices. The maps also reveal low ρ and high a′ in the CST. High a′ in the CST is clearest in subject 2 (labelled in panel g of Fig. 7). Although high a′ is less clear qualitatively in the CST of subject 1, the mean a′ in the midsagittal CST for subject 1 is 9.8 μm. That value is consistent with the mean a′ in other slices and comparable to the highest regional mean in the corpus callosum (region B3), where the mean a′ is 10.3 μm (Fig. 9). The plots in Fig. 9 reveal that a′ consistently overestimates the expected α. The histology data points in Fig. 9 come from singlesubject histograms so require caution in interpretation. Axon diameter histograms in the corpus callosum do vary among subjects (Aboitiz et al., 1992b). In particular, the proportion of large axons increases with age (Aboitiz et al., 1996). However, statistics of the axon diameter distribution from 20 subjects (Aboitiz et al., 1992a) suggest that the low–high–low trend in α is consistent in the normal population. The variation in axon diameter distribution within the normal population is not strong enough to explain the overestimation of α by a′. The overestimation more likely arises from lack of sensitivity to small axon diameters, as the simulation results in Fig. 4 predict. The index a′ is consistently lower in human subject 1 than subject 2. This may reflect age differences, as subject 2 is over 20 years older than subject 1, but we require a much larger cohort to test whether the difference reflects any genuine effect. In both the corpus callosum and the CST, f3′ b 0.1 in all but one or two isolated voxels. Fig. 6 (top right) suggests that f3′ up to 0.1 may arise from f3 up to 0.2 and Fig. 6 (bottom right) shows f3 = 0.2 may cause slight downward bias of ν. This may affect ρ slightly in localized regions, but is unlikely to influence the broad trend significantly across the midsagittal corpus callosum where f3′ is mostly zero.

1383

Fig. 9 shows that the broad trends across the corpus callosum reproduce quite well in subject 1, although individual voxels can show quite large variation, which is consistent with the human protocol simulations. Subject 2 shows greater between-scan variation, which arises in part from poorer alignment of the midsagittal slices, although greater physiological noise in subject 2 and more residual misalignment among the diffusion-weighted images may also contribute. Quality of fit Fig. 10 compares measurements with predictions from the fitted model in several voxels from different parts of the corpus callosum in scan 1 of human subject 1. The examples are typical of measurements and fits in the neighbourhood surrounding each voxel. The model captures the broad trends in the data but consistently underestimates the high b value measurements and the signal reduction is slightly slower than the model as we move from the perpendicular to parallel gradient direction. Although high b measurements with parallel gradient direction may hit the noise floor, the Rician noise model in the fitting minimizes the effect on the parameter estimates. Restriction in glial cells, nonparallel fibres and permeability more likely account for the minor departures. Parameter correlation Fig. 11 shows scatter plots and correlation coefficients for each pair of model parameter estimates from scan 1 of human subject 1. The figure also shows the correlation of each parameter estimate with various DTI indices: fractional anisotropy (FA) (Basser and Pierpaoli, 1996), the average of the two smaller eigenvalues D⊥ and the largest eigenvalue D||. Correlations and scatter plots from the other human scans are similar. The correlations we see correspond to observations in Barazany et al. (2009). Positive correlation coefficients show that FA increases most significantly with ν and ρ, which we expect, because extracellular water contributes less and has more tortuous perpendicular diffusion as ν and ρ increase. The FA decreases weakly with a′, which also makes sense as larger-diameter fibres accommodate greater perpendicular

Fig. 9. Regional averages of a′ in the corpus callosum from the human data sets (left) and monkey data sets (right). The sections are similar to those in Aboitiz et al. (1992a) and Lamantia and Rakic (1990), respectively; the schematic inset on the left illustrates the sections. Dotted lines join data points from the individual scans; solid lines join data points averaged over the two scans. The “Histology” data points show the statistic α in Eq. (4) computed from the histograms of axon diameters in Fig. 4 of Aboitiz et al. (1992a) (left) and Fig. 7 of Lamantia and Rakic (1990) (right). We multiply each diameter in the histology data by 1.5 to correct for shrinkage during histological preparation, as Aboitiz et al. (1992a) suggest. Some of the histograms in Aboitiz et al., (1992a) combine data from two regions so we repeat one point in G3 and B1 (“anterior midbody”) and B3 and I (“posterior body”); the precise location for their “genu” and “splenium” histograms is unclear so we use them to represent all subregions. Lamantia and Rakic (1990) provide full histograms for only three regions of the corpus callosum, but Fig. 8 in Lamantia and Rakic (1990) provides mean axon diameters in all regions. Here, we use a linear model to predict α from the mean in the regions without full histograms. The figure shows α from histograms with black markers and predictions from the model with grey markers.

1384

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

Fig. 10. Plots of the measurements in three example voxels at different locations in the corpus callosum in scan 1 of human subject 1. The top row shows a close up of the map of a′ overlaid on the FA in the midsagittal corpus callosum. Cyan squares highlight three example voxels: one in the genu (left) one in the midbody (middle) and one in the splenium (right). The plots on the bottom row show the measurements in each example voxel in the same left–right order. The solid lines show the predicted signals from the fitted model. Position on the x-axis indicates alignment of the gradient direction with the estimated fibre direction. The black dotted lines show the b = 0 measurements. All measurements are normalized by the single estimate of S⋆0 from the model fitting. The b-values in the legend have the unit seconds per squared millimeter (s mm− 2).

displacements. Similarly, D⊥ increases with a′, but decreases with ν′ and ρ. Decreasing FA and increasing D⊥ with f3 also correspond with intuition. The negative correlations of a′ and positive correlation of ν′ and ρ with D|| most likely arise from orientational heterogeneity in some voxels. Nonaligned fibres cause some restriction and hindrance to water mobility in the average parallel direction reducing apparent parallel diffusivity. Meanwhile, the fibre diameter appears larger because some parts of some fibres are oblique to the single direction the model assumes. Similarly, the positive correlation of ν′ and ρ with D|| may arise because homogeneously oriented fibres can pack to higher densities. The strong negative correlation of a′ with ρ is unsurprising as ρ is proportional to a′− 2 by definition. Monkey data Fig. 12 shows maps of a′ and ρ over the midsagittal slice of each acquisition from each monkey brain. Fig. 9 (right) compares the regional trend of mean a′ with estimates of α from histology data in Lamantia and Rakic (1990). The low–high–low a′ and high–low–high ρ trends are stronger and more consistent in both parameter maps from both monkeys than in the human results. Fig. 9 confirms better reproducibility than for the human data, which we expect from the simulations in Figs. 4 and 5 and because postmortem imaging avoids motion and misalignment. Overestimation of α is still apparent, but less dramatic than in the human results, which agrees with predictions from simulation experiments in Fig. 4. The trend in a′ follows that of the histology

data points in the genu and midbody. However, histology data in Lamantia and Rakic (1990) show large axons in the splenium, whereas we observe low a′ more in common with the human data. The histology data (Lamantia and Rakic, 1990) come from 8 adult rhesus monkeys (Macaca mulatta), so age and species differences with our juvenile vervets, as well as the sparsity of the histological sampling, may account for the discrepancy. As in the human results, f3′ is close to zero in all but a few isolated voxels. Apart from a few isolated voxels, f4′ lies consistently in the range 0.2 to 0.3. Fig. 13 compares fitted models with signals in example voxels both with and without the stationary component. The benefit of nonzero f4 is clear and allows the model to capture the parallel signal much more closely. Correlation plots (not shown) similar to those in Fig. 11 show relationships among the parameters similar to those in the human results. They show stronger positive correlation between D⊥ and a′ (r = 0.31) and stronger negative correlation between FA and a′ (r = −0.26). Both observations reflect greater sensitivity to the axon diameter distribution using the monkey protocol. The correlation of a′ with D|| is weaker (r = 0.09) than for the human data, which may reflect less orientational heterogeneity because of smaller voxel sizes relative to brain volume, but nonzero f4 may also have a significant influence. Discussion In summary, this paper proposes and tests orientationally invariant indices of axon diameter, a′, and density, ρ, from diffusion MRI.

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

1385

Fig. 11. Scatter plots of the model parameter estimates against each other and against DTI indices. The plots contain one point for each processed voxel in scan 1 of human subject 1. The annotation in each subfigure shows the correlation coefficient r between the plotted quantities.

Orientational invariance is important, because it enables combination of axon diameter and density information with tractography to study fibre composition along pathways with different and varying orientation. Previous techniques only work in regions with known and constant orientation, such as the spinal cord or the midsagittal corpus callosum. The discussion here concentrates on similar regions because histology data are available for comparison. Results in those structures demonstrate orientational invariance, however, because the protocols are not tuned for their specific orientations. The method provides indices in many other regions (Fig. 8), such as more lateral areas of the corpus callosum where orientation varies considerably. Simulation experiments verify that the indices obtained from in vivo and postmortem imaging protocols are sensitive to known differences in axon diameter distribution and density. Maps of the indices in live human brains show high–low–high ρ and low–high– low a′ anterior–posterior variation in the midsagittal corpus callosum, which is consistent with known trends in axon diameter and density from histology (Aboitiz et al., 1992a). They also reflect the low density and high diameter in the long-range connections in the CST known from histology (Graf von Keyserlink and Schramm, 1984). Postmortem monkey results broadly support the in vivo human results. They show similar trends with greater contrast, less intrasubject variability and index values closer to expectations. Simulations show that parameter estimates are less stable in vivo even without motion and misalignment. Higher gradient strength and lower diffusivity in

the postmortem experiment shifts the window of sensitivity to smaller axon diameters with greater overlap of the realistic range, which provides the extra contrast. Limitations An important limitation of the proposed technique is that the axon diameter index is only a single summary statistic of the axon diameter distribution. Previous works (Assaf et al., 2008; Barazany et al., 2009) estimate two parameters, which still provide only a coarse approximation, but nevertheless a more complete model of the whole distribution. Neuroanatomical studies typically focus on quite different summary statistics, such as the mean diameter or the number of fibres of, or over, a particular size. Such statistics are simpler to estimate from the two-parameter model. However, we anticipate that the contrast a′ does provide is still useful. The precise nature of the summary statistic that a′ represents remains unclear and characterization of its relationship both to the full axon diameter distribution and to the acquisition protocol requires further work. We hypothesize initially that the axon diameter index estimates α. Although α is not the true mean diameter, it is at least as useful a summary statistic as the true mean. The true mean is relatively insensitive to large axons because smaller axons are much greater in number, whereas α emphasizes the presence of large axons, which is an important feature of the

1386

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

our in vivo experiments because the threshold is higher than the realistic range of axon diameters. In fixed tissue, however, sensitivity to the largest axon diameters (10–20 μm) is difficult to obtain: lower diffusivity than in vivo means that diffusion times must be longer to observe restriction, but T2 is also lower than in vivo, which prevents long diffusion time measurements using PGSE. Various confounding effects potentially provide alternative explanations for the trends we observe in the brain data:

Fig. 12. Maps of a′ (left) and ρ (right) in the midsagittal slices of the monkey data sets. The top two rows show results from the two scans of monkey 1. The bottom two rows show the two scans from monkey 2.

distribution. However, the summary statistic we actually obtain overestimates α, particularly for human protocol. Previous studies, e.g., Lätt et al. (2007), observe that finite-pulse durations can lead to overestimation of pore sizes if the model assumes short gradient pulses. However, the model we use here does not assume short pulses, and simulations in Alexander (2008) show that the overestimation does not occur with single axon diameters. A more likely explanation is that a′ is preferentially sensitive to larger axon diameters, which it weights more strongly than α does. For any protocol, all axon diameters below some threshold cause negligible signal attenuation so are indistinguishable from one another. With the human protocol, for example, the greatest intra-axonal signal attenuation is less than 1% for diameters less than 5 μm, but around 10% for 10 μm diameters. Thus, the larger axons have stronger influence on the average signal attenuation causing a′ to weight them more heavily than α. The higher gradient strength in the monkey protocol extends sensitivity to smaller populations of axons and the 1%-threshold is around 2 μm. Barazany et al. (2009) also report overestimation of the axon diameter in comparison to histology. Although their overestimation may arise from shrinkage during histological preparation, as the authors suggest, nonuniform diameter sensitivity of their protocol may also contribute. We note further that, although a′ lacks sensitivity to small axon diameters, the intra-axonal volume fraction parameter ν does capture contributions from small axons accurately, as Fig. 5 demonstrates. Thus, a pathology of small axons does not go unnoticed, as we would expect ρ to reduce; we also expect some increase of a′, as the average intracellular signal attenuation increases. An upper threshold on the axon diameter also exists for any protocol, above which it does not provide sensitivity because diffusion times are too short to observe restriction. We do not see this effect in

1. The thickness of the corpus callosum in the midsagittal plane correlates with the expected trends in axon diameter and density, so partial volume with CSF is often higher in the midbody than the genu or splenium, which might produce the trends we observe artifactually. However, first, the sagittal slice orientation minimizes partial volume by fitting more than one voxel within the corpus callosum even at its thinnest part. Moreover, the threshold on S⋆0 removes voxels with significant partial volume. Second, simulations in Fig. 6 suggest that partial volume does not bias a′, especially for the low f3 voxels the S⋆0 threshold leaves behind. 2. Subvoxel heterogeneity of fibre orientation increases a′. The threshold on DT linearity and planarity still admits some fibre orientation heterogeneity. The curvature of fibres in the corpus callosum may account for higher a′ in the midbody. However, highresolution MRI and streamlines from tractography qualitatively suggest the curvature is slightly lower in the midbody than at the two ends. The midbody may still have greater isotropic heterogeneity of fibre orientation, i.e., splaying rather than fanning or bending, causing higher a′ than the genu or splenium. However, micrographs in, for example, Lamantia and Rakic (1990) and Barazany et al. (2009), show no obvious increase in eccentricity in axon cross section in the midbody compared to other regions of the corpus callosum. Moreover, the CST is likely to contain the most homogeneous fibre orientation and provides among the highest a′, which suggests splaying alone does not explain the trends we observe in a′. On the other hand, splaying and fanning may contribute to the change in pattern of a′ in the corpus callosum in more lateral slices in Fig. 8. 3. The model fitting procedure fixes the diffusivity parameters to biologically plausible values, which stabilizes the fitting procedure and hastens convergence. The precise setting of the diffusivity parameters does not affect the trends in the parameter maps a great deal, although it does affect the absolute values of a′ and ρ. The fixed diffusivities we use approximate the mean D|| from the lowest b-value HARDI shell in regions of coherent white matter and CSF. The fixed settings of dI also approximate the diffusivity of water at body (in vivo) and room (postmortem) temperature. Moreover, the fixed settings are very close to those that give minimum average fitting errors over the set of processed voxels, which supports the choices. However, since d∥ is fixed, trends in the true diffusivity could cause the trends we observe in a′ and ρ. Indeed, allowing the diffusivity to vary during fitting does produce trends in d′∥ in the corpus callosum similar to those in a′ and ρ and weakens (but does not remove) the trends in a′ and ρ. However, simulations (not shown) show correlation between d′∥ and a′ even when the true d∥ used in the simulation does not vary and that stable estimates of a require fixing d∥. Barazany et al. (2009) also fix d∥. It seems unlikely that the substructure of white matter varies enough among different regions to cause significant differences in d∥. In summary, we use a simple model to obtain measurements from a clearly more complex system. The agreement with independent measurements from histology is compelling. However, as with all model-based techniques, we must acknowledge that effects other than those we model might influence the trends we observe. The promising results we obtain here, however, motivate further simulation and

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

1387

Fig. 13. As Fig. 10 but for scan 1 of monkey 1. The middle row shows fits of the model with f4 = 0 and the bottom row shows fits with nonzero f4.

validation experiments to confirm their source and test the influence of the effects above. Applications Indices of axon diameter and density in live human brains offer a variety of applications. The diameter and density of axons varies with their function. Little is known about white matter fibre composition outside the major pathways and mapping these indices over the whole brain will provide new insight in neuroanatomy and neuroscience. The variation of fibre composition over the population and during development is largely unstudied because access is traditionally via painstaking postmortem histology and results have inaccuracies of their own, such as cell shrinkage. The postmortem technique also offers major advantages over histological techniques for assessing axon diameter and density because it provides estimates over extended brain regions and, importantly, does not destroy samples. Knowledge of axon diameter and density informs models of connectivity, since larger axons transfer information more quickly (Ritchie, 1982). In combination with tractography, indices of axon diameter and density potentially provide measures of connectivity that express the information carrying potential of each connection. Such measures of structural connectivity are valuable in neuroscience studies into functional connectivity using, for example, dynamic causal models (Stephan et al., 2009). Methods to recover structural connectivity from diffusion MRI tractography may also exploit

knowledge of the axon diameter and density in different fibres to resolve ambiguities that can occur at fibre crossings and other complex fibre configurations. The indices may also prove useful markers for specific diseases. In particular, multiple sclerosis attacks small axons preferentially (De Luca et al., 2004) and motor–neurone disease preferentially attacks large axons (Cluskey and Ramsden, 2001). Highley et al. (1999) show significant differences in fibre composition in the corpus callosum, assessed by postmortem histology, between healthy subjects and certain groups of patients with schizophrenia. Other diseases and disorders may also have characteristic microstructural abnormalities. Further work The results prompt many areas of further work at each stage of the methodological pipeline. Modelling The white matter model is appealing for its simplicity and we purposefully include the minimum number of components required to capture the main trends in the measured data. However, minor departures of the signals from the models are visible in Figs. 10 and 13 and some benefits may come from increasing complexity to capture more effects. The model might use a better tortuosity estimator (Fiermans et al., 2008) and predict the signal better with finite δ (Codd and Callaghan, 1999), but since the model is crude anyway, we may see little

1388

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

improvement in the final result from these minor refinements. We might include a distribution of axon diameters, as in Assaf et al. (2008), to increase sensitivity to smaller axons and provide a more complete description of the diameter distribution. Preliminary experiments (not shown) with our sparse data sets, however, fail to recover consistent estimates of the parameters of a gamma distribution on a. Single-parameter distribution models, such as Poisson, may prove useful. An important difference we observe between the in vivo and postmortem data is the apparent restriction in the parallel direction, which we only observe in the fixed tissue. Here we use the simplest possible model to capture the effect in the postmortem data. Future work will examine the effect more closely and compare models that explain the effect in different ways more thoroughly. The most significant limitation of the current model is the requirement for a tight distribution of axon orientations, since white matter fibres in the brain form a variety of complex configurations. One possibility is to keep the current model, but interpret a′ instead as a feature of the diameter distribution perpendicular to the mean fibre direction. However, more useful indices may come from models that account for nontrivial distributions of fibre orientations. Various techniques estimate complex subvoxel distributions of fibre orientations; see Seunarine and Alexander (2009) for a review. These techniques currently use only simple models of diffusion in fibres, but extend naturally to include axon diameter parameters, which may provide the key to relaxing the constraint on fibre orientation distribution. The current model also assumes impermeable axon walls. The assumption is reasonable for myelinated axons over the timescales of the measurements we make (Stanisz et al., 1997). However, some white matter structures contain higher densities of unmyelinated axons than the CST or corpus callosum. Lack of myelin increases permeability and models may need to incorporate exchange for accurate parameter estimates in such tracts. Permeability and density of unmyelinated axons would be useful measurements in themselves if they prove possible to estimate. Sampling Advances in experiment design may come from improving the set of PGSE measurements or exploring different kinds of measurement. One simple improvement to the current design is to allow the sizes of the HARDI shells to differ. The current optimization divides the 360 measurements into four shells of equal size. With this constraint, the optimization often finds exact or near repeats, as Figs. 1 and 2 demonstrate. The optimized protocols suggest that only three shells are necessary to estimate the model parameters. However, the relative importance of the three unique PGSE settings is not equal, so the fixed-shell-size optimization repeats the most important to provide extra signal. In fact, preliminary experiments suggest the benefits of variable shell sizes are minor. We note that, although the model has up to ten free parameters in total, three shells are sufficient because additional sensitivity comes from orientational variation of the measurements. In a similar way, good estimates of multitensor models can come from single-shell HARDI protocols (Tuch et al., 2002; Parker and Alexander, 2003; Behrens et al., 2003). The a priori model parameter settings in the optimization affect the optimized protocol and the precision and accuracy of parameter estimates. More informed choices of a priori settings may improve the sensitivity of the axon diameter and density indices to differences in fibre composition. For a particular experimental set up, the available gradient strength, free diffusivity and T2 determine the range of diameters for which sensitivity is feasible (Alexander, 2008). The choices of a priori diameters for both protocol optimizations here were chosen to cover that range. In fact, the optimized protocol is fairly insensitive to the number and specific values of a priori diameters so long as they are all within that range. The optimization

becomes unstable when the a priori set includes diameters outside the range because the Cramer–Rao lower bound becomes much larger. The settings in the optimization of the human protocol represent only the largest axons. In fact, the 40-μm setting is unrealistically large, and its inclusion, which is mainly for consistency with earlier work (Alexander, 2008), does affect the optimized protocol. In simulation, however, it has only a minor effect on the indices. Further work can certainly refine the choice of a priori settings, although we do not expect this alone to produce dramatic improvements. The human imaging protocol is at the limit of tolerable times for live subjects and too long for most patient studies. Here, we choose the maximum imaging time purposefully to obtain the best data possible to demonstrate the new technique. One obvious saving is to reduce angular resolution, which is currently very high (N = 90). Such a saving would allow extension of the image volume to the whole brain providing indices in other major tracts such as the visual pathways. For some applications, full orientational invariance may be unnecessary. The experiment design optimization adapts naturally to exploit prior knowledge of fibre orientation. Such a procedure could design protocols for imaging, say, just the midsagittal corpus callosum with much lower acquisition demands. More significant improvements may come from replacing the standard PGSE sequence with other diffusion-sensitive sequences, such as oscillating gradient (Callaghan and Stepisnik, 1995; Does et al., 2003) or double wave-vector (Komlosh et al., 2008; Koch and Finsterbusch, 2008; Ozarslan and Basser, 2007) experiments, which may be more sensitive to microstructural parameters. Stimulated echo sequences (Merboldt et al., 1991) enable longer diffusion times, which should provide better sensitivity to larger axon diameters in the post-mortem protocol. The experiment design optimization extends naturally to protocols using other pulse sequences. For example, Clayden et al. (2009) demonstrate the technique for the dual spin-echo sequence (Reese et al., 2003), which reduces Eddy-current effects compared to PGSE in vivo and may enhance sensitivity slightly. Combination of the experiment design with these other pulse sequences should allow the a priori range of axon diameters to extend to include smaller diameters. This should provide protocols with sensitivity to wider ranges and provide more discriminative axon diameter indices. Fitting Future work could improve the current model fitting procedures. The MCMC step increases the likelihood of finding global minima and provides robustness to spurious fitting results from direct search, but is slow. The data-smoothing step also improves stability, although the trends we observe in the a′ and ρ maps do not depend on it (Alexander et al., 2009). More sophisticated smoothing may improve qualitative appreciation of those trends further. Validation We limit validation here to simulations and comparison with sparse histology data in the literature. More direct comparison with histology is possible with our fixed monkey brains, which will be revealing particularly outside the CST and midsagittal corpus callosum. However, we plan further imaging studies before sectioning these valuable samples. Acknowledgments The EPSRC supports DCA and MGH with grants EP/E056938/1 and EP/E007748. MP is supported by the National Science and Engineering Research Council of Canada (NSERC). TD is supported by the Lundbeck Foundation. The future and emerging technologies (FET) program of the EU FP7 framework fund the CONNECT consortium www.brain-connect.eu, which also supports this work. The authors thank Yaniv Assaf for seed code and Lee Wright for artwork.

D.C. Alexander et al. / NeuroImage 52 (2010) 1374–1389

References Aboitiz, F., Scheibel, A.B., Fisher, R.S., Zaidel, E., 1992a. Fiber composition of the human corpus callosum. Brain Res. 598, 143–153. Aboitiz, F., Scheibel, A.B., Fisher, R.S., Zaidel, E., 1992b. Individual differences in brain asymmetries and fiber composition in the human corpus callosum. Brain Res. 598, 154–161. Aboitiz, F., Rodriguez, E., Olivares, R., Zaidel, E., 1996. Age-related changes in fiber composition of the human corpus callosum. Neuroreport 7, 1761–1764. Alexander, D.C., 2008. A general framework for experiment design in diffusion MRI and its application in measuring direct tissue-microstructure features. Magn. Reson. Med. 60, 439–448. Alexander, D.C., Hubbard, P.L., Hall, M.G., Moore, E.A., Ptito, M., Parker, G.J.M., Dyrby, T.B., 2009. Orientationally invariant axon-size and density weighted MRI. Proc. ISMRM, Honolulu. ISMRM, p. 357. Assaf, Y., Basser, P.J., 2005. Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. NeuroImage 30, 1100–1111. Assaf, Y., Freidlin, R.Z., Rohde, G.K., Basser, P.J., 2004. New modeling and experimental framework to characterize hindered and restricted water diffusion in brain white matter. Magn. Reson. Med. 52, 965–978. Assaf, Y., Blumenfeld-Katzir, T., Yovel, Y., Basser, P.J., 2008. AxCaliber: a method for measuring axon diameter distribution from diffusion MRI. Magn. Reson. Med. 59, 1347–1354. Barazany, D., Basser, P.J., Assaf, Y., 2009. In-vivo measurement of the axon diameter distribution in the corpus callosum of a rat brain. Brain 132, 1210–1220. Basser, P.J., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative diffusion tensor MRI. J. Magn. Reson., Ser. B 111, 209–219. Basser, P.J., Matiello, J., Le Bihan, D., 1994. MR diffusion tensor spectroscopy and imaging. Biophys. J. 66, 259–267. Behrens, T.E.J., Woolrich, M.W., Jenkinson, M., Johansen-Berg, H., Nunes, R.G., Clare, S., Matthews, P.M., Brady, J.M., Smith, S.M., 2003. Characterization and propagation of uncertainty in diffusion-weighted MR imaging. Magn. Reson. Med. 50, 1077–1088. Callaghan, P.T., Stepisnik, J., 1995. Frequency domain analysis of spin motion using modulated-gradient NMR. J. Magn. Reson. A 117, 118–122. Clayden, J.D., Nagy, Z., Hall, M.G., Clark, C.A., Alexander, D.C., 2009. Active imaging with dual spin-echo diffusion MRI. In Proc. IPMI, Williamsburg. Springer, pp. 264–275, LNCS 5636. Cluskey, S., Ramsden, D.B., 2001. Mechanisms of neurodegeneration in amyotrophic lateral sclerosis. J. Clin. Pathol. Mol. Pathol. 54, 386–392. Codd, S.L., Callaghan, P.T., 1999. Spin echo analysis of restricted diffusion under generalized gradient waveforms: planar, cylindrical and spherical pores with wall relaxivity. J. Magn. Reson. 137, 358–372. Cook, P.A., Bai, Y., Nedjati-Gilani, S., Seunarine, K.K., Hall, M.G., Parker, G.J.M., Alexander, D.C., 2006. Camino: open-source diffusion–MRI reconstruction and processing. Proc. 14th Annual Meeting of the ISMRM, Berlin. De Luca, G.C., Ebers, G.C., Esiri, M.M., 2004. Axonal loss in multiple sclerosis: a pathological survey of the corticospinal and sensory tracts. Brain 127, 1009–1018. Does, M.D., Parsons, E.C., Gore, J.C., 2003. Oscillating gradient measurements of water diffusion in normal and globally ischemic rat brain. Magn. Reson. Med. 49, 206–215. Dyrby, T.B., Baaré, W.F.C., Alexander, D.C., Jelsing, J., Garde, E., Søgaard, L.V., 2010. An exvivo imaging pipeline for producing high quality and high resolution diffusion weighted imaging datasets. Hum. Brain Mapp. doi:10.1002/hbm.21043. Fiermans, E., De Deene, Y., Delputte, S., Ozdemir, M.S., D'Asseler, Y., Vlassenbroeck, J., Deblaere, K., Acten, E., Lemahieu, I., 2008. Simulation and experimental verification of the diffusion in an anisotropic fiber phantom. J. Magn. Reson. 190, 189–199. Graf von Keyserlink, D., Schramm, U., 1984. Diameter of axons and thickness of myelin sheaths of the pyramidal tract fibres in the adult human medullary pyramid. Anat. Anz. 157, 97–111. Hall, M.G., Alexander, D.C., 2009. Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI. IEEE Trans. Med. Imaging 28, 1354–1364.

1389

Highley, J.R., Esiri, M.M., McDonald, B., Cortina-Borja, M., Herron, B.M., Crow, T.J., 1999. The size and fibre composition of the corpus callosum with respect to gender and schizophrenia: a post-mortem study. Brain 122, 99–110. Jenkinson, M., Bannister, P., Brady, M., Smith, S., 2002. Improved optimization for the robust and accurate linear registration and motion correction of brain images. NeuroImage 17, 825–841. Jones, D.K., Horsfield, M.A., Simmons, A., 1999. Optimal strategies for measuring diffusion in anisotropic systems by magnetic resonance imaging. Magn. Reson. Med. 42, 515–525. Koch, M.A., Finsterbusch, J., 2008. Compartment size estimation with double wave vector diffusion-weighted imaging. Magn. Reson. Med. 60, 90–101. Komlosh, M.E., Lizak, M.J., Horkay, F., Freidlin, R.Z., Basser, P.J., 2008. Observation of microscopic diffusion anisotropy in the spinal cord using double-pulsed gradient spin echo MRI. Magn. Reson. Med. 59, 803–809. Lamantia, A.-S., Rakic, P., 1990. Cytological and quantitative characteristics of four cerebral commissures in the rhesus monkey. J. Comp. Neurol. 291, 520–537. Lätt, J., Nilsson, M., Malmborg, C., Rosquist, H., Wirestam, R., Ståhlberg, F., Topgaard, D., Brockstedt, S., 2007. Accuracy of q-space related parameters in MRI: simulations and phantom measurements. IEEE Trans. Med. Imaging 26, 1437–1447. Merboldt, K.D., Hanicke, W., Frahm, J., 1991. Diffusion imaging using stimulated echoes. Magn. Reson. Med. 19, 233–239. Murday, J.S., Cotts, R.M., 1968. Self-diffusion coefficient of liquid lithium. J. Chem. Phys. 48, 4938–4945. Ozarslan, E., Basser, P.J., 2007. MR diffusion/diffraction phenomenon in multi-pulsefield-gradient experiments. J. Magn. Reson. 188, 285–294. Packer, K., Rees, C., 1972. Pulsed NMR studies of restricted diffusion: I. Droplet size distributions in emulsions. J. Colloid Interface Sci. 40, 206–218. Panagiotaki, E., Fonteijn, H., Siow, B., Hall, M.G., Price, A., Lythgoe, M.F., Alexander, D.C., 2009. Two-compartment models of the diffusion MR signal in brain white matter. In Proc. MICCAI, London. Springer, pp. 329–336. LNCS 5761. Parker, G.J.M., Alexander, D.C., 2003. Probabilistic Monte Carlo based mapping of cerebral connections utilising whole-brain crossing fibre information. Proc. IPMI, Ambleside. Springer, pp. 684–695. LNCS 2732. Ptito, M., 2003. Functions of the corpus callosum as derived from split chiasm studies in cats. In: Zaidel, E., Iacobini, M. (Eds.), The Parallel Brain. MIT Press, pp. 139–153. Reese, T.G., Heid, O., Weisskoff, R.M., Wedeen, V.J., 2003. Reduction of eddy-currentinduced distortion using a twice-refocussed spin echo. Magn. Reson. Med. 49, 177–182. Ritchie, J.M., 1982. On the relation between fibre diameter and conduction velocity in myelinated nerve fibres. Proc. R. Soc. Lond. B 217, 29–35. Seunarine, K.K., Alexander, D.C., 2009. Multiple fibers: beyond the diffusion tensor. In: Johansen-Berg, H., Behrens, T.E.J. (Eds.), Diffusion MRI: From Quantitative Measurement to In Vivo Neuroanatomy. Academic Press, pp. 56–74. Stanisz, G.J., Szafer, A., Wright, G.A., Henkelman, M., 1997. An analytical model of restricted diffusion in bovine optic nerve. Magn. Reson. Med. 37, 103–111. Stejskal, E.O., Tanner, T.E., 1965. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42, 288–292. Stephan, K.E., Tittgemeyer, M., Knösche, T.R., Moran, R.J., Friston, K.J., 2009. Tractography-based priors for dynamic causal models. NeuroImage 47, 1628–1638. Szafer, A., Zhong, J., Gore, J.C., 1995. Theoretical model for water diffusion in tissues. Magn. Reson. Med. 33, 697–712. Tuch, D.S., Reese, T.G., Wiegell, M.R., Makris, N., Belliveau, J.W., Wedeen, V.J., 2002. High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity. Magn. Reson. Med. 48, 577–582. Van Gelderen, P., DesPres, D., van Zijl, P.C.M., Moonen, C.T.W., 1994. Evaluation of restricted diffusion in cylinders. Phosphocreatine in rabbit leg muscle. J. Magn. Reson., Ser. B 103, 255–260. Waxman, S.G., Kocsis, J.D., Stys, P.K., 1995. The axon. Structure, Function and Pathophysiology. Oxford University Press, New York. Westin, C.-F., Maier, S.E., Khidhir, B., Everett, P., Jolesz, F.A., Kikinis, R., 1999. Image processing for diffusion tensor magnetic resonance imaging. Proc. MICCAI, Cambridge. Springer, pp. 441–452.