Optimal Configuration for Relaxation Times Estimation in Complex Spin Echo Imaging

Many pathologies can be identified by evaluating differences raised in the physical parameters of involved tissues. In a Magnetic Resonance Imaging (MRI) framework, spin-lattice T1 and spin-spin T2 relaxation time parameters play a major role in such an identification. In this manuscript, a theoretical study related to the evaluation of the achievable performances in the estimation of relaxation times in MRI is proposed. After a discussion about the considered acquisition model, an analysis on the ideal imaging acquisition parameters in the case of spin echo sequences, i.e., echo and repetition times, is conducted. In particular, the aim of the manuscript consists in providing an empirical rule for optimal imaging parameter identification with respect to the tissues under investigation. Theoretical results are validated on different datasets in order to show the effectiveness of the presented study and of the proposed methodology.


Introduction
Relaxation times define the rate of spin magnetic equilibrium recovery in nuclear magnetic resonance (NMR) [1,2]. For each tissue, several relaxation times can be defined. Besides, the main interest is in and T E combinations, it is possible to emphasize the effect of one tissue intrinsic parameter with respect to others, obtaining the well-known ρ-weighted, T 1 -weighted or T 2 -weighted images.
Given the previously mentioned motivations, we present a deeper analysis of the complex SE model considered in [18] extended to three parameters (i.e., ρ, T 1 and T 2 ). The analysis is conducted exploiting the Cramer-Rao lower bounds (CRLBs) [16]. Since CRLBs provide the best achievable performances in the unbiased estimation of one or more parameters, by minimizing them with respect to the MR scanner imaging parameters, it is possible to find the optimal acquisition configuration for the relaxation time estimation. Practically speaking, we look for the acquisition parameters that allow achieving lower relaxation time estimation errors. The result of the study is the introduction of a general empirical rule for determining the optimal (with respect to CRLBs) MRI scanner parameter configuration. In particular, the identification of these parameters in the case of several tissues has been conducted. The effectiveness of the theoretical results and of the empirical rule is validated and verified on different datasets.
The manuscript is organized as follows: in Section 2, the acquisition model for an MRI spin echo sequence is presented, and in Section 3, the achievable performances of the estimation are analyzed via the CRLBs. In Section 4, the CRLB-based empirical rule for the optimal acquisition parameter configuration is presented. Finally, validation on different datasets is presented in Section 5, and conclusions are drawn.

The Model
Let us consider an MRI acquisition system using a spin echo imaging sequence. The amplitude of the recorded complex signal after the image formation process, i.e., after the computation of the 2D Fourier transform, is related to the tissue parameters, ρ, T 1 and T 2 , via [2,21]: where T E and T R are the echo and repetition time, respectively, which are two imaging parameters that can be set in the MRI scanner, and θ = [ρ T 1 T 2 ] T is the vector containing the tissue parameters in which we are interested. Note that Equation (1) is valid in the case of a homogeneous imaged object. In the case of clinical data, the presence of different hydrogen environments within each voxel has to be taken into account. The acquisition model reported in Equation (1), which is a solution to Bloch equations, assuming that T E is short with respect to T R , is related to the noise-free case and does not take into account the dependency on the static magnetic field, B. Considering noise, in the complex domain, the model becomes: where n R and n I are the real and imaginary parts of the noise samples, which are distributed as independent circularly Gaussian variables [22], and φ represents the angle of the complex data [23,24]. Thus, the statistical distributions of the real and imaginary parts of the acquired signal are: where σ 2 is the variance of real and imaginary noise components. Due to the independence of the real and imaginary parts of noise, the joint statistical distribution of y R and y I is the product of the two probability density functions of Equation (3). Once N acquisitions with different T R and T E combinations have been recorded and collected in the data vector y = [y R , y I ], with y R = [y R (1), · · · y R (N )] and y I = [y I (1), · · · y I (N )], we can derive the likelihood function from the factorization of the Probability Density Functions (PDFs): Starting from the likelihood function of Equation (4), the CRLBs for θ are derived and analyzed in the following sections.

Cramer-Rao Lower Bounds Evaluations
In order to evaluate the performances of the optimal estimator for the model presented in Section 2, the Cramer-Rao lower bounds have to be computed. According to Statistical Estimation Theory [16], given an observation model, the accuracy of any estimator can be evaluated according to its mean and its variance. In particular, in order to be optimal, an estimator needs to have its mean equal to the value to be estimated (i.e., unbiased estimator) and to have the smallest possible variance. CRLBs represent the lower bound of the variance of any unbiased estimator, resulting an interesting and powerful tool for evaluating the achievable performances of a considered model. By computing the CRLBs for different configuration of the parameters involved in the acquisition model, it is possible to find the best parameter configuration, the one that provides the minimum values of CRLBs (i.e., the minimum achievable variances). Considering the vector parameter θ, the minimum variance that any unbiased estimator of parameter θ i can reach is provided by the i-th diagonal element of the inverse of matrix I [16]: with I being the Fisher matrix, which is equal to: where E [·] is the expected value operator. A closed form for the second order derivatives of Equation (6) has been derived and reported in the Appendix. The closed form greatly improves the computational accuracy of the CRLB evaluation and decreases the computational burden of the simulations reported in the following.
In order to experimentally compute the matrix of Equation (6), Monte Carlo simulations with 10 5 iterations have been considered for statistical average computation.
For the following evaluations, we considered a tissue, named A, with parameters θ = [ρ T 1 T 2 ] T = [2.5 1600 90] T . Note that within the paper, all relaxation times are expressed in milliseconds, while the proton density is in percentage. The following simulations are reported and analyzed in order to investigate CRLB dependency and behavior with respect to the signal-to-noise ratio (SNR), the number of acquisitions and the scanner acquisition parameters.

CRLB vs. SNR
Let us start by computing CRLBs varying the noise standard deviation (i.e., the SNR). Sixteen images have been considered, which refer to the all combinations of four T R and four T E values equally spaced in the intervals [500, 3500] ms and [80, 350] ms, respectively. Note that the lower T E value has been set according to the minimum echo time for the SE sequence accepted by the Philips Achieva 3.0 T, the MR scanner we worked on, while the maximum value of T R has been set in order to limit the global acquisition time. The CRLBs in the case of different SNRs are shown in Figure 1. As expected, the square root of the CRLBs related to all considered parameters decreases with respect to SNR growth, i.e., high SNRs positively affect the estimator performances. In the considered range of SNRs, no saturation appears. Very similar behaviors are obtained varying T R and T E combinations. Globally, it can be stated that SNR linearly affects CRLBs, so in the following the results of each simulation can be easily extended to any SNR configuration.

CRLB vs. the Number of Acquisitions
A second case study has been conducted in order to evaluate the advantage of increasing the number of acquisitions. Two vectors of T R and T E , of a length of N R and N E , respectively, have been generated in the [500, 3500] ms (for T R ) and [80, 350] ms (for T E ) intervals. The square root of CRLBs, i.e., the minimum achievable standard deviations, are reported in Figure 2 for ρ, T 1 and T 2 , respectively, for different N R and N E combinations. The noise variance has been fixed in order to obtain an SNR of 16 dB for the image with the lowest signal intensity (i.e., T R = 500 ms and T E = 80 ms). It can be noted that the number of T R values mainly affects the achievable performances with respect to T 1 estimation, while CRLBs of ρ and T 2 are dependent on the number of both T E and T R values, with a higher dependency on echo times. The results confirm the strict connections between T R and T 1 and also between T E and T 2 , as expected from the exponential terms of the SE signal model (Equation (1)). However, it is interesting to stress how the CRLB of ρ is very tightly related to T E values rather than to T R ones.

CRLB vs. T R and T E Values
As a further case study, an evaluation of CRLBs with respect to T R and T E values with a fixed number of acquisition has been performed. Four acquisitions have been considered, corresponding to all the combinations of CRLBs have been computed while varying T R (1) and T E (2) and considering T R (2) = 3, 500 ms and T E (1) = 80 ms, again in the case of tissue A parameters. Results are reported in Figures 3. Figure 4a shows that the ρ estimation would prefer low T R (1) and high T E (2) values. The behaviors of CRLBs for T 1 and T 2 differ remarkably from Figure 4a, as it can be noticed that the estimation of T 1 is almost unresponsive with respect to T E (2) values, as far as T 2 estimation with respect to T R (1). In particular, for the estimation of T 1 , the ideal T R (1) is as low as possible, while the ideal T E (2) for the estimation of T 2 is between 150 and 250 ms. For this experiment, a second dataset has also been considered: the same simulation has been conducted in the case of a second tissue, named B, with parameters θ = [ρ T 1 T 2 ] T = [2.8 1800 60] T , in order to know if the results of Figure 3 are always valid or if they are highly dependent on the considered tissue. The results are reported in Figure 4. It can be noticed that the lower regions remain in the same position, although being increased in value, but for CRLBs of ρ and T 2 , the ideal T E (2) range reduces to [130,180] ms. This is mainly due to the lower T 2 value of tissue B with respect to tissue A. Thus, it can be concluded that the general trend is confirmed, although the position of the global minimum is strictly related to the considered tissue. These two simulations show that the choice of optimal parameters is strictly dependent on the relaxation times of the imaged tissues. In the next section, we investigate the possibility of finding a rule for ideal imaging parameter identification.

Optimal Parameter Configuration
After the evaluation of ρ, T 1 and T 2 CRLB behaviors, an analysis dedicated to the computation of optimal T R and T E combinations is presented. In the following, it will be shown that a proper imaging configuration can greatly improve the performances with respect to such a choice. In particular, the aim of this section is to identify the ideal imaging parameters with respect to imaged tissues.
Let us show how the optimal imaging parameters can be determined. Initially, we have focused on the minimization of T 2 CRLB, which consist in finding T E values that minimize the element (3, 3) of the inverse of the Fisher matrix, I(θ), of Equation (6) for different spin-spin relaxation times T 2 . The optimization has been performed by searching the three optimal T E values in the [82, 350] ms range for a fixed value of T R . The evaluation has been done varying the tissue T 2 relaxation times in the [20,200] ms range, obtaining the results shown in Figure 5. Some considerations can be drawn: • there is no T E value combination that is simultaneously ideal for tissues with different spin-spin relaxation times. As a consequence, we can only find the T E combination that is ideal for a specific tissue; • by analyzing Figure 5, it can be noticed that the lowest T E value of the ideal configuration is always equal to the lower bound of the considered variability range, which, in our case, was fixed to 82 ms. As stated before, this value is the minimum echo time for the SE sequence accepted by the Philips Achieva 3.0 T, the MR scanner we worked on; • the two higher T E values, which are the red and the green lines of Figure 5, practically coincide.
This can be explained considering that we are interested in the estimation of relaxation times, i.e., of decay rates. In order to achieve such a goal, it is crucial that the measurement of the signal decrease, i.e., the ratio of the signal acquired in two different echo times. Therefore, instead of values T E , it is only important the difference between them. A third echo time, T E (3), equal to T E (2), allows us to compute twice the signal decay, which is the quantity in which we are interested; • the red and the green lines of Figure 5 show a clear trend: their values grow linearly when increasing T 2 . In particular, we found that the distance with the blue line (i.e., lowest T E , 82 ms) is a little bit bigger than the value of the considered spin-spin relaxation time, T 2 . For example, in the case of T 2 = 100 ms, the ideal echo times were T E = [83, 197, 205] ms; the last two values are approximately 110% of (T E (1) + T 2 ). By considering the other simulated T 2 values, we found that this coefficient is 110% ± 10%. Within this range, the CRLB of T 2 can be considered constant. From these simulations, we can derive an empirical rule for the optimal T E selection: the lower one should be fixed to the minimum value accepted by the MR scanner, while the other values should to be set in the range of 100%-120% of the value of (T E (1) + T 2 ), considering the T 2 of the tissue in which we are interested.
A similar evaluation has been conducted for the minimization of T 1 CRLB varying MRI scanner repetition times T R , with a fixed value of T E ; the results are shown in Figure 6. The higher T R value is fixed to the right edge of the considered variability range, which we set equal to 4, 000 ms. The intermediate and low T R s have similar values, which, starting from 500 ms in the case of tissue with T 1 = 700 ms, grow almost linearly up to 1, 400 ms for tissues with higher T 1 (about 3, 000 ms). It is hard to determine an empirical rule in this case; anyway, we can say that a choice of around 1, 000 ms for T R (1) and T R (2) will fit a wide class of tissues, i.e., those with 1, 200 < T 1 < 2, 000 ms. Figure 6. T R values that minimize the CRLB of T 1 for tissues with different spin-lattice relaxation times, T 1 . Three values have been considered: the red line is for the highest T R value, the blue line for the lowest one and green for the intermediate one.
Concluding this section, one more evaluation has been conducted. Instead of optimizing T R and T E values separately, a joint minimization has been done. Nine acquisitions have been considered, related to three repetition and three echo times. Among the three values, the lower and the higher have been fixed to the search range bounds, so only the intermediate T E and T R values were variable. Results are shown in Figure 7, respectively. It is evident from the figure that T E values can be considered independent from T 1 , as far as T R from T 2 , proving the correctness of the separate optimization of the echo and repetition times. In particular, from Figure 7a, we can state that tissues with equal T 2 , but very different T 1 values share the same three optimal echo times for T 2 estimation, and vice versa. That said, the behaviors of Figures 5 and 6, i.e., the minimization, one parameter at a time, are confirmed.
In order to easily apply the obtained results, the ideal acquisition parameters for different tissues have been computed exploiting CRLB minimization in the case of a 1.5 T and a 3 T MRI scanner. The results are shown in Tables 1 and 2 for T 1 and T 2 , respectively. According to the results reported in Figure 7, the minimizations have been computed independently for spin-lattice and spin-spin relaxation time estimation. Tissue relaxation times have been simulated according to reference values present in the literature [25], which are reported in Table 3.
In Table 4, optimal echo times in the case of gray matter T 2 estimation for different minimum T E are reported. It can be noticed that the lower optimal echo time is always the minimum and that the empirical rule is confirmed.  Note that the usefulness of a proper T R and T E selection, besides the lower estimation variance, consists also in reducing the acquisition time. In order to make such an advantage evident, Table 5 reports the achievable performance in the case of 16 images (4 T R and 4 T E values) when moving from equally spaced to optimized acquisition parameters. In particular, the last column of Table 5 shows that 12 acquisitions, with properly chosen parameters, can lead to better results with respect to 16 equally spaced images, while definitely reducing the global acquisition time.

Numerical Experiments
Within this section, some numerical results are shown in order to validate the advantage of the optimal selection of the imaging parameters according to the previously reported theoretical studies. A tissue with parameters [ρ T 1 T 2 ] = [5.5 775 44.5] has been considered. Three noisy datasets (SNR = 30 dB) have been simulated, each one composed of four acquisitions. The parameters of Dataset 1 have been chosen according to the results of Figures 5 and 6 in order to optimize the estimation for the considered tissue. Datasets 2 and 3 have been generated with non-ideal parameters. The dataset characteristics are summarized in Table 6. To asses and validate the CRLB studies, the estimation of the relaxation times has been implemented via Monte Carlo simulation. In particular, a maximum likelihood estimator (MLE) has been set up in the complex domain. Considering that the noise is circularly Gaussian distributed, MLE corresponds to a non-linear least squares (NLLS) estimator [18]. It is important to note that the previously reported theoretical studies about the optimal selection of the imaging parameters are valid for any unbiased estimator, since CRLBs are related only to the acquisition model. Among different estimators, NLLS has been chosen thanks to its low computational times and complexity. We recall that the choice of the optimal estimator is not the aim of this paper. The NLLS estimator for the ρ, T 1 and T 2 parameters has been implemented on the three datasets of Table 6. A quantitative analysis of the results, in terms of estimation means and variances, has been reported in Table 7. By analyzing it, it is possible to infer that the estimator means are very close, while variances significantly vary from one dataset to the other. In particular, the smallest variances are obtained in the case of Dataset 1. This fully agrees with the theoretical studies reported in Section 4; as a matter of fact, Dataset 1 has been generated by using the previously developed optimal T E and T R parameter selection for the considered relaxation times. It is evident that choosing a non-ideal imaging parameters configuration can lead to very inaccurate results. For example, the T 2 estimator variance of Dataset 3 is approximately six times larger than the one of Dataset 1. In order to visualize such results, the normalized standard deviations of ρ, T 1 and T 2 in the case of Datasets 1, 2 and 3 are reported in Figure 8. Figure 8. Square root of the CRLB for proton density (blue), spin-lattice (T 1 ) relaxation time (green) and spin-spin (T 2 ) relaxation time (red) for the dataset with different acquisition parameters. CRLBs values have been normalized for the parameter values in order to be plotted in the same graph.
The higher achievable accuracy in the case of optimal imaging parameters selection can also be inferred from the empirical probability density functions of the estimators, reported in Figure 9. In each image, the blue, the green and the red curves refer to Datasets 1, 2 and 3 of Table 6. As expected, most of the presented estimators follow a Gaussian distribution, with a different width. Blue curves obtained using Dataset 1, characterized by the optimal T R /T E values for the simulated tissue, are always the narrowest (smallest variances). Moving to curves obtained from Datasets 2 and 3, the estimation error becomes larger. Moreover, in the case of the ρ estimator, the results start showing a bias in the case of Dataset 3, i.e., the one with the worst acquisition parameters, and the empirical PDF does not look like a Gaussian function any more.  Table 8. Taking into account the developed procedure, Dataset 1 parameters represent the ideal configuration for the first tissue, while Dataset 2 is the ideal for the second one.
The empirical probability density functions for T 1 and T 2 estimators have been computed for both datasets and are shown in Figures 10, respectively. Once again, The results validate the theoretical study of Section 4. Estimation based on Dataset 1 (red line) shows lower variance in the case of spinal cord, i.e., the tissue with the lowest relaxation times (the left peaks of Figures 10). Considering gray matter, Dataset 2 (blue line) -based estimation gives better results, although the improvement of the T 1 estimator is not pronounced. Once again, the result highlights the need of properly tuning the acquisition parameters. Table 8. Acquisition parameters in the case of two tissues; the datasets are composed of four images.

Conclusions
Within this paper, an analysis on the spin echo signal model in MR imaging has been addressed. In particular, Cramer-Rao lower bounds for relaxation time estimation in the case of a complex Gaussian acquisition model have been evaluated. Several CRLB-based evaluations have been presented in order to investigate the possibility of finding the optimal, in terms of reconstruction accuracy, imaging parameter configuration for the estimation of T 1 and T 2 maps. According to these theoretical studies, an empirical rule together with the identification of the optimal imaging parameter combination (echo and repetition times) in case of different tissues (different T 1 and T 2 ) has been proposed. Moreover, the optimal acquisition parameters for several tissues have been computed for both 1.5 T and 3 T acquisitions. The theoretical results have been numerically validated on different datasets. It should be underlined that such optimal parameters are independent from the implemented estimators, as CRLBs only depend on the signal model. Once the data have been acquired, different estimators proposed in the literature can be applied. It is important to underline that the theoretical studies reported within the paper can be easily adapted to different imaging sequences.
In order to compute the expected value of Equation (7), Monte Carlo simulations have to be implemented.

Author Contributions
The author contributions are substantially equal. In particular, the main contribution of Vito Pascazio, Giampaolo Ferraioli and Fabio Baselice was the methodology development. Moreover, Giampaolo Ferraioli and Fabio Baselice were specifically focused on the numerical implementation. Alessandro Grassia worked both on code implementation and on software simulation tasks.