Computational MRI with Physics-based Constraints: Application to Multi-contrast and Quantitative Imaging

TL;DR

This paper introduces a physics-based modeling approach for MRI that combines explicit physical constraints with compressed sensing to improve multi-contrast and quantitative imaging, enabling better image quality and parameter recovery.

Contribution

It presents a novel framework integrating physics-based models with compressed sensing for MRI reconstruction and quantitative parameter mapping.

Findings

Enhanced image reconstruction quality from highly accelerated scans

Successful multi-contrast imaging using physics-based constraints

Accurate recovery of bio-physical parameters from MRI data

Abstract

Compressed sensing takes advantage of low-dimensional signal structure to reduce sampling requirements far below the Nyquist rate. In magnetic resonance imaging (MRI), this often takes the form of sparsity through wavelet transform, finite differences, and low rank extensions. Though powerful, these image priors are phenomenological in nature and do not account for the mechanism behind the image formation. On the other hand, MRI signal dynamics are governed by physical laws, which can be explicitly modeled and used as priors for reconstruction. {1}These explicit and implicit signal priors can be synergistically combined in an inverse problem framework to recover sharp, multi-contrast images from highly accelerated scans. Furthermore, the physics-based constraints provide a recipe for recovering quantitative, bio-physical parameters from the data. This article introduces physics-based modeling constraints in MRI and shows how they can be used in conjunction with compressed sensing for image reconstruction and quantitative imaging. We describe model-based quantitative MRI, as well as its linear subspace approximation. We also discuss approaches to selecting user-controllable scan parameters given knowledge of the physical model. We present several MRI applications that take advantage of this framework for the purpose of multi-contrast imaging and quantitative mapping.

Equations

$$M_z(t)$$

$$M_{xy}(t)$$

$$M_{xy}(\mathbf{r},t)=\rho (\mathbf{r})f_t\left(\mathbf{\theta }(\mathbf{r}),{\mathbf{u}}_{\tau }(\mathbf{r}){|}_{\tau =0}^t\right),$$

$$x_i(\mathbf{r})=M_{xy}(\mathbf{r},i{T}_s)=\rho (\mathbf{r})f_{i{T}_s}\left(\mathbf{\theta }(\mathbf{r}),{\mathbf{u}}_{\tau }(\mathbf{r}){|}_{\tau =0}^{i{T}_s}\right).$$

$$s(t)={\int }_{\mathbf{r}}{M}_{xy}(\mathbf{r},t)e^{-j2\pi \mathbf{k}(t{)}^{\top }\mathbf{r}}d\mathbf{r}+w(t),$$

$$\mathbf{y}=\mathbf{P}\mathbf{F}\mathbf{x}+\mathbf{w},$$

$$\underset{\mathbf{x}}{\mathrm{argmin}}\frac{1}{2}{\Vert \mathbf{y}-\mathbf{E}\mathbf{x}\Vert }_2^2+\lambda \mathcal{R}(\mathbf{x}),$$

$${\mathbf{E}}_i={\mathbf{P}}_i\mathbf{F}.$$

$$x_i(\mathbf{r})$$

$${\mathbf{y}}_i$$

$$\underset{\mathbf{\rho },\mathbf{\theta }}{\mathrm{minimize}}$$

$$x_i(\mathbf{r})=\rho (\mathbf{r})f_{i{T}_s}\left(\mathbf{\theta }(\mathbf{r}),\mathbf{u}(\mathbf{r})\right),$$

$$i=1,\dots ,T.$$

$$\mathbf{X}=\mathbf{\Phi }{\mathbf{\Phi }}^H\mathbf{X}.$$

$$\Vert \mathbf{X}-{\mathbf{\Phi }}_K{\mathbf{\Phi }}_K^H\mathbf{X}\Vert <ϵ,$$

$$\mathbf{\alpha }={\mathbf{\Phi }}_K^H\mathbf{x}.$$

$$\underset{\mathbf{x}}{\min }\frac{1}{2}{\Vert \mathbf{y}-\mathbf{E}\mathbf{x}\Vert }_2^2+\frac{\mu }{2}{\Vert \mathbf{x}-{\mathbf{\Phi }}_K{\mathbf{\Phi }}_K^H\mathbf{x}\Vert }_2^2+\lambda \mathcal{R}\left(\mathbf{x}\right),$$

$$\underset{\mathbf{x}}{\mathrm{minimize}}$$

$$\mathbf{x}={\mathbf{\Phi }}_K\mathbf{\alpha }.$$

$$\underset{\mathbf{\alpha }}{\min }\frac{1}{2}{\Vert \mathbf{y}-\mathbf{E}{\mathbf{\Phi }}_K\mathbf{\alpha }\Vert }_2^2+\lambda \mathcal{R}\left(\mathbf{\alpha }\right),$$

$$\mathbf{P}\mathbf{F}{\mathbf{\Phi }}_K$$

$$⟹{\mathbf{\Phi }}_K^H{\mathbf{F}}^H\mathbf{P}\mathbf{F}{\mathbf{\Phi }}_K$$

$$={\mathbf{F}}^H{\mathbf{\Psi }}_K\mathbf{F},$$

$$\underset{\rho (\mathbf{r}),\mathbf{\theta }(\mathbf{r})}{\mathrm{argmin}}\frac{1}{2}{\Vert \hat{\mathbf{x}}(\mathbf{r})-\rho (\mathbf{r})\mathbf{f}(\mathbf{\theta }(\mathbf{r}),\mathbf{u}(\mathbf{r}))\Vert }_2^2.$$

$$\underset{\rho ,\mathbf{\theta }}{\mathrm{argmin}}\frac{1}{2}{\Vert \hat{\mathbf{\alpha }}(\mathbf{r})-{\mathbf{\Phi }}_K^H\rho \mathbf{f}(\mathbf{\theta },\mathbf{u}(\mathbf{r}))\Vert }_2^2.$$

$$\mathbf{y}=\mathbf{f}\left(\mathbf{\theta },\mathbf{u}\right)+\mathbf{w}.$$

Full text

These explicit and implicit signal priors can be synergistically combined in an inverse problem framework to recover sharp, multi-contrast images from highly accelerated scans. Furthermore, the physics-based constraints provide a recipe for recovering quantitative, bio-physical parameters from the data. This article introduces physics-based modeling constraints in MRI and shows how they can be used in conjunction with compressed sensing for image reconstruction and quantitative imaging. We describe model-based quantitative MRI, as well as its linear subspace approximation. We also discuss approaches to selecting user-controllable scan parameters given knowledge of the physical model. We present several MRI applications that take advantage of this framework for the purpose of multi-contrast imaging and quantitative mapping.

I Introduction

MRI is a flexible and rich imaging modality with a broad range of applications including visualizing soft tissue contrast, capturing motion, and tracking functional and behavioral dynamics. However, long scan times remain a major limitation. A typical MRI exam consists of several scans, each several minutes long; in comparison, a full-body computed tomography exam takes a few seconds. To lower scan time, concessions are often made in the acquisition process, leading to reduced resolution, image blurring, and lower signal to noise ratio (SNR). Since data collection in MRI consists of Fourier sampling, these tradeoffs can be understood from a signal processing perspective: scan time and SNR increase with the number of Fourier measurements collected, and sampling theory dictates resolution and field of view (FOV).

A recent trend for faster scanning is to subsample the data space while leveraging computational methods for reconstruction. One established technique that has become part of routine clinical practice is parallel imaging [1], which takes advantage of multiple receive coils placed around the object to acquire data in parallel. Another avenue that is especially suited to MRI is compressed sensing (CS) [2]. In CS MRI, data are acquired below the Nyquist rate using an incoherent measurement scheme such as pseudo-random sampling. A non-linear reconstruction then leverages image structure such as sparsity and low rank to recover the image as if they were sampled at the Nyquist rate.

Though extremely powerful, common image priors in CS fundamentally leverage the natural image statistics, yet often neglect to leverage the constraints of the imaging physics that create the underlying image in the first place. Natural image statistics are phenomenological in nature and do not account for how each voxel (3D pixel) obtained its particular value. Modeling the signal dynamics can help elucidate additional structure in the image, as well as remove unwanted artifacts that are not related to under-sampling. For example, Figure 1 shows a high-resolution “gold-standard” image of a volunteer’s foot obtained with a spin-echo sequence over 12 minutes. In comparison, the middle image simulates a 48-second fast spin-echo (FSE) acquisition [3], one of the most common sequences used clinically. As the name suggests, FSE is considerably faster, but leads to tissue-dependent image blurring due to signal decay during the acquisition. Even when combined with CS, blurring due to FSE persists, because the sparsity does not address the mechanism behind the blurring. As we will show in this article, by modeling and exploiting structure in the signal dynamics as well as leveraging CS, it is possible to mitigate the blurring and obtain sharp images [4], shown by the right-most image in Figure 1.

In this article, we provide an overview of computational MRI methods that incorporate physics-based constraints and show how they can be used to accelerate acquisitions, eliminate artifacts, and extract quantitative bio-physical information from the scan. To simplify the exposition, we limit our scope to spin density and relaxation effects, as these are the most commonly used contrast mechanisms in clinical MRI.

In general, the concepts can be applied to other intrinsic tissue parameters (e.g. diffusion, spectroscopy, chemical exchange, and others), so long as an appropriate physical signal model is included. Importantly, many of these approaches have proceeded beyond the “bench-top” testing phase, and are actively used in clinical practice.

The growing interest in characterizing the underlying bio-physical tissue properties follows the larger trend toward personalized, quantitative medicine. Compared to conventional imaging, quantitative MRI can potentially better aid in identifying abnormal tissue, evaluating patients in longitudinal studies, discovering novel biomarkers, and more.

II MRI Physics

II-A Spin Dynamics

The MRI system can be approximated by a dynamical system based on the Bloch equations [5], introduced by Felix Bloch and illustrated conceptually in Figure 2. An aggregate of spins in each voxel creates a net magnetization that is initially aligned with the scanner’s main magnetic field, and evolves based on intrinsic biophysical tissue parameters and user control inputs consisting of radio-frequency (RF) pulses and magnetic field gradients. RF pulses act to rotate the magnetization vector away from the main (longitudinal) field direction and toward the transverse plane. The strength and duration of the RF pulse determine the degree of rotation, or flip angle, experienced. The transverse component of the magnetization is sensed by nearby receive coils through Faraday’s Law of Induction. We will denote the transverse magnetization as $M_{xy}$ and the longitudinal magnetization as $M_{z}$.

Although many tissue parameters influence the signal evolution, here we limit our scope to relaxation parameters, which are the most common contrast mechanism used in MRI. Relaxation is a fundamental component of nuclear magnetic resonance and dictates the rate that magnetization returns to the equilibrium state, effectively “resetting” the MRI system. Longitudinal magnetization exponentially recovers to its initial state with time constant $T_{1}$, and transverse magnetization exponentially decays to zero with time constant $T_{2}$.

As Figure 3 shows, relaxation can be understood through an analogy with a toilet111Introduced by Al Macovski in the 2009 ISMRM Lauterbur Lecture., where the water in the tank represents longitudinal magnetization, water in the bowl represents transverse magnetization, and the toilet flush is an RF excitation. When the toilet is flushed (excitation), water transfers from the tank to the bowl, producing a detectable signal. As the toilet bowl drains ($T_{2}$ relaxation), the tank refills with water ($T_{1}$ recovery). Successive flushes will transfer new water from the tank to the bowl. Neglecting other system effects, the magnetization evolution after a 90∘ RF pulse is given by

[TABLE]

where $M_{z}(t)$ and $M_{xy}(t)$ are the longitudinal and transverse magnetization components at time $t$ after the excitation, respectively, and have initial conditions of $M_{z}(0)=0$ and $M_{xy}(0)=1$ immediately following the RF pulse.

In general, the transverse magnetization distribution in space and time can be described by

[TABLE]

where $\bm{r}\in\real{3}$ represents the 3D spatial position, $\rho(\bm{r})$ is the intrinsic amount of magnetization at position $\bm{r}$ (called proton-density, or PD), and $f_{t}(\cdot)$ is a spatio-temporal signal evolution that depends on both the biophysical tissue parameters given by $\bm{\theta}(\bm{r})$ and the full history of user-controllable pulse sequence and scanner parameters222Although the pulse sequence parameters may not have a spatial component, the scanner hardware may introduce a spatial dependence, e.g. due to spatial inhomogeneities. given by $\bm{u}_{\tau}(\bm{r})$, $0\leq\tau\leq t$.

Although there are numerous tissue parameters and user controls that influence the magnetization, we will focus on the case of FSE imaging, in which the magnetization is sampled at $T$ equispaced intervals with spacing $T_{s}$. In this case, the image of interest at position $\bm{r}$ and state $i$ is the sampled transverse magnetization and is given by

[TABLE]

For FSE, the magnetization is primarily sensitive to relaxation parameters, i.e., $\bm{\theta}=\left(T_{1},T_{2}\right)$, and to RF refocusing flip angles, i.e. $\bm{u}_{iT_{s}}(\bm{r})=\mathrm{RF}_{i}(\bm{r})$, where $\mathrm{RF}_{i}$ represents the flip angle of the $i^{\mathrm{th}}$ RF pulse.

Although the RF pulse is not prescribed for each position, in practice the flip angles smoothly vary spatially due to RF transmit field inhomogeneity effects that impart varying levels of RF power across the imaging volume. We represent the vector of magnetization points at the echo times by $\bm{f}\left(\bm{\theta}(\bm{r}),\bm{u}(\bm{r})\right)\in\mathbb{C}^{T}$, where the $i^{\mathrm{th}}$ component of $\bm{f}$ is equal to $f_{iT_{s}}$ (i.e. the transverse magnetization signal).

Based on the sequence timing and RF flip angle inputs, different types of image contrasts can be created. Figure 4 shows four common FSE image contrasts primarily due to PD, $T_{1}$, and $T_{2}$, and created by using different sequence parameters in independent scans.

In particular, $\bm{f}\left(\bm{\theta}(\bm{r}),\bm{u}(\bm{r})\right)$ can be modeled by solutions to the Bloch equations [5], which are differential equations that describe the magnetization evolution as a function of time333In fact the Bloch equations themselves are also phenomenological!. The Bloch equations allow us to calculate the magnetization (signal) evolution for individual spins in a spatial region, given pulse sequence inputs. For example, the relaxation behavior given by (missing) 1 and (missing) 2 is the solution to the Bloch equations when no time-varying fields are present. Many extensions have been introduced to model additional contrast mechanisms, including diffusion and chemical exchange, though these are beyond our scope.

II-B Simulating Spins

Simulating the magnetization evolutions of spins can be a powerful tool for developing and evaluating physics constrained reconstruction methods. For example, in Section V, we will show how to use simulated signal evolutions to construct an approximate linear subspace for reconstruction. The Bloch equations are linear differential equations, and their numerical simulation can be efficiently computed under reasonable assumptions and approximations.

In particular, when the RF pulses are discretized, spin simulation using the Bloch equations reduces to successive applications of rotation followed by relaxation to the magnetization vector $\bm{M}=(M_{xy},M_{z})$. Figure 5 shows an illustration of the simulation process. The overall signal evolution is non-linear but differentiable with respect to the relaxation and system parameters.

In real imaging applications, multiple resonant frequencies often appear within each voxel due to local magnetic field inhomogeneities.

That is, the magnetic field experienced by spins within a voxel can be slightly different. For a realistic simulation, the Bloch equations have to be solved for each resonant frequency separately that appears within one particular voxel, and then summed together to form the resulting signal. This can be computationally demanding.

An alternative to the Bloch equations is the extended phase graph (EPG) formalism, which simulates the signal evolution on a voxel level [6]. EPG uniformly discretizes the resonant frequencies within one voxel. The discretization allows EPG to efficiently keep track of signal evolutions across multiple resonant frequencies using the Fourier series. Similar to the Bloch equations, spin simulation using the EPG involves successive applications of rotation followed by relaxation to the underlying magnetization representation. The overall signal evolution using EPG is also non-linear but differentiable.

III MRI Sampling

III-A Spatial Encoding

The received MRI signal represents spatial frequencies of the transverse magnetization distribution in space, and the linear relationship is described by the integral

[TABLE]

where $s(t)\in\mathbb{C}$ is the acquired signal at time $t$, $\bm{k}(t)\in\real{3}$ is a trajectory through the 3D frequency space, and $w(t)$ is complex-valued white Gaussian noise. The symbol $(\cdot)^{\top}$ denotes the transpose operation and $j=\sqrt{-1}.$ Because the spatial frequency wave-number is typically denoted as $\bm{k}$, MRI acquisitions are often described as sampling in k-space [7].

Throughout this article we consider a discrete Fourier approximation of (missing) 5. Though clinical MRI systems come standard with parallel imaging receive arrays [1], we will omit their discussion for brevity, with the understanding that they can be flexibly incorporated into the sampling model. We first consider the acquisition model when the magnetization stays constant over time, i.e. $M_{xy}(\bm{r},t)=M_{xy}(\bm{r})$ with no relaxation effects. Then, given an underlying image $\bm{x}\in\mathbb{C}^{N}$, which consists of the transverse magnetization of all $N$ voxels, the full forward model is represented in matrix form as

[TABLE]

where $\bm{y}\in\mathbb{C}^{M}$ are the acquired k-space measurements, $\bm{F}\in\mathbb{C}^{N\times N}$ is a discrete Fourier transform operator, $\bm{P}\in\mathbb{C}^{M\times N}$ is a sampling operator that selects the acquired k-space measurements, and $\bm{w}\in\mathbb{C}^{M}$ is the noise. The encoding operator is succintly represented as $\bm{E}=\bm{P}\bm{F}$.

Since the scan time is directly proportional to the number of measurements, we are typically interested in solving problems for the case where $M<N$. Compressed sensing offers an avenue for targeting this regime by exploiting low-dimensional structure in the image representation [2]. A common inverse problem approach to CS MRI is a regularized least-squares optimization given by

[TABLE]

where $\mathcal{R}$ is a sparsity-promoting regularization, e.g. $\ell_{1}$ norm of the wavelet coefficients, and $\lambda>0$ is a regularization term.

III-B Multi-dimensional Extensions

Although initial work on CS MRI focused on sparsity of static, anatomical images using spatial wavelet transforms and total variation [2], many extensions have been proposed to handle additional imaging dimensions, including joint sparsity, low rank, and their variants [8, 9, 10]. In particular, the linear forward model can be extended to represent additional image states, $\bm{x}=\begin{bmatrix}\bm{x}_{1}&\cdots&\bm{x}_{T}\end{bmatrix}$, where $\bm{x}_{i}\in\mathbb{C}^{N}$ is the image at the $i^{\mathrm{th}}$ image state and $T$ is the number of states. The forward model, illustrated in Figure 6 for signal relaxation, includes different encoding operators for each state $i$ based on the user-specified sampling patterns:

[TABLE]

This concept can be used to include additional dimensions representing signal relaxation [4], cardiac and respiratory motion [2, 11, 12], and many others [10]. Selecting the $T$ sampling patterns in a way that maintains compatibility with parallel imaging and CS is an active research area [13, 14].

IV Model-based MRI

The signal evolution model developed in the previous sections describes how the signal received by the MRI scanner is formed. Combining with spatial encoding, the physics based forward model considers the following non-linear evolution:

[TABLE]

The most explicit use of the physics model in a reconstruction is to directly solve for the tissue parameters from the raw k-space measurements. This is in contrast to first reconstructing a time series of images followed by a parameter fit. It is possible to formulate this model-based inversion even in the case of undersampled k-space [15, 16]. This can be written as a non-linear, non-convex least squares objective in which we aim to solve for ($\bm{\rho}$, $\bm{\theta}$):

[TABLE]

In the MRI literature, equation (missing) 10 is often referred to as quantitative MRI with a model-based reconstruction, as the physical model is incorporated into the objective function.

Compared to the conventional image series reconstruction followed by a fit, solving for the parameters directly serves as the ultimate dimensionality reduction: instead of computing a set of images, the aim is to recover only the intrinsic information represented by $\bm{\rho}$ and $\bm{\theta}$. The main benefit of this approach is that the problem size is significantly reduced. For example, for mapping PD and $T_{2}$, solving a problem of the form (missing) 7 consists of $TN$ unknown variables, while the model-based representation (missing) 10 has only $2N$ unknowns. In addition, sparsity-promoting regularization penalties can be applied directly to the parameter maps. Since the data fidelity term is differentiable, the overall problem can be optimized using first-order or second-order methods.

Additional system-related parameters such as RF field inhomogeneity and off-resonance can also be modeled and used as variables to be solved for as part of the reconstruction. Explicit incorporation of these systematic deviations can make the reconstruction more robust to imperfections. With the parameters in hand, synthetic contrast-weighted images could be potentially generated by evaluating (missing) 3 with specific scan and sequence parameters, a technique known as Synthetic MR [15, 17].

However, even though the individual terms are differentiable, the main downside is that the reconstruction problem is highly non-convex, thereby complicating the optimization. The non-convexity results in increased computation time and dependence on initialization. The problem requires a good initial guess of $\bm{\rho}$ and $\bm{\theta}$ in order to converge to a reasonable estimate. Another drawback is that mismatches between the model and the true acquisition can lead to error propagation in the estimated parameter maps.

For example, radiologists have described the presence of flow artifacts, white-noise artifacts, and other artifacts in synthetic FLAIR contrast images [17], likely due to unmodeled effects such as flow. Partial voluming is also known to impair the estimation accuracy [15]. Errors in sampling trajectories due to eddy currents and gradient delays can also manifest as blurring and streaking artifacts. These effects can potentially be reduced by expanding the signal model, e.g. to incorporate scanner non-idealities. The pulse sequence can also include navigator components to aid in the estimation.

V Relaxing the Multi-contrast Model

Instead of fully incorporating MRI physics into the reconstruction, it is possible to use relaxed constraints that are more amenable to optimization [18]. The natural dynamics of the MR signal are constrained by the Bloch equations, implying the existence of a low-dimensional manifold, which is non-linear in general. As Figure 7 illustrates, a linear subspace approximation to the manifold with higher dimensionality could provide a compromise between representation simplicity and size. This may be attractive for a few reasons; namely, to maintain convexity and computational efficiency in the optimization, to decouple the reconstruction from the quantitative fitting, and to reduce propagation of model error. In addition, many applications do not require quantitative parameters, and instead rely on high-quality contrast-weighted images. This is the case for nearly all clinical diagnostic imaging. Unfortunately, contrast-weighted images derived from parameter maps are susceptible to error propagation due to unmodeled components, e.g. from partial voluming and flow effects, and have seen limited use clinically [17].

V-A Subspace Constraint

Continuing with the FSE sequence as a guiding example, Figure 8 shows the signal evolutions for an FSE simulation with a particular flip angle schedule. Despite differences in relaxation parameters, the signal evolutions for different tissues follow similar trends. This correlation implies low-dimensional structure; namely, the signal evolutions of different tissues form a low-dimensional subspace [19, 20, 21, 4, 22].

Many approaches can be taken to design the subspace. Here we focus on simulating a set of training signals derived from EPG simulation. We assume a prior distribution $p(\bm{\theta})$ is known and draw $L$ samples from the distribution to create the training signals. The prior distribution can be taken from known literature values, or from a conventional mapping procedure focused on a particular anatomy.

Consider a data matrix $\bm{X}\in\mathbb{C}^{T\times L}$ consisting of an ensemble of $L$ signal evolutions sampled at $T$ echo times. Each column in $\bm{X}$ represents the signal evolution of a spin population with a particular $\bm{\theta}\in p(\bm{\theta})$. Let $\bm{\Phi}\in\mathbb{C}^{T\times T}$ be an orthonormal temporal basis, i.e.

[TABLE]

The goal is to design $\bm{\Phi}=\begin{bmatrix}\bm{\varphi}_{1}&\cdots&\bm{\varphi}_{T}\end{bmatrix}$ and a $K-$dimensional subspace, $\bm{\Phi}_{K}=\operatorname{span}\{\bm{\varphi}_{1},\dots,\bm{\varphi}_{K}\}$, such that

[TABLE]

where $\epsilon$ is a modeling error tolerance. The choice of norm in (missing) 12 will affect the chosen subspace and can be used to capture average, worst-case, and other error metrics.

When the Frobenius norm is used, the solution corresponds to the truncated singular value decomposition [19, 20, 4], i.e. principal component analysis (PCA), where $\bm{\Phi}_{K}$ consists of the left singular vectors of $\bm{X}$ corresponding to the $K$ largest singular values. The principal component images, corresponding to the principal component vectors, are given by

[TABLE]

Based on (missing) 12, we also have $\bm{x}\approx\bm{\Phi}_{K}\bm{\alpha}$.

V-B Reconstruction

The subspace relationship can be incorporated into the reconstruction as additional prior knowledge in the form of regularization [22]:

[TABLE]

where $\mu>0$ controls the degree of subspace modeling. Alternatively, the subspace can be used as an explicit constraint, i.e.,

[TABLE]

Although the analysis form described by (missing) 14 is a more faithful representation, the reconstruction still requires solving for $TN$ parameters in $\bm{x}$, in addition to introducing a hyper-parameter $\mu$. In contrast, the synthesis form used in (missing) 15 introduces explicit model error but is significantly more computationally efficient.

When using the hard constraint, we can now solve for the subspace coefficient images directly [4],

[TABLE]

followed by back-projection: $\hat{\bm{x}}=\bm{\Phi}_{K}\bm{\alpha}$. This is a significant dimensionality reduction! Instead of solving for $TN$ variables, we only need to solve for $KN$ variables. In addition, when the normal equations are used in an iterative optimization (as is the case for many first-order iterative algorithms), we can take advantage of the commutativity of the subspace operator and the Fourier transform, as the former only operates on the parametric dimension and the latter only operates on the spatial coordinates [4]:

[TABLE]

where $\bm{\Psi}_{K}\in\mathbb{C}^{KM\times KM}$ is a block-wise diagonal operator with a $K\times K$ symmetric block for each spatial frequency point. In other words, the optimization does not have to perform any computation in the ambient space, and the complexity grows with $K$, independent of $T$.

Compared to model-based quantitative mapping, the subspace-constrained forward model is convex and easier to solve. In addition, the subspace is less sensitive to model error, as it does not strictly impose a specific physical model. As an example, a voxel with partial voluming will contain a linear combination of signal evolutions, $\bm{x}(\bm{r})=a_{1}\bm{x}^{(1)}(\bm{r})+a_{2}\bm{x}^{(2)}(\bm{r})$, comprising different tissue parameters $\bm{\theta}^{(1)}$ and $\bm{\theta}^{(2)}$, which is inconsistent with (missing) 3. In contrast, if $\bm{x}^{(1)}$ and $\bm{x}^{(2)}$ are separately represented by the subspace, then so is their combination.

The main drawback to the subspace formulation is that the problem is not reduced to its intrinsic dimension governed by $\bm{\rho}$ and $\bm{\theta}$. However, in practice a subspace size of $K<5$ is practical for applications even when the ambient dimension is in the hundreds [4, 23]. A second drawback is that the subspace can represent points off the manifold, that are not physically meaningful. Thus, data inconsistencies can manifest as inaccurate images after back-projections.

V-C Quantitative Mapping

A straightforward approach to quantitative mapping from reconstructed multi-contrast images is to perform a voxel-wise non-linear least squares fit based on the signal model. Given the reconstructed image $\hat{\bm{x}}$ at voxel $\bm{r}$, this amounts to solving

[TABLE]

As with Section missing IV, it is also possible to solve for RF field inhomogeneity by additionally solving for the non-negative scalar $\eta$ that multiplies the flip angles: $\bm{u}\rightarrow\eta\bm{u}$. This formulation also covers dictionary-based methods which have grown in recent popularity with the advent of MR Fingerprinting [24]. Rather than solving a continuous non-linear least squares problem, dictionary-based fitting utilizes a grid search across the parameter space and is equivalent to matched filtering. Optimization-based methods can be used as well.

When the reconstruction strictly enforces a subspace, there is a known model error between the reconstruction $\hat{\bm{x}}$ and the signal evolution $\bm{f}$, and we do not expect them to match even in the absence of noise and under-sampling artifacts. To improve parameter estimation when using a subspace-constrained reconstruction, we can solve for the parameters directly in the subspace:

[TABLE]

Since $\bm{f}$ is differentiable with respect to $\bm{\theta}$ and $\bm{u}$, first-order and second-order solvers can be used. Figure 9 shows an example of solving (missing) 21 following a $T_{2}$ Shuffling reconstruction444Example code available on https://eecs.berkeley.edu/~mlustig/Software.html [4], where $K=3$ and $\bm{\theta}=(|\rho|,\angle\rho,T_{2})$. The non-linear least squares was solved using the Trust Region Reflective algorithm included in the Python Scipy package, and the Jacobians were calculated using the adjoint states method [25].

When solving under-determined inverse problems, it is important to recognize the bias introduced through modeling assumptions and regularization. In particular, the optimization problem (missing) 15 introduces two forms of bias: model error due to the subspace constraint, and error due to the regularization. The subspace constraint leads to a straightforward tradeoff between model error and noise amplification: noise standard deviation increases with $\sqrt{K}$, where $K$ is the subspace size [4]. A small subspace will also reduce sensitivity in the parameter space, and manifests as reduced contrast and parameter-dependent noise amplification, as all voxels are pushed to the same curve fit [20]. Bias due to regularization can also be an issue in reconstruction and parameter fitting, but its impact depends on the specific regularization used. With non-linear regularization such as $\ell_{1}$ commonly used in CS, the tradeoff can be harder to quantify. Blocking and smoothing artifacts are common when using wavelet regularization and total variation due to loss of high-frequency content [2]. Low-rank regularization and can lead to washed-out contrast due to loss of spectral information, and its multi-scale variants can lead to blocking artifacts due to their translation-variant structure [9]. In general, under-regularization can cause residual incoherent under-sampling artifacts and noise amplification, while over-regularization can lead to blurring [15, 18, 20, 21, 16, 22, 11, 4, 23, 12].

VI Choosing Scan Parameters

In addition to incorporating MRI physics in the reconstruction, using knowledge of the physics to optimize the MRI pulse sequence could potentially improve the quality of the acquired raw data. Designing an MRI pulse sequence consists of two main components: (i) designing the pulse waveform and timing to guide the signal evolution, and (ii) designing the spatial encoding to appropriately sample k-space. Initial work focused on designing sequences with simple goals in mind such as maximizing contrast and limiting blurring. A recent trend in computational MRI is to use optimal control to motivate scan parameter selection [25], and incoherent sampling properties to motivate k-space sampling [2]. Optimal experiment design has also been applied to spatial encoding under a synthesis sparsity model [14]. Importantly, both spatial encoding and sequence parameters should be designed with a specific reconstruction in mind.

VI-A Optimizing Sequence Parameters

Many approaches to optimizing sequence parameter selection involve the use of a Cramer-Rao lower bound (CRLB), which imposes a lower bound on the variance of an unbiased estimator. To simplify the optimization, the spatial encoding is decoupled from the sequence evolution, and a 1D experiment is used for optimization:

[TABLE]

Then, the Fisher information matrix, $\bm{I}(\bm{\theta};\bm{u})$, which measures the sensitivity of $\bm{\theta}$ captured by $\bm{y}$ given sequence parameters $\bm{u}$, can be used to lower bound the variance of an unbiased estimator. The sequence parameters can then be optimized to minimize this lower bound on the variance of the biophysical tissue parameter estimates. This approach has been used in [26] to develop a general framework to optimize sequence parameters from combinations of pulse sequences for precise estimation of $T_{1}$ and $T_{2}$ parameters jointly. It has also been used to determine sequence parameters for MR Fingerprinting to obtain maximal SNR while not violating physics-based MRI constraints [27].

Another optimal control design based approach uses the EPG formalism to develop a model for the signal evolution in terms of various sequence and tissue parameters [25]. An optimization problem is then formulated according to the goal, e.g. maximizing signal intensity or minimizing RF power, and solved using the adjoint states method. This approach can flexibly incorporate other objective functions, e.g. to minimize the CRLB, similar to methods discussed above. Figure 10 shows an example of optimizing the RF flip angles in an FSE experiment given RF power constraints [28]. The objective function maximized the component in the Fisher Information matrix corresponding to the variance with respect to $T_{2}$ estimation. Constant flip angles that achieve the same RF power limit are shown in comparison.

Data-driven methods for sequence parameter optimization are also emerging. With these methods, a certain set of sequence parameters are chosen, the resulting time-evolution of the signal for these parameters is either simulated or determined experimentally, and finally an image is reconstructed. The sequence parameters are then updated according to the reconstruction error. Compared to the optimal experiment based approaches, these methods have the advantage of incorporating spatial encoding and reconstruction into the evaluation.

When choosing sequence parameters, it is important to understand the impact of system imperfections on the optimized sequence. These imperfections can often be incorporated into and mitigated in the optimization framework, though they may greatly increase the computational complexity. For example, to reduce sensitivity to RF transmit field inhomogeneity, one may optimize the average (or min-max) CRLB across spins experiencing a distribution of RF transmit field inhomogeneities by simulating each spin independently. This quickly increases in complexity as additional components are considered.

VI-B Optimizing Sampling Pattern

Conventional methods to determine CS sampling trajectories are based on exploiting the incoherence conditions for CS by the use of a transform point spread function (TPSF), which determines the leakage of one transform coefficient to another transform coefficient due to subsampling [2]. A sampling pattern can be determined by a Monte-Carlo design procedure where a variable density pattern is obtained by randomly drawing indices using a probability density function, the TPSF of the obtained trajectory is calculated, and the procedure is repeated, choosing the pattern with the lowest peak interference. The reconstruction is performed according to (missing) 7 with a sparsity-promoting regularizer on the transform coefficients. Thus, incoherence in the transform domain is a good metric for sampling pattern quality.

Several data-driven approaches have been developed in order to improve on the current incoherence-based methods. These data-driven approaches use a learning-based framework to find the best sampling pattern for a set of training signals. In particular, methods that optimize a sampling pattern for a specific reconstruction rule and anatomy [14] as well as methods that jointly optimize a sampling pattern and a reconstruction strategy have been developed [29]. In this way, physics-based constraints in the reconstruction are implicitly propagated through to the sampling pattern optimization.

Although various control design based and data-driven strategies have been developed to optimize one of scan parameters, sampling trajectories and reconstructions, the ultimate goal of developing methods that optimize parameter, trajectory and reconstruction in conjunction remain an area for future work.

VII Summary

Since the advent of MRI, physics-based knowledge has been used to recover and understand the signal dynamics and tissue parameters. Many early works have leveraged physical constraints to mitigate image artifacts due to systematic errors and derive quantitative maps. However, long scan time and the lack of sophisticated reconstruction algorithms have prevented the clinical adoption of these techniques.

In the last decade, compressed sensing and other computational imaging approaches have transformed the landscape of what is possible with MRI. Scan time can be appreciably reduced by leveraging natural image statistics in the reconstruction. Many advanced reconstruction algorithms have been developed for this express purpose. Using these numerical tools and combining sparsity-based modeling, it is now feasible to run physics-constrained computational MRI methods in the clinic with reasonable scan and reconstruction times [30]. By incorporating the physical dynamics due to both tissue-specific and scanner-specific parameters, acquisitions and reconstructions can be designed in tandem to work across a broad patient population in a robust manner.

The methods introduced in this article provide a framework for modeling the MR dynamics, modifying the acquisition to account for the signal evolutions, and incorporating them into the reconstruction. Although we limited our focus to only modeling dynamics due to relaxation effects and RF pulses, there are many other contrast mechanisms and scanner controls that can be accounted for. On the tissue characterization side, these include diffusion, water-fat imaging, susceptibility, spectroscopy, pharmacokinetics, magnetization transfer, chemical exchange, and blood flow. On the imaging system side, these include field inhomogeneity, eddy current effects, gradient delays, temperature, and more. Each of these components can be modeled jointly, but will greatly increase the dimensionality of the problem. Model-based quantitative imaging reduces the inverse problem to recovering the intrinsic parameters, but careful modeling and simulation must be used to avoid artifacts due to model error. In contrast, subspace constraints and other low-dimensional representations can be used to flexibly capture the dynamics without making strong assumptions.

Many exciting modeling techniques are emerging in the signal processing community, and are great candidates at improving physical modeling in MRI. On the other hand, purely data-driven approaches have grown in popularity, in the hopes of learning the signal characteristics from real data.

With the growing trend of applying deep learning to inverse problems, there is great promise to incorporating additional physics-based constraints directly in the learning in order to restrict the feasible solution space, reduce the dependence on large training data sets, and model effects not described by the simplified physics [31]. Like compressed sensing, the rapid empirical progress in deep learning-based imaging should be suitably counter-balanced with theoretical guarantees to guide its use in clinical settings. Combining physical modeling with data-driven learning is an active area of research in the signal processing community in general, and in the MRI field in particular.

Acknowledgements

This work was supported by NIH grant R01EB009690, Sloan Research Fellowship, Bakar Fellowship, and research support from GE Healthcare. The authors thank Volkert Roeloffs, Zhiyong Zhang, and Karthik Gopalan for valuable comments.