Fang Cao

Olivier Commowick

Elise Bannier

Christian Barillot

Simultaneous Estimation of T 1 , T 2 and B 1 Maps From Relaxometry MR Sequences

Interest in quantitative MRI and relaxometry imaging is rapidly increasing to enable the discovery of new MRI disease imaging biomarkers. While DESPOT1 is a robust method for rapid whole-brain voxel-wise mapping of the longitudinal relaxation time (T1), the approach is inherently sensitive to inaccuracies in the transmitted flip angles, defined by the B1 inhomogeneity field, which become more severe at high field strengths (e.g., 3T). We propose a new approach for simultaneously mapping the B1 field, M0 (proton density), T1 and T2 relaxation times based on regular fast T1 and T2 relaxometry sequences. The new method is based on the intrinsic correlation between the T1 and T2 relaxometry sequences to jointly estimate all maps. It requires no additional sequence for the B1 correction. We evaluated our proposed algorithm on simulated and in-vivo data at 3T, demonstrating its improved accuracy with respect to regular separate estimation methods.

1I n t r o d u c t i o n

Quantitative MRI is becoming more and more important to study the brain micro-structure and to provide calibrated measures of brain tissue properties from MRI that are crucial to obtain MR imaging biomarkers. The gold standard methods for longitudinal and transverse relaxation times (T 1 and T 2 )m a p p i n g s require long acquisition times that are often not applicable for clinical use. Many methods have been proposed to speed up the scanning procedures.

The driven-equilibrium single-pulse observation of T 1 (DESPOT1) [START_REF] Deoni | Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state[END_REF] and the Carr-Purcell-Meiboom-Gill (CPMG) multi-echo sequence for measuring T 2 [START_REF] Hennig | Multiecho imaging sequences with low refocusing flip angles[END_REF] have gained significant popularity over the past decades due to their superior time efficiency, allowing fast T 1 and T 2 mapping with high resolution and large spatial coverage. However, both DESPOT1 and CPMG methods suffer from a strong dependence on an accurate knowledge of the Flip Angle (FA), which leads to significant errors in the presence of the non uniform radiofrequency excitation field B 1 [START_REF] Collins | SAR and B1 field distributions in a heterogeneous human head model within a birdcage coil[END_REF]. This is particularly problematic at high field strengths, where B 1 variations can be over 20% and in situ FA correction is necessary [START_REF] Collins | SAR and B1 field distributions in a heterogeneous human head model within a birdcage coil[END_REF][START_REF] Dowell | Fast, accurate, and precise mapping of the RF field in vivo using the 180 degrees signal null[END_REF]. It is therefore commonly accepted that DESPOT1 and CPMG at 3T or higher field strengths need to be combined with an appropriate B 1 correction method [START_REF] Parker | Accurate multislice gradient echo T1 measurement in the presence of non-ideal RF pulse shape and RF field nonuniformity[END_REF][START_REF] Lebel | Transverse relaxometry with stimulated echo compensation[END_REF], especially for T 1 estimation [START_REF] Deoni | High-resolution T1 mapping of the brain at 3T with driven equilibrium single pulse observation of T1 with high-speed incorporation of RF field inhomogeneities (DESPOT1-HIFI)[END_REF]. Different prospective methods have been presented to compensate for B 1 inhomogeneities by acquiring additional or new sequences. Parker et al. proposed a methodology for estimating T 1 using a two-point technique with a standard multislice gradient echo 2F . C a o et. al.

sequence [START_REF] Parker | Accurate multislice gradient echo T1 measurement in the presence of non-ideal RF pulse shape and RF field nonuniformity[END_REF]. Treier compensated B 1 field inhomogeneities in T 1 mapping by performing an additional measurement using an optimized fast B 1 mapping technique [START_REF] Treier | Optimized and combined T1 and B1 mapping technique for fast and accurate T1 quantification in contrast-enhanced abdominal MRI[END_REF]. Deoni presented a method that combines the DESPOT1 data with at least one inversion-prepared SPGR data to obtain an unique solution for the B 1 map [START_REF] Deoni | High-resolution T1 mapping of the brain at 3T with driven equilibrium single pulse observation of T1 with high-speed incorporation of RF field inhomogeneities (DESPOT1-HIFI)[END_REF]. Dowell et al. proposed a B 1 mapping sequence using spoiled gradient echo (SPGR) and 180 • signal null [START_REF] Dowell | Fast, accurate, and precise mapping of the RF field in vivo using the 180 degrees signal null[END_REF]. Yarnykh developed a fast B 1 mapping sequence that consists of two identical RF pulses followed by two delays of different durations [START_REF] Yarnykh | Actual flip-angle imaging in the pulsed steady state: A method for rapid three-dimensional mapping of the transmitted radiofrequency field[END_REF].

However, since B 1 map acquisition takes up precious scanning time and most retrospective studies do not include B 1 field mapping sequences, we wish to keep the data acquisition protocol of DESPOT1 and CPMG multi-echo sequences. We therefore introduce a new combined algorithm to estimate T 1 and T 2 relaxation times, and B 1 inhomogeneities without requiring additional B 1 acquisitions. The estimation of all maps is formulated as an energy minimization problem on the whole image without any prior knowledge of the inner structure of the subject.

In Sect. 2, we review the regular estimation methods and propose our new algorithm for joint estimation with B 1 correction. We validate our method using synthetic phantom and in vivo data sets in Sect. 3. The experimental results are given in Sect. 4 to demonstrate the effectiveness of the proposed method.

2M e t h o d s

Regular T 1 and T 2 Estimation

Following the DESPOT1 algorithm [START_REF] Deoni | Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state[END_REF], traditional T 1 estimation uses two SPGR sequences. The MR signals measured at two flip angles α i are modeled as

I † αi (T 1 ,M 0,T1 )= M 0,T1 (1 -exp(-TR/T 1 )) sin α i 1 -exp(-TR/T 1 ) cos α i ,i=1, 2 (1)
where T 1 is the longitudinal relaxation time, M 0,T1 is proportional to the proton density M 0 and TR is the repetition time. T 1 and M 0,T1 are optimized through a least squares problem to fit I † αi (T 1 ,M 0,T1 ) to the observed signals S αi .S u c hap r o b l e mh a sa na n a l y t i c a ls o l u t i o nf o rM 0,T1 and T 1 [START_REF] Deoni | Rapid combined T1 and T2 mapping using gradient recalled acquisition in the steady state[END_REF]. Then, the traditional T 2 estimation is performed by minimizing P n j=1 (S TEj -I † TEj ) 2 [START_REF] Tofts | Quantitative MRI of the Brain: Measuring Changes Caused by Disease[END_REF] where S TEj (j =1 ,...,n, n being the echo train length) are the acquired signals measured at TE j ,t h ej-th echo time in the CPMG sequence, and I † TEj is the corresponding simulated signal of the j-th echo. We use the Extended Phase Graph (EPG) algorithm [START_REF] Hennig | Multiecho imaging sequences with low refocusing flip angles[END_REF][START_REF] Lebel | Transverse relaxometry with stimulated echo compensation[END_REF] to model the simulated signals. This algorithm tracks the multiple coherence pathways of spins after consecutive periods that model precession, relaxation and refocusing, and has been used to calculate the echo decay curves in CPMG sequence yielding precise monoexponential T 2 quantification [START_REF] Lebel | Transverse relaxometry with stimulated echo compensation[END_REF][START_REF] Prasloski | Applications of stimulated echo correction to multicomponent T2 analysis[END_REF][START_REF] Layton | Modelling and estimation of multicomponent T2 distributions[END_REF]. Rather than computing the time evolution of magnetization vectors, the EPG algorithm uses the magnetization phase state vector

F =[ F 1 ,F ⇤ 1 ,Z 1 ,F 2 ,F ⇤ 2 ,Z 2 , ••• ,F n ,F ⇤ n ,Z n ]
T to describe the spin system at a given time. The states F and F ⇤ refer to the dephasing and rephasing of the spins in the transverse plane. The states Z refer to the spins along the longitudinal axis that maintain their phase. For the j-th echo, F j is calculated by

F j = E(T 1 ,T 2 )TR j (α)E(T 1 ,T 2 )F j-1 ,j=1, ••• ,n (2)
where F 0 is proportional to the proton density M 0 . The matrix R j represents the effects of RF refocusing, depending on the nominal flip angle α of the refocusing pulse for each echo train. The matrix E represents T 1 and T 2 relaxation between successive echoes. The matrix T represents the phase transitions for each echo train. The amplitude I † TEj of the j-th echo is given by the first element of F j . More details on the EPG algorithm can be found in [START_REF] Lebel | Transverse relaxometry with stimulated echo compensation[END_REF][START_REF] Prasloski | Applications of stimulated echo correction to multicomponent T2 analysis[END_REF][START_REF] Layton | Modelling and estimation of multicomponent T2 distributions[END_REF].

Simultaneous Relaxation Times and B 1 Field Estimation

The traditional DESPOT1 algorithm and T 2 relaxometry estimation provide a solution for T 1 , T 2 and M 0 estimation. However, this solution lacks accuracy in the presence of inhomogeneities in the B 1 transmission field. In this study, we use a combined model of DESPOT1 and EPG algorithms to estimate T 1 , T 2 , M 0 and B 1 simultaneously. This joint model, thanks to the intrinsic dependence between the T 1 and T 2 relaxometry sequences, makes it possible to fit both T 1 and T 2 relaxometry signal decays (intensities) at the same time as B 1 .W ef o r m u l a t e the estimation as a single energy problem depending on four parameters per voxel (T 1 ,T 2 ,M 0 ,B 1 ) plus one global parameter k over the image. Our approach amounts to minimizing the following energy function over the entire image space V :

E = λE r (B 1 )+ Z V 2 X i=1 (S αi -kI αi (T 1 ,M 0 ,B 1 )) 2 + n X j=1 S TEj -I TEj (T 1 ,T 2 ,M 0 ,B 1 ) 2
(3) where B 1 represents the field inhomogeneities under the assumption that the transmission field is a linear system. The coefficient k is a voxel-independent scale factor, accounting for the global change of M 0 between T 1 and T 2 relaxometry sequences. This assumption of a stationary k factor is reasonable when considering that the time between T 1 and T 2 relaxometry sequences acquisitions is short. There is therefore no large variation arising from transmitter and receiver effects between these two relaxometry sequences and the k values originated from different voxels are approximately constant in the whole brain.

This energy function is made of three terms. The first one in the integral depends on the T 1 relaxometry sequences with a modified signal equation I αi with α i accounting for B 1 inhomogeneities. This is done by replacing the flip angles α in T 1 and T 2 relaxometry signals (Sect. 2.1) by α = B 1 α nomi where α nomi is the nominal value of the flip angle. The second term in the integral is related to the T 2 relaxometry sequences. Again, we utilize the EPG expression for the signal equation I TEj at the j-th echo, with the flip angle α accounting for B 1 .F i n a l l y ,t h et e r mE r introduces an L 2 regularization of the B 1 field to account for the fact that those inhomogeneities vary slowly in space.

4F . C a o et. al.

Proposed Algorithm: Parameters Estimation

It is worth noting that the energy function E is convex with regard to k, M 0 and T 2 respectively. To simplify the optimization process, we estimate the parameters of this function using an alternated minimization process, as follows:

1. Initialize T 1 , T 2 and M 0 from regular estimation (Sect. 2.1). B 1 is set to 1. 2. Minimize E alternately (fixing other parameters) with respect to k, T 2 , M 0 , B 1 and T 1 . 3. Regularize B 1 with a Gaussian kernel. 4. Go to step 2 until convergence has been reached.

Estimating k amounts to solving a least squares problem between T 1 relaxometry signals and the simulated ones from current values of the parameters. However, to account for outliers in the estimation caused by voxels outside of the brain (vessels, bones), we estimate k using Tukey's Bisquare M-estimator [START_REF] Tukey | Exploratory Data Analysis[END_REF]. Similarly, M 0 has an analytical solution at each voxel. Estimating the other parameters is more complicated as derivatives are not easily computed. We therefore perform the alternated minimizations using the BOBYQA algorithm for bounded optimization without gradient [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF] subject to constraints extracted from the literature: M 0 ≥ 0, 0 ≤ T 1 ≤ 5000 ms, 0 ≤ T 2 ≤ 1000 ms and 0.2 ≤ B 1 ≤ 2 [START_REF] Gelman | MR imaging of human brain at 3T: preliminary report on transverse relaxation rates and relation to estimated iron content[END_REF][START_REF] Stanisz | T1, T2 relaxation and magnetization transfer in tissue at 3T[END_REF][START_REF] Dowell | Fast, accurate, and precise mapping of the RF field in vivo using the 180 degrees signal null[END_REF]. We regularize B 1 as a separate smoothing step to account for the B 1 inhomogeneities in space.

3E x p e r i m e n t a l D e s i g n f o r V a l i d a t i o n

The validation process included synthetic phantom and in vivo measurements. Both synthetic phantom and in vivo tests used the following parameters to simulate/acquire the relaxometry sequences: (1) for the T 1 relaxometry: voxel size: 1.

Synthetic Phantom

We generated a synthetic phantom to quantitatively evaluate the results of T 1 and T 2 estimations. The phantom was made of three types of tissues organized in concentric spheres, which are cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM). All reference values are taken as similar to values in the human brain1 (Tab. 1). To demonstrate the capability to retrieve realistic B 1 maps, B 1 was set to different values: smaller than 1 for GM and higher than 1 for CSF and WM. Since B 1 is generally high at the center of the scene [START_REF] Dowell | Fast, accurate, and precise mapping of the RF field in vivo using the 180 degrees signal null[END_REF], we set B 1 for CSF to be the largest among all tissues. The reference value of k is set to be 7. The reference relaxation maps (see Fig. 1) were used to feed into an MRI simulator2 in order to synthesize T 1 and T 2 relaxometry sequences following signal equations from Sect. 2.1 with α accounting for B 1 . Then we added Rician noise to the simulated relaxometry sequences (SNR= 20dB). The estimation was run with the regular DESPOT1 and EPG estimations (Sect. 2.1) and with the proposed method. The results were then compared to the reference values to quantitatively evaluate our algorithm.

In Vivo Study

In vivo validation was performed on a high field MRI scanner (3T Siemens Verio VB17), which has significant effects of B 1 field inhomogeneities. Wholebrain MR images were acquired on 13 healthy subjects (6 male, 7 female, mean age= 29.2 ± 17.8 y.o.). To comp ensate for b etween-scans sub ject motion, a sixparameter rigid-body registration was carried out for each subject based on normalized mutual information. We acquired T1-w SE and T2-w/PD-w MR sequences for each subject and considered these real acquisitions as references for in vivo validation. Then the estimated maps (with our method and with the methods described in Sect. 2.1) were used to feed into the MRI simulator 2 to generate T1-w, T2-w and PD-w images with the same sequences and parameters as the real acquisitions. Note that for all simulations, we applied bias field correction of the M 0 map as it is sensitive to the coil sensitivity. This correction is based on a unified segmentation [START_REF] Ashburner | Unified segmentation[END_REF]. To quantitatively evaluate the performance of the algorithms, we use the squared correlation coefficient R,computed between the real and the simulated images, as an indicator of the estimation accuracy. Assuming the simulations are perfect, the closer the R values are to 1, the better the map estimations.

4R e s u l t s

Synthetic Phantom

We present in Tab. 1 the comparison b etween the reference and estimated values for T 1 , T 2 , B 1 and M 0 using the two evaluated methods. It can be seen that the DESPOT1 algorithm is sensitive to B 1 inhomogeneities due to its inability to estimate those values. The estimated T 1 and M 0 therefore vary significantly from the reference values. However, when we applied the proposed algorithm, the estimated B 1 , T 1 , T 2 and M 0 are in close agreement with the reference values. The reference values are always included within one standard deviation of the estimated values using the proposed method.

In Vivo Results

Fig. 2 presents the maps estimated using the regular method and the proposed algorithm. This figure confirms visually that the regular algorithm is sensitive to B 1 inhomogeneities, especially seen on T 1 maps where the center values are much higher for the same tissues. On the contrary, the proposed method removes the influence of B 1 inhomogeneities and the obtained maps are much more uniform. Fig. 3 presents the T1-w images from the real acquisition and from the simulation using the regular method and the proposed algorithm. The results of T2-w and PD-w simulations were computed for quantitative evaluation but are not presented as all simulations look very similar to the real acquisitions. On this figure, it is clear that the regular method's inability to estimate B 1 results in degraded simulations for T1-w sequences. These observations are confirmed by the R coefficients in Tab. 2. All values are significantly increased (paired t-test, p<10 -6 for each of the simulated sequences) between the regular and proposed methods. Moreover, the correlation coefficients for T1-w images increase over 50% using the proposed method (0.653 compared to 0.327).

5C o n c l u s i o n

We have prop osed a new approach for simultaneous mapping of B 1 inhomogeneity field, T 1 and T 2 relaxation times, and proton density M 0 . The method is based on the intrinsic correlation between the T 1 and T 2 relaxometry sequences. We use a combined model of DESPOT1 and EPG algorithms to estimate T 1 , T 2 , M 0 and (e) B1 (f) T1 (g) T2 (h) M0 Fig. 2. Estimated maps on one healthy subject. Row 1 shows the maps using the regular method, row 2 using the proposed method. From left to right for each row: B1, T1, T2 and M0 maps respectively. The same window level is set for each column. Our experiments on simulated data demonstrated that the proposed method is able to accurately estimate T 1 , T 2 , M 0 and B 1 maps. Experiments on in vivo data showed high correlation (up to 50 % increase) between the weighted acquisitions and the simulated sequences, which confirms the great potential of the proposed method to handle clinical applications where quantitative MRI has been shown to be highly relevant (MS, stroke, pediatrics, ...).

3 × 1 . 3 ×

 313 3 mm 3 , TR = 15 ms, TE =1 .54 ms, FA: 5 • and 30 • ;( 2 ) for the T 2 relaxometry: voxel size: 1.3 × 1.3 × 3 mm 3 , TR = 4530 ms, TE = [13.8, 27.6, 41.4, 55.2, 69.0, 82.8, 96.6] ms.

Fig. 1 .

 1 Fig. 1. The synthetic phantom. From left to right: B1, T1, T2 and M0 reference maps.

  6F . C a o et. al.

Fig. 3 .Table 2 .

 32 Fig. 3. Simulated T1-w images and the corresponding real acquisition on one healthy subject. From left to right: real acquisition, simulated images using the regular and the proposed methods respectively. T1-w Ave ( R) T1-w Std(R) T2-w Ave ( R) T2-w Std(R) PD-w Ave ( R) PD-w Std(R) Regular 0.327 0.092 0.894 0.018 0.729 0.051 Proposed 0.653 0.107 0.923 0.016 0.834 0.026 Table 2. Ave r a g e c o r r e l a t i o n c o e ffi c i e nt R between the simulated images and the real acquisitions on 13 healthy subjects, for T1-w, T2-w and PD-w images. Row 1 shows R for the regular method, row 2 for the proposed algorithm. Differences are significant (paired t-test, p<10 -6 ).

B 1

 1 simultaneously. This combination requires no additional sequence for the B 1 correction, making use of the traditional T 1 and T 2 relaxometry sequences only. 8F . C a o et. al.

Table 1 .

 1 Tissue Method Ave(B 1 )Std(B 1 ) Ave(T 1 )Std(T 1 ) Ave(T 2 )Std(T 2 ) Ave(M 0 )Std(M 0 ) Statistical values of B1, T1 (ms), T2 (ms) and M0 on the synthetic phantom. Ave a n d S t d a r e t h e m e a n a n d s t a n d a r d d e v i a t i o n . T h e e s t i m a t e d k is 6.94.

		Reference	1.10	0.00 500.00	0.00 70.00	0.00	80.00	0.00
	WM	Regular	-	-608.06	0.71 71.31	1.97	77.55	1.59
		Proposed	1.10	0.00 495.97 15.71 69.84	1.69 80.32	1.13
		Reference	0.80	0.00 830.00	0.00 80.00	0.00	90.00	0.00
	GM	Regular	-	-526.88	0.87 87.26	2.47	79.10	1.47
		Proposed	0.80	0.00 823.50 60.30 79.91	2.57 90.22	3.12
		Reference	1.30	0.00 2500.00	0.00 330.00	0.00 100.00	0.00
	CSF	Regular	-	-4364.33 449.69 476.14 48.61	75.95	1.12
		Proposed	1.30	0.00 2485.96 193.76 329.73 27.13 100.50	3.33

The T1, T2 and M0 values are given by http://www.bic.mni.mcgill.ca/brainweb/ MICCAI Workshop -IntellMR

SimuBloch v0.3 at http://vip.creatis.insa-lyon.fr