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

1Magnetic susceptibility is one of the following major categories of tissue contrast mechanisms in proton MRI: 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. 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 and the echo time (. This phase may be used to further attenuate the signal for enhancing T2* image contrast, which is called susceptibility weighted imaging –. 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. 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).

2QSM was contemplated at the early days of MRI. However, the inversion from field to susceptibility is ill-posed,: There are zeroes in the kernel connecting the susceptibility distribution and the field, and a simple kernel division causes large errors that present as streaking artifacts in the reconstructed susceptibility map,. Regularization or conditioning is necessary to select a unique solution for a given field –. Fortunately, MRI provides plenty of information on tissue anatomical structures. This information can serve as a prior in Bayesian regularization to overcome this ill-posed inverse problem, generating a reasonably accurate susceptibility map –. Various regularizations have since been developed,–, making QSM a feasible tool for the MRI community.

3As 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.

4In the MRI main field, a molecule in tissue gains a magnetic moment through its electron cloud polarization. Correspondingly, a tissue with volume magnetic susceptibility gains magnetization (see TISSUE 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.

5According to Maxwell’s equations in vacuum, a magnetic dipole moment at a source location generates a magnetic field at an observation location ( is a unit vector along),

6Here, the inverse-cube of the distance term characterizes the spatial extent of the dipole field. The delta-function term can be understood from the field of a current loop with a fixed magnetic moment and a radius (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) observed by a water proton is the sum of contributions from all surrounding susceptibility sources [their distribution defined by magnetization], excluding that from the proton’s own location:

7The 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 that has a uniform interior field from the sum of the fields of all sources according to Eq. 1. In Eq. 2,. (Here [expression] – 1 if expression is true and 0 otherwise.)

8Eq. 2 relates the field at to the magnetization distribution over the whole space. This can be expressed in a differential form that relates the field at to the magnetization located at:

9Eq. 3 can be Fourier transformed into. For notational convenience, we use lower case for-space quantities and upper case for corresponding-space quantities (except constant), the z-component of is referred to as the dipole kernel with Fourier transform, and the z-component of is noted as. For scalar susceptibility, Eq. 2 becomes with Fourier transform.

10On 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 induces a shielding magnetic field (Fig. 2), or chemical shift (referenced to water), that alters the field experienced by protons within the molecule: 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 –.

11The 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). These complicated effects are characterized as relaxation enhancement –. Susceptibility estimation from MRI signal magnitude affected by relaxation is prone to large errors. Fortunately, the phase of a water proton spin linearly increases with the field. Using multiple radio frequency (RF) coils, with theth coil element having a complex coil sensitivity function and acquisition noise, the k-space signal measured in coil at time is:

12Here, and is the proton transverse magnetization [is proportional to proton density and is much smaller than the electronic magnetization]. Typically, is homogeneous within a voxel and is much longer than the readout duration, and the data sampling gradient is much larger than. Then, the complex image can be reconstructed using the Moore-Penrose pseudo-inverse for both full -space sampling and parallel imaging. Approximating the voxel sensitivity function as a box function, the signal detected at in a voxel centered at with size (volume) is the summation over signals from all proton spins in the voxel,

13Here, the integration is evaluated by assuming small phase inhomogeneity in a voxel: The quantities in Eq. 5 are defined as: is 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 times the proton density-weighted average of the field in the voxel, and the T2* decay rate differs from by 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.

14The three-dimensional (3D) GRE sequence is currently the most frequently used to acquire data for QSM. Multiple echoes (number of echoes – 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 mm3) 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, or zero TE 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. Resolution should be as high as possible to minimize averaging effects in Eq. 5. Parallel imaging and spiral sampling trajectories, 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), 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.

15The 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: as the voxel center, as the field value at, and as the average of the susceptibilities in the source voxel positioned at. This digitization can be analytically derived for spherically shaped voxels: The field’s spherical mean value equals the field value at the sphere center; 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 sinc-interpolation or zero-padding may be used to reduce the digitization error. Accordingly, Eq. 2 will be used to model the MRI signal with voxels indexed by.

16Eq. 5 models the MRI signal as a product of an exponential factor that describes a phase linear in time and a complex amplitude that contains a constant phase and an amplitude decay with time. Then, the total field can be estimated from MRI signals at multiple TEs as a nonlinear least-squares problem per voxel,,:

17In the case of large SNR, Eq. 6 can be approximated as a linear regression of the signal phases at various, allowing for a closed form solution. In general, Eq. 6 can be solved iteratively and efficiently using a gradient search method, such as the Gauss-Newton or Levenberg-Marquardt method, with initialized as the linear rate of phase evolution between the first two echoes. Eq. 6 may also be solved rapidly using autoregression. Interestingly, although Eq. 5 requires the estimation of complex coil sensitivity from calibration data,, the field can be estimated without full knowledge of complex coil sensitivities from the phase difference at two TEs – or even at a single echo.

18The field estimated from Eq. 6 can contain many artificial jumps from voxel to voxel because is periodic in and can only define up to a period such as. Consequently, outside is aliased or wrapped into, resulting in abrupt artificial jumps of ( is an integer) in the output phase, which is called the wrapped phase denoted by. To estimate, these phase wraps have to be compensated by adding as needed:

19is 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, typically via path finding (including noncontinuous path, residual balancing, or branch cut) 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,; 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. Phase unwrapping may be avoided for field data estimated from direct water saturation lineshape. A frequently used method is the path-finding type, image-quality-guided region growth, which works robustly and rapidly on MRI data using SNR as the image-quality guidance,,.

20After 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. 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.

21The 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).

22The 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. Magnetic susceptibility sources exist outside, 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 field-to-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 generate the tissue field. Background susceptibility sources generate the background field for. Together, they form a total field in tissue. is tissue susceptibility when.

23The second fundamental challenge is the zero cone surface in the dipole kernel in-space or in-space: with respect to the main magnetic field (magic angle) (Fig. 4). The dipole kernel in-space is nearly flat (with all derivatives vanishing) at the ends of the cone neighborhood, implying that the standard analysis cannot be used to define in. Noise in the data would allow many possible susceptibility solutions that differ by an arbitrary amount in for a given field map, causing a substantial dipole kernel null space. Therefore, the magnetic field-to-susceptibility inverse problem is ill-posed and lacks a unique solution.

24To 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 high-pass 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,,, 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 or the harmonic property in the sophisticated harmonic artifact reduction on phase data (SHARP) method. These priors are systematically organized in the following as solutions to Maxwell’s Equations (Eq. 3).

25By definition, the background field has no source in, whereas the local tissue field has sources inside (Eq. 3). From electrodynamics or partial differential equation, for a finite domain, a unique solution to Laplace’s equation Eq. 3 can be obtained according to the values at the boundary,. In MRI, the tissue field is typically weak. At, the background field may be approximated as the total field. This is the Laplacian boundary value (LBV) method (Fig. 5) to estimate background field.

26Eq. 8 can be solved numerically as a system of linear equations by expressing the Laplacian with a difference operator:. Here, denotes the nearest neighbor average operator. It is advantageous in SNR to express Eq. 8 using the spherical mean value (SMV) operator,,:. The system matrix is diagonally dominant (with diagonal elements), and the Jacobian iterative method can be used to give the (n+1)th iteration solution. With an initialization that satisfies the boundary value condition in Eq. 8, the background field can be obtained numerically by repeatedly applying, which is the iterative spherical mean value (iSMV) method. 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.

27When 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),

28Here, requires is smaller by a border layer than the available tissue region for computing. Eq. 9 can be solved as:, using truncated singular value decomposition (TSVD) (with a carefully chosen truncation value) in-space, and, using for possible denoising benefit. This is the sophisticated harmonic artifact reduction for phase data (SHARP). in the border layer may be recovered using a continuity assumption as in harmonic (background) phase removal using the Laplacian operator [HARPERELLA ]. Another regularization is the minimal norm of the tissue field (implicitly modifying inverse Laplacian),:

29Here, for any field (squared L2 norm). This method is called regularization-enabled SHARP (RESHARP). Eq. 2 can be expressed in a Lagrangian form.

30There is another approach to estimate the background field based on the equivalent source or charge simulation method. The background field can be represented by fields of dipoles outside that are approximately orthogonal to the fields of dipoles inside except near. Then the background field can be estimated by all possible PDFs,,:

31Here, noise-weighting is the phase SNR, which is assumed to be large by linearizing the signal model Eq. 5. Note that in tissue is unique, although in backgroundis not. This PDF method may provide a slight advantage in dealing with noise in by avoiding the inversion of the Laplacian in Eq. 8 through Eq. 10, and by extending to the nonlinear form, as in Eq. 6. However, the orthogonality between fields of dipoles breaks down near the boundary, making PDF prone to over-fitting errors near the border.

32Because 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 in Eq. 8, making RESHARP similar to LBV. The minimization in Eq. 10 is similar to that in Eq. 11, without noise weighting, 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.

33To 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 [and we can extend to the general SNR case using the complex signal as in Eq. 6 ]. Then Eq. 5, after background field removal, is reduced to a linear field-to-susceptibility inverse problem (Eq. 2 with noise),

34Arbitrary susceptibility values in the dipole zero cone neighborhood are allowed in Eq. 12. A regularization can be used to specify susceptibility values in. 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,. This method is known as the calculation of susceptibility using multiple orientation sampling (COSMOS), which delivers an exact reconstruction by fully sampling the susceptibility,. Unfortunately, acquiring multiple scans of patients in different orientations is not clinically acceptable. We should focus on the regularization approach to QSM using single orientation data.

35Regularizations 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-space approach, exemplified by TSVD based on matrix computation, and 2) the Bayesian approach based on optimization, also known as the-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,,,,,. Additionally, noise will always propagate into the final solution.

36Closed-form solutions form a class of noniterative-space methods. One example is TSVD: with, which is similar to a Tikhonov-regularized minimal norm (MN) solution,. More commonly used in QSM is a TSVD variant called truncated -space division (TKD):

37Truncation obviously leads to an underestimation of the susceptibility values in, and consequently produces streaking artifacts along the magic angle in the susceptibility map. This error in the TKD method is: the first part is the regularization error from points in where the kernel has been modified by truncation; the second part is the noise error from all data points in k-space, but interestingly is also dominated by points in and near. The underestimation in Eq. 13 may be compensated by a scale factor. The streaking in Eq. 13 may be reduced using a high threshold [, susceptibility – field convolving with a kernel ] or using iterative filtering (defined by in-space) of signals outside high-susceptibility structures in-space defined by a mask [iterative Susceptibility Weigted Imaging and susceptibility Mapping (iSWIM), ]. A variant of MN is a Tikhonov-regularized minimal gradient norm [:].

38Bayesian regularizations form a class of iterative optimization methods. To allow optimal treatment of regularization error, prior information is regarded as a probability distribution function (pdf) typically expressed as: is a tunable regularization parameter and 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. For Gaussian noise with pdf ( as in Eq. 11) in the estimated field, the posterior probability is, the maximization of which is the MAP solution:

39The 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,, 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 –.

40There are connections between Bayesian optimization and noniterative-space methods. If in Eq. 14, then L2 norm-based Tikhonov regularization leads to and leads to; both are noniterative-space methods. Because MRI phase noise requires spatially varying weighting, noniterative-space methods may suffer from noise errors. A wide range of forms for in Eq. 14 have been reported, including piece-wise constant susceptibility,, smooth susceptibility or susceptibility gradient, sparse susceptibility gradient or wavelet,,, and morphological consistency of the susceptibility map,,,, 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 are calculated firstly (and properly) according to Krylov sequence, and later iterations fill structured noise in that cause streaking artifacts,, (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-space center,. For a solver of Eq. 14 that includes use of LSQR, such as in the Gauss-Newton method, its final solution may have. Eq. 14 may only determine the susceptibility up to a constant [similar to the Neumann boundary value problem for Laplace’s equation ]. 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.

41Available anatomic information in a specific imaging situation defines a physical prior for morphology-enabled dipole inversion (MEDI),,. 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. This minimizes streaking artifacts common to the dipole kernel nullspace. From Eq. 14, one MEDI implementation to reconstruct QSM can be formulated as,: Here, is evaluated in Fourier space by,for anatomic image and threshold ( is chosen such that approximately 30% of voxels are labeled as edges), and(L1 norm). Because the dipole kernel (Fig. 4), similar to a HPF, preserves susceptibility edges in the tissue field image and accordingly in the T2*-weighted (T2*w) magnitude image, the GRE magnitude and phase images can be used as the anatomic images. The nonlinear Eq. 15 can be solved using the Gauss-Newton method without explicit formation of the memory costly inverse Hessian matrix. 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.

42New 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. Multicenter trials are needed to establish a consensus on application-specific priors, which leads to the topic of the next section.

43To 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; TKD; iSWIM; MEDI; compressed sensing compensated (CSC) inversion; homogeneity-enabled incremental dipole inversion (HEIDI); and total variation using split Bregman (TVSB). 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 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.

44Rigorous 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,. 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).

45The 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.

46In the accuracy assessment (Table 1), MEDI and TVSB achieved the best slope and coefficient of determination (R2) 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 R2 values. A possible cause may be voxels of WM with susceptibility anisotropy. The reconstruction speed of the-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.

47QSM 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, 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.

48QSM can easily differentiate diamagnetic calcification from paramagnetic materials such as hemosiderin,. Both calcification and chronic hemorrhage appear hypointense on GRE magnitude images and may be undetectable on conventional T1- and T2-weighted SE imaging,. GRE phase imaging has long been recognized for its ability to identify diamagnetic calcifications,, but there is no study demonstrating its diagnostic accuracy. As such, CT is widely used for detecting calcifications despite its use of ionizing radiation. QSM may replace CT for detecting calcification in neurocysticercosis (Fig. 8) and tumors –. 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. QSM can also be used to measure the loss of myelin, another important diamagnetic biomaterial.

49During blood degradation in hemorrhage (Fig. 1b), susceptibility progressively increases from oxyhemoglobin (diamagnetic) to deoxyhemoglobin (paramagnetic with 4 unpaired electrons in Fe2+), methemoglobin (strongly paramagnetic with 5 unpaired electrons in Fe3+), and hemosiderin (super paramagnetic with possible magnetic domain formation or ferromagnetic) –. GRE is more sensitive than CT at detecting intracerebral hemorrhage,. However, the T2* hypointensity in GRE suffers from blooming artifacts that are highly dependent on imaging parameters. Reliably estimating the hematoma volume is critical for managing hemorrhagic stroke patients,, but it is difficult to do on GRE. QSM removes these blooming artifacts,, and can be used as a universal measurement of microbleeds (Fig. 9) and hematoma volume in GRE MRI,,.

50Susceptibility values in QSM can be converted to the venous deoxyhemoglobin concentration according to deoxyhemoglobin’s molar susceptibility–10765 ppb,–, allowing quantitative fMRI. The tissue metabolic rate of oxygen consumption is regarded as a fundamental biomarker for assessing viability of aerobic tissue such as those in the brain, heart, and kidney. Measurements of regional blood flow by quantitative perfusion, such as the arterial spin labeling and by QSM, can be used to map according to oxygen mass conservation –:. For tissues with nonheme iron such as ferritin, susceptibility contributions can be corrected using iso-metabolism manipulation.

51Iron overload can generate biologically toxic reactive oxygen species, causing oxidative stress and damaging macromolecules including proteins, lipids, and DNA,. Excessive iron deposition in specific regions of the brain has been observed in many neurodegenerative diseases –, including Parkinson’s disease, Alzheimer’s disease, amyotrophic lateral sclerosis, Huntington’s disease, Friedreich’s ataxia, and multiple sclerosis –. Consider the example of Parkinson’s disease (PD). Several studies demonstrate an increase of iron deposits in the substantia nigra in PD –, perhaps increasing with disease progression. The abnormal increase in nigral iron generates reactive oxygen species, possibly causing nigrostriatal dopamine neuron degeneration, and alpha-synuclein aggregation,. MRI R2* (1/T2*) has been used to measure tissue iron content,–, but R2* in a voxel reflects the field variance within that voxel (see definition after Eq. 5). The latter depends on the background field, surrounding tissue susceptibilities, and imaging parameters including field strength, voxel size, and TE. The MRI phase has also been used to measure brain iron in PD,–, but phase is a convolution of tissue susceptibility in space (Eq. 2). QSM overcomes the problems of R2* and phase images, enabling reliable iron quantification when there is no other substantial susceptibility contributor,,,,–. Perhaps for this reason, QSM has been widely used to measure brain susceptibility,.

52QSM can improve visualization of the target in deep brain stimulation (DBS), the surgical implementation of stimulating electrodes in the subthalamic nuclei (STN), or the globus pallidus pars interna (GPi) for treating neurological disorders including PD –. The anatomical accuracy of electrode lead placment is critical for a successful surgical outcome –. Intraoperative CT, typically used to guide DBS, has poor contrast for the basal ganglia structures,,–. MRI provides better tissue contrast than CT, but the visualization contrast for STN and GPi is still poor in standard T2-weighted imaging, and high field MRI using T2*w imaging has been sought for DBS –. QSM can be used to remove the blurring present in T2*w and improve contrast-to-noise-ratio by for the visualization of the STN,, (Fig. 0).

53QSM can be applied to measure the biodistribution of highly paramagnetic contrast agents (CA), providing an effective tool for quantifying CA concentration [CA] in MRI,. The quantitative study of the phase change observed in contrast-enhanced MRA was an initial motivation to formulate the field-to-susceptibility inverse problem. Absolute determination of [CA] according to T1/T2 enhancement effects requires calibration and is highly susceptible to flip angle errors –. CAs with limited access to water demonstrate the well-known T1/T2 relaxation quench,,–; relaxation enhancement requires CA binding with bound water, which, in turn, exchanges chemically with surrounding bulk water (CA↔bound H2O↔bulk H2O),– (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 map can be obtained using QSM and fast imaging, from which quantitative perfusion map may be generated, (Fig. 11).

54GRE phase images have been used to study iron distribution and demyelination in multiple sclerosis (MS) lesions –. Iron distribution has been reported to be abnormally high in both the basal ganglia and lesions in MS patients,,,, and may vary with lesion age and inflammatory status,,. QSM can be used to measure susceptibility changes in both lesions and nonlesion tissues in MS brains,. A recent QSM study of magnetic susceptibilities of MS lesions 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, chronic nonenhancing) (Figs.12 and 13). This MS lesion susceptibility time course is consistent with the no-phase-variation on the 2.5y follow-up of nonenhancing MS lesions seen in another study, and with the rapid increase of off-resonance frequency observed in acute enhancing lesions in a third study. 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.

55QSM 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 (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.

56QSM 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.

57The 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, such as fat 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,. The spatial smoothness of the tissue magnetic field can be used for fat–water separation –:

58Note 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 higher-order terms in the evaluation of the exponential and include contributions from the neighboring voxels using an accurate voxel sensitivity function,. These complications lead to a voxel signal phase that varies nonlinearly with TE.

59There is growing interest in modeling tissue microstructure using MRI,–. Subvoxel structures may be characterized as gradients and higher-order spatial derivatives in spin density and magnetic field. 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,,. 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 within the voxel. An example compartment is the cylinder or generalized Lorentz model. When water exchange among compartments is small, the signal model is a simple extension of Eq. 5,

60Here, for a given subvoxel compartment model, Maxwell’s Equations can be used to determine the field’s dependence on subvoxel structures such as their orientations and underlying molecular susceptibility anisotropies,,. With a sufficient number of measurements, the compartmental susceptibility may be estimated from the MRI signals:

61The 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). 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,. The increased number of variables in STI requires acquisitions at many orientations, which may be reduced by using prior information obtained from diffusion tensor imaging,. Similar to Eq. 17, subvoxel structures may be incorporated into the MRI signal equation, introducing phase nonlinear in time and other complexities. The most interesting biomaterial demonstrating susceptibility anisotropy may be myelin,, and the assessment of myelin integrity using MRI remains an important unmet clinical need.

62Magnetic 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.