High resolution in-vivo MR-STAT using a matrix-free and parallelized reconstruction algorithm

TL;DR

This paper introduces a novel matrix-free, parallelized reconstruction algorithm for high-resolution MR-STAT imaging, enabling accurate quantitative parameter maps from short scans without relying on FFT, demonstrated through simulations, phantoms, and in-vivo brain experiments.

Contribution

It presents a new inexact Gauss-Newton based algorithm that efficiently reconstructs high-resolution quantitative maps in MR-STAT without matrix operations, suitable for parallel computing.

Findings

Achieved 1mm resolution in simulations, phantoms, and in-vivo brain scans.

Reconstructed T1 and T2 values match standard measurements and literature.

Conventional MR Fingerprinting fails with the used pulse sequences.

Abstract

MR-STAT is a recently proposed framework that allows the reconstruction of multiple quantitative parameter maps from a single short scan by performing spatial localisation and parameter estimation on the time domain data simultaneously, without relying on the FFT. To do this at high-resolution, specialized algorithms are required to solve the underlying large-scale non-linear optimisation problem. We propose a matrix-free and parallelized inexact Gauss-Newton based reconstruction algorithm for this purpose. The proposed algorithm is implemented on a high performance computing cluster and is demonstrated to be able to generate high-resolution ($1mm \times 1mm$ in-plane resolution) quantitative parameter maps in simulation, phantom and in-vivo brain experiments. Reconstructed $T_1$ and $T_2$ values for the gel phantoms are in agreement with results from gold standard measurements and for the in-vivo experiments the quantitative values show good agreement with literature values. In all experiments short pulse sequences with robust Cartesian sampling are used for which conventional MR Fingerprinting reconstructions are shown to fail.

Tables

Table 1. Table 1: On-disk sizes of MR-STAT matrices for N p = 4 subscript 𝑁 𝑝 4 N_{p}=4 and N t = N = 4 × N v subscript 𝑁 𝑡 𝑁 4 subscript 𝑁 𝑣 N_{t}=N=4\times N_{v} for an increasing number of voxels N v subscript 𝑁 𝑣 N_{v} . The memory sizes are computed as 2 ⋅ N v ⋅ N t ⋅ 8 ⋅ 2 subscript 𝑁 𝑣 subscript 𝑁 𝑡 8 2\cdot N_{v}\cdot N_{t}\cdot 8 bytes ( 𝐌 𝐌 \mathbf{M} ), 2 ⋅ N ⋅ N t ⋅ 8 ⋅ 2 𝑁 subscript 𝑁 𝑡 8 2\cdot N\cdot N_{t}\cdot 8 bytes ( 𝐉 𝐉 \mathbf{J} ) and N 2 ⋅ 8 ⋅ superscript 𝑁 2 8 N^{2}\cdot 8 bytes ( 𝐇 𝐇 \mathbf{H} ) respectively. The factors of 2 come from the real and imaginary components and the factor of 8 represents the bytes necessary to store 64-bit floating point numbers.

Image size	Voxels $\left(N_v\right)$	$𝐌$	$𝐉$	$𝐇$
$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}64}\times 64$	$4.096$	1 GB	4 GB	2 GB
$128\times 128$	$16.384$	16 GB	64 GB	32 GB
$256\times 256$	$65.536$	256 GB	1.024 GB	512 GB
$512\times 512$	$262.144$	4.096 GB	16.384 GB	8.192 GB

Table 2. Table 2: Acquisition parameters for in-vivo MR-STAT brain scans.

Acquisition parameter	Subjects 1 and 2	Subject 3
Field strength	$1.5$ T	$1.5$ T
In-plane resolution	$1$ mm $\times$ $1$ mm	$1$ mm $\times$ $1$ mm
Field-of-view	$224$ mm $\times$ $224$ mm	$224$ mm $\times$ $224$ mm
Slice thickness	$5$ mm	$3$ mm
TR	$7.6$ ms	$7.9$ ms
TE	$3.8$ ms	$3.95$ ms
Readout bandwidth	$85.6$ kHz	$85.6$ kHz
Pulse duration	$0.76$ ms	$0.81$ ms
Scan time	$13.6$ s	$14.15$ s

Table 3. Table 3: Three different error metrics (NRMSE, HFEN, MAPE) computed for the MRF (Cartesian, Radial and Spiral) and MR-STAT (Cartesian) reconstructions on the numerical brain phantom. No noise was added to the data in these reconstructions. The MR-STAT reconstructions result in the lowest errors because the reconstructions do not suffer from undersampling artefacts and because there are no discretization errors due to a finite dictionary.

Parameter	Metric	Units	MRF Cartesian	MRF Radial	MRF Spiral	MR-STAT Cartesian
$T_1$	NRMSE	[a.u.]	0.2302	0.0432	0.0110	0.0025
	MAPE	[%]	65.6	8.5	2.5	0.4
	HFEN	[a.u.]	18.1	15.6	15.7	15.7
$T_2$	NRMSE	[a.u.]	0.2486	0.0756	0.0492	0.0048
	MAPE	[%]	56.1	8.1	4.9	0.9
	HFEN	[a.u.]	3.2	2.5	2.4	2.5
$\rho$	NRMSE	[a.u.]	1.0979	0.2626	0.1129	0.0830
	MAPE	[%]	15.8	4.8	2.7	1.8
	HFEN	[a.u.]	8.5	5.6	5.5	5.4

Table 4. Table 4: Mean T 1 subscript 𝑇 1 T_{1} and T 2 subscript 𝑇 2 T_{2} values and standard deviation in white- and gray matter regions for each of the three in-vivo brain scans.

Tissue type	Subject	$T_1$	$T_2$
Frontal white matter	1	$505\pm 48$ ms	$53.3\pm 4.0$ ms
	2	$542\pm 48$ ms	$57.4\pm 3.8$ ms
	3	$519\pm 54$ ms	$56.1\pm 4.3$ ms
Putamen (gray matter)	1	$874\pm 64$ ms	$74.8\pm 4.4$ ms
	2	$956\pm 66$ ms	$80.2\pm 4.5$ ms
	3	$895\pm 107$ ms	$78.4\pm 7.0$ ms

Equations

$$s(t)={\int }_V\rho (\mathbf{r})m(\mathbf{r},\mathbf{\theta }(\mathbf{r}),t)d\mathbf{r}.$$

$$s(t)=\sum _{j=1}^{N_v}{\rho }_j{m}_j({\mathbf{\theta }}_j,t){\Delta }_V.$$

$${\mathbf{m}}_j:=\left(m_j({\mathbf{\theta }}_j,t_1),\dots ,m_j({\mathbf{\theta }}_j,t_{N_t})\right)\in {\mathbb{C}}_t^N$$

$$\text{s}=\sum _{j=1}^{N_v}{\rho }_j{\mathbf{m}}_j$$

$$M_{i,j}:=m_j({\mathbf{\theta }}_j,t_i),$$

$$\mathbf{\rho }=\left({\rho }_1,\dots ,{\rho }_{N_v}\right),$$

$$\mathbf{s}=\mathbf{M}\mathbf{\rho }.$$

$$\mathbf{r}(\mathbf{\alpha })=\mathbf{d}-\mathbf{s}(\mathbf{\alpha })$$

$$f(\mathbf{\alpha })=\frac{1}{2}\Vert \mathbf{r}(\mathbf{\alpha }){\Vert }_2^2$$

$${\mathbf{\alpha }}^{\ast }={\text{argmin}}_{\mathbf{\alpha }}f(\mathbf{\alpha }),$$

$$\mathbf{g}=\nabla f=\Re ({\mathbf{J}}^H\mathbf{r}),$$

$$\mathbf{J}(\mathbf{\alpha }):=\left[\frac{\partial \mathbf{r}}{\partial {\alpha }_1}\dots \frac{\partial \mathbf{r}}{\partial {\alpha }_N}\right],$$

$$\mathbf{H}(\mathbf{\alpha }):={\left[\frac{{\partial }^{2}f}{\partial {\alpha }_i\partial {\alpha }_j}\right]}_{i,j=1}^{N_p}.$$

$$\mathbf{s}=\mathbf{M}\mathbf{\rho }=\sum _{j=1}^{N_v}{\rho }_j{\mathbf{m}}_j.$$

$$\text{g}=\Re ({\text{J}}^H\text{r}).$$

$$\mathbf{J}=\left[\frac{\partial \mathbf{r}}{\partial {\alpha }_1}\dots \frac{\partial \mathbf{r}}{\partial {\alpha }_N}\right],$$

$$\text{g}=\left[\begin{array}{c}\Re \left(⟨\frac{\partial \mathbf{r}}{\partial {\alpha }_1},\mathbf{r}⟩\right)\\ ⋮\\ \Re \left(⟨\frac{\partial \mathbf{r}}{\partial {\alpha }_N},\mathbf{r}⟩\right)\end{array}\right].$$

$$\frac{\partial \mathbf{r}}{\partial {\alpha }_i}=-\frac{\partial \left({\sum }_{r=1}^{N_v}{\rho }_r{\mathbf{m}}_r\right)}{\partial {\alpha }_i}=-\frac{\partial \left({\rho }_j{\mathbf{m}}_j\right)}{\partial {\alpha }_i}.$$

$$\Re \left({\mathbf{J}}^H\mathbf{J}\right)\mathbf{p}=-\mathbf{g}$$

$${\mathbf{x}}^{MRF}={\text{argmin}}_{\mathbf{x}}\Vert {\mathbf{d}}^{MRF}-{\mathbf{E}}^{MRF}\mathbf{x}{\Vert }_2^2$$

Full text

Keywords: Quantitative MRI, MR-STAT, MR Fingerprinting, Large-scale non-linear optimization, parallel computing

Accepted for publication in NMR in Biomedicine on 2019-12-05

Running title: High resolution in-vivo MR-STAT

Word count: 5777

Corresponding Author: Oscar van der Heide, University Medical Center, Heidelberglaan 100, 3508 GA Utrecht, The Netherlands. E-mail: o.vanderheide@umcutrecht.nl.

1 Center for Image Sciences, University Medical Center Utrecht, Utrecht, The Netherlands

Abbreviations

bSSFP

Balanced Steady-State Free Precession

CSF

Cerebrospinal Fluid

HFEN

High-Frequency Error Norm

L-BFGS

Limited-Memory Broyden-Fletcher-Goldfard-Shanno

FFT

Fast Fourier Transform

MAPE

Mean Absolute Percentage Error

MRF

Magnetic Resonance Fingerprinting

MR-STAT

Magnetic Resonance Spin Tomography in Time-Domain

NRMSE

Normalized Root Mean Square Error

qMRI

Quantitative Magnetic Resonance Imaging

SVD

Singular Value Decomposition

TE

Echo time

TR

Repetition time

VARPRO

Variable Projection

Introduction

Conventional magnetic resonance imaging (“MRI”) methods rely on the Fourier-Transform relationship between signal and local magnetization value for spatial encoding. Tissue differentiation is possible in the resulting qualitative images because different tissue types have distinct MR-related biophysical properties like $T_{1}$ and $T_{2}$ relaxation times. Quantitative MRI (“qMRI”) methods aim to estimate MR-related biophysical properties like $T_{1}$ and $T_{2}$ relaxation times. Quantitative images could provide additional diagnostic value and are more suited for the purpose of multi-center studies and computer-aided diagnosis [1, 2]. The most straightforward and robust choices for $T_{1}$ and $T_{2}$ mapping sequences, i.e. single echo (inversion recovery) spin echo sequences have prohibitively long scan times. Over time, a multitude of alternative pulse sequences have been developed that reduce acquisition times to clinically acceptable levels [3, 4, 5, 6]. In recent years acquisition times have been reduced even further with advanced reconstruction techniques that include more information of the underlying physical processes in the reconstructions [7], add a-priori knowledge in the form of sparsity or low-rank constraints [8] and/or allow estimation of multiple parameter maps simultaneously [9, 10]. A prime example is MR Fingerprinting (“MRF”, [11]). In MRF a transient-state pulse sequence with quasi-random components is used and many highly undersampled k-spaces are acquired in a single, short acquisition. The Fast Fourier Transform (“FFT”) is applied on each k-space to generate many snapshot images. Then, on a voxel-per-voxel basis, the measured fingerprints are matched to a precomputed Bloch-equation based dictionary to obtain the quantitative parameters. Through this novel combination of transient-state sequences with a pattern recognition step, MRF has been able to drastically reduce qMRI acquisition times.

MR-STAT [12] is a recently proposed qMRI framework that, similarly to MRF, aims to estimate multiple parameter maps from a single short scan simultaneously. However, instead of performing FFTs in a separate step for spatial localisation of signal, parameter maps are fitted directly to the measured time-domain signal using a Bloch-equation based volumetric signal model. That is, a single large-scale non-linear optimisation problem is numerically solved in which the spatial localisation and parameter estimation are performed simultaneously. In addition, instead of using a dictionary matching procedure, in MR-STAT gradient-based iterative methods are used to solve the optimisation problem. Compared to MRF the MR-STAT approach results in different trade-offs made in the reconstruction. Since the FFT is no longer explicitly used to transform back and forth between image space and frequency space, the spatial gradient encoding is entangled directly into the MR-STAT signal model. With this approach, data from different readouts of a transient state pulse sequence can be naturally combined into a single reconstruction process. There is no reliance on dictionary compression [13] or compressed sensing [14] techniques to suppress aliasing artefacts. As will be demonstrated, MR-STAT allows for the reconstruction of high quality parameter maps from very short scans even when using standard and experimentally reliable Cartesian sampling strategies.

Solving the non-linear optimisation problem that results from using the volumetric signal model in MR-STAT does introduce new computational challenges. As will be discussed, the computational and memory requirements scale quadratically with the resolution and parallelizing the computations is non-trivial because the FFT is not used to spatially decouple the unknowns. In Sbrizzi et al. [12], to alleviate the computational challenges at high resolution, a 1D FFT along the readout direction was still employed to decouple the problem in one direction in space, resulting in many smaller and independent 1D subproblems to be solved. This hybrid approach only partly benefits from the above mentioned advantages of using a volumetric signal model, e.g., dynamical behaviour during readouts cannot be taken into account. Furthermore, it can only be used with Cartesian sampling strategies. Thirdly, if the technique is applied to 3D acquisitions, each of the resulting 2D subproblems will itself be a large-scale problem. Therefore, to unlock the full potential of MR-STAT, a specialized reconstruction algorithm is required that:

does not require storage of large model matrices (i.e. is matrix-free), 2. 2.

is suitable for a parallel computing implementation to reduce computation times, 3. 3.

is extensible to non-Cartesian sampling strategies.

In the current work we present a reconstruction algorithm for MR-STAT based on an inexact Gauss-Newton method (see [15] and Algorithm 7.2 in [16]) that satisfies the above requirements. For partial derivative computations that are needed in the reconstruction we propose to use exact algorithmic differentiation. With the new reconstruction algorithm we demonstrate the potential of the MR-STAT framework through simulation studies, phantom experiments and by reconstructing high-resolution in-vivo brain maps. Although in principle the reconstruction algorithm can be used with non-Cartesian sampling, in all experiments we will use Cartesian sampling patterns. The reason is that Cartesian sequences - which are used in the vast majority of clinical exams - are challenging to work with in the context of conventional MRF [17] whereas with MR-STAT parameter maps can be reconstructed successfully even when using very short acquisitions in the order of seconds per slice.

Theory

In this section we first review the MR-STAT framework as presented by Sbrizzi et al. [12]. Then we discuss the computational challenges resulting from the large scale reconstruction problem and we propose techniques to deal with these challenges.

MR-STAT Framework

The time evolution of a single spin isochromat $\mathbf{m}=\left(m_{x},m_{y},m_{z}\right)$ with spatial coordinates $\mathbf{r}=(x,y,z)$ and tissue properties $\boldsymbol{\theta}=\left(T_{1},T_{2},\ldots\right)$ is governed by the Bloch equations.

Let $m=m_{x}+im_{y}$ be the transverse component of the magnetization in the rotating frame. The demodulated time-domain signal $s$ is equal to the volume integral of the transverse magnetisation of all spins within the field of view $V$, weighted by their effective proton spin densities $\rho$. For the purpose of this work, $\rho$ includes also the amplitude of the receive coil sensitivity and the transceive phase, thus $\rho$ is a complex quantity, i.e. $\rho=\rho_{x}+i\rho_{y}$. In short:

[TABLE]

After discretization of the field of view $V$ into $N_{v}$ voxels, each having volume $\Delta_{V}$, equation (1) becomes

[TABLE]

Here $m_{j}$ is the magnetization in voxel $j$, which can be computed by numerical integration of the Bloch equations.

Let $N_{t}$ be the total number of signal samples and let $t_{1},\ldots,t_{N_{t}}$ denote the sampling times. Define the magnetization vector $\mathbf{m}_{j}$ in voxel $j$ as

[TABLE]

and the signal vector $\mathbf{s}\in\mathbb{C}^{N_{t}}$ as

[TABLE]

Note that if we introduce the magnetization matrix $\textbf{M}\in\mathbb{C}^{N_{t}\times N_{v}}$,

[TABLE]

and proton density vector $\boldsymbol{\rho}\in\mathbb{C}^{N_{v}}$

[TABLE]

then s can be written as

[TABLE]

Let $N_{p}$ denote the number of distinct parameters per voxel (including real and imaginary parts of the proton density). Then ${s}$ depends on $N:=N_{v}\times N_{p}$ different parameters. All parameters are concatenated into a single vector $\boldsymbol{\alpha}\in\mathbb{R}^{N}$ in such a way that indices $\left\{j+kN_{v}\mid k=0\ldots,N_{p}-1\right\}$ denote the parameters associated with voxel $j$.

Now, given a vector of measured time-domain samples $\mathbf{d}\in\mathbb{C}^{N_{s}}$, define the residual vector $\textbf{r}\in\mathbb{C}^{N_{t}}$ as

[TABLE]

and define the non-linear least-squares objective function $f:\mathbb{R}^{N}\rightarrow\mathbb{R}$ as

[TABLE]

The parameter maps $\boldsymbol{\alpha}^{*}$ are obtained by numerically solving

[TABLE]

subject to physical constraints represented by the Bloch equations and realistically attainable intervals for the parameters.

Computational Challenges

Note that (10) is a non-linear optimization problem that requires iterative algorithms to be solved. At each iteration, the signal $\mathbf{s}=\mathbf{M}\boldsymbol{\rho}$ needs to be computed and that requires the Bloch equations to be integrated for each voxel. In addition, the gradient of $f$ (i.e. the vector of partial derivatives of $f$ with respect to each of the parameters) needs to be computed. From the least-squares structure of the problem it follows that the gradient can be expressed as

[TABLE]

where $\mathbf{J}\in\mathbb{C}^{N_{t}\times N_{v}}$ is the Jacobian matrix defined as

[TABLE]

$\mathbf{J}^{H}$ is the Hermitian transpose of $\mathbf{J}$ and $\Re$ is the real-part operator.

A gradient-descent type algorithm could be used to minimize (10) but it may result in poor convergence (see Chapter 3 of Nocedal and Wright [16]). Second-order methods (i.e. Newton methods) typically lead to better convergence. At each iteration, these methods require the inversion of a linear system involving (an approximation to) the Hessian matrix $\mathbf{H}\in\mathbb{R}^{N\times N}$, which includes curvature information and is defined as

[TABLE]

A second-order MR-STAT reconstruction algorithm would follow the steps as outlined in Algorithm 1:

Using Algorithm 1 for MR-STAT poses several practical challenges due to the large scale of the problem.

First of all, to estimate $N$ parameters, the number of sample points $N_{t}$ will in general be in the order of $N=N_{p}\times N_{v}$ as well. Assuming $N_{t}\approx N$, it follows that $\mathbf{M}$ will be of size $N_{t}\times N_{v}\approx(N_{p}N_{v})\times N_{v}$ (complex entries) and $\mathbf{J}$ will be of size $2N_{t}\times N\approx(N_{p}N_{v})\times(N_{p}N_{v})$ (complex entries). Since $\mathbf{H}$ will be of size $(N_{p}N_{v})\times(N_{p}N_{v})$ as well, it follows that all three matrices scale with $N_{v}^{2}$. In Table 1, the required computer memory to store matrices of these sizes is reported for various values of $N_{v}$ for the case $N_{p}=4$. It can be seen that, even for 2D acquisitions, it will be infeasible to store these matrices in memory for clinically relevant resolutions.

Secondly, the actual time needed to compute the entries of $\mathbf{M}$, $\mathbf{J}$ and $\mathbf{H}$ scales with $N_{v}^{2}$ as well. When using a regular desktop computer the reconstruction times quickly become too long to make MR-STAT useful in clinical practice.

Fortunately, as will be detailed in the next section, Algorithm 1 only requires matrix-vector products with the matrices $\mathbf{M},\mathbf{J}$ and (approximations to) $\mathbf{H}$. These matrix-vector product can be computed without having to store the full matrices in memory. Moreover, the computation of the matrix-vector products can be efficiently distributed among multiple computing cores on a high performance computing cluster, reducing the MR-STAT computation times to acceptable levels for off-line reconstructions.

Solution Strategies

Computing the time-domain signal s

In the first step of Algorithm 1 we need to compute $\mathbf{r}=\mathbf{d}-\mathbf{s}$ for the current estimate of the parameters $\boldsymbol{\alpha}$. Recall that

[TABLE]

Since the time evolution of the magnetization in each voxel is assumed to be independent from other voxels, the $\mathbf{m}_{j}$ can be computed independently from each other. In particular, storage of the matrix $\mathbf{M}$ is not required for computing $\mathbf{s}$, see Algorithm 2.

Note that Algorithm 2 only requires the allocation of two vectors of length $N_{t}$, which is feasible on modern computing architectures for both 2D and 3D acquisitions. The computation of $\mathbf{s}$ can then be parallelized using $N_{c}$ computing cores by following the procedure outlined in Algorithm 3 (see also [18, 19]).

The communication requirements for this parallelized algorithm can be summarized as follows:

• •

To distribute the parameters, the master process sends $N/N_{c}$ parameters to each of the $N_{c}$ slaves.

• •

To receive the local signals from the slaves, each slave sends a vector of length $N_{t}\approx N$ to the master process.

Computing the gradient g

To compute $\textbf{g}=\nabla f$ for the current estimate of the parameters $\boldsymbol{\alpha}$, recall that

[TABLE]

Since $\mathbf{J}$ is defined as

[TABLE]

it follows that

[TABLE]

To compute the $\frac{\partial\textbf{r}}{\partial\alpha_{i}}$, again note that the magnetization in different voxels is assumed to evolve independently. Hence if $\alpha_{i}$ is a parameter associated with voxel $j$ (i.e. $j=i\mod N_{v}$), it follows that

[TABLE]

Using Algorithm 4 only requires storage of one vector of length $N$ for the output and - in principle - one complex vector of length $N_{t}$ to store the intermediate partial derivative vector. In practice we will compute the $N_{p}$ partial derivatives for each voxel simultaneously so that $N_{p}$ complex vectors of length $N_{t}$ are stored simultaneously. For a high-resolution 2D scan this only requires limited memory (in the order of tens of megabytes).

Utilizing Algorithm 4, the computation of $\mathbf{g}$ can be performed in parallel as outlined in Algorithm 5. The parallelization schemes for both the signal and gradient computations are visualized in Figure 1.

Communication requirements for the parallel gradient computation can be summarized as follows:

• •

To distribute the parameters, the master process sends $N/N_{c}$ parameters to each of the $N_{c}$ slaves.

• •

To distribute the residual vector the master process sends a vector of length $N_{t}$ to each slave.

• •

To receive the local gradients from the slaves, each slave sends a vector of length $N/N_{c}$ to the master process.

Note that for both algorithms 3 and 5, the communication requirements scale linearly with the number of parameters $N$ for a fixed number of cores $N_{c}$. Since $N=N_{v}\times N_{p}$, it follows that the communication requirements scale linearly with $N_{v}$ as well. As discussed in the previous section the computational requirements scale quadratically with $N_{v}$. Therefore we hypothesize that, as long as $N_{c}\ll N_{v}$, the communication overhead is negligible compared to reduction in computation times achieved by dividing the computation load over $N_{c}$ computing cores. That is, we expect the total computation time to decrease linearly with the number of cores available under this assumption. This hypothesis is confirmed in Figure 2 in the results section.

Incorporating Curvature Information

Given the ability to compute the gradient $\mathbf{g}$ using the matrix-free, parallelized algorithm from the previous subsection, in principle the so called limited-memory Broyden–Fletcher–Goldfarb–Shanno (“L-BFGS”, [20]) method can be applied to obtain the update step $\mathbf{p}$ at each iteration. The L-BFGS method approximates the inverse of the Hessian matrix using a limited number of gradient vectors from previous iterations. However, in practice the L-BFGS method was observed to result in poor convergence.

Alternatively, since we are dealing with a least-squares problem, a Gauss-Newton method might be used in which the Hessian matrix $\mathbf{H}$ in Algorithm 1 is approximated by $\Re(\mathbf{J}^{H}\mathbf{J})$ and

[TABLE]

is solved to obtain update steps $\mathbf{p}$. Note that the matrix $\Re(\mathbf{J}^{H}\mathbf{J})$ is of the same size as the Hessian matrix itself and thus, in principle, cannot be stored into computer memory. If, however, we use iterative techniques (e.g. a Conjugate Gradient method) to solve the linear system $\Re(\mathbf{J}^{H}\mathbf{J})\mathbf{p}=-\Re(\mathbf{g})$, we only need matrix-vector products with $\Re(\mathbf{J}^{H}\mathbf{J})$. In the previous subsection it was outlined how matrix-vector products of the form $\mathbf{J}^{H}\mathbf{v}$ may be computed in a matrix-free, parallelized fashion. Similar techniques can be applied to matrix-vector products of the form $\mathbf{J}\mathbf{v}$. Hence matrix-vector products of the form $\Re(\mathbf{J}^{H}\mathbf{J})\mathbf{v}$ can be computed in a matrix-free, parallelized fashion by first computing $\mathbf{y}=\mathbf{J}\mathbf{v}$ and subsequently computing $\Re(\mathbf{J}^{H}\mathbf{y})$. With this technique, the linear system in equation (19) can be solved numerically even for large scale problems. In practice it will not be necessary to solve equation (19) to high precision and the number of iterations in this inner loop can be limited, resulting in an inexact Gauss-Newton method (see [15] and Algorithm 7.2 in [16]) as outlined in Algorithm 6

Methods

The matrix-free, parallelized MR-STAT reconstruction algorithms was tested on both simulated and experimentally acquired data.

Pulse Sequence

In all test cases, a transient-state 2D balanced gradient-echo pulse sequence similar to the pulse sequence in Sbrizzi et al [12] was used. Throughout the whole sequence the TR was fixed and TE was set to TR/2. A linear, Cartesian sampling strategy was employed together with time-varying flip angles that change according to a smoothly varying pattern. We refer to Supporting Information S1 for more details on the sampling trajectory and flip angle pattern. The phase of the RF pulse alternated between [math] and $180$ degrees. Changing the flip angles prevents the spins from reaching a steady-state and by following a smoothly varying pattern the spin-echo behaviour of bSSFP sequences [21] is preserved to a large extent. This spin-echo like behaviour is needed for proper $T_{2}$ estimation and at the same time it also effectively eliminates sensitivity to $\Delta B_{0}$ within a certain passband of off-resonances [22]. An added benefit of the smoothly changing flip angle train is the improved convexity of the minimization landscape [23].

Each RF pulse has a Gaussian envelope and at the start of the pulse sequence a non-selective inversion pulse is played out for enhanced $T_{1}$ encoding. The pulse sequence was implemented on a 1.5 T clinical MR system (Ingenia, Philips Healthcare, Best, The Netherlands).

Reconstructions

All reconstruction code was written in the open-source Julia programming language [24]. To compute the MR-signal for a given set of parameters, an optimized Bloch-equation solver was implemented which also takes into account also the slice profile [25]. To compute exact partial derivatives algorithmic differentiation in forward mode [26] was implemented. We refer to the Supporting Information S2 for more details.

The inexact Gauss-Newton method was implemented using a trust-region framework (following [15] and Algorithm 7.2 in [16]). In order to facilitate bound constraints on the parameters, reflection at feasible boundaries was incorporated [27]. For the L-BFGS method, an implementation from the Optim.jl package [28] was used. The reconstruction algorithm was implemented on a high performance computing cluster which consists of multiple Intel Xeon Gold 6148 nodes with 40 cores each, on which the CentOS Linux 7 (Core) operating system is installed.

For all experiments, $T_{1},T_{2}$ and $\rho$ (complex) maps are reconstructed. For the data obtained with clinical MR systems we also reconstruct $|B_{1}^{+}|$ to take into account transmit field inhomogeneities. The off-resonance $\Delta B_{0}$ was set to zero and thus it was not reconstructed because of the flat spectral response of the balanced sequence within the passband. The non-linear parameters were initialized as follows: $T_{1}=1000$ ms, $T_{2}=100$ ms, $|B_{1}^{+}|=1\,$ a.u. and $\Delta B_{0}=0$ Hz. In previous work [12] the Variable Projection method (“VARPRO”, [29]) was utilized to separate the linear parameters (i.e. proton density) from the non-linear parameters. The VARPRO method in principle requires computing (through SVD or QR decomposition) and storing an orthogonal basis for the matrix $\mathbf{M}$. For the matrix sizes in the current work that would be computationally infeasible and it is non-trivial to extend the VARPRO technique to matrix-free methods. Therefore, in the current work we treat the proton densities as non-linear parameters. We only make use of the linearity to provide an initial guess for the proton densities. That is, given the initial guess for the non-linear parameters, the (complex) proton density was initialized as the least squares solution to the linear system $\mathbf{M}(\boldsymbol{\alpha}_{0})\boldsymbol{\rho}=\mathbf{d}$ obtained using a linear solver (LSQR). Based on the resulting initial guess for the proton density, a mask was drawn to exclude regions with no significant proton density from subsequent simulations.

In all reconstructions, logarithmic scaling is applied to both $T_{1}$ and $T_{2}$ parameters. The variable substitution brings both variables in a similar range and it thus improves convergence of the algorithm.

The reconstruction code will be made available online at https://gitlab.com/mrutrecht/mrstat after acceptance of the manuscript for publication.

Numerical Brain Simulation

Signal from a numerical brain phantom [30] with a field-of-view of $192$ mm $\times$ $192$ mm and voxel size of $1$ mm $\times$ $1$ mm was simulated using the transient-state pulse sequence. A total number of $8\times 192=1536$ readouts were simulated (each phase encoding line was acquired eight times but note that for each readout line the flip angle and thus state of the spins is different) with a TR of $7.88$ ms and a TE of $3.94$ ms. The total sequence duration was $12.1$ s.

Reconstructions were performed using 64 cores. The number of outer and inner iterations for the inexact Gauss-Newton method were limited to 40 and 20, respectively.

For comparison purposes, we also performed MRF reconstructions on signal from the numerical brain phantom using the Cartesian trajectory, as well as signal from radial and spiral trajectories for which MRF is known to work well. In all three cases the same flip angle train, TE and TR were used. For the radial case, $k_{max}$ was extended by a factor of $√2$ and each readout the spoke was rotated by the golden angle. For the spiral acquisition a variable density spiral was generated [31, 32] that would require 24 interleaves to fully sample the inner region of k-space and 48 interleaves for the outer region of k-space. The spiral was rotated by the golden angle each readout. Data $\mathbf{d}^{MRF}$ was then simulated by applying a forward operator, consisting of the (non-uniform) FFT [33] and an undersampling operator, on fingerprints simulated using the numerical brain phantom. To perform the MRF reconstructions, a dictionary consisting of 100 $T_{1}$ values from $0.1$ s to $5.0$ s in increments of $4$ % and 100 $T_{2}$ values from $0.01$ s to $2.0$ s in increments of $5.5$ % was generated. The $T_{1}$ and $T_{2}$ values of the phantom were not used in generating the dictionary. The dictionary was compressed in the time direction to rank 5 [34] using the SVD. For all trajectories (linear) low-rank forward operators $\boldsymbol{E}^{MRF}$ were formed that consisted of the expansion of low-rank coefficients to the full time-series, a nuFFT operator, and a sampling operator compression [13]. Low-rank snapshot images $\mathbf{x}^{MRF}$ were reconstructed from the undersampled data $\mathbf{d}^{MRF}$ by solving the linear system

[TABLE]

with LSQR (similar to e.g. [35] and low-rank inversion in [13]). Finally, dictionary matching with the compressed dictionary was performed on $\mathbf{x}^{MRF}$ to obtain the parameter estimates.

To further study the effect of noise on the MR-STAT reconstruction algorithm, additional reconstructions were performed where complex Gaussian noise was added to the simulated signal such that $\|\text{signal}\|_{2}/\|\text{noise}\|_{2}=50,25$ and $10$.

Gel Phantom Experiment

Signal from a 2D transverse slice of six gadolinium-doped gel phantoms (TO5, Eurospin II test system, Scotland) was collected on the 1.5 T MR system using the manufacturer’s thirteen-channel receive headcoil. In total $8\times 96=768$ readout lines were acquired with a spatial resolution of $1$ mm $\times$ $1$ mm $\times$ 5 mm and a field-of-view of $96$ mm $\times$ $96$ mm. The TR and TE were $7.4$ ms and $3.7$ ms, respectively, resulting in a total acquisition time of $5.7$ s. For reproducibility purposes the MR-STAT scan was repeated four times with full relaxation in between the different scans.

Parameters that describe the pulse sequence were exported from the scanner and subsequently loaded into Matlab [36]. The measured signals from different receive channels were compressed into a single signal by applying the principal component analysis and choosing the principle mode [37].

Reconstructions of the parameter maps were performed using the inexact Gauss-Newton method on the computing cluster using 32 cores. The number of inner iterations was limited to fifteen whereas the number of outer iterations was limited to ten.

To assess correctness of the $T_{1}$ and $T_{2}$ maps reconstructed with MR-STAT, data was also acquired using gold standard methods in the form of an inversion-recovery single spin-echo protocol with inversion times of [50, 100, 150, 350, 550, 850, 1250] ms for $T_{1}$ mapping and a single echo spin-echo protocol with echo times of [8, 28, 48, 88, 138, 188] ms for $T_{2}$ mapping.

In-vivo experiments

Using the 1.5 T clinical MR system we also acquired signal from 2D transverse slices of the brain in three healthy volunteers. Each volunteer gave written informed consent. For each acquisition, a total of $8\times 192=1536$ readout lines were acquired with acquisition parameters as reported in Table 2. After masking approximately 25000 voxels remain for which quantitative parameters are estimated. The MR-STAT reconstructions were performed with 64 cores using the reconstruction settings as for the gel phantom experiment.

To demonstrate the effect of accelerated acquisitions, we also performed reconstructions using time-domain data corresponding to the first 896 TRs and the first 448 TRs from one of the subjects. The corresponding acquisition times were $6.8$ s and $3.4$ s respectively.

One of the in-vivo brain datasets was also used to test the effectiveness of the parallelization scheme. Individual matrix-vector products of the form $\mathbf{M}\mathbf{v}$ and $\mathbf{J}^{H}\mathbf{v}$ were computed and timed for $5,10,20,40,60,80,100,120$ and $240$ cores respectively.

Results

Parallelization

In Figure 2 the time required to compute matrix-vector products of the form $\mathbf{M}\mathbf{v}$ and $\mathbf{J}^{H}\mathbf{v}$ for the one of the in-vivo datasets is shown for an increasing number of computing cores $N_{c}$. Initially we observe a linear decrease in computation times, however this linear decrease flattens beyond approximately 64 cores. This effect can be explained by the increase in communication overhead when using more cores and increased competition between cores for shared resources like memory bandwidth and cache memory. Although the linear decrease flattens beyond 64 cores, a decrease in computation times is still observed even when going towards 240 cores. Because for MR-STAT reconstruction times are dominated by the computation of these matrix-vector products, the reconstruction times can thus be effectively reduced by the proposed parallel implementation.

Numerical Brain Phantom

The $T_{1},T_{2}$ and proton density maps reconstructed using (Cartesian) MR-STAT and Cartesian, radial and spiral MRF are shown in Figure 3 as well as the corresponding absolute relative error maps. It can be seen that the parameter maps reconstructed with either MR-STAT or spiral MRF are in excellent agreement with the ground truth. The radial MRF reconstructions show stronger residual streaking artefacts but in general the the estimated parameter values are close to the ground truth. For the Cartesian case the MRF reconstruction is unable to cope with the high undersampling (factor 192), resulting in severaly biased parameter maps.

To quantify the quality of the reconstructions, normalized root mean square errors (“NRMSE”), high-frequency error norms (“HFEN”, [38], with standard deviation of 1.5 pixels) and mean absolute relative errors (“MAPE”) were computed and are reported in Table 3. It can be seen that the MR-STAT reconstruction results in the lowest NRSME and MAPE for all three parameters. The HFEN for the radial and spiral MRF and Cartesian MR-STAT reconstructions are similar.

In Figure 4 convergence curves for MR-STAT with the inexact Gauss-Newton method for different SNR levels (50, 25 and 10) are shown as well as mean absolute percentage errors per iteration for $T_{1}$, $T_{2}$ and proton density. For the higher SNR cases the error values stabilize and the method converges after relatively few, e.g. ten, outer iterations. It can be seen that for the lowest SNR case, overfitting occurs after around six iterations. Based on these observations the number of outer iterations for the in-vivo case was chosen to be ten.

Gel Phantoms

In Figure 6, reconstructed $T_{1}$ and $T_{2}$ maps for the gel phantoms are shown and the mean $T_{1}$ and $T_{2}$ values per tube are compared to the gold standard measurements. It can be seen that the mean values are in excellent agreement. The mean values reported for the different repetitions of the MR-STAT scans are also in good agreement with each other (i.e. within standard deviations). In general the standard deviations for the reconstructed $T_{2}$ values is higher than for $T_{1}$ values, indicating a much stronger encoding of $T_{1}$ information into the signal which can be explained by the inversion pulse at the start of the sequence.

To reconstruct the parameter maps, only five iterations of the reconstruction algorithm were needed and the total reconstruction time was approximately nine minutes using 32 computing cores. In Figure 5 a logarithmic plot of the measured signal magnitude and the residual vector after the fifth iteration are displayed for one of the MR-STAT repetitions. Histograms of the measured noise and the residual vectors are also shown. It can be seen that the residual vector follows a zero-mean Gaussian distribution with standard deviation similar to the noise, indicating that the model used in MR-STAT is able to adequately describe the measured time-domain signal.

High-resolution 2D brain scan

In Figure 7, the reconstructed $T_{1},T_{2}$ and proton density (magnitude) maps for the in-vivo brain scans performed on the three volunteers are shown. The maps show clear contrast between white matter, gray matter and cerebrospinal fluid (“CSF”). The maps corresponding to subject 3 appear noisier compared to the maps corresponding to subjects 1 and 2, which can be explained by the differences in slice thickness in the acquisition ($3$ mm vs $5$ mm). Mean $T_{1}$ and $T_{2}$ values and standard deviations in regions of white- and gray matter are reported in Table 4. The mean values are generally in good agreement with values found in literature for 1.5 T experiments [39, 6, 40] although we do observe an underestimation compared to some other studies, especially in white matter. We expect the underestimation is related to magnetization transfer that is known to affect the signal of balanced gradient-echo sequences (in a way that depends on the used $TR$ and RF pulse duration) [41, 42]. The reconstruction time for each slice was approximately five hours using 64 cores.

In Figure 8 we show $T_{1},T_{2}$ and proton density (magnitude) for the same 2D brain slice but reconstructed using data corresponding to, respectively, $13.6$ s, $7.8$ s and $3.4$ s acquisitions. It can be seen that the maps corresponding to the 6.8 s acquisition are comparable to the maps corresponding to the $13.6$ s acquisition except that more noise is present. Depending on the application it might be more beneficial to repeat such a shorter sequence multiple times for noise averaging instead of scanning with the longer sequence. An added benefit of a shorter sequence duration is that the Bloch simulations are faster and thus reconstruction times are reduced by approximately the same factor with which the scantime is reduced. For the $3.4$ s acquisition the MR-STAT problem (Eq. 10) is underdetermined in the sense that the number of datapoints is less than the number of unknowns in the problem. As can be seen in the reconstructed maps, this mostly results in biases in the CSF values. Note that in none of the reconstructions parallel imaging or compressed sensing techniques were utilized.

Discussion & Conclusion

MR-STAT is a framework for obtaining multiple quantitative parameter maps by fitting directly to measured time-domain data obtained from one short scan. Rather than relying on the FFT for spatial localisation of signal in a separate step, the spatial localisation and parameter estimation are performed simultaneously by iteratively solving a single non-linear optimization problem using a signal model that explicitly includes the spatial encoding gradients. The inherent large scale of the problem brings along new challenges in terms of computer memory requirements and computation times that make it difficult to perform MR-STAT reconstructions at high resolutions. To address these issues, we have presented a parallel and matrix-free reconstruction algorithm for MR-STAT and demonstrated that it can be used to generate high-resolution quantitative parameter maps.

All MR-STAT experiments in the current work have been performed with linear, Cartesian sampling strategies. This sampling strategy offers important advantages in the form of robustness to hardware imperfections (e.g. eddy currents, especially for gradient-balanced sequences [43, 44]), less susceptibility to $\Delta B_{0}$ related blurring artefacts [45] and direct availability on clinical MR systems. Within the conventional MRF framework it is more challenging to work with Cartesian sampling strategies, as demonstrated using the simulation experiments. Studies that perform Cartesian MRF [46, 47] therefore typically acquire multiple readout lines per snapshot, resulting in much longer acquisition times compared to non-Cartesian MRF acquisitions. A formal explanation of why Cartesian acquisitions are less suitable for MRF is reported in Stolk et al. [17]. More advanced iterative MRF reconstructions [48, 49, 13, 14] might perform better with Cartesian sampling than the currently used MRF reconstructions (low-rank inversion followed by low-rank dictionary matching) and an in-depth comparison will be the subject of further studies. It should also be noted that neither the MR-STAT framework nor the currently proposed reconstruction algorithm are restricted to Cartesian sampling and further research is also aimed at incorporating non-Cartesian trajectories into MR-STAT.

An additional benefit of the volumetric signal model used in MR-STAT over FFT-based methods is that dynamic behaviour during the readouts (e.g. $T_{2}$-decay and $\Delta B_{0}$ induced phase accumulation) is taken into account. This may especially be beneficial for improving reconstructions based on acquisitions with long readouts (e.g. spiral readouts).

MR-STAT reconstructions are performed by solving a non-linear optimization problem using gradient-based iterative methods. No pre-computed dictionary is used. Compared to a dictionary-matching approaches there are no discretization errors and the reconstruction procedure is also flexible with respect to changes in sequence parameters (e.g. no rebuilding of a dictionary required when scan settings change). A downside of using iterative reconstruction algorithms to solve non-linear optimization problems is the risk of landing in a local minimum. In practice, with the currently used pulse sequence with smoothly changing flip angles and initial guess of the parameters, we have not encountered issues with local minima [23].

Whereas with MRF the addition of new parameters results in an exponential increase in dictionary size (and thus also an exponential increase in dictionary generation and matching time), with MR-STAT additional parameters can be added at a quadratic increase in computation time. The quadratic increase can be explained as follows. The total number of parameters to be reconstructed $N$ increases linearly with the number of parameters per voxel ($N=N_{p}N_{v}$). Since the minimum number of time points $N_{t}$ that needs to be acquired - and thus simulated - is in the order of $N$, the computation time per Bloch simulation increases linearly as well. In addition, the number of partial derivative computations that needs to be performed per voxel also increases linearly. That is, both the number of rows and columns of the Jacobian $\mathbf{J}$ increase linearly, resulting in approximately a quadratic increase in computation time. In this respect we do note that, although currently $\Delta B_{0}$ maps are not reconstructed (because the employed bSSFP sequence used in this work is designed not to be sensitive to $\Delta B_{0}$ within the “passband”), it is part of our all Bloch simulations and partial derivative computations. In addition, for the MR-STAT experiments described in the manuscript we used pulse sequences such that $N_{t}\approx 2N$ so that the problem remains overdetermined when an additional parameter is reconstructed. Therefore, assuming a pulse sequence is used that has sufficient $\Delta B_{0}$ encoding [10, 50], we do not expect to see an increase in computation times when reconstructing $\Delta B_{0}$ as an additional parameter.

For the phantom experiment we observed that the noise levels was reached for the residuals. However, this was not observed for the in-vivo case as certain effects are still accounted for in the model. Examples include patient motion, blood flow, magnetization transfer and diffusion effects.

A limitation of the proposed method is that at this moment reconstruction times are still long for high-resolution scans, especially when compared to the dictionary matching procedures used in MRF. Even when employing a high performance computing cluster, reconstruction times in the order of hours for a single 2D brain slice. Although possible from a memory point-of-view, 3D reconstructions will take too long for practical purposes with the current reconstruction setup. The main bottleneck in the reconstructions is formed by the partial derivative computations needed to solve equation (19). Further research is aimed at performing these computations on GPU architectures [51, 52], reducing the computational effort through algorithmic improvements [53] and through the use of surrogate models [54]. Together with (cloud) computing resources becoming cheaper and more accessible over time, we believe it is possible to accelerate the computations to an extent that MR-STAT becomes applicable in clinical settings.

Further research is also aimed at reduction of acquisiton time and improving precision and accuracy of the MR-STAT parameter maps by incorporating parallel imaging [55], compressed sensing and through sequence optimization.

The main aim of the MR-STAT project is to explore possibilities to achieve very short acquisition times beyond what is possible with FFT-based frameworks. Although the MR-STAT framework in principle allows for much flexibility in the data acquisition process (e.g. non-Cartesian acquisitions), in the current work we have opted for Cartesian sampling patterns because of their robustness to hardware imperfections and because they clearly exemplify the benefits of skipping the FFT step (i.e. no introduction of artificial aliasing noise through application of the FFT on undersampled k-spaces). An additional benefit is the direct availability of such sequences on clinical MR systems. In the current work we used constantly varying flip angle trains, however, as shown in Supporting Information S3, MR-STAT could even be used with Cartesian bSSFP sequences with a fixed flip angle per k-space that require little to no pulse programming for their implementation.

Acknowledgements

This work was funded by the Dutch Technology Foundation, grant #14125.

Legend of Supporting Information

Supporting Information S1: MR-STAT pulse sequence 2. 2.

Supporting Information S2: Algorithmic Differentiation of Bloch equation solver 3. 3.

Supporting Information S3: MR-STAT with Cartesian sampling and constant flip angle per k-space