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

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 fast Fourier transform (FFT). To do this at high resolution, specialized algorithms are required to solve the underlying large-scale nonlinear 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 (1 mm  ×  1 mm 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 MR fingerprinting reconstructions are shown to fail.

balanced steady‐state free precession

cerebrospinal fluid

fast Fourier transform

high‐frequency error norms

limited‐memory Broyden–Fletcher–Goldfard–Shanno

mean absolute percentage error

magnetic resonance fingerprinting

magnetic resonance spin tomography in the time domain

normalized root‐mean‐square errors

non‐uniform fast Fourier transform

quantitative magnetic resonance imaging

singular value decomposition

echo time

repetition time

variable projection method.

INTRODUCTION

Conventional magnetic resonance imaging (MRI) methods rely on the Fourier transform relationship between the signal and the 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 and relaxation times. Quantitative magnetic resonance imaging (qMRI) methods aim to estimate MR‐related biophysical properties like and relaxation times. Quantitative images could provide additional diagnostic value and are more suited for the purpose of multicenter studies and computer‐aided diagnosis.1, 2 The most straightforward and robust choices for and mapping sequences, that is, 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 about 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 magnetic resonance fingerprinting (MRF).11 In MRF, a transient‐state pulse sequence with quasi‐random components is used and many highly undersampled ‐spaces are acquired in a single, short acquisition. A fast Fourier transform (FFT) is applied to each ‐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 reduce qMRI acquisition times drastically.

Magnetic Resonance Spin TomogrAphy in Time‐domain (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 nonlinear optimisation problem is solved numerically, where 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 with MRF, the MR‐STAT approach results in different trade‐offs made in the reconstruction. Since the FFT is no longer used explicitly 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 combined naturally into a single reconstruction process. There is no reliance on dictionary compression13 or compressed sensing14 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 nonlinear 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 nontrivial, because the FFT is not used to decouple the unknowns spatially. 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 benefits only partly from the abovementioned advantages of using a volumetric signal model, for example, 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);

is suitable for a parallel computing implementation to reduce computation times; and

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 Steihaug15 and Algorithm 7.2 in Nocedal and Wright16) 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 of the order of seconds per slice.

Theory

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

MR‐STAT framework

The time evolution of a single‐spin isochromat with spatial coordinates and tissue properties is governed by the Bloch equations.

Let be the transverse component of the magnetization in the rotating frame. The demodulated time‐domain signal is equal to the volume integral of the transverse magnetisation of all spins within the field of view , weighted by their effective proton spin densities . For the purpose of this work, also includes the amplitude of the receive coil sensitivity and the transceive phase, thus is a complex quantity, that is, . In short:

After discretization of the field of view into voxels, each having volume , Equation (1) becomes

Here, is the magnetization in voxel , which can be computed by numerical integration of the Bloch equations.

Let be the total number of signal samples and let denote the sampling times. Define the magnetization vector in voxel as and the signal vector as

Note that if we introduce the magnetization matrix , and proton density vector , then can be written as

Let denote the number of distinct parameters per voxel (including real and imaginary parts of the proton density). Then depends on different parameters. All parameters are concatenated into a single vector in such a way that indices denote the parameters associated with voxel .

Now, given a vector of measured time‐domain samples , define the residual vector as and define the nonlinear least‐squares objective function as

The parameter maps are obtained by solving numerically subject to the physical constraints represented by the Bloch equations and realistically attainable intervals for the parameters.

Computational challenges

Note that Equation (10) is a nonlinear optimization problem that requires iterative algorithms to be solved. At each iteration, the signal needs to be computed and that requires the Bloch equations to be integrated for each voxel. In addition, the gradient of (i.e., the vector of partial derivatives of 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 where is the Jacobian matrix defined as

is the Hermitian transpose of , and is the real‐part operator.

A gradient‐descent type algorithm could be used to minimize Equation (10), but it may result in poor convergence (see chapter 3 of Nocedal and Wright16). 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 , which includes curvature information and is defined as

A second‐order MR‐STAT reconstruction algorithm would follow the steps 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 parameters, the number of sample points will in general be of the order of as well. Assuming , it follows that will be of size (complex entries) and will be of size (complex entries). Since will be of size as well, it follows that all three matrices scale with . In Table 1, the required computer memory to store matrices of these sizes is reported for various values of for the case . It can be seen that, even for 2D acquisitions, it will be infeasible to store these matrices in memory for clinically relevant resolutions.

Table 1

On‐disk sizes of MR‐STAT matrices for and for an increasing number of voxels . The memory sizes are computed as bytes ( ), bytes ( ), and bytes ( ), 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 (N_v)	M	J	H
64×64	4.096	1 GB	4 GB	2 GB
128×128	16.384	16 GB	64 GB	32 GB
256×256	65.536	256 GB	1.024 GB	512 GB
512×512	262.144	4.096 GB	16.384 GB	8.192 GB

Secondly, the actual time needed to compute the entries of , , and scales with 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, 1 only requires matrix–vector products with the matrices , and (approximations to) . These matrix–vector products can be computed without having to store the full matrices in memory. Moreover, the computation of the matrix–vector products can be distributed efficiently among multiple computing cores on a high‐performance computing cluster, reducing the MR‐STAT computation times to acceptable levels for offline reconstructions.

Solution strategies

Computing the time‐domain signal

In the first step of Algorithm 1, we need to compute for the current estimate of the parameters . Recall that

Since the time evolution of the magnetization in each voxel is assumed to be independent of other voxels, the can be computed independently from each other. In particular, storage of the matrix is not required for computing ; see Algorithm 2.

Note that Algorithm 2 only requires the allocation of two vectors of length , which is feasible on modern computing architectures for both 2D and 3D acquisitions. The computation of can then be parallelized using computing cores by following the procedure outlined in Algorithm 3 (see also Stöcker et al.18 and Liu et al.19).

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

to distribute the parameters, the master process sends parameters to each of the slaves;

to receive the local signals from the slaves, each slave sends a vector of length to the master process.

Computing the gradient

To compute for the current estimate of the parameters , recall that

Since is defined as it follows that

To compute the , again note that the magnetization in different voxels is assumed to evolve independently. Hence, if is a parameter associated with voxel (i.e., ), it follows that

Using Algorithm 4 only requires storage of one vector of length for the output and—in principle—one complex vector of length to store the intermediate partial derivative vector. In practice, we will compute the partial derivatives for each voxel simultaneously, so that complex vectors of length 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 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.

Figure 1

Visualization of the (matrix‐free) algorithms to compute the signal [left] and the gradient [right] in a parallelized fashion

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

to distribute the parameters, the master process sends parameters to each of the slaves;

to distribute the residual vector, the master process sends a vector of length to each slave;

to receive the local gradients from the slaves, each slave sends a vector of length to the master process.

Note that, for both Algorithms 3 and 5, the communication requirements scale linearly with the number of parameters for a fixed number of cores . Since , it follows that the communication requirements scale linearly with as well. As discussed in the previous section, the computational requirements scale quadratically with . Therefore we hypothesize that, as long as , the communication overhead is negligible compared with reduction in computation times achieved by dividing the computation load over 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.

Figure 2

Time needed to compute matrix–vector products of the form and for different numbers of cores used on a high‐performance computing cluster

Incorporating curvature information

Given the ability to compute the gradient using the matrix‐free, parallelized algorithm from the previous subsection, in principle the so‐called limited‐memory Broyden–Fletcher–Goldfarb–Shanno (L‐BFGS20) method can be applied to obtain the update step 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 in Algorithm 1 is approximated by and is solved to obtain update steps . Note that the matrix is of the same size as the Hessian matrix itself and thus, in principle, cannot be stored in computer memory. If, however, we use iterative techniques (e.g., a conjugate gradient method) to solve the linear system , we only need matrix–vector products with . In the previous subsection we outlined how matrix–vector products of the form may be computed in a matrix‐free, parallelized fashion. Similar techniques can be applied to matrix–vector products of the form . Hence matrix–vector products of the form can be computed in a matrix‐free, parallelized fashion by first computing and subsequently computing . 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 Steihaug15 and Algorithm 7.2 in Nocedal and Wright16) 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 0 and . 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 sequences21 is preserved to a large extent. This spin‐echo‐like behaviour is needed for proper estimation and at the same time it also effectively eliminates sensitivity to 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 encoding. The pulse sequence was implemented on a 1.5T 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 the slice profile.25 To compute exact partial derivatives, algorithmic differentiation in forward mode26 was implemented. We refer to Supporting Information S2 for more details.

The inexact Gauss–Newton method was implemented using a trust‐region framework (following Steihaug15 and Algorithm 7.2 in Nocedal and Wright16). 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 package28 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, , and (complex) maps are reconstructed. For the data obtained with clinical MR systems, we also reconstruct to take transmit field inhomogeneities into account. The off‐resonance was set to zero and thus it was not reconstructed, because of the flat spectral response of the balanced sequence within the passband. The nonlinear parameters were initialized as follows:  ms,  ms, a.u. and Hz. In previous work,12 the variable projection method (VARPRO29) was utilized to separate the linear parameters (i.e., proton density) from the nonlinear parameters. The VARPRO method in principle requires computation (through SVD or QR decomposition) and storage of an orthogonal basis for the matrix . For the matrix sizes in the current work, that would be computationally infeasible and it is nontrivial to extend the VARPRO technique to matrix‐free methods. Therefore, in the current work we treat the proton densities as nonlinear 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 nonlinear parameters, the (complex) proton density was initialized as the least‐squares solution to the linear system 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 and parameters. The variable substitution brings both variables into a similar range and 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 phantom30 with a field of view of 192 mm 192 mm and voxel size of 1 mm 1 mm was simulated using the transient‐state pulse sequence. A total number of 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 was 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, was extended by a factor of and each readout the spoke was rotated by the golden angle. For the spiral acquisition, a variable‐density spiral was generated31, 32 that would require 24 interleaves to sample the inner region of ‐space fully and 48 interleaves for the outer region of ‐space. The spiral was rotated by the golden angle each readout. Data point was then simulated by applying a forward operator, consisting of the (non‐uniform) FFT33 and an undersampling operator, on fingerprints simulated using the numerical brain phantom. To perform the MRF reconstructions, a dictionary consisting of 100 values from 0.1–5.0 s in increments of 4% and 100 values from 0.01–2.0 s in increments of 5.5% was generated. The and values of the phantom were not used in generating the dictionary. The dictionary was compressed in the time direction to rank 534 using the SVD. For all trajectories, (linear) low‐rank forward operators 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 were reconstructed from the undersampled data by solving the linear system with LSQR (similar to e.g. Zhao et al.35 and low‐rank inversion in Assländer et al. 13). Finally, dictionary matching with the compressed dictionary was performed on to obtain the parameter estimates.

To study the effect of noise on the MR‐STAT reconstruction algorithm further, additional reconstructions were performed, where complex Gaussian noise was added to the simulated signal such that , 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.5T MR system using the manufacturer's 13‐channel receive headcoil. In total, readout lines were acquired with a spatial resolution of 1 mm 1 mm 5 mm and a field of view of 96 mm 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 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 15, whereas the number of outer iterations was limited to ten.

To assess correctness of the and maps reconstructed with MR‐STAT, data were 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 mapping and a single‐echo spin‐echo protocol with echo times of [8, 28, 48, 88, 138, 188] ms for mapping.

In‐vivo experiments

Using the 1.5T 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 readout lines were acquired, with acquisition parameters as reported in Table 2. After masking approximately, 25 000 voxels remain for which quantitative parameters are estimated. The MR‐STAT reconstructions were performed with 64 cores using the same reconstruction settings as for the gel phantom experiment.

Table 2

Acquisition parameters for in vivo MR‐STAT brain scans

Acquisition parameter	Subjects 1 and 2	Subject 3
Field strength	1.5T	1.5T
In‐plane resolution	1 mm × 1 mm	1 mm × 1 mm
Field of view	224 mm × 224 mm	224 mm × 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

To demonstrate the effect of accelerated acquisitions, we also performed reconstructions using time‐domain data corresponding to the first 896 TR and the first 448 TR from one of the subjects. The corresponding acquisition times were 6.8 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 and were computed and timed for , and 240 cores, respectively.

RESULTS

Parallelization

In Figure 2, the time required to compute matrix–vector products of the form and for the one of the in vivo datasets is shown for an increasing number of computing cores . 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 reduced effectively by the proposed parallel implementation.

Numerical brain phantom

The , 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 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 severely biased parameter maps.

Figure 3

First column: ground‐truth , and proton density maps for the numerical brain phantom. Second, third and fourth columns: reconstructed parameter maps and relative error maps for MRF with linear Cartesian, golden angle radial, and golden angle spiral trajectories, respectively. Fifth column: reconstructed parameter maps and relative error maps for MR‐STAT using a linear, Cartesian sampling trajectory. The MRF spiral and MR‐STAT reconstructions both show excellent agreement with the ground‐truth values. The radial MRF reconstructions show residual aliasing artefacts and the Cartesian MRF reconstruction is heavily biased

To quantify the quality of the reconstructions, normalized root‐mean‐square errors (NRMSE), high‐frequency error norms (HFEN38, with standard deviation of 1.5 pixels) and mean absolute percentage 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.

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 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
ρ	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

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 , , and proton density. For the higher SNR cases, the error values stabilize and the method converges after relatively few, for example, 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.

Figure 4

Top left: convergence curves for the inexact Gauss–Newton MR‐STAT method applied to data generated from the numerical brain phantom with different noise levels (SNR 50, 25, 10). In all cases, the value of the cost function converges to the value expected based on the noise level. Top right and bottom row: mean absolute percentage errors for , and proton density (magnitude) maps per iteration of the inexact Gauss–Newton method for different noise levels

Gel phantoms

In Figure 5, reconstructed and maps for the gel phantoms are shown and the mean and values per tube are compared with 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 values are higher than for the values, indicating a much stronger encoding of information in the signal, which can be explained by the inversion pulse at the start of the sequence.

Figure 5

Top row: and maps reconstructed with MR‐STAT from gel‐phantom data. Middle and bottom rows: comparison of mean and values obtained with MR‐STAT and gold standard methods for each of the six gel phantom tubes. For MR‐STAT, the acquisition has been repeated four times

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 6, 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 describe the measured time‐domain signal adequately.

Figure 6

Top row: logarithmic plot of the magnitude of the measured time‐domain data obtained from the gel phantoms and the magnitude of the residual vector entries after the fifth iteration of the inexact Gauss–Newton method. Bottom left: histogram of noise values (real and imaginary values concatenated). The noise was measured using the receive channels right before the actual acquisition and was subjected to the same preprocessing steps as the data used in the reconstruction (e.g., compression to a single channel using SVD). Bottom right: histogram of the residual vector entries (real and imaginary values concatenated) after the fifth iteration of the inexact Gauss–Newton method

High‐resolution 2D brain scan

In Figure 7, the reconstructed , 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 with the maps corresponding to subjects 1 and 2, which can be explained by the differences in slice thickness in the acquisition (3 mm versus 5 mm). Mean and 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 the literature for 1.5T experiments,6, 39, 40 although we do observe an underestimation compared with some other studies, especially in white matter. We expect the underestimation is related to magnetization transfer, which is known to affect the signal of balanced gradient‐echo sequences (in a way that depends on the TR and RF pulse duration used).41, 42 The reconstruction time for each slice was approximately five hours using 64 cores.

Figure 7

, and proton density (magnitude) maps reconstructed with MR‐STAT from in vivo brain data obtained at 1.5T (Philips, Ingenia) from multiple healthy volunteers. The in‐plane resolution was  mm2 for all three subjects. For subjects 1 and 2, the acquisition time was 13.6 s and the slice thickness was 5 mm. For subject 3, the acquisition time was 14.15 s and the slice thickness was 3 mm

Table 4

Mean and 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±48 ms	53.3±4.0 ms
	2	542±48 ms	57.4±3.8 ms
	3	519±54 ms	56.1±4.3 ms
Putamen (gray matter)	1	874±64 ms	74.8±4.4 ms
	2	956±66 ms	80.2±4.5 ms
	3	895±107 ms	78.4±7.0 ms

In Figure 8, we show , and proton density (magnitude) for the same 2D brain slice but reconstructed using data corresponding to, respectively, 13.6, 7.8, 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 by which the scan time is reduced. For the 3.4 s acquisition, the MR‐STAT problem (Equation (10)) is underdetermined in the sense that the number of data points is fewer 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 were parallel imaging or compressed sensing techniques utilized.

Figure 8

In vivo , and proton density (magnitude) maps at 1 mm 1 mm in‐plane resolution reconstructed with MR‐STAT based on acquisitions of, respectively, 13.6, 6.8, and 3.4 s on a 1.5T MR system (Philips, Ingenia)

DISCUSSION AND 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 solving a single nonlinear optimization problem iteratively 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, which make it difficult to perform MR‐STAT reconstructions at high resolution. 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 sequences43, 44), less susceptibility to 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 MRF46, 47 therefore typically acquire multiple readout lines per snapshot, resulting in much longer acquisition times compared with 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 reconstructions13, 14, 48, 49 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 is 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 readouts (e.g., decay and induced phase accumulation) is taken into account. This may be especially beneficial for improving reconstructions based on acquisitions with long readouts (e.g., spiral readouts).

MR‐STAT reconstructions are performed by solving a nonlinear optimization problem using gradient‐based iterative methods. No precomputed dictionary is used. Compared with 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 nonlinear 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 increases linearly with the number of parameters per voxel ( ). Since the minimum number of time points that needs to be acquired—and thus simulated—is of the order of , 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, the number of both rows and columns of the Jacobian increases linearly, resulting in approximately a quadratic increase in computation time. In this respect, we do note that, although currently maps are not reconstructed (because the bSSFP sequence used in this work is designed not to be sensitive to within the passband), is considered as parameter in all Bloch simulations and partial derivative computations. In addition, for the MR‐STAT experiments described in the article, we used pulse sequences such that , so that the problem remains overdetermined when an additional parameter is reconstructed. Therefore, assuming a pulse sequence is used that has sufficient encoding,10, 50 we do not expect to see an increase in computation times when reconstructing as an additional parameter.

For the phantom experiment, we observed that the noise level 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 with the dictionary matching procedures used in MRF. Even when employing a high‐performance computing cluster, reconstruction times are of 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 improvements53 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 such an extent that MR‐STAT becomes applicable in clinical settings.

Further research is also aimed at reducing the acquisition time and improving the precision and accuracy of 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 to undersampled ‐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 ‐space that require little to no pulse programming for their implementation.

supporting_information.pdf

Click here for additional data file.