Ph.D. candidates, graduates and scientists who are interested in understanding the mathematical and numerical foundations of clinically applicable (single-orientation) QSM techniques and challenges associated with their translation to the clinic.This course will discuss current approaches of data acquisition and reconstruction for quantitative susceptibility mapping (QSM) in the clinical setting. Upon completion of the course, participants will be able to list strategies for accelerated data acquisition, understand the numerical basics of clinically applicable algorithms, be able to explain the strengths and limitations of different solution strategies, and employ the relationship between unit dipole response, magnetic susceptibility distribution and MRI phase to implement their own inversion algorithm.

The requirement of additional measurement time is likely one of the major reasons why the translation of advanced imaging techniques is not a story of particular success. Adding another (lengthy) pulse sequence to an existing neuroimaging protocol naturally prolongs the total exam time and reduces the throughput. This could create a financial burden for the radiology department, because advanced MRI sequences are often not reimbursed directly by health insurances. It also implies that patients generally have to wait longer to receive an MRI. However, apart from these workflow-related issues and even more important, the patients’ comfort decreases over time, increasing the risk of motion artifacts and the likelihood that the patient terminates the exam prematurely. In particular, multi-orientation QSM techniques [9, 11], which rely on multiple scans and re-positioning of the patient, do not seem viable for the clinical setting.

Another matter of concern for the translation is the time required for the reconstruction of a susceptibility map after the data has been acquired. Although new technical means promise substantial room for acceleration [4], reconstruction time with current research-oriented software packages can range in the order of tens of minutes to hours. While long reconstruction times are not a matter of concern in the research setting, it may be required that susceptibility maps are reconstructed and available for inspection within a few minutes after the exam, depending on the urgency of a diagnosis or departmental workflow constraints.

In this lecture, we will take a look at recent progress toward fast data acquisition and susceptibility map reconstruction that will ultimately set the foundation for a successful translation of QSM to the clinic.

Similar to Susceptibility Weighted Imaging (SWI), QSM relies on one of the simplest pulse sequences, a gradient recalled echo (GRE) sequence. The rapid adoption of QSM in the research field (reflected by a relatively high number of published studies within a few years from its introduction) can be attributed partly to the established place of SWI (and hence the GRE sequence) in clinical routine neuroimaging protocols. Both SWI and QSM rely on the accrual of signal phase after the radio-frequency excitation pulse in the GRE sequence. Consequently, raw (unfiltered) phase images acquired with a conventional clinical SWI protocol may in principle be converted to susceptibility maps. However, due to the quantitative interpretation of the phase, for QSM usually more isotropic voxels are preferred than the highly anisotropic voxels typically used in SWI (0.5x0.5x2mm³).

One of the major downsides of GRE imaging for SWI and QSM that still hampers the clinical establishment is the relatively long acquisition time; between 8 and 15 minutes. A reason for this is that phase contrast $$$\Delta \varphi$$$ is a linear function of time, $$\Delta \varphi = -\gamma \cdot \Delta B \cdot \mathrm{TE},$$ where $$$\mathrm{TE}$$$ is the echo time, $$$\Delta B$$$ is the background-corrected field map, and $$$\gamma$$$ is the proton gyromagnetic ratio. The optimal phase contrast-to-noise ratio (CNR) is obtained when $$$\mathrm{TE}=T_2^*$$$, which translates to optimal echo times between 20 and 50ms in the human brain at 3 Tesla. However, long TEs are intrinsically associated with longer repetition times (TR) in GRE imaging, because TR>TE. Since the total acquisition time ($$$T_\mathrm{aq}$$$) depends on both TR and the total number of phase encoding steps $$$N_\mathrm p$$$, the longer TR the higher is the total scan time. Assuming a TR of 30ms and a rather low isotropic 1mm spatial resolution, which would be achieved with a 256x256x60 matrix size ($$$N_\mathrm p = 256 \cdot 60$$$), the acquisition time is approximately 8 minutes. The degree to which image resolution and spatial coverage can be sacrificed in QSM is limited because of the quantitative analysis of the field perturbations [5]. Since phase contrast at different field strengths $$$B_0$$$ is comparable when $$$B_0 \cdot \mathrm{TE}=\mathrm{const}$$$ (see above), higher field strengths generally allow choosing a lower TE and, hence, TR. Consequently, a higher field strength is generally preferable and should be used if the choice exists. Apart from this several other acceleration strategies reduce the total number of repetitions:

Partial Fourier (PF): A simple means to reduce the number of phase encoding steps is to sample k-space asymmetrically. This technique is available on most platforms. The missing parts of k-space are commonly either filled with zeros or reconstructed with sophisticated algorithms. PF reconstruction commonly relies on the assumption that the image is real-valued or phase varies only slowly, which does not hold for GRE imaging with long echo times, and results in Gibbs-like artifacts (which are, however, relatively well behaved) and slightly reduced image sharpness. PF with a k-space coverage of 75-85% has been successfully applied in several QSM studies, reducing the acquisition time to 56-72%, if applied in both phase encoding directions.

Echo-shifting techniques employ pairs of crusher gradients to shift the signal to a later repetition, creating an effective TE longer than TR. The technique has been shown to allow a 50% decrease of measurement time without substantially decreasing image SNR [13].

Accelerated parallel imaging: Similar to PF, parallel imaging techniques reduce the number of phase encoding steps by symmetrically under-sampling k-space [15]. Parallel imaging is available on all modern MRI scanners and may be combined with PF imaging. Acceleration factors of 2-3 are commonly used in routine imaging. However, reconstructed images generally suffer from reduced signal-to-noise ratio (SNR), the technique requires the acquisition of patient-specific calibration data, and, at higher acceleration factors, residual fold-over artifacts are often discernible on the images.

Multi-shot trajectories: Segmented and multi-shot trajectories acquire multiple echoes per repetition allowing a substantially more efficient traversal of k-space than conventional Cartesian trajectories. QSM at 1mm isotropic resolution in only 10 seconds [12, 8] has recently been demonstrated with a segmented 3D-EPI sequence [14]. Given that 3D-EPI sequence packages are available for most MRI platforms it is surprising that this technique has not more widely been used for QSM or SWI. However, EPI sequences are intrinsically sensitive to field inhomogeneities (often referred to as “susceptibility artifacts”), requiring post hoc distortion correction. Even more efficient acquisition has recently been demonstrated with a multi-echo stack-of-spiral readout, allowing whole-brain QSM without substantially reduced SNR compared to a conventional GRE sequence at 1mm isotropic resolution in 2.5 minutes.

WAVE-CAIPI: WAVE-CAIPI [2] allows highly accelerated parallel imaging with minimal noise amplification by employing a special cork-screw trajectory. The technique has recently been demonstrated to enable acceleration factors as high as 15, allowing QSM at an isotropic resolution of 1.1mm in 90 seconds [3]. Unfortunately, WAVE-CAIPI is not yet available on clinical systems.

While the future of QSM may lie in advanced acquisition techniques such as multi-shot EPI and WAVE-CAIPI, these techniques are not yet available in many clinical imaging centers. Translation of QSM to the clinic will have to utilize acceleration techniques widely available, such as parallel imaging and PF sampling, which allow QSM at acquisition times in the range of 5-10 minutes.

The mathematical foundation of all QSM algorithms is the (forward) relation between magnetic susceptibility distribution and magnetic field perturbation, $$$\Delta B$$$ [18]. The Fourier convolution theorem allows a straightforward evaluation of the related convolution integral to calculate the field perturbation from a known susceptibility distribution. To this end, one simply has to Fourier transform both the 3D unit dipole response and the susceptibility distribution and to point-wise multiply their Fourier spectra:

$$ \mathrm{FT}\{\Delta B\} = \mathrm{FT}\{d\} \cdot \mathrm{FT}\{\chi\}, \mathrm{(1)}$$

where $$$d$$$ is the unit dipole response, $$$\chi$$$ is the magnetic susceptibility distribution, and $$$\mathrm{FT}$$$ denotes the Fourier transform. Using Fast Fourier Transform (FFT) algorithms on current computer hardware, Eq. (1) may be evaluated within seconds even for very large matrix sizes. QSM algorithms solve for the unknown $$$\chi$$$ based on measurements of $$$\Delta B$$$, which is often referred to as the inversion from field to susceptibility. Solution strategies vary in different aspects of their numerical implementation and mathematical properties of the obtained solution. However, all algorithms ultimately rely on the Fourier-domain relation in Eq. (1).

The major challenge of solving Eq. (1) is that for certain spatial frequencies the magnitude of $$$\mathrm{FT}\{d\}$$$ is relatively low or even zero, implying that corresponding frequencies of the susceptibility distribution are attenuated in the field perturbation. Since field perturbation measurements are affected by random (white) noise, it is practically impossible to reconstruct heavily attenuated frequency components of the susceptibility distribution without additional assumptions or a priori knowledge on the (true) susceptibility distribution. This issue is also referred to as ill-conditioning of the inverse problem.

Multi-angle QSM techniques rely on multiple field maps acquired with the object (e.g. the head in a brain examination) in different orientations relative to the main magnetic field. If the orientations are chosen properly, the critical attenuation affects different spatial frequency components in each measurement, allowing the reconstruction of the true susceptibility distribution solely based on data (and without other assumptions). However, the requirement for repositioning and associated prolonged acquisition times make it unlikely that such techniques will be translated to the clinics.

Reconstruction of susceptibility maps from single GRE scans requires that the ill-conditioning is addressed. In other words, over-fitting of measurement noise needs to be prevented, because it would otherwise result in excessive amplification of noise in the susceptibility map. This may be achieved by introducing additional information (or constraints) on the susceptibility distribution, a technique also referred to as regularization.

Direct (non-iterative) solvers obtain a solution of Eq. (1) by applying a fixed set of arithmetic operations to the measured field map. The solution is obtained in a fixed amount of time and is a unique and exact solution of the posed mathematical (regularized) problem. Using matrix-vector formalism, Eq. (1) can be written in a simplified form as

$$ \mathrm F \hat b = \mathrm D \mathrm F \hat \chi, \mathrm{(2)}$$

where $$$F$$$ is the Fourier transform matrix, $$$\hat b$$$ and $$$\hat \chi$$$ are vectors concatenating voxel values of $$$\Delta B$$$ and $$$\chi$$$, respectively, and $$$\mathrm D$$$ is a diagonal matrix constructed from the elements of $$$\mathrm{FT}\{d\}$$$. A least squares solution of Eq. (2) with Tikhonov regularization is given by [1]

$$\hat \chi = \mathrm{argmin}_\chi \frac{1}{2} ||\mathrm F^\mathrm H \mathrm D \mathrm F \hat \chi - \hat b||_2^2 + \frac{\beta}{2} ||\mathrm T \hat \chi||_2^2, \mathrm{(3)}$$ where $$${}^\mathrm H$$$ denotes the conjugate transpose and $$$\mathrm T$$$ is the Tikhonov matrix, which imposes implicit constraints on the susceptibility distribution. The discrete image gradient operator in 3 dimensions is commonly used for $$$\mathrm T$$$, imposing smoothness on $$$\chi$$$. The scalar $$$\beta$$$ determines the degree of smoothness of the susceptibility map ($$$\beta=0$$$ implies no regularization).

The direct solution of Eq. (3) is given by $$$\hat \chi = \mathrm F^\mathrm H \cdot \mathrm A \cdot \mathrm F \hat b$$$, where $$$\mathrm A=(\mathrm D^\mathrm H \mathrm D + \beta \cdot \mathrm E^\mathrm H \mathrm E)^{-1} \mathrm D^\mathrm H$$$ is a diagonal matrix that can be pre-calculated [1]. Consequently, a regularized solution may be obtained by 1) Fourier transforming the field map, 2) point-wise multiplying with a pre-calculated filter function given by matrix $$$\mathrm A$$$, and 3) inverse Fourier transforming. Again, using FFT algorithms, a direct solution may be obtained within seconds on current computation hardware. An explicit construction of the Fourier transform matrix $$$\mathrm F$$$ is not required.

However, regularization means that the obtained susceptibility distribution depends on the choices for$$$\mathrm T$$$ and $$$\beta$$$. For example, a minimum-norm solution is obtained if $$$\mathrm T$$$ is the identity matrix. Regularization may also be achieved by modifying $$$\mathrm A$$$ directly, e.g. setting to zero [25] or thresholding [21] elements that would otherwise amplify noise. It has been demonstrated that even the pre-processing steps of phase unwrapping and background field correction may be incorporated into the direct solution, yielding a one-in-all direct QSM algorithm [19]: $$$\hat \chi = \mathrm F^\mathrm H \mathrm A \mathrm L^{-1} \mathrm F \mathrm M \mathrm L’ \varphi$$$, where $$$\mathrm L’$$$ is a wrap-insensitive Laplace operator [16], $$$\mathrm M$$$ is a subject-specific binary brain mask defining the region of interest for background field removal and $$$\mathrm L$$$ is the Fourier domain Laplace kernel.

Direct solvers have successfully been applied in several studies and can produce susceptibility maps with relatively high quality, in particular at ultra-high field strengths where phase noise is low [6]. Their computational efficiency, simple implementation, and numerical robustness makes them promising candidates for translation to the clinic. In particular, since direct algorithms rely on very simple arithmetic operations, their implementation seems feasible on current scanner platforms. However, the performance of direct algorithms declines when phase suffers from imaging artifacts or signal void pathologies exist, e.g. in the presence of dense hemorrhagic brain lesions [23]. A reason for this is that direct algorithms do not allow the selective exclusion (or weighting) of voxels with potentially unreliable phase values. Consequently, incorrect phase measurements are always directly converted to incorrect susceptibility values.

Iterative algorithms successively approach a solution in multiple steps. The true mathematical solution of the posed problem is obtained only in the limit of an infinite number of iteration steps. Consequently, different from direct methods, the obtained solution also depends on the number of iteration steps. However, iterative solvers provide much higher flexibility regarding the utilization of a priori information. This involves information on the location of edges in susceptibility maps[10], local smoothness [20] or sparseness in transform domains [26, 7]. Similar to direct solvers, iterative solvers rely on the fast forward computation of the field given in Eq. (1). However, while direct solvers require only a small number of Fourier transforms and arithmetic operations, iterative solvers evaluate Eq. (1) several times in each iteration step, resulting in orders of magnitude higher computation times. Nevertheless, the inclusion of spatial priors results in substantially improved susceptibility maps, with image quality comparable to multi-orientation approaches [24]. Furthermore, iterative algorithms allow the inclusion of spatial weights to exclude regions with unreliable or unavailable field information [23].

The higher quality and the ability to correctly calculate the susceptibility of signal-void lesions [17] makes iterative solvers a more likely candidate for the translation to the hospital than direct solvers. However, the computational complexity of the involved algorithms renders implementation on the scanner platform difficult.

Given that QSM relies on a conventional GRE sequence, which is widely available, the collection of data for QSM is relatively straightforward. Two primary challenges for the translation of QSM to the clinic may be identified: Technical aspects of the implementation on the scanner platform and evidence of practical significance of QSM for the clinical decision making.

The processing associated to QSM is entirely different from most other MRI techniques, because it consists of several highly sophisticated steps and complex numerical algorithms that solve large physically-motivated differential equations. It is challenging to convert existing prototype QSM software developed in research labs into software for online-reconstruction on the scanner platform, because the research software is often developed in high-level programming languages (e.g. MATLAB) and relies on other freely available software packages. While next generations of MRI scanners will provide powerful reconstruction environments, the implementation of QSM on current and older systems is challenging also due to limited computational power of the reconstruction engines and the limited capabilities of existing software libraries on these platforms. In particular, QSM relies on a highly accurate definition of brain tissue that is usually obtained using tools like FSL BET [22]. Furthermore, since QSM algorithms are under active development, the robustness and reproducibility of the various algorithms has yet to be investigated in the clinical setting. The underlying assumptions of some of the complex algorithms involved in QSM may not hold in very abnormal patient brains, e.g., when large paramagnetic hemorrhages, signal void pathologies, or tumor tissue is present close to the brain’s surface. For example, automated mask generation required for the background field correction may fail in these cases, resulting in artifacts or even a total failure of downstream processing steps. In particular, it needs to be understood if the typical (non-local) artifacts of QSM may mistakenly interpreted as pathology leading to wrong clinical decisions. While QSM has been successfully applied to study a multitude of neurological diseases, most published research focusses on the added informatyion on tissue pathology provided by QSM. These are important steps toward a better understanding of disease mechanisms and may potentially lead to new imaging biomarkers. However, for a translation to the clinic, it is important to identify applications in which QSM provides information that changes the clinical decision.