Whole-brain susceptibility mapping at high field: A comparison of multiple- and single-orientation methods

NeuroImage 53 (2010) 515–525

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

Whole-brain susceptibility mapping at high field: A comparison of multiple- and single-orientation methods Sam Wharton, Richard Bowtell ⁎ Sir Peter Mansfield Magnetic Resonance Centre, School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK

a r t i c l e

i n f o

Article history: Received 21 April 2010 Revised 14 June 2010 Accepted 28 June 2010 Available online 6 July 2010 Keywords: Susceptibility mapping High field Whole brain

a b s t r a c t Optimisation and comparison of the performance of three different methods for calculating threedimensional susceptibility maps of the whole brain from gradient-echo (phase and modulus) image data acquired at 7 T is described. The methods studied are a multiple-orientation method in which image data acquired with the head at several different angles to the main field are combined and two methods which use data acquired at a single orientation: the first of these is based on exclusion of some k-space data from the calculation (through thresholding of the dipolar field kernel), while the second incorporates a regularisation method that is based on using information from the modulus images. The methods were initially optimised via analysis of data from a phantom containing different compartments of known susceptibility. As part of this work, a novel high-pass filtering methodology was introduced to remove background fields from field maps based on phase data. The optimised methods were successfully applied to high-resolution (0.7 mm isotropic) whole-brain modulus and phase data acquired in vivo from five healthy male subjects, 25–30 years of age. The multiple-orientation method yielded high quality susceptibility maps, out-performing the single-orientation methods. Venous blood vessels as well as the substantia nigra and globus pallidus brain regions showed particularly high positive susceptibility offsets relative to surrounding tissue, consistent with high deoxyhemoglobin and non-heme iron content, respectively. To compare the performance of the different methods, regions of interest were drawn in deep grey matter structures and in cortical grey and white matter. The threshold-based approach was fast and simple to use, but underestimated susceptibility differences and showed significant artefacts due to noise amplification in difficult regions of k-space. The regularised single-orientation method yielded contrast dependent on the choice of spatial priors, but demonstrated the potential to yield susceptibility maps of a similar quality to those calculated using data acquired at multiple orientations to the field. © 2010 Elsevier Inc. All rights reserved.

Introduction Phase images of the human brain acquired at high field using gradient-echo-based sequences show excellent contrast (Duyn et al., 2007). This is due to the small differences in magnetic susceptibility between tissues that perturb the local magnetic field, thus producing spatial variation of the NMR frequency (Haacke et al., 2005; Duyn et al., 2007; Deistung et al., 2008; Koopmans et al., 2008; Haacke et al., 2009; Schäfer et al., 2009). Differences in the relative abundance of myelin, non-heme iron, and deoxyhemoglobin are thought to be major contributors to these susceptibility offsets in the brain (Annese et al., 2004; Haacke et al., 2005; Duyn et al., 2007). A common usage of phase data is in susceptibility weighted imaging (SWI), where processed phase data are combined with standard modulus images to improve the visualisation of veins. This approach exploits the large phase shifts

⁎ Corresponding author. Fax: + 44 115 9515166. E-mail address: [email protected] (R. Bowtell). 1053-8119/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2010.06.070

associated with the high concentrations of paramagnetic deoxyhemoglobin in venous blood (Haacke et al., 2004; Deistung et al., 2008; Koopmans et al., 2008). Extraction of accurate anatomical information from phase data is, however, made difficult by the non-local relationship between the field perturbation that gives rise to the phase variation and the associated susceptibility distribution (Marques and Bowtell, 2005; Schäfer et al., 2009). Recently, the development of a fast, Fourier-based method for rapidly simulating the field shifts produced by arbitrary susceptibility distributions (Salomir et al., 2003; Marques and Bowtell, 2005) has led the way for the development of a new technique called susceptibility mapping (Kressler et al., 2010; de Rochefort et al., 2008; Liu et al., 2009; Shmueli et al., 2009; de Rochefort et al., 2010; Wharton et al., 2010). Susceptibility mapping involves inversion of phase data by one of a number of possible Fourier-based methods so as to yield quantitative susceptibility maps. Susceptibility mapping has many possible applications including accurate in vivo measurement of the concentration of contrast agents, as well as investigation of the relationship between iron content and the progression of neurodegenerative diseases such as Parkinson's disease and multiple sclerosis.

516

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

The main challenge in creating susceptibility maps from phase data is the ill-posed nature of the process of inversion. The reason for using Fourier-based methods is that the relationship between the field perturbation and the magnetic susceptibility distribution becomes simple and local in the Fourier domain. However, when inverting this relationship to obtain the susceptibility in terms of the field variation a problem arises, since the denominator of the resulting expression tends to zero over two conical surfaces in k-space. These surfaces lie at the magic angle to the direction of the main magnetic field and to its reflection. Without careful conditioning the resulting division-by-zero gives rise to artefacts in the calculated susceptibility distribution. Several conditioning strategies have been proposed for ameliorating this problem. These include combining data obtained with the susceptibility distribution orientated at multiple angles with respect to the main magnetic field (Liu et al., 2009; Wharton et al., 2010), regularisation of the inversion problem via the incorporation of prior knowledge of the form of the susceptibility distribution (Kressler et al., 2010; de Rochefort et al., 2010), and utilisation of a simple k-space threshold (Shmueli et al., 2009; Wharton et al., 2010). The first method requires multiple data acquisitions and will be referred to as the multiple-orientation (MO) method. The MO method has the obvious disadvantage of increased scanning time and is also dependent on the ability of subjects to rotate their heads between scans. Despite these drawbacks the MO method deals with noise robustly and produces consistently high quality susceptibility maps. The second and third methods require only a single acquisition and will be referred to as the regularised single-orientation (RSO) method, and the threshold-based single-orientation (TSO) method, respectively. The RSO and TSO methods, although requiring less scanning time and no head rotation, are more sensitive to noise than the MO method. The aims of the study described here were firstly, to optimise the performance of each susceptibility mapping method for the calculation of high quality quantitative susceptibility maps from whole-brain phase and modulus data acquired at 7 T, and secondly, to investigate by comparison the strengths and weaknesses of the different methods. The methods and their optimisations were validated by analysis of data from a phantom containing different compartments of known susceptibility, and then applied to brain data acquired in vivo at 7 T. Theory Susceptibility and field shifts in the Fourier domain In the Fourier domain, the field perturbation, ΔBz(r), produced by a susceptibility distribution, χ(r), that is exposed to a magnetic field, B0^z, is given by a simple, local expression Δ B˜ z ðkÞ = B0 χ˜ ðkÞ

1 k2 − z2 3 jkj

! = B0 χ˜ ðkÞ⋅CðkÞ;

ð1Þ

Δ B˜ z ðkÞ 1 : ⋅ B0 CðkÞ

Conditioning using a multiple-orientation (MO) method One method that can be used to stabilize this inversion problem is to acquire and combine data with χ˜ ðkÞ rotated to several different orientations with respect to the main magnetic field. This approach has been most fully explored by Liu et al. (2009) who introduced the COSMOS (calculation of susceptibility through multiple-orientation sampling) method and applied it to simple phantoms and excised tissue. More recently, we showed that rotating χ˜ ðkÞ about two axes, as opposed to the single axis rotations used in COSMOS, may be advantageous for susceptibility mapping when the angle of rotation is restricted to small values, as is the case for head imaging (Wharton et al., 2010). These multiple-orientation methods take advantage of the rotational dependency of C(k). For a rotation of the object about the x-axis by angle, θ, followed by a rotation about the y-axis by ϕ, the Fourier transform of the convolution kernel becomes 0

 2 1 k cosθ cos ϕ−k sin θ cos ϕ + k sin ϕ z y x B1 C C ðk; θ; ϕÞ = @ − A: 3 jkj2

ð2Þ

Although this deconvolution appears to be straightforward, problems arise because the denominator in Eq. (2) goes to zero when the angle between the k-vector and the field direction is equal to 54.7° or 125.3°. This means that χ˜ ðkÞ tends to infinity on two conical surfaces in k-space and there is significant amplification of any

ð3Þ

This equation indicates that the conical surfaces where C(k, θ, ϕ) → 0 rotate with the magnetic field. Accordingly, if several rotations are carried out and field maps recorded at each orientation it is possible to combine the data sets such that the whole of k-space can be sampled without including elements with C(k, θ, ϕ) values that are close to zero. An estimate of χ˜ ðkÞ which best satisfies 2

3

2

Δ B˜ z ðkÞ 1

3

C1 6 7 6 ˜2 7 6 C2 7 Δ Bz ðkÞ 7 6 7⋅ χ˜ = 6 6 7 4 ::: 5 7 6 4 ::: 5 CN N Δ B˜ z ðkÞ

ð4Þ

(where CN = C(k, θN, ϕN) is the Fourier transform of the convolution N kernel and Δ B˜ z ðkÞ is the Fourier transform of the measured field map at the Nth orientation) can then be identified using weighted, leastsquares fitting methods. Liu et al. (2009) found that an iterative real space method, driven by the LSQR algorithm for solving sparse linear equations (Christopher and Michael, 1982), gave a good approximation of χ. The goal of such an algorithm is to find χ that satisfies minχ



‖ O1 ‖ 22 + ‖ O2 ‖ 22 ::: + ‖ ON ‖ 22



ð5Þ

where ON = WN ðFT

where the overstrike “~” denotes a three-dimensional Fourier transform (3DFT), and C(k) represents the Fourier transform of the dipole convolution kernel that links susceptibility and field. Rearranging this expression yields χ˜ ðkÞ as a function of Δ B˜ z ðkÞ χ˜ ðkÞ =

noise or errors arising from the measured Δ B˜ z ðkÞ in the region adjacent to these surfaces. For this reason, the deconvolution can be classed as an ill-posed problem requiring careful conditioning.

−1

ðCN χ˜ Þ−ΔBN Þ

ð6Þ

and WN is the chosen weighting matrix for the Nth orientation. In each iteration, N field maps are simulated using Eq. (1) and compared to the measured field maps to form a residual, which is then weighted by WN. A sensible choice in setting the values of WN is to use the modulus signal amplitude, as the noise in the field measurements is proportional to the inverse of the signal to noise ratio in the modulus data. Conditioning using a regularised single-orientation (RSO) method Another method for stabilising the inversion problem is to incorporate prior knowledge of the susceptibility distribution to form a regularisation problem requiring only data acquired at a single

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

orientation. This approach has been most fully explored by de Rochefort et al. (2010), who introduced the use of edge information extracted from modulus data as effective priors. The iterative method which these authors used to find the best fitting solution is based on a least-squares, conjugate gradient algorithm (Hestenes and Stiefel, 1952). In this approach, the expression that χ must satisfy can be written as minχ



‖ O1 ‖ 22

2

+α

‖ Mouter χ‖ 22 + β2 ‖ Wg Gχ‖ 22



ð7Þ

where α and β are regularisation parameters, and Mouter is a mask which defines the region outside of the brain/phantom where the susceptibility is assumed to be zero. The final, regularisation term can be expanded as

‖ Wg Gχ‖ 22 = ‖ Wgx Gx χ‖ 22 + ‖ Wgy Gy χ‖ 22 + ‖ Wgz Gz χ‖ 22

ð8Þ

where Wg is made up of a set of weighting matrices linked to each Cartesian direction, and Gi = ∂i∂ is the gradient operator which is used to extract edge information from χ along the ith dimension. de Rochefort et al. (2010) showed that the inverse of the gradient of the modulus data is a good choice for Wg as it is rich in structural detail and its use allows the inversion to be carried out without the acquisition of additional scans. This approach was also used here. In each iteration, a field map was simulated, and the gradient along each dimension of χ measured. The effect of the gradient weighting terms is to allow the χ distribution to vary in regions where the modulus signal varies with position, but to force it to be smooth in regions where the modulus data shows little spatial variation. Conditioning using a threshold-based single-orientation (TSO) method Eq. (2) suggests that dividing Δ B˜ z ðkÞ by C(k) directly will yield χ˜ ðkÞ. However, as explained above any noise in the measured field map will be strongly amplified in the regions of k-space where C(k) → 0. A relatively simple method for conditioning the inversion process is to threshold the data so as to reduce or remove completely the contribution of these badly behaved regions. In particular, Shmueli et al. (2009) showed susceptibility maps with interpretable detail could be created by truncating 1/C(k) to an appropriately chosen maximum value. Likewise, Wharton et al. (2010) showed similar maps could be created by setting k-space data to zero in regions where 1/C(k) is larger than a chosen threshold value. With such threshold-based methods, a balance must be reached between reducing noise by restricting the contribution of badly behaved kspace regions, and preserving contrast and spatial resolution by including in as much k-space information as possible. Materials and methods Initial experiments were carried out on a spherical agar phantom of 10 cm diameter which contained a random arrangement of small cylindrical sections of agar that had been doped with an iron-oxidebased contrast agent. In separate experiments, the susceptibility of this agent was measured to be 0.15 ± 0.01 ppm/mM at 7 T. Small cylindrical sections containing two different concentrations of the contrast agent (0.5 and 1 mM) were used, yielding a test sample showing a structured variation of magnetic susceptibility that spanned similar values to those found in the human brain. Subsequent in vivo studies, involved scanning five, healthy, male subjects aged between 25 and 30 years. The study was approved by the local ethics committee. Data processing was carried out using a 64-bit Linux system with a 2 GHz Dual Core AMD processor and 8 GB of RAM. Data fitting and simulations were carried out in Matlab (The Mathworks Inc, Massachusetts).

517

Data acquisition Images were acquired on a 7 T scanner (Philips Medical Systems, Best, the Netherlands) using a 3D spoiled GE FLASH sequence with 0.7 mm isotropic resolution, 200 × 200 × 100 mm3 FOV with axial slab orientation,, TE/TR = 15/23 ms, EPI factor = 3, SENSE factor = 2, a flip angle of 12° and a scan time of 210 s per acquisition. The phantom was ˆ (i) rotated −10° about the x-axis imaged at four orientations to B0 z: from an arbitrary starting position, (ii) rotated + 10o about the x-axis, (iii) rotated +10° about the y-axis, and finally (iv) rotated − 10° about the y-axis. Similarly, for the in vivo experiments each subject ˆ (i) tipped was scanned with their head at four orientations to B0 z: back by ~ 10°, corresponding to a rotation about the right-left axis (x-axis) from a neutral head position, (ii) tipped forward by ~ 10°, (iii) tipped to the side by ~10o, corresponding to a rotation about the anterior-posterior axis (y-axis), and (iv) tipped to the opposite side by ~10°. Approximate angles are quoted, as the subjects rotated their heads by varying angles ranging from 5° to 15° in each direction. High pass filtering Before susceptibility mapping can be carried out, high-pass spatial filtering of the field maps must be performed to remove unwanted slowly varying, background fields due to air/tissue interfaces, such as the sinuses, and incorrectly set shimming fields. Fig. 1 shows the steps involved in filtering the phase data. Before filtering, the phase images were unwrapped using the Prelude tool in FSL (Jenkinson, 2003), and scaled by γB0TE × 10−6 to yield field maps in ppm, as shown in Fig. 1b. Next, a dipole fitting technique similar to that used by de Rochefort et al. (2010) was employed, with the addition of a low-order polynomial fit. This approach identifies a susceptibility distribution outside the region of interest which generates a field that best fits the measured field variation inside the region of interest. The filtering was carried out by solving a minimization problem analogous to that described for the RSO method above (Eq. (7)):  minχ;p

‖ W1 ðFT −1 ðC1 χ˜ Þ + Vp−ΔB1 Þ‖ 2 2

+α

2

‖ Mχ‖ 2 2

 ð9Þ

using an iterative conjugate gradient based algorithm (Hestenes and Stiefel, 1952). Here M is a mask that is only non-zero inside the region of interest (i.e. telling the algorithm where susceptibility, for the purposes of modelling background fields, is not allowed to exist); V is a matrix describing low-order polynomials in x, y, and z; p is a vector describing the polynomial coefficients associated with these orders; ΔB1 is the measured field map; and χ is the susceptibility distribution outside the region of interest. The polynomial was incorporated to account for the slowly varying fields produced by B0-inhomogeneity. The BET tool in FSL (Stephen, 2002) was applied to the modulus data (Fig. 1a) to create a mask of the brain, M (Fig. 1c). To decrease computation time the field map was masked by M and resampled to yield a lower resolution field map, ΔB1, with effective voxel size = 2.1 mm isotropic. The weighting matrix W1 was set equal to the modulus data, after being masked and resampled to match the resolution of ΔB1. M was also resampled to match ΔB1 and W1. To allow room for susceptibility offsets to be populated around the field map, ΔB1, W1, and M were zero padded to form matrices twice as large in each dimension as the originals. The matrix V was set up to describe polynomials up to 2nd order over the region of space defined by M, so that p contained 10 coefficients. Before the fitting could take place an appropriate value for the regularisation parameter, α, had to be found. To help simplify this choice for different acquisitions/subjects W1 was normalised by dividing each element of W1 by the mean value of W1 over the region of space defined by M. An acceptable value of α = 1000 was established by trial and error. This was the lowest value at which the masking term restricted the susceptibility fitted in the region

518

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

Fig. 1. The steps used in high-pass filtering of field maps. (a) Axial modulus image. (b) Unwrapped phase data. (c) Mask of brain region, M. (d) Fitted susceptibility, χfit. (e) Simulated dipole field, FT −1 ðC1 χ˜ fit Þ. (f) Fitted 2nd order polynomial, Vpfit. (g) Interpolated and cropped, fitted field map based on (e) + (f), ΔBfit. (h) High-pass filtered field map calculated from (b) to (g).

defined by M, to numerical values close to zero (b0.0001 ppm). It was found that after 2000 iterations of the fitting algorithm, taking about 5 h, acceptable values for χ and p were found to produce a fitted field, ΔBfit, given by ΔBfit = FT

−1

identical to that used in the RSO method (Eq. (7)), which was found to improve the stability of the fitting algorithm:  minχ

ðC1 χ˜ fit Þ + Vpfit

‖ O1 ‖ 2 + ‖ O2 ‖ 2 + ‖ O3 ‖ 2 + ‖ O4 ‖ 2 + α2 ‖ Mouter χ‖ 2 2

2

2

2

2

 :

ð11Þ

ð10Þ

where χfit is the fitted susceptibility in the outer region shown in Fig. 1d, FT −1 ðC1 χ˜ fit Þ is the simulated dipole field shown in Fig. 1e, pfit describes the fitted polynomial coefficients, and Vpfit is the fitted polynomial shown in Fig. 1f. ΔBfit was then interpolated back to 0.7 mm isotropic resolution and cropped back to a matrix size identical to the original ΔB1, see Fig. 1g. The fitted field map, ΔBfit, was subtracted from ΔB1to yield a high-pass filtered field map, ΔBfit as shown in Fig. 1h. Measuring orientation Before susceptibility maps can be calculated it is necessary to establish the exact orientation of each measured data set relative to B0. To do this, the modulus data sets were first co-registered using the FLIRT tool in FSL (Jenkinson and Smith, 2001). This process yields rotation matrices describing the transformation of each data set to a common space. Combining this information with the angle of each imaging stack to B0, recorded directly from the scanner, it was possible to characterise the orientations accurately in terms of angles θ and ϕ that can then be used to describe C(k, θ, ϕ) via Eq. (3). Susceptibility calculations For the phantom and in vivo data the MO method was applied to data acquired at four orientations. The minimisation described in Eq. (5) was slightly modified by the addition of a mask-related term,

The mask was calculated from modulus data acquired at one of the orientations and following the process used for the high-pass filtering, the weighting matrices were set equal to the modulus data for each orientation and mean normalised. In a similar method to that used in defining the high-pass phase filter an acceptable value for the regularisation parameter was established at α = 10. For the RSO method similar weighting and masking matrices were extracted from the modulus data and the mask regularisation parameter, α, was set to a value of 15. The next stage was to establish an acceptable value for the regularisation parameter, β (see Eq. (7)), and also for an extra noise-related parameter. The latter is used to clip the maximum value of Wg, as defined in Eq. (8), to a chosen threshold level (de Rochefort et al., 2010). This threshold defines what is considered to be a significant gradient in the modulus image. Preliminary experiments were carried out on the phantom data to determine the best values for the noise threshold and regularisation parameters. These experiments involved varying each parameter and inspecting the associated susceptibility map within carefully chosen ROIs. In contrast to the MO and RSO methods the TSO method requires less preparation, the only variable parameter being the k-space threshold value. Preliminary experiments, similar to those used for the RSO parameter selection, were carried out on the phantom data to determine the best value for this threshold. The first data set acquired (with the head/phantom tipped back by a rotation about the rightleft/x-axis of about 10°) was used for the susceptibility calculations in the single-orientation methods.

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

ROI selection It is important to note that the susceptibility mapping described here yields relative rather than absolute values of the susceptibility due to the arbitrary setting of the resonant frequency and the fact that χ˜ ðkÞ is not defined at k = 0. It is therefore necessary to work in terms of susceptibility differences between regions or compartments, which we write as Δχ. For the phantom, ROIs were drawn inside the doped cylinders, and also in adjacent regions of the surrounding agar allowing the difference in susceptibility between the doped and undoped agar compartments to be calculated. Several brain areas divided across lower (red nucleus [RN] and substantia nigra [SN]), middle (internal capsule [IC], globus pallidus [GP], putamen [PU], thalamus [TH], and caudate nucleus [CN]) and upper (grey matter [GM] regions in the cortex) sections were targeted in the in vivo data. The upper GM regions in the cortex were divided into those located in the frontal lobe [GMFL], anterior parietal lobe [GMAPL], and posterior parietal lobe [GMPPL]. ROI were drawn on transverse slices of susceptibility maps created using the three different methods and then combined so that only voxels common to all ROI were included in the final composite ROI. Fig. 2 shows a representative set of ROI drawn on axial susceptibility maps (created using the MO method) from each of the three sections. Although this figure only shows ROIs on the left side of the brain for lower and middle sections, the data analysis included the average of ROI drawn on 4–5 slices at each brain level on both the right and left sides. As discussed above, it is important to choose a reference susceptibility value for meaningful data analysis. Here, ROI were drawn in WM throughout the brain for each subject. The average susceptibility value measured in the WM was then subtracted from the values measured in each of the brain structures listed above to yield values of Δχ relative to WM. Results Phantom results Fig. 3 shows susceptibility values and associated errors calculated from maps produced in preliminary phantom experiments by using the RSO and TSO methods. Results are plotted for varying values of the noise threshold parameter and edge information regularisation parameter, β, for the RSO method, and for different values of the kspace threshold parameter for the TSO method. To simplify the plots, only the results obtained from the ROI in cylinders doped with a 1 mM concentration of the contrast agent are shown. As discussed above, these cylinders have an expected Δχ of 0.15 ± 0.01 ppm relative to the surrounding undoped agar. The Δχ values produced using the RSO method, for three different noise thresholds (0.5, 16, and 512) are plotted against log(β) in Fig. 3a, while Fig. 3b shows plots of the

519

associated averaged standard deviations (SD) within the ROI. These plots show that small values of β yield accurate Δχ values, but also produce a high standard deviation that indicates significant variation of the calculated susceptibility within what should be homogeneous compartments, while a high value of β yields a low SD, but also produces Δχ values that are significantly reduced in magnitude. As the noise threshold is increased, the transition between these two regimes occurs at higher values of β. The points which are circled on Fig. 3a and b show the parameter values (noise threshold = 16; β = 1) which were used in the later calculations. These values gave the best compromise between maintaining the correct average Δχ value whilst limiting its standard deviation. The TSO method showed a similar dependence on the k-space threshold parameter (Fig. 3c and d). An approximately linear reduction in Δχ was measured on increasing the threshold value, while the SD values decreased quickly at first and then levelled off for higher threshold settings. A k-space threshold value of 0.06 (assuming C(k) is scaled to have a range between 1/3 and − 2/3, see Eq. (1)) was chosen as an acceptable compromise, corresponding to the points circled on Fig. 3c and d. To test the validity of using the phantom data to optimise the TSO method for the brain studies, similar analysis was carried out on the in vivo data using the susceptibility values in all ROI. The results showed a similar dependence of susceptibility on the threshold parameter, supporting the use of a threshold of 0.06 for the in vivo data analysis. Fig. 4 shows representative coronal images of the phantom along with susceptibility maps calculated using the optimal parameter settings for the three different methods. Fig. 4a shows the meannormalised modulus data acquired at one orientation, which was used for weighting in the MO and RSO methods. Fig. 4b shows a corresponding field map which was calculated from spatially filtered phase data. Fig. 4c shows the susceptibility map calculated by applying the MO method to data acquired with the phantom at four different orientations to B0. No parameter selection was necessary for this method as the data from each orientation were weighted equally in the minimization problem given in Eq. (5). Calculation of the susceptibility map required 400 iterations, taking approximately 3 h of computation time. Fig. 4d shows the susceptibility map calculated using the TSO method with a k-space threshold value of 0.06. The total computation time for this method was about 2 min. The susceptibility map calculated using the RSO method with a noise threshold of 16 and β = 1 is shown in Fig. 4e. This calculation required around 200 iterations, taking ~ 2 h of computation time, to converge on a solution. The combined magnitude of the weighting matrices (W2gx + W2gy + 2 1/2 Wgz ) , which was used to impose edge information in the RSO method, is also shown in Fig. 4f. Table 1 lists the average Δχ values measured in cylinders containing 0.5 and 1 mM concentrations of contrast agent for the MO, RSO, and TSO methods.

Fig. 2. ROIs (coloured lines) used for data analysis superimposed on (a–c) transverse slices of a susceptibility map from a single subject created using the MO method. The regions chosen were: the substantia nigra (SN), red nucleus (RN), internal capsual (IC), globus pallidus (GP), putamen (PU), caudate nucleus (CN), thalamus (TH), GM in a posterior region of the parietal lobe (GMPPL), GM in a anterior section of the parietal lobe (GMAPL), and GM in the frontal lobe (GMFL).

520

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

Fig. 3. Measured susceptibility differences (a and c) and associated standard deviations (b and d) in the compartments of the agar phantom doped with 1 mM contrast agent, for the RSO (a and b) and TSO (c and d) susceptibility mapping methods, plotted against different values of relevant mapping parameters. The noise threshold and regularisation parameter, β, were varied in the RSO method, while the k-space threshold was altered in the TSO method.

Fig. 4. Representative coronal images acquired from the agar phantom and the corresponding susceptibility maps calculated using the optimised MO, RSO, and TSO methods. (a) Mean-normalised modulus data from a single orientation. (b) Filtered field map from a single orientation. (c) χ-Map calculated using the MO method. (d) χ-Map calculated using the optimised TSO method. (e) χ-Map calculated using the optimised RSO method. (f) Combined weighting matrix used to impose spatial priors based on the gradient of (a) in the RSO method (W2gx + W2gy + W2gz)1/2, with a noise threshold = 16.

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525 Table 1 Susceptibility differences measured for the two differently doped compartments of the agar phantom using the three different mapping methods. The standard deviation (SD) was calculated from the variance within each ROI.

0.5 mM 1 mM

MO

RSO

TSO

Δχ ± SD (ppm)

Δχ ± SD (ppm)

Δχ ± SD (ppm)

0.16 ± 0.01 0.076 ± 0.009

0.15 ± 0.01 0.067 ± 0.009

0.13 ± 0.03 0.07 ± 0.02

In vivo results Fig. 5 shows representative axial brain images drawn from the lower, middle and upper sections. The modulus data and field maps, upon which the susceptibility calculations are based, are shown in Fig. 5a, f and k, and b, g and l, respectively. χ-Maps calculated using the MO (Fig. 5c, h and m), RSO (Fig. 5d, i and n), and TSO (Fig. 5e, j and o) methods are also shown. For the RSO method the noise threshold and regularisation parameters were set to the optimum values derived from the phantom data (noise threshold= 16; β= 1). Similarly, for the TSO method the k-space threshold was set to the optimum value of 0.06. Table 2 lists the average Δχ values in each of the brain regions measured relative to the white matter reference region using the three different methods. ΔB values taken from the field maps are also detailed. The numbers quoted are the average over the five subjects along with the standard deviations of the means across subjects. The χ-map calculated using the MO method shows excellent structural detail, with high

521

Table 2 Measured susceptibility differences averaged over the five subjects using each mapping method, for each of the different brain regions. The standard deviation (SD) was calculated from the difference in the mean values over the five subjects. The measured field differences for each ROI are also listed for comparison.

SN RN IC GP PU CN TH GMPPL GMAPL GMFL

MO

RSO

TSO

Field map

Δχ ± SD (ppm)

Δχ ± SD (ppm)

Δχ ± SD (ppm)

ΔB ± SD (ppm)

0.17 ± 0.02 0.14 ± 0.02 − 0.01 ± 0.01 0.19 ± 0.02 0.10 ± 0.01 0.09 ± 0.02 0.045 ± 0.009 0.043 ± 0.007 0.053 ± 0.003 0.040 ± 0.006

0.16 ± 0.03 0.13 ± 0.02 − 0.014 ± 0.009 0.19 ± 0.02 0.09 ± 0.01 0.09 ± 0.01 0.05 ± 0.01 0.05 ± 0.01 0.052 ± 0.009 0.05 ± 0.01

0.11 ± 0.01 0.09 ± 0.01 − 0.03 ± 0.006 0.14 ± 0.01 0.08 ± 0.01 0.08 ± 0.009 0.025 ± 0.006 0.04 ± 0.006 0.05 ± 0.01 0.04 ± 0.01

− 0.003 ± 0.010 − 0.006 ± 0.004 − 0.022 ± 0.002 0.004 ± 0.004 0.001 ± 0.005 0.006 ± 0.004 0.005 ± 0.003 0.006 ± 0.001 0.004 ± 0.003 0.002 ± 0.001

contrast in areas associated with large susceptibility offsets such as the GP and SN. Veins are also easily discernable as areas of positive susceptibility, consistent with a high deoxyhemoglobin content. The RSO method shows similar large scale detail to the MO method, with the GP, SN, and large veins having high positive Δχ values relative to the surrounding tissue. However, the finer, more subtle detail present in the MO method, such as the GM/WM contrast in the mid and upper brain regions, is attenuated. The χ-map calculated using the TSO method shows increased artefacts relative to the other methods, particularly around blood vessels.

Fig. 5. Transverse modulus images, field map data, and susceptibility maps for a single subject at three levels in the brain. (a, f and k) Modulus data. (b, g and l) Field map data. (c, h and m) Susceptibility maps calculated using MO method. (d, i and n) Susceptibility maps calculated using RSO method. (e, j and o) Susceptibility maps calculated using TSO method.

522

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

Analysis of the phantom experiments indicates that it is reasonable to use the results of the MO method as benchmark data when assessing the other susceptibility mapping methods. We therefore plotted the susceptibility differences measured in the different brain regions using each of the three mapping methods against the results produced using the MO method in Fig. 6. For comparison, differences measured directly from the field map are also shown. In each case a linear fit was carried out and is depicted as a dashed line on the plot. Table 3 details the linear fits, along with the correlation coefficient produced by comparing the RSO, TSO, and field map values to the results of the MO method. The RSO method yielded average susceptibility offsets in good agreement with the MO results. The TSO method however, underestimated the susceptibility values compared to the MO method. Discussion The results presented here suggest that susceptibility mapping at high field is a viable and robust technique. High quality susceptibility maps were produced from image data acquired in vivo using a multiple-orientation method and similar, but slightly lower quality maps, were generated using two different single-orientation methods. The mapping procedure benefited from the use of a novel high-pass filtering methodology (Fig. 1) which did an excellent job of removing remotely generated, slowly varying fields, while preserving the local field perturbations from which the susceptibility variation can be inferred. This filtering method is equally applicable to other imaging techniques, such as SWI, in which phase variation due to local field inhomogeneities is exploited. The experimental data from the dopedagar phantom allowed accurate optimisation of each susceptibility mapping method to be carried out. For the MO method, the phantom results showed that an accurate χ-map could be created from data measured at four different orientations produced by small rotations about two orthogonal axes (Fig. 4c). For the RSO method, the phantom data allowed the correct choice of regularisation parameters to be made (Fig. 4e and f). Similarly, the phantom data allowed an acceptable k-space threshold to be chosen for the TSO method that gave the best compromise between contrast and noise (Fig. 4d). Comparing the results produced from the phantom using the different methods, it is clear the MO method yielded the highest quality susceptibility map. This is free of noise-related artefacts and shows susceptibility values which are in good agreement with those

Fig. 6. Δχ values measured in each brain region using the RSO (blue), and TSO (red) methods plotted against the results for the MO method. For comparison the measured field differences are plotted against the MO results (green). The dashed lines represent least-squares linear fits. The regions associated with each Δχ value are labelled.

Table 3 Calculated linear fits and correlation coefficients of RSO, and TSO methods to the MO method. For comparison, the field map fit and correlation are also included.

RSO TSO Field Map

Linear fit to MO (ppm)

Correlation to MO data

0.94 × ΔχMO + 0.0043 0.74 × ΔχMO − 0.0008 0.036 × ΔχMO − 0.0035

0.99 0.97 0.27

expected from the known concentrations of contrast agent (Table 1). The TSO method slightly underestimated the susceptibility values and produced maps showing “streaking” artefacts, but which still contain readily interpretable structural information about the underlying susceptibility distribution. The RSO method, while not achieving quite the same quality of χ-map as the MO method, out-performed the TSO method in reducing artefacts and producing accurate average susceptibility values. This is a consequence of the RSO method's utilisation of the structural information present in the modulus data. The in vivo experiments showed that high quality whole-brain susceptibility mapping at 7 T is feasible, and allowed further characterisation of the strengths and weaknesses of the different mapping methods in the context of human brain imaging. The susceptibility maps created by the MO method show exquisite detail, with particularly high positive Δχ values in the GP, PU, CN, SN, and RN regions, as well as a negative offset in the IC, as shown in Fig. 5c, h and m and in Table 2. The rich structural detail and high contrast to noise present in the field map are improved upon, in the χ-map generated using the MO method. This is particularly emphasised in the upper brain slice, shown in Fig. 5m, where the χ-map shows clear boundaries between GM and WM in the parietal lobe. In contrast, the field map of the same region, shown in Fig. 5l, displays varying field shifts at the GM/WM interfaces which make boundaries harder to discern relative to those in the χ-map. Another impressive property of the χ-map produced using the MO method is the visualisation of veins. The characteristic dipole artefacts that occur in and around large blood vessels (Deistung et al., 2008) in field maps are completely eliminated in the susceptibility map which shows high, positive, and localised Δχ-offsets in the veins. This effect is particularly clear in the large vessels close to the SN/RN region in the lower brain slices shown in Fig. 5b and c. The susceptibility maps calculated using the RSO method display similar large scale structure to the MO χ-map with the GP and SN regions particularly visible. However, finer more subtle detail such as GM/WM boundaries in the cortices appear less well defined. The behaviour of the RSO method can be explained to some extent by considering the definition of compartments in the associated modulus data (Fig. 5a, f and k). Since the regularisation algorithm uses spatial priors based on the gradient of the modulus data, structures that are clearly defined in the modulus data, such as the GP shown in Fig. 5f and the cylindrical compartments of the phantom depicted in Fig. 4a, are likely to be better represented in the final χ-map, than less clearly defined structures, such as the GM and WM in the upper brain sections (Fig. 5k). To demonstrate the dependence of the RSO method on the choice of spatial priors, some simple simulations were carried out, the results of which are shown in Fig. 7. A digital phantom containing two small spheres 12 voxels in radius were set inside a larger sphere of radius 50 voxels, itself set inside a 1283 matrix of zeros. The small spheres were given susceptibly values of 1 ppm relative to a value of zero for the rest of the phantom. A field map was simulated using Eq. (1) and masked using the larger sphere and Gaussian noise with a standard deviation of 0.02 ppm was added (Fig. 7a). Modulus data were also created by giving the small spheres signal magnitudes of 0.2, and the larger sphere a magnitude of 1, Gaussian noise with standard deviation 0.01 was then added (Fig. 7e). The RSO method was then used to calculate a χ-map (Fig. 7b) utilising edge information from the modulus data to form the spatial priors (Fig. 7f). As would be expected

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

523

Fig. 7. Results of simulations demonstrating the dependence of the RSO method on the choice of spatial priors. (a) Simulated field map. (e) Simulated modulus data. (f) Composite spatial prior based on (e). (b) χ-Map calculated using (a) and (f). (g) Spatial prior with one sphere removed from the modulus data. (c) χ-Map calculated using (a) and (g). (h) Spatial prior based on modulus data with one sphere replaced with a cube. (d) χ-Map calculated using (a) and (h).

from such high quality spatial priors, the RSO method yields an excellent χ-map, retrieving accurate susceptibility values within the spheres. If one of the spheres is removed from the modulus data, the spatial prior (Fig. 7g) contains no edge information and the RSO method yields a smoothed susceptibility map in this region (Fig. 7c). This is analogous to the low contrast between GM and WM in the cortex in the in vivo modulus data sets (see Fig. 5k). We can also add erroneous modulus information by changing one of the spheres to a cube. Use of this incorrect spatial prior (Fig. 7h), yields a susceptibility map (Fig. 7d) reflecting the underlying cubic region in the modulus data. The results of this experiment suggest the performance of the RSO method is highly sensitive to the choice of spatial priors. Accordingly, the RSO method should be strongly dependent on the imaging parameters used to acquire the modulus data, such as TE, TR and flip angle. Investigating this dependence on scanning parameter selection was beyond the scope of this work, but forms an important area for future study. It should also be noted that there are no restrictions on the image contrast that is used to generate the spatial priors. Here, the modulus signal was used as it comes free with the phase information, but other image data could also easily be used (e.g. T1- or T2-weighted images). Recently, Yao et al. (2009) showed that R*2 measurements at high field are highly correlated with iron content, making R*2-maps a potentially rich source of susceptibility-based structural information for driving RSO calculations. Susceptibility maps calculated using the TSO method showed similar bulk structure to the MO χ-maps, but also displayed “streaking” artefacts especially around areas of the brain with large Δχ-offsets, such as veins (Fig. 5e, j and o). Veins which appear clearly with high, positive Δχ-offsets in the χ-maps produced using the MO and RSO methods are harder to visualise in the TSO-generated map. The artefacts around veins most likely result from the over-simplified treatment of noise in the TSO method. Both the MO and RSO methods employ weighting matrices to reduce the contribution of phase measurements in regions where the noise is high. In contrast, the TSO method does not use noise information and weights all phase measurements equally. Thus, in areas of high noise such as around veins, where the magnitude signal has significantly decayed due to the short T*2 of venous blood, the

susceptibility maps calculated using the TSO method show severe artefacts. The in vivo susceptibility measurements shown in Table 2 and Fig. 6 indicate that the TSO method underestimated Δχ by about 20% relative to the MO and RSO methods. This is most likely due to the unavoidable removal of k-space data that is implicit to the TSO method. The measured Δχ values for the TSO method in this study are higher than those found in recent work by Shmueli et al. (2009). However, the threshold parameter used here (0.06) was significantly lower than the equivalent truncation value used by Shmueli et al. (0.2), meaning that less k-space data were effectively nulled in our study. This comparison highlights the importance of keeping a fixed threshold parameter when carrying out studies of multiple subjects. In this study, the single-orientation methods utilised a single data set which was acquired in about 3 min, while the MO method utilised four data sets which took a total of 12 min to acquire. Field maps acquired at a single orientation over 12 min would obviously have a higher signal to noise ratio (SNR) and use of such maps might provide a fairer comparison with the MO data. However, we note that just increasing the SNR of the field maps cannot compensate for the illposed nature of the deconvolution problem. This point was discussed in detail by Wharton et al. (2010) who showed that averaging over multiple scans at a single orientation yielded little gain in the quality of the final χ-map created using a threshold-based method. Comparison of the field measurements to the Δχ values produced using the MO method, in particular the low correlation coefficient listed in Table 3, highlights the effects of the non-locality of the field perturbations arising from a given susceptibility distribution. This illustrates the fact that care must be taken when interpreting differences in phase or field measurements between brain regions. Although the correlation of field and susceptibility is lower than in other phase based studies (Xu et al., 2008; Shmueli et al., 2009), the use of different filtering approaches and ROI definitions in the different studies make direct comparison difficult. In general, the in vivo susceptibility maps produced in this work, showed large Δχ offsets and well resolved structural detail in the GP, SN, RN, IC, PU, and CN regardless of the mapping method utilised. This suggests that these brain regions have particularly large Δχ offsets

524

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525

Instead, a small number of subjects could be scanned at multiple orientations with the desired FOV and voxel size, in preliminary experiments. χ-Maps calculated by applying the MO method to these data could then be used to calibrate a single-orientation method, such as the RSO method. The calibration would involve choosing appropriate regularisation parameters, as well as the selection of scanning parameters such as TE, TR and flip angle which would yield modulus data showing the structural information required by the RSO method. The bulk of the subjects in the study could then be scanned at a single-orientation and χ-maps calculated using the calibrated RSO method. Conclusion

Fig. 8. Average susceptibility values calculated from the MO method plotted against estimated non-heme iron content (Hallgren and Sourander, 1958) for GMFL, TH, PU, CN, RN, SN, and GP brain regions.

relative to surrounding tissues. The positive offsets in the GP, SN, RN, PU, and CN are most likely due to the high concentration of non-heme iron known to exist in these brain areas (Hallgren and Sourander, 1958; Haacke et al., 2005). Fig. 8 shows the average susceptibility values calculated from the MO method plotted against the estimated nonheme iron content based on the work of Hallgren and Sourander (1958), for brain regions common to both studies. The excellent correlation (r = 0.96) between susceptibility and estimated iron content suggests that susceptibility differences in the human brain are closely linked to iron content (Shmueli et al., 2009). The slope of the plot is found to be 0.75 ± 0.1 ppm per mg iron per g tissue agreeing within error with the value of 0.72 ppm per mg iron per g tissue that has previously been quoted (Duyn et al., 2007). The negative Δχ-offset of the IC could be attributed to the effect of the large amount of myelin in the IC (Mori et al., 2006; Fukunga et al., 2010). As would be expected from Fig. 6, a plot of field rather than susceptibility against non-heme iron content produces a much poorer correlation (r = 0.43). The sensitivity of susceptibility mapping to iron and myelin rich regions could be particularly beneficial in clinical studies. Recently, for example, field shifts based on phase data measured in the GP, PU, and CN, of patients with multiple sclerosis were shown to be significantly different to those measured in healthy subjects (Hammond et al., 2008). Susceptibility mapping could yield further insight into these differences by eliminating confounding effects of the non-local relationship between field measurements and the underlying Δχoffsets. In addition, changes in the levels of non-heme iron present in the SN are thought to be closely linked with the progression of Parkinson's disease (Sofic et al., 1988). High resolution susceptibility mapping could form a useful tool for accurately characterising these changes. While this work focuses on 7 T data, Liu et al. (2010) recently applied a method that is similar to the RSO method utilised here to data acquired at 1.5 T and measured susceptibility values of 0.18 and 0.13 ppm in the GP and SN, respectively. These results, which are in good agreement with values calculated in our study (see Table 2), suggest that despite the coarser spatial resolution available at lower field strengths, acceptable bulk susceptibility values can be measured in reasonably sized brain structures, potentially making susceptibility mapping a viable tool at 1.5 and 3T. Finally, we speculate that the methodology used here to compare the different susceptibility mapping methods could also be beneficially utilised in studies involving longitudinal measurements and/ or large numbers of subjects in which the MO method may not be suitable due to the need for subjects to execute head rotations.

Three susceptibility mapping methods have been optimised for application to whole-brain modulus and phase data acquired at 7 T. These are a multiple-orientation method, a regularised single-orientation method, and a threshold-based single-orientation method. Comparison of these methods showed that the multiple-orientation method yielded high quality susceptibility maps, out-performing the single-orientation methods. The threshold-based single-orientation approach was fast and simple to use, but underestimated susceptibility values and showed significant noise-related artefacts. The regularised method yielded contrast dependent on the choice of spatial priors, but demonstrated the potential to yield single-orientation susceptibility maps of a similar quality to those calculated using data acquired at multiple orientations to the field. Acknowledgments The authors would like to thank Ali M. Al-Radaideh for his help in the scanning of human subjects. This work was funded by MRC Programme Grant G9900259 and an MRC Ph.D. studentship for S.W. References Annese, J., Pitiot, A., Dinov, I.D., Toga, A.W., 2004. A myelo-architectonic method for the structural classification of cortical areas. Neuroimage 21 (1), 15–26. Christopher, C.P., Michael, A.S., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Softw. 8 (1), 43–71. de Rochefort, L., Brown, R., Prince, M.R., Wang, Y., 2008. Quantitative MR susceptibility mapping using piece-wise constant regularized inversion of the magnetic field. Magn. Reson. Med. 60 (4), 1003–1009. de Rochefort, L., Liu, T., Kressler, B., Liu, J., Spincemaille, P., Lebon, V., Wu, J., Wang, Y., 2010. Quantitative susceptibility map reconstruction from MR phase data using bayesian regularization: validation and application to brain imaging. Magn. Reson. Med. 63 (1), 194–206. Deistung, A., Rauscher, A., Sedlacik, J., Stadler, J., Witoszynskyj, S., Reichenbach, J.R., 2008. Susceptibility weighted imaging at ultra high magnetic field strengths: theoretical considerations and experimental results. Magn. Reson. Med. 60 (5), 1155–1168. Duyn, J.H., van Gelderen, P., Li, T.Q., de Zwart, J.A., Koretsky, A.P., Fukunaga, M., 2007. High-field MRI of brain cortical substructure based on signal phase. Proc. Of The Natl Acad. Sci. 104 (28), 11796–11801. Fukunga, M., Li, T.-Q., van Gelderen, P., de Zwart, J.A., Shmueli, K., Yao, B., Lee, J., Maric, D., Aranova, M.A., Zhang, G., Leapman, R.D., Schenck, J.F., Merkle, H., Duyn, J.H., 2010. Layer-specific variation of iron content in cerebral cortex as a source of MRI contrast. PNAS 107 (8), 3834–3839. Haacke, E.M., Xu, Y.B., Cheng, Y.C.N., Reichenbach, J.R., 2004. Susceptibility weighted imaging (SWI). Magn. Reson. Med. 52 (3), 612–618. Haacke, E.M., Cheng, N.Y., House, M.J., Liu, Q., Neelavalli, J., Ogg, R.J., Khan, A., Ayaz, M., Kirsch, W., Obenaus, A., 2005. Imaging iron stores in the brain using magnetic resonance imaging. Magn. Reson. Imaging 23 (1), 1–25. Haacke, E.M., Mittal, S., Wu, Z., Neelavalli, J., Cheng, Y.C.N., 2009. Susceptibility-weighted imaging: technical aspects and clinical applications, part 1. Am. J. Neuroradiol. 30 (1), 19–30. Hallgren, B., Sourander, P., 1958. The effect of age on the non-haemin iron in the human brain. J. Neurochem. 3 (1), 41–51. Hammond, K.E., Metcalf, M., Carvajal, L., Okuda, D.T., Srinivasan, R., Vigneron, D., Nelson, S.J., Pelletier, D., 2008. Quantitative in vivo magnetic resonance imaging of multiple sclerosis at 7 Tesla with sensitivity to iron. Ann. Neurol. 64 (6), 707–713. Hestenes, M.R., Stiefel, E., 1952. Methods of conjugate gradients for solving linear systems. J. Res. Natl Bur. Stand. 49 (6), 409. Jenkinson, M., 2003. Fast, automated, N-dimensional phase-unwrapping algorithm. Magn. Reson. Med. 49 (1), 193–197.

S. Wharton, R. Bowtell / NeuroImage 53 (2010) 515–525 Jenkinson, M., Smith, S., 2001. A global optimisation method for robust affine registration of brain images. Med. Image Anal. 5 (2), 143–156. Koopmans, P.J., Manniesing, R., Niessen, W.J., Viergever, M.A., Barth, M., 2008. MR venography of the human brain using susceptibility weighted imaging at very high field strength. Magn. Reson. Mater. Phys. Biol. Med. 21 (1–2), 149–158. Kressler, B., de Rochefort, L., Liu, T., Spincemaille, P., Jiang, Q., Wang, Y., 2010. Nonlinear Regularization for Per Voxel Estimation of Magnetic Susceptibility Distributions From MRI Field Maps. Medical Imaging, IEEE Transactions on 29 (2), 273–281. Liu, T., Spincemaille, P., de Rochefort, L., Kressler, B., Wang, Y., 2009. Calculation of susceptibility through multiple orientation sampling (COSMOS): a method for conditioning the inverse problem from measured magnetic field map to susceptibility source image in MRI. Magn. Reson. Med. 61 (1), 196–204. Liu, T., de Rochefort, L., Khalidov, I., Prince, M., Wang, Y., 2010. Quantitative susceptibility mapping by regulating the field to source inverse problem with a sparse prior derived from the Maxwell Equation: validation and application to brain. 18th Annual Meeting of the ISMRM. Stockholm, p. 705. Marques, J.P., Bowtell, R., 2005. Application of a fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility. Concepts Magn. Reson. Part B Mag. Reson. Eng. 25B (1), 65–78. Mori, S., Wakana, S., Nagae-Poetscher, L.M., van Zijl, P.C.M., 2006. MRI atlas of human white matter. AJNR Am. J. Neuroradiol. 27 (6), 1384–1385.

525

Salomir, R., De Senneville, B.D., Moonen, C.T.W., 2003. A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts Magn. Reson. Part B Magn. Reson. Eng. 19B (1), 26–34. Schäfer, A., Wharton, S., Gowland, P., Bowtell, R., 2009. Using magnetic field simulation to study susceptibility-related phase contrast in gradient echo MRI. Neuroimage 48 (1), 126–137. Shmueli, K., de Zwart, J.A., van Gelderen, P., Li, T., Dodd, S.J., Duyn, J.H., 2009. Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data. Magn. Reson. Med. 62 (6), 1510–1522. Sofic, E., Riederer, P., Heinsen, H., Beckmann, H., Reynolds, G.P., Hebenstreit, G., Youdim, M.B.H., 1988. Increased iron (III) and total iron content in post mortem substantia nigra of parkinsonian brain. J. Neural Transm. 74 (3), 199–205. Stephen, M.S., 2002. Fast robust automated brain extraction. Hum. Brain Mapp. 17 (3), 143–155. Wharton, S., Schäfer, A., Bowtell, R., 2010. Susceptibility mapping in the human brain using threshold based k-space division. Magn. Reson. Med. 63, 1292–1304. Xu, X., Wang, Q., Zhang, M., 2008. Age, gender, and hemispheric differences in iron deposition in the human brain: an in vivo MRI study. Neuroimage 40 (1), 35–42. Yao, B., Li, T.-Q., Gelderen, P., Shmueli, K., de Zwart, J.A., Duyn, J.H., 2009. Susceptibility contrast in high field MRI of human brain as a function of tissue iron content. Neuroimage 44 (4), 1259–1266.