Model‐based super‐resolution reconstruction of T2 maps

High‐resolution isotropic T2 mapping of the human brain with multi‐echo spin‐echo (MESE) acquisitions is challenging. When using a 2D sequence, the resolution is limited by the slice thickness. If used as a 3D acquisition, specific absorption rate limits are easily exceeded due to the high power deposition of nonselective refocusing pulses. A method to reconstruct 1‐mm3 isotropic T2 maps is proposed based on multiple 2D MESE acquisitions. Data were undersampled (10‐fold) to compensate for the prolonged scan time stemming from the super‐resolution acquisition.


| INTRODUCTION
Spin-spin relaxation characterized by its relaxation time T 2 is 1 of the 2 principal relaxation mechanisms in MRI. Its quantification, typically referred to as T 2 mapping, provides important information about the tissue of interest. Previous studies have demonstrated the importance of T 2 mapping to study various neurological conditions such as stroke, 1 epilepsy, 2 multiple sclerosis, 3 and tumor detection. 4 Conventionally, T 2 is measured by sequentially acquiring several spin-echo (SE) images, each with a different echo time (TE) and subsequently fitting a mono-exponential decay. This is commonly acknowledged as a gold standard, despite residual diffusion effects affecting the T 2 quantification. Nonetheless, the long acquisition times (TAs) of around 50 min for whole-brain coverage in 2D is impractical for clinical applications. As an alternative approach, the multiple-echo SE (MESE) sequence in the Carr-Purcell-Meiboom-Gill condition 5 uses subsequent refocusing pulses to acquire multiple echoes for each excitation, reducing the total acquisition length. However, the quantitative accuracy of the sequence is compromised by imperfect refocusing. The imperfect refocusing results in the formation of stimulated (secondary) echoes that disrupt the T 2 decay of the primary SEs. 6 The effect of the stimulated echoes can be reduced by ignoring the first echo (because this pure SE biases the fitting of the subsequent stimulated-echo-contaminated data) while fitting the relaxation curve. 7 However, it has been shown that skipping echo approaches still yield highly variable results with errors that depend on flip angle, T 2 , and echo train length. 8 Therefore, more complex signal models are required to accurately estimate T 2 . 6 The ability to accurately and precisely map T 2 at high resolution (e.g., 1-mm 3 isotropic voxel size) in large volumes may help to improve the quantification of small focal changes such as multiple sclerosis lesions or brain areas causing focal seizure onset in epilepsy. However, using thin slices in a 2D MESE acquisition is challenging due to the reduced signal-to-noise ratio (SNR) and increasing difficulties to achieve a good slice profile. Furthermore, radio-frequency pulses with long duration are required to excite thin slices leading to longer echo spacings (ΔTE). 9 The increased ΔTEs aggravate and inhibit an accurate quantification of short T 2 values. In addition, true T 2 -weighting is difficult to obtain in reasonable TAs with 3D acquisition methods because a long repetition time (TR) is required to fully recover the magnetization to equilibrium between excitations. Because all spins in the field of view (FOV) of a 3D acquisition are excited by every pulse, the recovery time cannot be used for interleaved slice sampling, rendering the sequence less efficient. 10 Furthermore, a 3D acquisition is also limited by specific absorption rate safety constraints. A large amount of power is deposited if multiple nonselective 180° pulses are used which can easily exceed the specific absorption rate limits. 11 A fundamental consideration in any MRI experiment is to optimally balance image resolution, SNR, and TA. Various methods have been published to accelerate quantitative mapping methods, for example applying model-based reconstructions, 7,12-14 low-rank approaches, 15,16 or sparsity constraints 17,18 to highly undersampled acquisitions. However, these methods experience a low SNR and use low spatial resolution to counteract this effect. On the other hand, it has been shown that super-resolution (SR) reconstruction provides a better tradeoff between TA, spatial resolution, and SNR. 9,10 The earlier SR methods 19,20 focused on the improvement of the in-plane resolution of MR images. To achieve this in-plane resolution improvement, several images with a subpixel shifted FOV in the in-plane directions were acquired. However, the subpixel shift of the FOV in the in-plane directions correspond to a linear phase modulation in the k-space and does not acquire new frequency information. 21 Thus, the apparent improvement in the in-plane resolution in this case is due to an improvement in the SNR.
Later, most SR methods focused on decreasing the slice thickness and reaching voxel isotropy by exploiting the aliasing as a result of downsampling in the slice-select direction. 22 The resolution is thereby enhanced by acquiring multiple low-resolution images with a shift of the FOV in the slice direction, 10 3 orthogonal slice orientations 22 or rotated slice orientations, 23 and subsequently combining the images by solving a nonlinear inverse problem. In quantitative MRI, SR reconstruction benefits from combining the parametric model with the SR model. This has been shown in T 1 mapping 24 where the relaxation model was combined with the SR model allowing the direct estimation of a high-resolution T 1 map from low-resolution images. However, the acquisition of multiple low-resolution images results in long scan times (TA > 20 min), hence limiting its use for clinical applications.
Both the SR and model-based reconstructions are based on solving an inverse problem using iterative optimization methods. The present work investigates the combination of SR and model-based reconstruction to exploit the respective advantages (i.e., high-resolution and fast TA) in the application of quantitative T 2 mapping. The proposed method is based on multiple 2D MESE acquisitions, which were highly undersampled to compensate for prolonged scan time. The method was tested on a phantom and 4 healthy volunteers. An early version of this framework was presented at the Annual Meeting & Exhibition of the ISMRM in 2017. 25

| THEORY
SR reconstruction is a method to obtain a high-resolution image from a series of low-resolution images, where each low-resolution has a different FOV or orientation. Each FOV or orientation can be expressed as a different geometric transform T j (with j = 1, …, J, and J the number of transformations) from the high-resolution image to the low-resolution image. The resolution is enhanced because the different FOVs or orientations contain complementary information. The reconstructed high-resolution image on the other hand benefits from the high SNR of the lowresolution images, which are typically acquired with a high in-plane resolution and a low through-plane resolution, i.e., thick slices. Let x n (with n = 1, …, N, and N the number of SEs) represent a high-resolution T 2 -weighted series of vectorized images and y j,n,c (with c = 1, …, C, and C the number of coils) be the undersampled low-resolution k-space measurements acquired with a MESE sequence. The acquisition of the undersampled low-resolution y j,n,c , can be modeled by: where T j is the geometric transformation representing a rotation or translation of the FOV (Figure 1), ↓ is a downsampling operator that maps the high-resolution grid to the low-resolution grid, S j,c are the complex coil sensitivities calculated separately for each transformation, F is the forward discrete Fourier operator, and P is a binary undersampling mask. According to Gudbjartsson and Patz, 26 the noise η j,n,c in MRI can be assumed to be additive, white, and Gaussian when the SNR > 3, which is what we assume here. Using pertinent SR techniques, the highresolution image x n can be estimated by solving the inverse problem of Equation 1: where the choice of the least-squares criterion is motivated by the assumption that the noise is Gaussian distributed. The greater the number of acquisitions, J, with complementary information of the same object, the better the problem becomes conditioned. However, the problem in Equation 2 is often still ill-posed due to the downsampling operator, because several different intensity combinations in the high-resolution image can lead to the same intensity in the low-resolution image.
To further improve the conditioning of the problem, we can incorporate a signal model into the cost function. Given the task of T 2 quantification, it is natural to use the physical constraints imposed across the sequence of images, x n to further regularize the optimization. The signal model in its simplest form is typically described with a mono-exponential decay depending on the relaxation time map T 2 and proton density map M 0 : where the above equation represents the pixel-wise scalar nonlinear function for M 0 and T 2 and t n being the TE. In this work, the first echo is ignored to alleviate the effect of stimulated echoes. 27 However, even with ignoring the first echo, the standard exponential fitting will result in an overestimation of T 2 . 8 A more complex model that includes effects such as B + 1 , slice profiles, multiple T 2 compartments, diffusion and magnetization (1) y j,n,c = PF S j,c ↓ T j x n + j,n,c (2) F I G U R E 1 Illustration of the linear shifting (Left) and rotation (Right) experiment. The dashed lines within the boxes indicates the slice encoding direction, the frequency encoding direction is from head to toe. For the linear shifting each low-resolution scan (colors) is shifted by a known subpixel distance along the slice encoding direction (which in this case is sagittal). For illustration purposes, the shifting of the FOV (colored boxes) has been exaggerated. For rotation, each of the low-resolution scans are rotated (45°) across the common encoding direction. The red box represents sagittal section, green box represents coronal section, the blue and yellow boxes represent diagonal sections transfer is required to accurately fit the data. Because the present work aims at a proof of concept for SR for T 2 mapping, the common approach of ignoring the first echo is used due to its simplicity and easy implementation.
In the next step of the algorithm, the most intuitive step might be to substitute the high-resolution image x n in Equation 2 with the signal-model (Equation 3) and solving directly for M 0 and T 2 , resulting in the following cost-function: However, minimizing Equation 4 often requires special techniques, such as gradient scaling and repeated restarts of the optimization algorithm to achieve a fast convergence. 7,27 Therefore, in this work, we rather formulate the minimization problem as: where the first term ensures data-consistency of the highresolution image with the acquired data as in Equation 1 and an additional term ensures model-consistency, i.e., forces the signal intensities to decay exponentially across the different echoes N. In order to balance between the data-and model-consistency, a regularization parameter λ is introduced.
The above cost function in Equation 5 can be minimized using a split algorithm that estimates quantitative T 2 by minimizing data-and model-consistency terms alternately. 20 By fixing T 2 and M 0 , solving Equation 5 with respect to x n amounts to a standard linear least-squares problem with a closed form solution. On the other hand, for solving Equation 5 fixing variable x n corresponds to fitting a mono-exponential decay onto x n , intrinsically estimating T 2 and M 0 .

| Numerical simulation
Numerical noiseless T 2 and PD maps were generated from a segmentation of gray matter (GM), white matter (WM), and cerebrospinal fluid in a single axial slice (1 × 1 mm 2 resolution) of a numerical phantom. 28 For the 3 main tissues, the following T 2 values were used: 0.1 s for GM, 0.06 s for WM, and 2 s for cerebrospinal fluid. 29 From these maps, T 2 -weighted high-resolution images x n were simulated, each with equidistant TEs (TE = 10 …, 160 ms). The T 2 -weighted images were downsampled to the desired lowresolution using Equation 1 and different transformations T j . These different shifted images were then downsampled to low-resolution. Complex coil sensitivities for each of these transformations were simulated using the Parallel MRI Noisy Phantom Simulator. 30 These simulations exploit the fact that the reconstruction is separated into multiple 2D problems along the read-out direction (i.e., head-feet direction in this case). Therefore, for the phantom, a sagittal acquisition with phase-encoding in anterior-posterior and slice-encoding in left-right directions can be simulated. For the number of repetitions, the lowresolution single slice of the numerical phantom was shifted in the slice encoding direction (left-right). For the rotation, the FOV was rotated around the read-out (i.e., head-feet) axis for each repetition. Therefore, the applied transformations were not in-plane because the slice plane is sagittal/coronal but the numerical phantom is axial.
Three different types of transformations were tested: linear (where the images were shifted along the slice direction), 25 orthogonal (images were rotated to 2 orthogonal positions 0 and 90 degrees), and diagonal (2 orthogonal rotations plus 2 rotations at 45 and -45 degrees). An illustration of the different transformations is shown in Figure 1. The images were transformed to k-space and then 10-fold undersampled with a mixed parallel imaging and block undersampling pattern described in Hilbert et al. 31 The high-resolution images were reconstructed with the proposed SR model-based reconstruction. Root mean square error (RMSE) was calculated between ground truth and reconstructed T 2 maps for all the orientations.
In principle, we assume that adding more rotations (e.g., orthogonal with 2 rotations versus diagonal with 4 rotations) will increase the accuracy of the SR reconstruction but at the expense of increased TA. To that end, various numbers of rotations (2, 4, and 5) were tested in combination with different acceleration factors (6-fold, 10-fold, 14-fold) to determine a trade-off between accuracy and TA. Approximate TAs were calculated for all the simulations with the following parameters: 1 × 1 mm 2 resolution, matrix size 210 × 53, and TR = 5 s, 2 concatenations. The RMSE between the ground truth and the estimated T 2 map was calculated to quantify the reconstruction quality.

| Image acquisition
The proposed approach was tested on a standardized multipurpose phantom (E 38 19 195 K2130, Siemens) with 5 compartments of different concentration of MnCl2·4H2O in distilled water and 4 healthy subjects. Permission from the Institutional Review Board was obtained for all the in vivo imaging studies and written informed consent for the study and its publication was obtained from all participants prior to the experiments. The datasets were acquired with a standard 20-channel head/neck coil using a 10-fold undersampled GRAPPATINI 31 prototype sequence at 3T (MAGNETOM Skyra, Siemens Healthcare, Erlangen, Germany). The sequence implements a block undersampling pattern, 7 which samples the k-space in blocks (or segments) that are shifted across the echoes. A classic parallel imaging scheme with 2-fold acceleration is used additionally where only every other (phase-encoding) line inside each block is acquired. For the first dataset, 50 sagittal slices with 4 mm thickness and 1 × 1 mm 2 in-plane resolution were acquired. The FOV was moved by 1 mm in the slice-encoding direction with 4 shifts implementing the linear transformation. The total time of acquisition for this dataset was 16 min. For the second dataset, 60 slices each in 2 orthogonal orientations (sagittal and coronal) and 2 diagonal orientations (sagittal > coronal 45° and −45°) with the same slice thickness and in-plane resolution were acquired. The TA for this dataset was 18 min. Additional data were acquired in the same setup (i.e., orthogonal and diagonal) with 14-fold acceleration, which reduced the TA to 11 min. Of note, the dataset with orthogonal orientation was not acquired separately, but rather the 2 diagonal orientations were removed during reconstruction to create a third dataset. Spectral fat suppression was enabled, and other relevant acquisition parameters were: TR = 5.4 s, 16 echoes, ΔTE = 10 ms. For reference T 2 values, a fully sampled MESE dataset with 16 echoes and 29 slices and a resolution of 1.1 × 1.1 × 4 mm 3 was acquired. In addition, for the phantom, a conventional single-slice single-echo SE sequence with 16 TEs (10, 20, …, 160 ms) was acquired as reference. It should be noted that the single-echo SE sequence is not a gold standard because it is still affected by diffusion, but is commonly used as a reference acquisition. 32,33 Detailed overview of the acquisition parameters of these datasets can be found in Table 1.

| Image reconstruction
The reconstruction was implemented using Matlab (MATLAB2017a, The Mathworks Inc., Natick, MA). First, GRAPPA 34 was used to fill the missing lines in each block of the acquired k-space according to Hilbert et al. 31 Using the k-space samples available at different TEs, composite fully sampled images were reconstructed and upsampled to the high-resolution grid. Subsequently, an initial guess of a T 2 and M 0 map was estimated by nonlinear least-square fitting of these magnitude high-resolution images. This estimate of the maps was then used in the algorithm for SR T 2 estimation by alternately solving Equation 5 fixing T 2 and M 0 first as described in the theory section. Complex coil sensitivities were approximated for each of the transformations by summing k-spaces across echoes. The k-space is then transformed into coil images by inverse Fourier transform and then divided by the sum of squares to obtain coil sensitivities. 35 The regularization parameter was heuristically set to λ = 1, establishing equal contributions from model-and data-consistency in the cost function. The reconstruction was performed iteratively until no improvement was achieved in the previous 2 consecutive iterations: where X (i) was the reconstructed signal series at i th iteration and ε = 1e-04 . In addition, the maximum number of iterations was fixed to 30 after it was experimentally tested at what iteration the algorithm typically convergences (results not shown). The reconstruction scheme is illustrated in Figure 2.

| Validation
To evaluate the performance of the proposed SR modelbased method, T 2 maps were qualitatively compared to T 2 maps from 2 other methods using the numerical phantom and 1 in vivo dataset. The 2 approaches were as follows: • Low-resolution model-based + SR reconstruction: It is pertinent to assess whether minimizing the cost function jointly contributes to the better reconstruction. For this purpose, both the data consistency and model consistency were enforced 1 after the other rather in a joint cost function. First, low-resolution T 2 maps were reconstructed from individual low-resolution orientations using model-based reconstruction, followed by a SR algorithm to upsample to a high-resolution T 2 map. • SR reconstruction only: As described in the theory section, the signal model is used as a regularizer to better condition the minimization problem. To assess the significance of the signal model, the cost function was minimized without the model consistency term, which is equivalent to classical SR approach.
To ascertain the accuracy of the calculated T 2 values, a region of interest (ROI) analysis was performed on the phantom and in vivo T 2 maps for both accelerated and fully sampled datasets. T 2 values from different compartments of the phantom were compared with the T 2 values obtained from the conventional nonlinear least-squares fitting of the fully sampled MESE and gold standard SE data. The same comparison was made for the healthy brain where ROIs were drawn in the frontal WM, deep GM (putamen and caudate nucleus), and corpus callosum. ROI labeling and segmentation was conducted using ITK-SNAP 36 (www.itksn ap.org). The segmentation was performed in the native space of the fully sampled T 2 w images. Subsequently, the fully sampled T 2 w image was pairwise registered rigidly onto the T 2 maps of the other datasets using Elastix. 37 The resulting transformation was then applied to the label map of the ROIs to provide a similar segmentation among datasets. The datasets of the 4 healthy volunteers who were scanned with both a 10-and 14-fold

| Numerical phantom
T 2 maps and images obtained from the SR model-based reconstruction in the numerical phantom using different orientations (linear, orthogonal, and diagonal) are shown in Figure 3 in comparison to the ground truth. The diagonal rotations were better able to resolve the brain structures as compared to the linear and orthogonal datasets. This can be better appreciated in the zoomed images of the same area that demonstrates the diagonal rotation results in better delineation of the small structures in brain (Figure 3, bottom). The comparison between the 3 orientations in reconstructed T 2 -weighted images also shows sharper edges in the highresolution images when using the diagonal orientation. The trade-off between the number of rotations and the acceleration is demonstrated in Figure 4. The numerical simulation demonstrated that adding more rotations results in reduced error around the edges for all the acceleration factors (from left to right). At the same time, an increased acceleration factor influences the T 2 accuracy, i.e., an increased difference from the ground truth (top to bottom). Considering the prolonged TA with the number of rotations, the best tradeoff is a 10-fold accelerated scan with 4 rotations (RMSE = 7.8 ms, TA = 14 min). Five rotations and 6-fold acceleration showed the least error (RMSE = 3.2 ms) but required an TA of 29 min. The 14-fold acceleration with 5 rotations had an TA of 10 min but an RMSE of 12.5 ms. The details of the RMSE for all the accelerations and rotations are provided in Supporting Information Table S1, which is available online. Figure 5 shows the reconstructed T 2 maps and T 2 -weighted images from 1 low-resolution acquisition, and the SR modelbased reconstructed maps for linear, orthogonal, and diagonal orientations. The low-resolution T 2 maps and weighted images ( Figure 5A) show partial-volume effects within the compartments of the phantom and blurring of the fine structures in the brain. Estimating T 2 maps with SR model-based reconstruction enhances the spatial resolution, reducing the partial-volume effect. However, with linear shifting, there is F I G U R E 3 SR T 2 maps (upper row), reconstructed T 2 -weighted images (middle row), and the zoomed-in T 2 maps (bottom row) for still blurring of the edges around the compartments of the phantom and brain structures ( Figure 5B). For orthogonal rotations, the compartment edges in the phantoms appear sharper in both orthogonal directions ( Figure 5C, white arrows), whereas the diagonal orientation improves the sharpness further, especially in the diagonal direction ( Figure 5D, black arrows). The same can be observed for the human brain where diagonal rotation improved the spatial resolution and delineation of the brain structures. Figure 6 shows the axial and coronal views of the reconstructed T 2 -weighted images from a 10-fold accelerated low-resolution dataset with diagonal orientations and the fully sampled MESE data. It is evident that the SR-T 2 mapping allows satisfactory 3D visualization due to isotropic resolution as compared to conventional MESE, which inherently has low through-plane resolution. This can be appreciated even more from the zooms shown in Figure 6. A transversal zoom on the caudate nucleus head and the putamen is shown where the interface between different structures is well delineated in both SR-T 2 -weighted and MESE-T 2 -weighted images. However, for the coronal zoom, the interface between GM and WM is better defined for SR-T 2 -weighted image.

| Validation
A comparison between the SR model-based reconstruction and the SR-only and low-resolution model-based + SR reconstruction is shown in Figure 7. For numerical phantom, nRMSEs showed that the proposed approach outperforms the SR-only and model-based + SR reconstruction. The nRMSE for SR only reconstruction was 0.50, for modelbased + SR reconstruction was 0.23 and for proposed reconstruction was 0.12. The SR-only reconstruction showed visible artifacts due to the undersampling as no prior information from the signal model (mono-exponential decay) was incorporated in the reconstruction. The low-resolution model-based + SR reconstruction showed improvement in the reconstruction; however, blurring around the edges is evident in both numerical phantom and in vivo data. In comparison, SR model-based reconstruction demonstrated that combining SR and model knowledge jointly in 1 cost function improves the reconstruction.

| ROI analysis
The ROI analysis of the phantom compartments revealed that at shorter T 2 , 10-fold accelerated data were comparable F I G U R E 4 The gain in the quality of SR-T 2 reconstruction as a function of number of rotations and TA in a numerical phantom. The corresponding TAs are mentioned with the T 2 maps. Adding more rotations results in reduced error around the edges for all acceleration factors (from left to right). At the same time, the acceleration factor affects the quality of the reconstruction (top to bottom), especially for high T 2 values in the cerebral spinal fluid. Considering the increase in TA with the number of rotations, it can be deduced that the optimal configuration is 10-fold with 4 rotations with the fully sampled MESE ( Figure 8A). For compartment 1, the mean T 2 value for 10-fold accelerated data was 31.1.0 ± 1.6 ms, for 14-fold it was 29.25 ± 1.7 ms, whereas the T 2 value for the MESE was 25.41 ± 1.5 ms and for SE was 18.22 ± 0.5 ms. For compartment 2, the mean T 2 value for 10-fold, 14-fold, MESE and the SE was 29.19 ± 0.8 ms,   overestimation is well known and caused by the stimulated echoes formed in the MESE sequence. 6 For the in vivo data, the T 2 values in different regions of the brain showed good agreement with the fully sampled MESE acquisition ( Figure 9B). For frontal WM, the mean T 2 value for the 10-fold acceleration was 80.0 ± 1.3 ms, for 14-fold it was 76.5 ± 5.1, and for MESE fully sampled it was 83.0 ± 2.1 ms. For putamen, the mean T 2 value for 10-fold, 14-fold, and MESE fully sampled was 69.8 ± 1.2, 66.2 ± 1.3, and 73.0 ± 3.1 ms, respectively. In the corpus callosum (MESE T2 = 94 ± 2 ms), the 14-fold acceleration showed increased difference (80 ± 2.4 ms) as compared to 10-fold (87 ± 1.8 ms). The T 2 value for the caudate nucleus with the 10-fold acceleration (74.8 ± 2.0 ms) was closest to the MESE acquisition (75 ±1.8), whereas the T 2 at 14-fold acceleration was 67 ± 2.8 ms. Comparing the mean T 2 values for 10-fold accelerated data from all the 4 ROIs across volunteers demonstrated consistent values with some natural variation ( Figure 9C). The average T 2 value across all volunteers were 70.8 ± 8.5 ms for the frontal WM, 81.3 ± 10.5 ms for the putamen, 79.8 ± 9.5 ms for the corpus callosum, and 79.66 ± 13.6 ms for the caudate nucleus.

| DISCUSSION
The proposed work aims to address the challenges associated with high-resolution T 2 mapping with the widely used MESE sequence. The method combines SR and model-based reconstruction for T 2 estimation into 1 integrated approach, enabling the direct estimation of an isotropic high-resolution T 2 map from a set of anisotropic, undersampled, low-resolution k-spaces.
Numerical simulations demonstrated that rotating the slice orientations around the phase-encoding results in better image quality of the SR reconstruction than linear sub-pixel shifts in the slice-encoding direction. These results were confirmed by the phantom and the in vivo acquisitions where the rotation resulted in much sharper images. Low-resolution images with different slice orientations appear to better sample different spatial information than shifting the low-resolution images by subpixel distances along the slice selection direction. 9 This is because acquiring low-frequency information in the slice-selective direction limits the additional information available in the subpixel shift. In contrast, the different orientations populate the 3D k-space much better (similar F I G U R E 9 A, Different ROIs drawn in axial, coronal, and sagittal sections of a model-based SR T2w (TE = 80 ms) image. B, The bar chart represents the mean T 2 values and standard deviation for accelerated and fully sampled dataset in 4 different regions of brain. C, The mean and standard deviation T 2 values from model-based SR reconstruction (diagonal orientations, 10-fold acceleration) found in all subjects grouped by brain structures, each bar representing a subject to PROPELLER 38 sequences). The optimal number of slice orientations ensures that the low-resolution images contain enough spatial information to reconstruct the high-resolution image without extensively prolonging the TA. The comparison between number of rotations and acceleration factors showed that 4 slice orientations (0°, 90°, 45°, −45°) with a 10-fold acceleration was the optimum in this application with an TA of 18 min.
The advantage of combining SR with model-based reconstruction has already been demonstrated for T 1 mapping. 24 This work was performed in the image domain without k-space undersampling, whereas the proposed approach addresses a more complex problem by reconstructing quantitative maps from undersampled k-space data. It is important to ascertain whether model-based SR-T2 mapping performs better than the sequential application of model-based and SR reconstructions. The comparison of the proposed approach to the SR reconstruction without prior model and to the low-resolution model-based reconstruction followed by SR reconstruction demonstrated that incorporating the model information in the reconstruction yields improved results.
Comparing a T 2 map estimated with model-based SR to a T 2 map estimated with a conventional nonlinear least-square fitting on a fully sampled MESE dataset proved that T2 values are comparable between the 2 methods. For the phantom, the compartments with longer T 2 had a higher bias as compared to shorter T 2 compartments, which may be attributed to an imperfect signal model. One contributing factor, as already reported in the literature, 6 is that theT 2 values derived from MESE-type acquisitions are overestimated in comparison to single-echo SE acquisitions due to stimulated-echo signal contributions, which mostly stem from an imperfect slice profile of the refocusing pulses. This can be addressed by using a different signal model which accounts for slice profile shape and B1 inhomogeneity, e.g., by using an analytical 6,27 or numerical 33,39,40 signal model. The current model of the proposed method does not account for the slice profile. It is assumed that incorporating it in the reconstruction will not only allow addressing stimulated-echo issues but also lead to a more realistic downsampling/upsampling operator between high-resolution and low-resolution images. Because this more accurate model will potentially further improve the spatial resolution, the slice profile is subject of future investigation.
An alternative formulation of the reconstruction problem proposed here can be as follows: the signal in a low-resolution voxel may consist of different tissue types, e.g. at a WM/GM boundary. Neglecting microstructural T 2 components, this partial volume effect yields a multi-exponential signal decay. The fitting of a mono-exponential decay in a high-resolution voxel corresponds to a multi-exponential fitting in the lowresolution image space due to the downsampling operator.
The resulting ill-posed problem becomes better conditioned because the different orientations provide selective highresolution information, which allows us to separate out the multi-exponential components within a voxel.
Currently, the optimal-quality 10-fold accelerated acquisition lasts 18 min. Despite the high acceleration, the TA is too long to be acceptable for clinical routine or clinical studies. A higher acceleration factor (14-fold) was tested which reduced the TA to 11 min but decreased the sharpness of brain structures. Nevertheless, this long TA is still too long for clinical applications where a typical sequence requires 3 min. Various other techniques can be applied to further accelerate the TA, such as simultaneous multi-slice. 41,42 Moreover, the same sampling pattern was used in each slice and orientation. Varying the sampling across slices and orientations may improve the reconstructed T 2 maps and will be further investigated.
SR approaches require several acquisitions with different orientations, and the reconstruction method assumes that each voxel corresponds to the same anatomical location in all the orientations. Subject motion and image distortions (e.g. due to imperfect gradient performance) violate this assumption and may not perfectly align the multiple orientations for model fitting. This model violation will lead to image artifacts and wrong T 2 values which is a limitation of the proposed method. In the frame of this study, motion correction schemes were not explored, but could be used to facilitate a better alignment of different orientations and improve the reconstruction. The area of motion correction and compensation has been well studied specifically for SR reconstruction in fetal imaging. 43,44 Future work will aim to incorporate these strategies to further improve the robustness of the SR-T 2 mapping. In addition, incorporating a spatial regularization of the T 2 maps will potentially further improve the reconstruction.

| CONCLUSIONS
We propose a technique that uses four 10-fold undersampled low-resolution MESE 2D acquisitions to iteratively reconstruct a high-resolution T 2 map. The proposed technique enables high-resolution 1 mm 3 isotropic whole-brain T 2 mapping in 18 min. The proposed technique may allow the assessment of T 2 values in small brain structures valuable for the search of imaging biomarkers in the future.