Quantitative susceptibility mapping (QSM): Decoding MRI data for a tissue magnetic biomarker

In MRI, the main magnetic field polarizes the electron cloud of a molecule, generating a chemical shift for observer protons within the molecule and a magnetic susceptibility inhomogeneity field for observer protons outside the molecule. The number of water protons surrounding a molecule for detecting its magnetic susceptibility is vastly greater than the number of protons within the molecule for detecting its chemical shift. However, the study of tissue magnetic susceptibility has been hindered by poor molecular specificities of hitherto used methods based on MRI signal phase and T2* contrast, which depend convolutedly on surrounding susceptibility sources. Deconvolution of the MRI signal phase can determine tissue susceptibility but is challenged by the lack of MRI signal in the background and by the zeroes in the dipole kernel. Recently, physically meaningful regularizations, including the Bayesian approach, have been developed to enable accurate quantitative susceptibility mapping (QSM) for studying iron distribution, metabolic oxygen consumption, blood degradation, calcification, demyelination, and other pathophysiological susceptibility changes, as well as contrast agent biodistribution in MRI. This paper attempts to summarize the basic physical concepts and essential algorithmic steps in QSM, to describe clinical and technical issues under active development, and to provide references, codes, and testing data for readers interested in QSM. Magn Reson Med 73:82–101, 2015. © 2014 The Authors. Magnetic Resonance in Medicine Published by Wiley Periodicals, Inc. on behalf of International Society of Medicine in Resonance. This is an open access article under the terms of the Creative commons Attribution License, which permits use, distribution, and reproduction in any medium, provided the original work is properly cited.


INTRODUCTION
Magnetic susceptibility is one of the following major categories of tissue contrast mechanisms in proton MRI (1): 1) spin thermal relaxation in a voxel of water; 2) water motion, including diffusion, perfusion, flow and tissue deformation; and 3) molecular electron cloud polarization by the main magnetic field B 0 . A polarized molecule generates its own magnetic field, which is known as a chemical-shift shielding field for observer protons inside the molecule and as a magnetic-susceptibility inhomogeneity field for observer protons outside the molecule. This field adds phase accumulation and consequently causes intravoxel dephasing or magnitude T2* decay in the commonly available gradient echo (GRE) MRI. Therefore, noninvasive MRI is well suited for investigating the magnetic susceptibility of tissue. The GRE phase is equal to the magnetic field multiplied by the gyromagnetic ratio g and the echo time (TEÞ. This phase may be used to further attenuate the signal for enhancing T2* image contrast, which is called susceptibility weighted imaging (2)(3)(4). However, the field at an observer location is the sum of contributions from all surrounding magnetic susceptibility sources, with each contribution dependent on the source-observer distance and orientation (5). Consequently, the phase or T2* contrast does not exclusively depict the local tissue magnetic property but is a weighted summation of the magnetic properties of the surrounding tissue, reflecting only the "shadow" of the surrounding susceptibility sources. For example, the phase and T2* contrast of tissues with weak susceptibility may primarily come from nearby air-tissue interfaces, across which there are large susceptibility changes. To uncover local tissue magnetic properties, the field has to be deconvolved, which is referred to as quantitative susceptibility mapping (QSM).
As indicated by a PubMed search of QSM papers (6 in 2011; 18 in 2012; and 37 in 2013), there is rapidly growing interest in developing techniques for QSM data acquisition and processing, and in developing clinical and scientific applications ranging from iron distribution and metabolic consumption of oxygen to myelin in white matter (WM) tracts (Fig. 1). This review tries to serve these interests by summarizing the basic physical concepts in QSM, outlining the fundamental algorithmic steps in QSM, organizing the available MATLAB (MathWorks, Natick, MA) codes for QSM algorithms, and surveying the clinical and technical QSM issues that are under active development.
Preparation: Estimating Susceptibility-Generated Field from its Effects on MRI Signal In the MRI main field B 0 , a molecule in tissue gains a magnetic moment p through its electron cloud polarization. Correspondingly, a tissue with volume magnetic susceptibility x gains magnetization m r ð Þ % x r ð ÞB 0 =m 0 (see TIS-SUE MAGNETISM in supporting information for a brief summary of the molecular physics). Tissue magnetization generates its own magnetic field that affects MRI signal.
Here we review the mathematical relationships that link magnetization, field, and MRI signal, based on which the field can be estimated from the MRI signal.

Magnetic Dipole Field and Field Observed by a Proton in Tissue
According to Maxwell's equations in vacuum, a magnetic dipole moment p at a source location r 0 generates a magnetic field b r ð Þ at an observation location r (5) (ê rr 0 is a unit vector along r À r 0 ), b r ð Þ ¼ m 0 4p 3 p Áê rr 0 ð Þ e rr 0 À p jr À r 0 j 3 j r 0 6 ¼r þ 2m 0 3 pd r À r 0 ð Þ: [1] Here, the inverse-cube of the distance term characterizes the spatial extent of the dipole field. The deltafunction term can be understood from the field of a current loop with a fixed magnetic moment and a radius ! 0 (Fig. 2). In water MRI, the delta-function term is dropped; the probability of the polarized electron cloud penetrating into the space of the observer water protons is negligible. Thus, the field (scaled to B 0 ) observed by a water proton is the sum of contributions from all surrounding susceptibility sources [their distribution defined by magnetization m r ð Þ], excluding that from the proton's own location: r 0 6 ¼r 3 m r 0 ð ÞÁê rr 0 ð Þ e rr 0 À m r 0 ð Þ jr À r 0 j 3 B 0 d 3 r 0 ¼ ðd Ã xÞ r ð Þ: [2] The exclusion of the observer point in the integration of Eq. [2] represents the Cauchy principal value, commonly known as the Lorentz correction. This is conventionally, but unnecessarily, interpreted as subtracting a sphere of magnetization m that has a uniform interior field 2m 0 m=3 from the sum of the fields of all sources according to Eq. 1 (5). In Eq. [2], d r; r 0 ð Þ¼ 3mðr 0 ÞÁê rr 0 ð Þ e rr 0 Àmðr 0 Þ 4p jrÀr 0 j 3 if expression is true and 0 otherwise.) Eq. [2] relates the field at r to the magnetization distribution over the whole space. This can be expressed in a differential form that relates the field at r to the magnetization located at r (31): Eq. [3] can be Fourier Þ . For notational convenience, we use lower case for r-space quantities and upper case for corresponding k-space quantities (except On the other hand, the electron cloud of a molecule does penetrate into and interact with observer protons within the molecule. Consequently, electron cloud polarization by B 0 induces a shielding magnetic field À sðrÞB 0 (Fig. 2) (32), or chemical shift (referenced to water), that alters the field experienced by protons within the molecule: b r ð Þ ¼ ðd Ã xÞ r ð Þ À sðrÞB 0 : Both the magnetic susceptibility (observed by a large number of water protons outside the molecule) and the chemical shift (observed by protons inside the molecule) reflect the same molecular electron-cloud polarization (33)(34)(35).

Field Effects on MRI Signal
The magnetic field of a polarized molecule may affect the MRI signal magnitude through a chemical exchange between free water and water bound to the molecule (inner sphere relaxation) and through a free water diffusion in the field (outer sphere relaxation) (36). These complicated effects are characterized as relaxation enhancement (37)(38)(39)(40)(41). Susceptibility estimation from MRI signal magnitude affected by relaxation is prone to large errors (31). Fortunately, the phase of a water proton spin linearly increases with the field. Using multiple radio frequency (RF) coils, with the ath coil element having a complex coil sensitivity function c a ðrÞ and acquisition noise N a , the k-space signal measured in coil a at time t is (1): Here, v 0 ¼ gB 0 and mðrÞ is the proton transverse magnetization [mðrÞ is proportional to proton density and is much smaller than m r ð Þ the electronic magnetization]. Typically, T2 is homogeneous within a voxel and is much longer than the readout duration (e Àt k ð Þ=T2ðrÞ % e ÀTE=T2 ), and the data sampling gradient 2p g @k @t is much larger than rbðrÞ (e Àiv0b r ð ÞtðkÞ % e Àiv0b r ð ÞTE ).
Then, the complex image mðrÞe Àiv0b r ð ÞTE can be reconstructed using the Moore-Penrose pseudo-inverse (1) for both full k -space sampling (42) and parallel imaging (43). Approximating the voxel sensitivity function (44) as a box function, the signal detected at TE in a voxel centered at r with size Dr (volume DV ) is the summation over signals from all proton spins in the voxel, During metabolic consumption of oxygen in the brain, heart, and kidney, weakly diamagnetic oxyhemoglobin releases O 2 and becomes strongly paramagnetic deoxyhemoglobin. Whereas the 3d electron orbits of Fe 2þ in deoxyhemoglobin may be approximated as an isolated iron atom with four unpaired electrons (right), the intramolecular interaction between the porphyrin ring and Fe 2þ in oxyhemoglobin (ligand interaction) splits the Fe atom's 3d-orbit into two levels, e g and t 2g , with all six electrons paired in the three t 2g orbits. (b) Blood degradation in hemorrhage. Following the onset of a hemorrhage, a small fraction of red blood cells (RBCs) may be endocytosed by microglia/macrophages. The majority of RBCs undergo cell lysis and hemoglobin (Hb) degradation from deoxyhemoglobin into methemoglobin (Fe 3þ ) and hemosiderin (possible magnetic domain). Modeled after: Lancet Neurol 2012;11:720-731. (c) Susceptibility sources in the human brain. Major susceptibility sources in (i) the brain include myelin and ferritin. The white matter tracts in the brain consist of myelinated nerve fibers. (ii) Zoomed view of the box in (i) showing axon (yellow) and myelin sheath (purple). Myelin consists of several layers of lipid bilayer. (iii) Zoomed view of the box in (ii) showing a lipid bilayer and an individual lipid. (iv) Ferritin in a cross-section. Ferritin consists of a peptide spherical shell 2-nm thick with a 8-nm diameter cavity. Fe 2þ enters through a four-fold symmetric channel, is stored as Fe 3þ oxide mineral, and is released as Fe 2þ through a three-fold symmetric channel. There are five unpaired 3d electrons in Fe 3þ , generating strong paramagnetism.
T2 Ã e Àib r ð Þv0TE þ nðrÞ: [5] Here, the integration is evaluated by assuming small phase inhomogeneity in a voxel: e Àiv0ððbðr 0 ÞÀbðrÞÞTE $ 1 À iv 0 ððbðr 0 Þ À bðrÞÞTE À ðv 0 ððbðr 0 Þ À bðrÞÞTEÞ 2 =2 þ oð3Þ: The quantities in Eq. [5] are defined as: the noise in the reconstructed image. Common imaging situations involve weak tissue susceptibilities, short TEs (and small echo spacing in multiple echo sampling), and small voxels. Then, the phase of the signal in a voxel is approximately equal to gTE times the proton density-weighted average of the field in the voxel, and the T2* decay rate differs from 1 T2 by 1 2 g 2 TE times the proton density-weighted field variance in a voxel. As a result, MRI phase signals allow us to measure the magnetic field generated by tissue susceptibility.

Pulse Sequences Sensitizing Susceptibility
The three-dimensional (3D) GRE sequence is currently the most frequently used to acquire data for QSM. Multiple echoes (number of echoes N e ¼ 1-12 having been reported in literature) can be used with short TE for detecting strong susceptibility and long TEs for weak susceptibilities. The GRE sequence on a typical 3T MRI system may allow a high-resolution acquisition ($0.5 Â 0.5 Â 2.0 mm 3 ) with minimal first TE $5 msec, uniform TE spacing $5 msec, readout bandwidth $300 Hz/pixel, last TE $45 msec, and TR $50 msec (whole brain scan time $5-10 min). For imaging regions with very strong susceptibility variations, such as those containing the bone-tissue and air-tissue interfaces, ultra-short (45,46) or zero TE (47) and small TE increments may be used to achieve adequate signal-to-noise ratio (SNR). Flow compensation can be used to account for phases caused by motion (48). Resolution should be as high as possible to minimize averaging effects in Eq. [5]. Parallel imaging and spiral sampling trajectories (49,50) can be used to reduce scan time. Susceptibility effects may also be detected with asymmetric spin echo (SE) sequences. It is also possible to measure the magnetic field in a voxel through direct saturation of the water protons at a range of off resonance frequencies and Lorentzian lineshape fitting (water saturation shift referencing) (51), similar to continuous-wave NMR spectroscopy. Whereas GRE sequences generally are faster than SE sequences, there may be applications that benefit from a combination of GRE and SE acquisitions.

Digitized Dipole Field: Field and Susceptibility Averaged over a Voxel and Localized at the Voxel Center
The field averaged over a voxel in Eq. [5] requires digitization of Eq. [2]. The assumption that the distributions of the observation protons and susceptibility sources are fairly uniform within a voxel (physical smoothness) may allow for digital interpretation of Eq. [2]: r as the voxel center, b r ð Þ as the field value at r, and x r ð Þ as the average of the susceptibilities in the source voxel positioned at r. This digitization can be analytically derived for spherically shaped voxels: The field's spherical mean value equals the field value at the sphere center (5); the field outside of a uniformly magnetized sphere is equivalent to the field of a dipole defined by the sphere's total moment placed at the sphere's center; and the field inside the sphere is zero (with the Lorentz correction by Eq. [2]). The physical smoothness requirement is reasonably valid for high resolution MRI, and sincinterpolation or zero-padding (44) may be used to reduce the digitization error. Accordingly, Eq. [2] will be used to model the MRI signal with voxels indexed by r.

QSM
Step 0: Estimate the Total Field from MRI Data Eq. [5] models the MRI signal s r; TE ð Þ as a product of an exponential factor e Àib r ð Þv0TEj that describes a phase linear in time and a complex amplitude a r; TE ð Þ that contains a constant phase and an amplitude decay with time. Then, the total field bðrÞ can be estimated from MRI signals at multiple TEs as a nonlinear least-squares problem per voxel (20,21,27) js r; TE j À Á À a r; TE j À Á e Àib r ð Þv0TEj j 2 : [6] In the case of large SNR, Eq. [6] can be approximated as a linear regression of the signal phases at various TE j , allowing for a closed form solution (21). In general, Eq. [6] can be solved iteratively and efficiently using a gradient The electron cloud of a molecule polarized by B 0 generates magnetic shielding or chemical shift for the observer proton in the molecule and a susceptibility field (in dipole pattern) for the observer proton outside the molecule.
search method, such as the Gauss-Newton (20) or Levenberg-Marquardt method (52), with bðrÞv 0 initialized as the linear rate of phase evolution between the first two echoes (N e > 1). Eq. [6] may also be solved rapidly using autoregression (53). Interestingly, although Eq. [5] requires the estimation of complex coil sensitivity from calibration data (42,43), the field can be estimated without full knowledge of complex coil sensitivities from the phase difference at two TEs (54)(55)(56)(57) or even at a single echo (58).
The field estimated from Eq. [6] can contain many artificial jumps from voxel to voxel because e Àiw is periodic in w and can only define w up to a period such as ½Àp; pÞ. Consequently, b r ð Þv 0 TE j outside ½Àp; pÞ is aliased or wrapped into ½Àp; pÞ, resulting in abrupt artificial jumps of J Á 2p (J is an integer) in the output phase, which is called the wrapped phase denoted by w w . To estimate b r ð Þ, these phase wraps have to be compensated by adding J Á 2p as needed: J r ð Þ is determined by the physical consideration that the unwrapped phase is spatially continuous. Unfortunately, it is difficult to impose continuity in discrete space because discretization is an undersampled approximation of the continuous space, and noise-ubiquitous in MRI data-causes jumps in regions of low SNR. There are many algorithms to impose continuity or smoothness in phase unwrapping (59), typically via path finding (including noncontinuous path, residual balancing, or branch cut) (60) or global error minimization (least squares) (Fig. 3). Phase can also be unwrapped rapidly by the Fourier spectral solution of the Laplacian of the phase (23,61); however, second order derivatives in the Laplacian may not allow large phase changes to occur in regions with big susceptibility variations, such as near veins (62). Phase unwrapping may be avoided for field data estimated from direct water saturation lineshape (51). A frequently used method is the path-finding type, image-quality-guided region growth (63), which works robustly and rapidly on MRI data using SNR as the image-quality guidance (17,64,65).
After phase unwrapping, a few bad points with fields not consistent with Eq. [2] may remain due to turbulent flow or idiopathic causes, and their effects can be reduced using a consistency check during QSM iteration (20). As the phase measurement in MRI is relative to RF carrier frequency, the field can only be determined up to a constant uniform field, which may be removed by background field removal in the next section.

QSM Algorithm: Formulating and Solving the Field-to-Susceptibility Inverse Problem
The goal of QSM is to determine tissue susceptibility from the field (Eq. [6]) or the complex MRI signal (Eq. [5]) using Eq. [2], which connects the magnetic field and the tissue susceptibility. There are two fundamental challenges that are imbedded in Eq. [2] and Eq. [5], as outlined below, and we describe two QSM steps to address them correspondingly: 1) background field removal and 2) dipole inversion (Fig. 3).

Two Fundamental Challenges in QSM
The first fundamental challenge is the lack of MRI signal in regions with susceptibility sources. MRI signal in Eq. [5] can only be detected in the region with water or the tissue of interest (V). Magnetic susceptibility sources exist outside V ðR 3 nVÞ, which is the background. However, there is no MRI signal in the background. If we regard each voxel in MRI as a field detector, then the number of detectors is less than the number of sources, making the fieldto-source inverse problem ill-posed. Therefore, this lack of MRI signal in the background forms the first fundamental challenge for QSM. Susceptibility sources in tissue (V) generate the tissue field b t ðrÞ ¼ m 0 R r 0 2X r 0 6 ¼r d r À r 0 ð Þx r 0 ð Þd 3 r 0 for r 2 V. Background susceptibility sources generate the background field ð Þ is tissue susceptibility when r 2 V. The second fundamental challenge is the zero cone surface in the dipole kernel in r-space or in kspace: u ¼ 654:7 0 with respect to the main magnetic Step 0: Generate magnitude image and field with unwrapping.
Step 1: Remove background field and generate tissue field.
Step 2: Generate QSM from tissue field and magnitude image.
field (magic angle) ( Fig. 4) (14). The dipole kernel in kspace is nearly flat (with all derivatives vanishing) at the ends of the cone neighborhood G e ¼ k; jD k ð Þj e f g : . . ., implying that the standard analysis cannot be used to define X k ð Þ in G e (66). Noise in the data would allow many possible susceptibility solutions that differ by an arbitrary amount in G e for a given field map B k ð Þ, causing a substantial dipole kernel null space. Therefore, the magnetic field-to-susceptibility inverse problem is illposed and lacks a unique solution (21).

QSM Step 1: Background Field Removal
To address the first QSM challenge of the lack of MRI signal in the background regions with susceptibility sources, prior knowledge is needed. The background field may be regarded as slow varying and can be removed by highpass filtering (HPF). Unfortunately, HPF erroneously removes the low spatial frequency components of the tissue field and fails to remove the high spatial frequency components of the background field near the tissue border (3,67,68), causing substantial errors in QSM. Better priors on the separation between the background and tissue fields are required for accurate determination of the tissue field, such as approximate orthogonality in the projection onto dipole fields (PDF) method (69) or the harmonic property in the sophisticated harmonic artifact reduction on phase data (SHARP) method (70). These priors are systematically organized in the following as solutions to Maxwell's Equations (Eq. [3]).
By definition, the background field b b has no source in V, whereas the local tissue field b t has sources m 0 ðr 2 m=3 Àr r Á m ð ÞÞ inside V (Eq. [3]). From electrodynamics (5) or partial differential equation (71), for a finite domain V, a unique solution to Laplace's equation Eq. [3] can be obtained according to the values at the boundary @V (5,72). In MRI, the tissue field is typically weak. At @V, the background field may be approximated as the total field. This is the Laplacian boundary value (LBV) method ( Fig. 5) to estimate background field (73).
Eq. [8] can be solved numerically as a system of linear equations by expressing the Laplacian with a difference operator (74): Here, o NN denotes the nearest neighbor average operator. It is advantageous in SNR to express Eq. [8] using the spherical mean value (SMV) operator o SMV (5,70,75): The system matrix is diagonally dominant (with diagonal elements ¼ À1), and the Jacobian iterative method can be used (74) ð Þ ¼ bðrÞ that satisfies the boundary value condition in Eq. [8], the background field can be obtained numerically by repeatedly applying o SMV , which is the iterative spherical mean value (iSMV) method (76). However, iSMV is slow to converge. Eq. [8] can be efficiently solved using the full multi-grid solver that first finds the solution on a coarse grid and then successively refines the solution on finer grids (74). (v, vi) Susceptibility in image space obtained by truncated k-space division with the threshold t ¼ 0 and 0.1. As a consequence of the dipole kernel zero behavior in the cone surface neighborhood G t , there is substantial noise propagation from the field measurements into the susceptibility estimate (40), as illustrated in an example of reconstruction by direct division (v and vi). A little noise added in the phase map (peak SNR ¼ 20) leads to a totally corrupted susceptibility image that bears no physical resemblance to the true susceptibility source.
When the boundary condition is not known, a partial differential equation in difference form may be regarded as an ill-posed problem, which can be solved by regularization. One regularization is to impose spectral truncation when evaluating the inverse Laplacian (therefore altering the inverse Laplacian), Here, ð Þ 2 , using truncated singular value decomposition (TSVD) (with a carefully chosen truncation value l) in k-space, and r 2 b % o SMV À 1 ð Þ b, using o SMV for possible denoising benefit. This is the sophisticated harmonic artifact reduction for phase data (SHARP) (70). r 2 b in the border layer V À V 0 may be recovered using a continuity assumption as in harmonic (background) phase removal using the Laplacian operator [HARPERELLA (62)]. Another regularization is the minimal norm of the tissue field (implicitly modifying inverse Laplacian) (62,77): Here, jjf jj 2 2 ¼ P r2X jf ðrÞj 2 for any field f ðrÞ (squared L2 norm). This method is called regularization-enabled SHARP (RESHARP) (77). Eq. [10] can be expressed in a 2 . There is another approach to estimate the background field based on the equivalent source or charge simulation method (78). The background field can be represented by fields of dipoles outside V that are approximately orthogonal to the fields of dipoles inside V except near @V. Then the background field can be estimated by all possible PDFs (17,22,69): Here, noise-weighting w ¼ P j TE j a r; TE j À Á is the phase SNR, which is assumed to be large by linearizing the signal model Eq. [5]. Note that b b in tissue V is unique, although x b in background R 3 nV is not. This PDF method may provide a slight advantage in dealing with noise in b by avoiding the inversion of the Laplacian in Eq. [8] through Eq. [10], and by extending to the nonlinear form, as in Eq. [6] (20). However, the orthogonality between fields of dipoles breaks down near the boundary, making PDF prone to over-fitting errors near the border.
Because background and tissue fields are approximately orthogonal, the minimal norm of the tissue field in Eq. [10] is similar to the dominance of the background field at @V in Eq. [8], making RESHARP similar to LBV. The minimization in Eq. [10] is similar to that in Eq. [11], without noise weighting w, making RESHARP similar to PDF. Therefore, all four methods (LBV, SHARP, RESHARP, and PDF) based on Maxwell's equations in Figure 5 provide similar performance, whereas there are very strong values near the tissue boundary in the HPF processed tissue field. The assumption or regularization in any method may contain errors, which propagate through into errors in the final reconstructed tissue susceptibility, a challenge for future research.

QSM Step 2: Field-to-Susceptibility Inversion
To address the second fundamental challenge in QSM, the ill-posedness caused by the dipole kernel zeroes, prior knowledge is again needed to uniquely identify the susceptibility component in the dipole kernel null space. For simplicity, we consider the high SNR case for which the phase noise is approximately Gaussian with variance / j1=aj 2 [and we can extend to the general SNR case using the complex signal as in Eq. [6] (20)]. Then Eq. [5], after background field removal, is reduced to a linear field-tosusceptibility inverse problem (Eq. [2] with noise), Arbitrary susceptibility values in the dipole zero cone neighborhood G e are allowed in Eq. [12]. A regularization can be used to specify susceptibility values in G e . Alternatively, "missing data" about the susceptibility in one orientation can be recovered by reorienting the subject in a fixed magnet and resampling the MRI signal (14,79). This method is known as the calculation of susceptibility using multiple orientation sampling (COSMOS), which delivers an exact reconstruction by fully sampling the susceptibility (14,22). Unfortunately, acquiring FIG. 5. Background field removal using various algorithms. The tissue fields in a healthy volunteer estimated using HPF, LBV, SHARP, RESHARP, and PDF methods, respectively, demonstrate similar tissue patterns but with slight and different accents. HPF, high-pass filtering; LBV, Laplacian boundary value; PDF, projection onto dipole fields; RESHARP, regularization enabled SHARP; SHARP, sophisticated harmonic artifact reduction on phase data. multiple scans of patients in different orientations is not clinically acceptable. We should focus on the regularization approach to QSM using single orientation data.
Regularizations can be expressed in various mathematical forms including TSVD and assumptions of a smooth, sparse, or piece-wise smooth solution. Consequently, there are too many QSM methods to be characterized by a unifying framework. For the purpose of illustrating concepts, we outline two important classes of QSM methods: 1) the closed-form k-space approach, exemplified by TSVD based on matrix computation (80), and 2) the Bayesian approach based on optimization (52), also known as the r-space approach. Specific algorithm formulas, codes, and results are summarized in the next section on experimental validation. All of these solutions may suffer from model errors when the mathematical properties are not consistent with the physical reality (7,10,11,17,21,81). Additionally, noise will always propagate into the final solution.
Closed-form solutions form a class of noniterative kspace methods. One example is TSVD (82): 21,22). More commonly used in QSM is a TSVD variant called truncated k -space division (TKD) (10): Truncation obviously leads to an underestimation of the susceptibility values in G l , and consequently produces streaking artifacts along the magic angle in the susceptibility map. This error in the TKD method is ð Þj; l ð Þ N : the first part is the regularization error from points in G l where the kernel has been modified by truncation; the second part is the noise error from all data points in kspace, but interestingly is also dominated by points in and near G l . The underestimation in Eq. [13] may be compensated by a scale factor (26). The streaking in Eq. [13] may be reduced using a high threshold [l ¼ 2=3, susceptibility ¼ field convolving with a kernel (26)] or using iterative filtering (defined by G l in k-space) of signals outside high-susceptibility structures in r-space defined by a mask [iterative Susceptibility Weigted Imaging and susceptibility Mapping (iSWIM), (30)]. A variant of MN is a Tikhonov-regularized minimal gradient norm [(29): Bayesian regularizations form a class of iterative optimization methods. To allow optimal treatment of regularization error (83), prior information is regarded as a probability distribution function (pdf) typically expressed as / e ÀlRðxÞ=2 : l is a tunable regularization parameter and RðxÞ is a functional of the susceptibility map. Noise is also characterized by its pdf. Then, the optimal solution with minimal total error from regularization and noise is the maximum a posteriori (MAP) estimate (84). For Gaussian noise with pdf / e ÀjwðrÞn r ð Þj 2 =2 (w as in Eq. [11]) in the estimated field, the posterior proba-bility is e ÀlRðxÞ=2 Q r e ÀjwðrÞðb r ð ÞÀd r ð ÞÃx r ð Þj 2 =2 , the maximization of which is the MAP solution: x r ð Þ ¼ argmin The first term in Eq. [14] is the weighted data fidelity term, which contains noise; the second term is the regularization term, which contains the regularization error. The regularization parameter, l, may be chosen to provide a minimal total error in a given imaging situation, and is typically chosen such that the regularization error is approximately equal to the expected noise level, a criterion known as the discrepancy principle (85)(86)(87).
There are connections between Bayesian optimization and noniterative k-space methods. If w ¼ 1 in Eq. [14], then L2 norm-based Tikhonov regularization R x ð Þ ¼ jjxjj 2 2 leads to X MN k ð Þ and R x ð Þ ¼ jjrxjj 2 2 leads to X MGN k ð Þ; both are noniterative k-space methods. Because MRI phase noise requires spatially varying weighting w, noniterative k-space methods may suffer from noise errors (88). A wide range of forms for R x ð Þ in Eq. [14] have been reported, including piece-wise constant susceptibility (11,89), smooth susceptibility or susceptibility gradient (21), sparse susceptibility gradient or wavelet (21,27,90), and morphological consistency of the susceptibility map (17,19,21,25), some of which are detailed in the next section. We should note the property of the conjugate gradient method that typically is used to minimize the data fidelity term in Eq. [14] as in LSQR: The solution is initialized as zero by default, points outside the zero cone neighborhood G e are calculated firstly (and properly) according to Krylov sequence, and later iterations fill structured noise in G e that cause streaking artifacts (16,52,91) (Fig. 6). Use of a small iteration number (n ¼ 5) may be regarded as an implicit regularization for a solution with moderate streaking and zero value at the k-space center X 0 ð Þ ¼ 0 (16,88). For a solver of Eq. [14] that includes use of LSQR, such as in the Gauss-Newton method, its final solution may have X 0 ð Þ ¼ 0. Eq. [14] may only determine the susceptibility up to a constant [similar to the Neumann boundary value problem for Laplace's equation (5)]. Specific solvers may introduce implicit regularization in their output. A reference to water, such as the cerebrospinal fluid in the ventricles, may be used for consistent display of susceptibility values in QSM.
Available anatomic information in a specific imaging situation defines a physical prior for morphologyenabled dipole inversion (MEDI) (19,83,92). The edges in the desired susceptibility map are likely to be colocated with edges in known anatomical images, because they reflect the same anatomy. The colocalization of edges may be expressed as the sparsity of susceptibility edges outside the known edge locations using the L0-norm or the more manageable L1-norm (93). This minimizes streaking artifacts common to the dipole kernel nullspace. From Eq. [14], one MEDI implementation to reconstruct QSM can be formulated as (17,19): Here, d Ã x is evaluated in Fourier space by for anatomic image I and threshold t (t is chosen such that approximately 30% of voxels are labeled as edges), and jjMrxjj 1 ¼ P j;r2V jMr j xðrÞj(L1 norm). Because the dipole kernel (Fig. 4), similar to a HPF, preserves susceptibility edges in the tissue field image (5) and accordingly in the T2*weighted (T2*w) magnitude image, the GRE magnitude (17) and phase images (25) can be used as the anatomic images. The nonlinear Eq. [15] can be solved using the Gauss-Newton method (52) without explicit formation of the memory costly inverse Hessian matrix (19). The data fidelity term in Eq. [15] can be generalized to Gaussian noise in complex data (e.g., Eq. [6]) and can be solved using a procedure identical to Eq. 15 (20).
New priors continue to be proposed; the search for the best prior is ongoing. An optimal prior may be specific to the imaging application. There has been a preliminary attempt to compare several priors (88). Multicenter trials are needed to establish a consensus on applicationspecific priors, which leads to the topic of the next section.

QSM Source Codes and Experimental Validations
To enable readers to try QSM algorithms on their own, here we organize the available MATLAB (MathWorks) codes for certain QSM algorithms, along with testing data (http://weill.cornell.edu/mri/pages/qsm.html). We tried to consider all the QSM methods that were published before December 1, 2013, classifying similar methods into the same category. We asked the first or corresponding authors to share MATLAB (MathWorks) codes of their methods, proofread our implementations, and comment on results of testing data, and we thank them for their valuable feedback. The page limitation forced us to select one (perhaps the most widely used) of multiple algorithms published by each group, leading to the following seven methods in Table 1: TSVD (22); TKD (10); iSWIM (30); MEDI (19); compressed sensing compensated (CSC) inversion (27); homogeneity-enabled incremental dipole inversion (HEIDI) (25); and total variation using split Bregman (TVSB) (28). The first two methods do not require anatomic prior information. The third method, iSWIM, incorporates an anatomic prior iteratively into the k-space approach. The last four methods are based on the Bayesian approach and are listed in chronologic order. The Bayesian methods all aim to minimize a function consisting of a data fidelity term in the L2 norm (measuring noise power) and a regularization term in the L1 norm or total variation (promoting All the calculations were performed on a PC equipped with Intel V R Core i7-3770k CPU @ 3.5 GHz and 32 GB of memory, except (1) was performed on a personal laptop with Intel Core i5-M2450 CPU @ 2.5 GHz and 8 GB of memory, and (2) was performed on a PC with Intel Core i5-2320 CPU @ 3.00 GHz and 16 GB of memory. CPU, central processing unit; COSMOS, calculation of susceptibility using multiple orientation sampling; CSC, compressed sensing compensated; GB, gigabytes; GHz, gigahertz; HEIDI, homogeneity-enabled incremental dipole inversion; iSWIM, iterative susceptibility weighted imaging and susceptibility mapping; MEDI, morphology-enabled dipole inversion; PC, personal computer; TKD, truncated k -space division; TSVD, truncated singular value decomposition; TVSB, total variation using split Bregman.
FIG. 6. Evolution of susceptibility solutions in conjugate gradient. Susceptibility images in k-space (left column) and r-space after the first, third, fifth, 10th, and 100th iterations using the conjugate gradient solver demonstrate that the none-cone region in k-space converges quickly in the first few iterations; and the later iterations mainly contribute to signal in the cone region that manifests as streaking artifacts in the sagittal view and noise in the axial view in r-space. sparsity). Whereas MEDI and HEIDI attempt to sparsify the edge difference with known anatomical priors, CSC promotes image sparsity in the wavelet domain, and TVSB hugely accelerates reconstruction speed by dropping noise whitening. Rigorous experimental validation with a reference standard is required to assess the accuracy of any quantitative technique. This is particularly necessary for a technique involving regularization such as the QSM algorithms in Table 1. We performed three validation experiments: 1) numerical simulation with known truth, 2) MRI data of gadolinium phantom with known susceptibilities, and 3) in vivo brain MRI with susceptibilities assessed by COSMOS (92,94). These validations were imperfect, particularly because COSMOS contains errors from noise, orientation registration, and WM susceptibility anisotropy. However, they can serve as a starting point for readers to experience various QSM methods. Details of the experimental setup are described in Validation data in the supporting information. Although evaluations here are not intended to be conclusive, they allow readers to assess various QSM methods by their qualities (by examining the streaking artifacts), accuracies (by voxel-based linear regression with ground truth or reference), and computing costs (by taking the median running time of 5 consecutive runs).
The numerical simulation (Fig. S1a,b in the supporting information) demonstrated that all methods yielded satisfactory image quality, with minimal streaking in the sagittal view. The phantom experiment (Fig. S2a-c in the supporting information) demonstrated that constraining solutions' energies at the cone region alone was not sufficient to suppress streaking (TSVD, TKD). The iSWIM method reduced streaking by iterative filtering, and the most significant improvement was observed when spatial constraints were incorporated during dipole inversion (MEDI, CSC, HEIDI, TVSB). The in vivo brain MRI (Fig. 7) demonstrated that all methods successfully generated QSMs. The major structures containing high levels of ferritin (basal ganglia and nuclei in deep gray matter) or deoxyhemoglobin (veins) are shown with high paramagnetic values on QSM, and WM tracts are shown with diamagnetic values on QSM. Rapid streaking signal variations not on COSMOS were observed in TSVD and TKD, likely artifacts originating from veins. The iSWIM method reduced but did not eliminate these artifacts. TVSB overblurred compared to COSMOS; one of the causes may be its lack of noise whitening. MEDI, CSC, and HEIDI yielded QSM images similar to COSMOS.
In the accuracy assessment (Table 1), MEDI and TVSB achieved the best slope and coefficient of determination (R 2 ) in both numerical simulation and phantom experiments, presumably because the piece-wise constant nature of the susceptibility distribution matched the assumptions in MEDI and TVSB. In the human brain, although MEDI provided the highest slope, none of the methods provided adequate R 2 values. A possible cause may be voxels of WM with susceptibility anisotropy. The reconstruction speed of the k-space methods were much faster than that of the iterative Bayesian methods (Table 1). However, the split Bregman implementation in TVSB showed promise of fast online reconstruction for the Bayesian methods.

Clinical Applications Under Development
QSM has become sufficiently accurate for measuring the strong susceptibilities of biomaterials, including iron distribution (ferritin), in the deep brain nuclei and basal ganglia; deoxyhemoglobin in the veins; blood degradation products (hemosiderin in late stage); calcification (hydroxylapatite crystals); and exogenous species such as gadolinium. Clinical applications of QSM are being developed to probe neurodegenerative and inflammatory diseases, to assess hemorrhage, to measure metabolic consumption of oxygen, and to guide and monitor therapy. QSM can also remove blooming artifacts in traditional T2*w imaging (95), providing an accurate definition of the distribution of magnetic biomaterials in MRI. In this brief survey, we focus on neurological applications, although applications outside the brain are also promising (96).

Diamagnetic Biomaterial-Based Applications
QSM can easily differentiate diamagnetic calcification from paramagnetic materials such as hemosiderin (89,97). Both calcification and chronic hemorrhage appear hypointense on GRE magnitude images and may be undetectable on conventional T1-and T2-weighted SE imaging (98,99). GRE phase imaging has long been recognized for its ability to identify diamagnetic calcifications (100,101), but there is no study demonstrating its diagnostic accuracy. As such, CT is widely used for detecting calcifications despite its use of ionizing radiation (102). QSM may replace CT for detecting calcification in neurocysticercosis (Fig. 8) and tumors (103)(104)(105)(106). A clinical study demonstrated that QSM is superior to phase imaging and has a very high sensitivity (90%) and specificity (95%) for the detection of intracranial hemorrhage and calcification (97). QSM can also be used to measure the loss of myelin (107), another important diamagnetic biomaterial.
Susceptibility values in QSM can be converted to the venous deoxyhemoglobin concentration ½dHb according to deoxyhemoglobin's molar susceptibility x dHb ¼10765 ppb (48,(118)(119)(120), allowing quantitative fMRI (117). The tissue metabolic rate of oxygen consumption (MRO 2 ) is regarded as a fundamental biomarker for assessing viability of aerobic tissue such as those in the brain, heart, and kidney (121). Measurements of regional blood flow f by quantitative perfusion, such as the arterial spin labeling (122) and dHb ½ by QSM, can be used to map MRO 2 according to oxygen mass conservation (123)(124)(125)(126): For tissues with nonheme iron such as ferritin, susceptibility contributions can be corrected using iso-metabolism manipulation (127).

Paramagnetic Contrast Agent Biodistribution Quantification-Based Application
QSM can be applied to measure the biodistribution of highly paramagnetic contrast agents (CA), providing an effective tool for quantifying CA concentration [CA] in MRI (188,189). The quantitative study of the phase change observed in contrast-enhanced MRA (190) was an initial motivation to formulate the field-to-susceptibility inverse problem (31). Absolute determination of [CA] according to T1/T2 enhancement effects requires calibration and is highly susceptible to flip angle errors (191)(192)(193). CAs with limited access to water demonstrate the well-known T1/T2 relaxation quench (191,192,(194)(195)(196)(197); relaxation enhancement requires CA binding with bound water, which, in turn, exchanges chemically with surrounding bulk water (CA$bound H 2 O$bulk H 2 O) (36,(198)(199)(200) (see Supporting Figure S4).
[CA] has no well-defined relationship with R2*. QSM overcomes these problems of mapping [CA] in T1/T2/T2* imaging and may be useful for determining the biodistribution of targeted CAs in molecular MRI. A high temporal-spatial gadolinium concentration ½Gd map can be obtained using QSM and fast imaging (201), from which quantitative perfusion map may be generated (50,202) (Fig. 11).

Mixed Diamagnetic and Paramagnetic Applications
GRE phase images have been used to study iron distribution and demyelination in multiple sclerosis (MS) lesions (203)(204)(205)(206)(207). Iron distribution has been reported to be abnormally high in both the basal ganglia and lesions in MS patients (137,203,204,208,209) and may vary with lesion age and inflammatory status (204,210,211). QSM can be used to measure susceptibility changes in both lesions and nonlesion tissues in MS brains (212,213). A recent QSM study of magnetic susceptibilities of MS lesions (214) demonstrates that MS lesion susceptibilities start at the level of normal appearing WM (NAWM) (age ¼ 0y, acute enhancing), quickly increase (within 0.5y) above that of NAWM, remain almost constant for a period (0-4y, intermediately aged and nonenhancing), and then decay gradually back to that of NAWM (> 7y, FIG. 11. Quantitative susceptibility mapping for quantifying paramagnetic contrast agents. In an in vivo dynamic gadolinium (Gd) enhancement study of the brain, time-resolved Quantitative susceptibility mapping (QSM) was developed using spiral readout and temporal resolution acceleration with constrained evolution reconstruction (TRACER) complex image reconstruction. The difference image divided Gd molar susceptibility generates time-resolved Gd concentration map. Source: Xu et al, MRM 2014, epub.
FIG. 12. Quantitative susceptibility mapping for quantifying a mixture of paramagnetic and diamagnetic biomaterials. (i) Acute enhancing lesions in a 32-year-old male with RRMS at two time points: initial study (T1w þ c1 and QSM1, 1st row) and 3-month follow-up (T1w þ c2 and QSM2, 2nd row) (T1w þ c ¼ T1-weighted imaging with Gd). Lesions appear in the right frontal WM (white arrows) and in the lcc (black arrows). Both lesions are enhancing (arrows) on T1w þ c1 and isointense (white boxes) on QSM1. The lesions changed on follow up to nonenhancing on T1w þ c2 and hyperintense on QSM2 (arrows). The right frontal WM matter lesion shrunk between QSM1 and QSM2. The lcc lesion (black arrow) recovered to normal appearing on T2w and T1w (not shown), T1w þ c. (ii) Nonenhancing lesions (33y, f, RRMS) on QSM at two time points (2nd row was 6 months later). All QSM lesions at time point 1 were estimated to be 1.2y using prior MRIs. All lesions (arrows) are QSM hyperintense on both time points with similar values. (iii) Chronic nonenhancing lesions (50y f RRMS) on QSM and T2W. Two lesions over 10 years old were detected (white arrows). They appear isointense on both QSM (white boxes, only initial study shown). lcc, left cingulate cortex; QSM, quantitative susceptibility mapping; RRMS, relapsing-remitting multiple sclerosis; WM, white matter. Source: Chen et al, Radiology 2014;271: 183-192. chronic nonenhancing) (Figs.12 and 13). This MS lesion susceptibility time course is consistent with the nophase-variation on the 2.5y follow-up of nonenhancing MS lesions seen in another study (215), and with the rapid increase of off-resonance frequency observed in acute enhancing lesions in a third study (216). Investigations are ongoing to connect susceptibility time course and MS cellular activities. QSM may constitute an important new biomarker for the inflammatory and neurodegenerative activities in MS.
Initial Results in Aorta, Breast, Extremity, and Kidney QSM applications beyond the brain are also under active development (Fig. 14). The susceptibility values from phase data of the aortic arch during a Gd bolus passage may provide quantitative contrast-enhanced MRA (31) (Supp. Fig. S5). QSM is feasible for applications in other body parts including the breast, extremity, and abdomen (liver and kidney) for studying hemorrhage, metabolic oxygen consumption, mineral distribution, and contrast agent kinetics (96).

QSM of Tissue Complexity: Multiple Species and Microstructures
QSM techniques have started to proliferate, an indication that QSM is a fertile field for innovation. This review so far has focused on modeling a voxel of tissue with a scalar susceptibility. Here we briefly survey investigations to model MRI signal with tissue complexity: multiple species of different chemical shifts, subvoxel structures, and molecular structures.
Nonlinear Phase Behavior of Multiple Spectral Species, Long TE, Large Voxel The signal model in Eq. [5] may be regarded as a first order (linear temporal phase evolution) approximation, which may be good enough for many brain applications. For imaging other body parts, there may be signal contribution from proton sources other than free water w r ð Þ, such as fat f r ð Þ with chemical shift (characterized by a constant frequency offset $ $ À3.4 ppm). Eq. [14] can be generalized to account for the chemical shift effects on signal phase (217,218). The spatial smoothness of the tissue magnetic field can be used for fat-water separation (219-221):  Note that the signal phase from a voxel is now nonlinear in its temporal evolution. Initial results in solving Eq. [16] are very encouraging, promising to extend QSM to body parts with fat (Fig. 14). The approximation in Eq. [5] works very well in most imaging situations but may break down in the presence of unusually strong susceptibility sources, long TEs, large voxels, or a combination of these factors. We may need to include higherorder terms in the evaluation of the exponential and include contributions from the neighboring voxels using an accurate voxel sensitivity function (44,222). These complications lead to a voxel signal phase that varies nonlinearly with TE.

Signal Behavior with Subvoxel Structure
There is growing interest in modeling tissue microstructure using MRI (215,(223)(224)(225)(226)(227)(228)(229)(230)(231). Subvoxel structures may be characterized as gradients and higher-order spatial derivatives in spin density and magnetic field (232). These violations of the smoothness assumption in digitizing Eq. [2] result in voxel signal phases with nonlinearly temporal evolutions. More useful models may include specific geometries for the underlying tissue microstructure such as solid or hollow cylinders for capillaries, fibers, and other linear microstructure (Supp. Fig. S6), also leading to phase nonlinear in time (225,227,233). Microstructures may be considered as static and observer water as undergoing rapid random motion. The ergodic hypothesis may be assumed: the sum over the observer water path becomes the sum over the ensemble distribution that is proportional to spin density. Then, voxel signal may be modeled as the sum of the contributions from water protons inside magnetic microstructures or compartments (V c ) within the voxel. An example compartment is the cylinder or generalized Lorentz model (225). When water exchange among compartments is small, the signal model is a simple extension of Eq. [5], s t ð Þ $ X c X r2Vc mðrÞe ÀibcðrÞv0t [17] Here, for a given subvoxel compartment model, Maxwell's Equations can be used to determine the field's dependence on subvoxel structures (b c ðrÞ) such as their orientations and underlying molecular susceptibility anisotropies (225,227,233). With a sufficient number of measurements, the compartmental susceptibility may be estimated from the MRI signals: mðrÞe ÀibcðrÞv0TEj jj 2 2 þ Rð x c f gÞ [18] Susceptibility Tensor The diamagnetic susceptibilities of anisotropic molecules (Supp. Fig. S6) must be described by recognizing the susceptibility in Eq. [2] as a tensor. If all types of anisotropic molecules are sufficiently smoothly distributed in the space, and the spatial dispersion of phase accruals is sufficiently small in a voxel-as assumed in Eq. [5]then the corresponding digital form of Eq. [2] with tensor susceptibility can be used, forming the foundation for susceptibility tensor imaging (STI) (234). Group symmetry theory suggests that susceptibility anisotropy can only be observed in a voxel if and only if anisotropic molecules are arranged orderly on a macroscopic scale (235,236). The increased number of variables in STI requires acquisitions at many orientations (237), which may be reduced by using prior information obtained from diffusion tensor imaging (235,236). Similar to Eq. [17], subvoxel structures may be incorporated into the MRI signal equation, introducing phase nonlinear in time and other complexities (233). The most interesting biomaterial demonstrating susceptibility anisotropy may be myelin (107,238), and the assessment of myelin integrity using MRI remains an important unmet clinical need.

CONCLUSION
Magnetic susceptibility directly reflects the molecular electron cloud behavior in the main magnetic field. Tissue susceptibility effects can be readily sensitized in MRI, for example using the widely available GRE sequence. Maxwell's equations and the MRI signal equation can be used to quantitatively model the relationship between MRI signal and tissue susceptibility. Regularization is necessary to obtain a unique solution for determining the tissue susceptibility map from the acquired MRI signal, which is an ill-posed problem due to the lack of MRI signal in the background and the zeroes in the dipole kernel. The current status of QSM is very encouraging. The first order solution of QSM can be robustly obtained using physically meaningful regularizations, including the Bayesian approach. QSM has promising clinical and scientific applications that involve large susceptibility changes by hemoglobin, ferritin, calcification, and contrast agents. The investigations of higher order solutions have also been initiated, including studies of important magnetic anisotropies and tissue microstructures.