Generalized min-max bound-based MRI pulse sequence design framework for wide-range T1 relaxometry: A case study on the tissue specific imaging sequence

This paper proposes a new design strategy for optimizing MRI pulse sequences for T1 relaxometry. The design strategy optimizes the pulse sequence parameters to minimize the maximum variance of unbiased T1 estimates over a range of T1 values using the Cramér-Rao bound. In contrast to prior sequences optimized for a single nominal T1 value, the optimized sequence using our bound-based strategy achieves improved precision and accuracy for a broad range of T1 estimates within a clinically feasible scan time. The optimization combines the downhill simplex method with a simulated annealing process. To show the effectiveness of the proposed strategy, we optimize the tissue specific imaging (TSI) sequence. Preliminary Monte Carlo simulations demonstrate that the optimized TSI sequence yields improved precision and accuracy over the popular driven-equilibrium single-pulse observation of T1 (DESPOT1) approach for normal brain tissues (estimated T1 700–2000 ms at 3.0T). The relative mean estimation error (MSE) for T1 estimation is less than 1.7% using the optimized TSI sequence, as opposed to less than 7.0% using DESPOT1 for normal brain tissues. The optimized TSI sequence achieves good stability by keeping the MSE under 7.0% over larger T1 values corresponding to different lesion tissues and the cerebrospinal fluid (up to 5000 ms). The T1 estimation accuracy using the new pulse sequence also shows improvement, which is more pronounced in low SNR scenarios.


Introduction
Quantitative estimation of longitudinal relaxation time, termed T 1 relaxometry, offers a very useful approach in magnetic resonance imaging (MRI) for enhancing tissue contrast, improving tissue characterization, and evaluating neuro-degenerative pathologies, etc. [1][2][3]. A T 1 map shows more accurately the brain tissue characteristics when compared against more commonly used contrast-based qualitative approaches, such as T 1 /T 2 -weighted a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 images. T 1 relaxometry enables better tissue classification and holds promise for improving detection of early-stage tissue degeneration, as well as characterization of advanced tissue destruction [3]. Therefore, a fast, accurate, and precise T 1 relaxometry technique can potentially be applicable in a variety of neuro-degenerative disorders, including multiple sclerosis (MS) [4], Alzheimer's disease [5] and Parkinson's disease [6], as well as in assisting imageguided surgeries.
T 1 relaxometry techniques estimate T 1 on a voxel-to-voxel basis using the magnetic signals acquired with specific MRI pulse sequences, which differ in their pulse times and flip angles. The pulse sequence parameters directly affect the acquired MR signals and thereby the performance of T 1 relaxometry techniques. Estimating T 1 by sampling the T 1 relaxation signals has a five decade history [7][8][9][10]. The conventional methods are based on inversion recovery (IR) [7] and saturation recovery (SR) [8] due to their relatively large signal dynamic ranges. However, the clinical applications of IR and SR sequences are severely hampered by the considerable time required for the partial recovery of the longitudinal magnetization. To achieve a reasonable scan time, Look and Locker proposed a "one-shot" method, which samples multiple points along the T 1 relaxation curve by continuously tipping the longitudinal magnetization with small flip angle pulses [9]. More recently, Deoni et al. presented a fast and high-resolution T 1 mapping approach, coined driven-equilibrium single-pulse observation of T 1 (DESPOT1) [10]. It uses a pair of spoiled gradient recalled echo (SPGR) images with different flip angles and achieves a shorter total scan time than other methods. However, researchers found that the T 1 estimates using DESPOT1 are generally biased because the data processing does not properly account for noise [11]. This bias can be significant, such that T 1 values are overestimated by as much as 10% to 20% for clinical SPGR images with T 1 of 800-1600 ms [11]. Moreover, DESPOT1 is limited in its ability to estimate large T 1 values corresponding to relatively advanced lesions and in its sensitivity to pulse flip angle perturbations [10].
Neuro-degenerative diseases like MS exhibit a wide range of damage to the brain's white matter. Consequently, T 1 relaxometry sequences must meet three requirements in order to be clinically useful. First, the technique needs to provide precise and accurate measurements in the range of T 1 s corresponding to normal white matter (WM, 800-900 ms at 3T) so as to identify early changes from the 'normal' condition. Second, the technique needs to be stable for measuring values up to those corresponding to cerebrospinal fluid (CSF, 4500-5500 ms at 3T) in order to be able to characterize the degree of damage in more advanced lesions. Third, the technique needs to do so within clinically acceptable scan times (of the order of 10-15 minutes) while being reasonably stable to variations of the radiofrequency (B 1 ) field. Currently, no existing technique fulfills all three requirements. These technological limitations preclude the use of T 1 relaxometry as a tool that may address both early white matter degeneration as well as lesion evolution, which are both crucial markers for monitoring the performance of neuroprotective drugs in treating dementia and progressive MS. This paper proposes a new pulse sequence design framework for optimizing the T 1 relaxometry performance for a broad range of T 1 values. There are two major components of the proposed framework that differ from previous ones. The first is to use the Cramér-Rao bound (CRB) to design the pulse sequence. The CRB provides a lower bound on the variance of any unbiased T 1 estimate. This bound-based design allows us to predict the performance of an MRI sequence based on its sequence parameters, and then to optimize the performance of that sequence by adjusting these parameters. Our CRB derivation leads to a geometric interpretation, which brings into sharper focus all of the factors that control the precision of T 1 relaxometry sequences. This differs from previous sequence design approaches that merely focus on increasing the dynamic range and signal-to-noise ratio (SNR) to improve T 1 mapping performance [10,12,13]. The second novel component of our approach is to design the pulse sequence by optimizing its performance over a broad range of T 1 values. Previous sequence design algorithms optimized the signal dynamic range or sensitivity for a single nominal T 1 value, or a small range of T 1 values [10][11][12][13]. When the tissue relaxation times fell outside the narrow range of nominal values, the relaxometry performance degraded rapidly. In contrast, our framework employs a min-max optimization strategy, which designs the pulse sequence to minimize the maximum estimate variance over a broad range of T 1 values. This guarantees the overall optimality of the resulting sequence for all T 1 values within the target region spanning thousands of milliseconds.
To demonstrate the effectiveness of the proposed sequence design framework, we optimize the tissue specific imaging (TSI) sequence for T 1 relaxometry. The proposed technique is sufficiently general that it can optimize nearly any T 1 relaxometry sequence. Optimizing the TSI sequence provides a powerful example of the benefits of bound-based min-max optimization because it yields three high-contrast images of WM, gray matter (GM) and CSF, which are employable as anatomical references in clinical settings. The TSI sequence is a relatively new imaging sequence and has been successfully applied for characterization of MS lesions [14,15]. The TSI pulse sequence includes three imaging pulses followed EPI acquisitions, which are interleaved with two inversion pulses in each pulse repetition period. This sequence was originally designed to acquire brain images for each of the two categories of brain tissues (WM and GM), as well as the CSF, with optimal contrast. We address the potential applications of TSI for T 1 relaxometry and optimize its sequence parameters to improve the T 1 estimate precision and accuracy while maintaining its total scan time. Another motivation for applying the TSI-type sequences for T 1 relaxometry is their improved precision relative to DESPOT1 over the T 1 range corresponding to normal brain tissues, and, more importantly, their improved stability for larger T 1 values corresponding to tissues such as advanced MS lesions [16,17].

MR signal model
T 1 relaxometry approaches apply the RF pulse sequences repeatedly to initialize the magnetization preparation. Denote the MR signal generated by a tissue voxel at time t i as x i = M 0 h i (T 1 ), i = 1, 2, . . ., n, where M 0 is the equilibrium longitudinal magnetization and h i (T 1 ) is the signal weighting factor at time t i . The acquired signal at time t i follows with additive uncorrelated noise w i for each acquisition. The vector s = [s 1 , . . ., s n ] T characterizes the acquired signals at different times t 1 , . . ., t n , where (Á) T denotes vector transpose. Eq (1) is known as a universal MR signal acquisition model assuming additive noise [25]. Assuming knowledge of the pulse sequence parameters (sequence repetition time TR, pulse times t i and flip angles α i ), the signal weighting factor as a function of T 1 after the nth pulse can be derived from the Bloch equations in Eq (2), where M eq z is the steady state magnetization. To derive M eq z , we assume that in each TR, the first imaging pulse is applied at the initial time t = 0 and there is a delay of (TR − t n ) after the nth pulse of one sequence and before the first pulse of the next sequence. h n ðT 1 Þ ¼ sin a n 1 þ The MR signals are most commonly acquired through quadrature detector channels, each of which typically suffer from independent additive zero mean, white and Gaussian noise. Thus, the noise in the reconstructed magnitude MR images should follow a Rician distribution [18]. Several references in the MRI literature use Rician distributed random noise in their signal models [19,20]. When the signal-to-noise ratio (SNR) is high enough, the Rician distribution converges to the Gaussian distribution [21,22]. For the convenience of the CRB evaluation, we assume the SNR is sufficiently high to exploit the approximation of the noise as additive zero mean, white and Gaussian noise in the signal model in Eq (1). As will be shown later, this Gaussian model assumption leads to a readily-derived closed-form expression to geometrically interpret the CRB, which provides insight into the factors controlling the T 1 estimate precision.

Cramér-Rao bound on joint M 0 and T 1 estimation
The CRB has been used as a quantitative tool for optimizing experimental MR protocols [20,23,24] and also evaluating the precision of specific relaxometry sequences [25][26][27]. The analytic expression for the CRB on the T 1 estimate, when jointly estimating M 0 and T 1 , is the foundation for the sequence design framework. Letting the parameter vector θ = [M 0 , T 1 ] T , the covariance matrix CðŷÞ of the unbiased estimatorŷ satisfies where I(θ) is the 2 × 2 Fisher information matrix (FIM) [28]. As a consequence, the variance of any unbiased estimation of T 1 is bounded from below by the (2,2) entry of matrix I −1 (θ). This quantity is also known as the CRB of the T 1 estimate The diagonal elements of the FIM represent the measured signals' sensitivity to the parameters in θ [28]. For joint estimation of M 0 and T 1 , the derivation of the FIM is straightforward following Eqs (6) and (7) IðyÞ ¼ E ½r y ln pðs; yÞ½r y ln pðs; yÞ where E(Á) takes the expectation over the acquired signal s. Eq (6) has entries ðh i ðT 1 ÞÞ 2 ; where σ is the standard deviation of the additive noise w. Note that the dependence of the FIM on the pulse sequence parameters TR, t 1 , . . ., t N and α 1 , . . ., α N is suppressed here for brevity.
Defining SNR = M 0 /σ yields the closed-form expression for the CRB on T 1 estimate To maximize the precision of the T 1 estimate, we should choose the sequence parameters to minimize Eq (8). The apparent complexity of Eq (8) can be simplified through a linear space interpretation [29]. To achieve this, define the signal weighting vector as h(T 1 ) = [h 1 (T 1 ), h 2 (T 1 ), . . ., h n (T 1 )] T , and the sensitivity vector @h/@T 1 as the derivative of h with respect to T 1 . Moreover, define ϕ as the principal angle between the vectors h and @h/@T 1 in the linear space (see Fig 1). Note that the principal angle ϕ between two vectors x and y in the linear space is defined as 0 ¼ arccos hx;yi kxkÁkyk , where hÁi is the inner product operator and kÁk is the Euclidean norm (or length) of a vector. Therefore, Eq (8) can be rewritten as h is the signal weighting vector containing the measured signals at all acquisition times. @h/@T 1 is the sensitivity vector, calculated as the derivative of the signal weighting vector with respect to T 1 . Conceptually, increasing the norm of the sensitivity term k@h/@T 1 k will increase the impact of small changes in T 1 on the acquired signals. The orthogonality term sin ϕ is a consequence of the joint estimation of T 1 and M 0 . The observed signal's sensitivity for M 0 is h, while that of T 1 is @h/@T 1 . The more orthogonal these vectors are, the easier it becomes to ascribe changes in the observed signal to M 0 or T 1 unambiguously. doi:10.1371/journal.pone.0172573.g001 Eq (9) makes clear that three factors control the CRB and thus there are three methods to improve the T 1 estimate precision. The first method is to increase the SNR, which is consistent with previous sequence designs seeking to increase the signal dynamic range or reduce the noise [10,12]. These approaches would improve T 1 estimate stability, but are not always the most effective methods of reducing variance. The second term @h/@T 1 in Eq (9) is the sensitivity of T 1 , which describes how sensitive the signal model is to T 1 variation. Increasing the norm of the sensitivity k@h/@T 1 k will increase the impact of small changes in T 1 on the overall signal weighting vector h. For example, a large sensitivity indicates that a small T 1 variation causes a large signal fluctuation. In this case, T 1 can be estimated from the signal measurements more precisely. In contrast, zero sensitivity indicates that the signal will remain constant regardless of the T 1 variation. In this case, the T 1 value can never be estimated from the measurements. Increasing the magnitude of the sensitivity vector can be an effective method of improving T 1 estimate precision. The third component in the CRB expression is the orthogonality term sin ϕ, which results from jointly estimating T 1 and M 0 . The observed signals' sensitivity for M 0 is h, while that of T 1 is @h/@T 1 . The more orthogonal these vectors are, the easier it becomes to ascribe changes in the observed signal to M 0 or T 1 unambiguously. In terms of T 1 estimation, increasing the orthogonality of the signal weighting vector and the sensitivity vector improves the T 1 estimate precision. The second and third approaches are notable as they can be exploited without requiring the costly hardware improvements usually needed to increase signal strength or decrease measurement noise. To the best of our knowledge, no prior MRI sequence design strategy has exploited these mechanisms simultaneously and analyzed these approaches explicitly for improving the precision of T 1 relaxometry.

Sequence optimization methods
As Eq (8) shows, the CRB is a function of the true T 1 value. Current publications suggest that the variation of T 1 within normal and diseased tissues scales with the mean of T 1 . For example, clinical T 1 measurements at 1.5T using histograms of normal-appearing white matter (NAWM) showed T 1 of 792 ms ± 36 for patients with secondary progressive MS. This represents a relative T 1 variation of ± 4.5% around its mean. However, histograms for the cortical normal-appearing gray matter (NAGM) showed T 1 of 1355 ms ± 62 for patients with secondary progressive MS, which also represents a relative T 1 variation of ± 4.5% around its mean [30,31]. In this example, the same stage of MS development corresponds to the same relative T 1 error for different tissue types. Therefore, rather than directly using the CRB, this paper uses a relative error as the metric for optimizing the pulse sequence parameters. The relative error is defined as the square root of the CRB normalized by the true T 1 value. The sequence parameters are optimized to minimize the maximum relative error over a broad range of T 1 values. This guarantees the relative error to be reasonably robust to the T 1 range of interest.
Fig 2 illustrates the structure of the TSI pulse sequence [14]. In each TR, there are three imaging pulses each followed by EPI acquisitions and interleaved by two inversion pulses. The sequence parameters were derived from a simulated annealing optimization process by maximizing the contrast-to-noise ratio (CNR) in the combined tissue-specific images. Using 3D sensitivity encoded (SENSE) EPI for data acquisition, the sequence achieves a 1.15 mm isotropic resolution in a FOV of 220 × 165 × 110 mm 3 within a scan time of 10 minutes. There are eight pulse parameters to optimize in this sequence: times for the two inversion pulses t 2 , t 4 , times for the second and third imaging pulse t 3 , t 5 , flip angles for the three imaging pulses α 1 , α 3 , α 5 and the sequence repetition time TR. The optimization process assumes the T 1 interval of interest to be 700-2000 ms for normal brain tissues at 3.0T [32]. However, the performance of the optimized pulse sequence is evaluated over a broader T 1 range of 700-5000 ms. This is to observe the sequence's robustness to larger T 1 variations that characterize more advanced lesions and the CSF region [3]. To achieve optimal sequence parameters that are practically applicable, we constrain the optimization process following the discussions in [14]. Specifically, the maximum allowed TR is set as 6 seconds to limit the total scan time to under 10 min. The inter-pulse interval needs at least 100 ms to allow enough time for the EPI acquisition. The flip angles of all imaging pulses are less than 90˚for the convenience of practical implementation. More explicitly, the optimal pulse parameters are achieved by minimizing C max (t, α) under the following constraints ðt; αÞ opt ¼ arg min t;α C max ðt; αÞ; TR À t 5 ! 100 ms; Due to the large number of parameters to optimize, classical optimization methods such as exhaustive search can be too computationally burdensome and gradient search can fail to In each TR period, there are three imaging pulses (dark gray) followed by EPI acquisitions and interleaved by two inversion pulses (light gray) (After Fig 1 from [14]). The three imaging pulses are characterized by their flip angles α 1 , α 3 , α 5 . The dashed lines indicate the times t when each pulse is applied. The first imaging pulse is applied at the beginning of each TR (t 1 = 0). There are 8 pulse parameters to optimize for in the TSI sequence: times for the two inversion pulses t 2 , t 4 , times for the second and third imaging pulses t 3 , t 5 , flip angles of the three imaging pulses α 1 , α 3 , α 5  converge to the globally optimal solution. To find the globally optimal sequence parameters, we adopt a hybrid of the Nelder-Mead downhill simplex method [33] and the simulated annealing approach [34]. The downhill simplex method finds the local minimum by expanding, contracting, and reflecting the simplex constructed by a group of different initial points in the N-dimensional hyperspace. With simulated annealing, the algorithm conditionally accepts uphill movement (leading to worse performance) during the optimization process and therefore improves the probability of finding the global optimum. This hybrid approach [35] has been proven successful in optimizing brain tissue contrast [14] and acquisition schemes for quantitative magnetization transfer MRI [36].
The optimization process is initialized with a TR of 5500 ms, evenly spread-out pulse times within TR, and a flip angle of 45˚for each imaging pulse. For each temperature T, the simplex routine iterates 2000 times to calculate the cost function in Eq (10). For each iteration, the relative error is calculated over the T 1 range 700-2000 ms, where the worst cases of the relative errors are compared and arranged in an ascending order for all simplex vertices. The relative error is randomly perturbed by a quantity T Á log(p), where p is a random variable uniformly distributed within [0, 1], to allow for conditional acceptance of uphill movement [34]. The temperature T and the simplex size D decrease according to the annealing schedule Tðn þ 1Þ ¼ 0:985TðnÞ The initial temperature is selected empirically as 0.8 and the initial simplex is constructed with a simplex size of 100 ms for the pulse times and 30˚for the flip angles. The overall simulated annealing algorithm runs for 1000 iterations to locate a hopefully global optimum. During the optimization process, a prohibitive penalty of 10 3 is added to the cost function whenever the new pulse sequence violates the constraints in Eq (11). Repeating the optimization process from several different initial sequence vectors allowed the pulse sequence to converge to a relatively stable and hopefully globally optimal CRB.

T 1 estimation methods
In general, there are two criteria to evaluate the performance of an estimator given an acquired signal s: accuracy and precision [28]. The accuracy of a T 1 estimatorT 1 ¼ f ðsÞ is measured by the bias which is the difference between the mean of the T 1 estimator EðT 1 Þ and the true T 1 . The smaller the bias is, the more accurate the estimatorT 1 is. The precision is determined by the variance ofT 1 The smaller the variance is, the more precise the estimatorT 1 is. An ideal T 1 estimatorT 1 is unbiased, with its variance achieving the CRB of T 1 . For the signal model in Eq (1), the least squared estimator (LSE) of T 1 for joint M 0 and T 1 estimation is asymptotically optimal. For infinitely high SNR, this estimator has its expected value achieve the true T 1 and its variance achieve the CRB(T 1 ) [28]. Mathematically, the LSE of θ = [M 0 , T 1 ] minimizes the squared error where x(θ) describes the noise-free signal. Specifically, when the signal model shows linearity in one parameter M 0 and non-linearity in T 1 , the LSE of T 1 can be calculated in a more computationally efficient way by minimizing another version of squared error [28] following where I is the identity matrix and h is the signal weighting vector defined in Eq (2). Since the squared error in Eq (16) depends only on one parameter T 1 , it is easy to find the T 1 value which minimizes the squared error using a grid search of T 1 in an appropriate range.

Numerical simulation methods
The accuracy and precision of T 1 estimates are evaluated in Monte Carlo simulations. The simulations consider four different T 1 estimation approaches: nonlinear least square estimation (NLSE) using the optimized TSI sequence, NLSE using the original TSI sequence [14], linear LSE using the SPGR sequence (also coined DESPOT1 [10]), and NLSE using the SPGR sequence [11]. The observed magnetic signals were simulated for an equilibrium longitudinal magnetization value of M 0 = 3000 a.u and a range of T 1 within 700-5000 ms covering almost all brain tissues and CSF at 3.0T [32]. The different structures of the TSI and SPGR sequences require different approaches to simulate the observed noise-free signals. For the TSI signals, the steady state magnetizations in Eq (3) are first evaluated based on different choices of TSI sequence parameters. The signal weighting factor in Eq (2) is then evaluated at different pulse times and scaled by M 0 as the raw noise-free magnetic signals. In contrast, the raw SPGR signals are simulated by evaluating with E 1 ¼ e À TR=T 1 at a given constant TR and varying flip angles [10]. For both TSI and SPGR, the simulated magnetic signals are then generated by adding varying levels of white Gaussian noise.
To accurately simulate and compare the performance of TSI and SPGR sequences on the same MRI machine with the same magnet and sensing coils, the SNR levels must be adjusted between the two sequences to account for the differences in their physical scan parameters. The simulated SNR is calculated as the ratio between M 0 and the noise standard deviation σ. To enable a fair comparison, we assume the B 0 strength, the voxel size, and the number of data measurements to be the same between TSI and SPGR. This assumption leaves out three factors for calibration: receiver bandwidth BW, T Ã 2 relaxation time decay, and sensitivity encoded (SENSE) EPI for TSI. A larger BW incorporates more noise in the acquired signals and therefore decreases the SNR by a relative factor of ffiffiffiffiffiffiffi ffi BW p [37]. A longer echo time T E in gradient echo imaging decreases the acquired signal amplitude and therefore decreases the SNR through e À T E =T Ã 2 . For SPGR, the chosen pulse parameters are TR = 7.8 ms, T E = 2.4 ms and receiver BW of ±31.3KHZ [38]. For TSI, the chosen pulse parameters are TR = 6 s, T E = 35 ms and readout time of 2.048 μs/sample (equivalent to receiver BW of ±244.1KHZ) [14]. The BW of TSI is eight times larger than that of SPGR, resulting in a SNR for TSI of ffiffi ffi 8 p lower than SPGR. Assuming an average T Ã 2 value of 48.9 ms for brain parenchyma [39], the relative T Ã 2 decay gives TSI a lower SNR than SPGR by a factor of e À ðTE TSI À TE SPGR Þ=T Ã 2 ¼ e À ð35À 2:4Þ=48:9 % 0:5. The SENSE EPI rate of 2 will further decrease the SNR of TSI by a factor of ffiffi ffi 2 p due to fewer phase encoding (PE) steps during signal acquisitions [37]. Combining all these factors, the simulations of SPGR signals must have a SNR level 8 times greater than TSI to match the simulated performances with practical experiments.

Pulse sequence optimization results
The main result of applying the proposed sequence design framework is that the optimized TSI pulse parameters achieve improved precision and accuracy over both the original TSI sequence [14] and the DESPOT1 sequence [10]. Table 1 shows the optimized TSI sequence parameters (TSI new ) obtained by the optimization algorithm in the Methods section, compared against the original TSI sequence from [14]. The original TSI sequence was obtained assuming nominal tissue T 1 values for WM of 800 ms, GM of 1550 ms and CSF of 3700 ms. In contrast, the optimization process assumes a range of T 1 values 700-2000 ms for the normal brain tissues. Table 1 shows that the parameters of TSI new differ from the TSI original in the pulse times and flip angles. For TSI new , all pulses occur in the first 3655 ms over a TR of 6 seconds, leaving a relatively longer time for the longitudinal magnetization to relax after the third imaging pulse. In contrast, the pulses in the original TSI sequence are more spread out over TR. Moreover, the flip angles of TSI new increase across the sequence, while the flip angles of TSI original first decrease then increase.

T 1 estimation performances: Precision and accuracy
To evaluate and compare the T 1 estimate precision and accuracy, we designed two Monte Carlo simulation experiments involving four estimators: NLSE with the new TSI sequence, NLSE with the original TSI sequence, DESPOT1, and NLSE with the SPGR sequence used by DESPOT1. The precision is compared in terms of the relative mean estimation error: the standard deviation of the estimated T 1 over the true T 1 . Five thousand Monte Carlo trials are repeated for each set of pulse parameters for each value of T 1 . Equivalent SNR levels are calibrated as SNR = 125 for the TSI sequences and SNR = 1000 for the SPGR sequence. For both TSI methods and SPGR methods, the relative mean estimation errors are compared against their own theoretical lower bounds, or relative errors, calculated as the square root of the CRB over the true T 1 . Fig 3 compares the T 1 estimate precision in terms of the relative mean estimation errors for the four different estimators. For the evaluated T 1 range of 700-5000 ms, the TSI sequences show overall improvement over the SPGR sequence. The new TSI sequence achieves mean estimation error less than 1.7% for normal brain tissues (T 1 700-2000 ms). The new TSI sequence also maintains the best robustness to T 1 variation, with mean estimation error less than 6.5% when T 1 reaches 5000 ms. In the SPGR family, DESPOT1 and SPGR NLSE provide similar errors: less than 7.0% for T 1 of 700-2000 ms (agreeing with the findings in [10,11]) and less than 10.5% when T 1 reaches 5000 ms. Both SPGR NLSE and TSI new NLSE achieve their theoretical lower bounds for all tested T 1 values. This implies that the CRB provides a reliable prediction of the precision performance for different pulse sequences. Generalized bound-based MRI pulse sequence design for T 1 relaxometry A second Monte Carlo experiment compares the accuracy of T 1 estimates among the four approaches over a broad range of SNRs. The accuracy is measured in terms of the relative bias : jðMean ofT 1 À True T 1 Þ=True T 1 j. This experiment uses a nominal T 1 value of 1500 ms, which is a typical T 1 value for brain GM at 3.0T. Five thousand Monte Carlo trials are repeated for each set of pulse parameters for each SNR level. Again, to calibrate for the SNR equivalence between the TSI sequences and the SPGR sequence, simulated SNRs are selected with 5 SNR 60 for TSI sequences and 40 SNR 480 for both DESPOT1 and SPGR NLSE. Fig 4 compares the accuracy of T 1 estimates in terms of the relative bias for the four different estimation approaches. We can see that the T 1 estimates using the SPGR methods are biased, with the biases getting more pronounced as SNR decreases. Specifically, at the highest tested SNR level of 480, DESPOT1 has a relative bias of 0.58% and SPGR NLSE of 0.66%. At the lowest tested SNR level of 40, DESPOT1 reaches a relative bias of 53.03% and SPGR NLSE reaches 22.87%. For the same SPGR data, using NLSE to estimate the T 1 value has a lower bias than using the linear data fitting process in DESPOT1. In general, SPGR signals have a relatively high SNR due to short TR. For example, the SNR ranges from 100-200 in brain tissues for the clinical whole-brain SPGR data acquired at 1.5T with a single channel receiver coil (with TR = 8 ms and flip angles of 2˚, 3˚, 14˚, 17˚) [11]. As shown in Fig 4, in this SNR region, the SPGR approaches have a relative bias between 4-16%. Note that for a nominal T 1 of 1500 ms, this range of relative bias corresponds to an absolute bias of 60-240 ms. Bias of this

estimates' precision for four different approaches: NLSE with the new TSI sequence (blue) against its theoretical lower bound (red solid), NLSE with the original TSI sequence (cyan), DESPOT1 (green), and NLSE with the SPGR sequence (black) against its theoretical lower bound (magenta dashed).
The precision is measured in terms of the relative mean estimation error, calculated as the standard deviation of T 1 estimates normalized by the true T 1 . The theoretical lower bound of the relative mean estimation error is the relative error, calculated as the square root of the CRB on T 1 estimates normalized by the true T 1 . SNR levels equalizing for both receiver bandwidths and echo times are calibrated as SNR = 125 for the TSI sequences and SNR = 1000 for the SPGR sequence. The new TSI sequence achieves the lowest mean estimation error and therefore highest precision for tested T 1 values. doi:10.1371/journal.pone.0172573.g003 Generalized bound-based MRI pulse sequence design for T 1 relaxometry range is significant and can severely ambiguate the detection of T 1 changes in NAGM, which was recently shown to be associated with cortical lesions and cognitive dysfunction for patients with long-standing MS [40,41]. The average T 1 values for normal frontal GM, per the studies in [42], are in the range of 1322ms ± 34 at 3.0T. A T 1 estimate bias on the order of 60-240 ms would imply tissue abnormalities, which actually results from the T 1 relaxometry approach inaccuracies. This finding agrees with the results in [11] that the T 1 values for brain tissues estimated using DESPOT1 can be overestimated by 10-20% in the clinical SPGR images.
In contrast, using the TSI sequences for T 1 relaxometry virtually eliminates the T 1 estimate bias for a broad range of simulated SNRs. At the highest tested SNR level of 60, the original TSI sequence achieves a relative bias of 0.56% and the new TSI sequence of 0.05%. Even at the lowest tested SNR level of 5, the original TSI sequence produces a relative bias of 8.17%, and the new TSI sequence of 4.19%. The clinical SNR of 100-200 for the SPGR sequences translates to an SNR of 12.5-25 after compensating for the SNR equivalence calibration. Fig 4 shows in the SNR range between 12.5-25, the T 1 estimates using the new TSI sequence have a relative bias between 0.25-0.88%. Again, for a nominal T 1 of 1500 ms, this range of relative bias corresponds to an absolute bias of 4-13 ms. Bias of this range may or may not have any clinical significance, given that recent studies showing for NAWM and NAGM, the relevant early T 1 changes are on the order of 10-20 ms for MS patients [30,31]. However, for clinically realistic SNR levels, the new TSI sequence improves the T 1 estimation accuracy by a factor of 16 times over DESPOT1 for the tested nominal T 1 of 1500 ms. This improvement would greatly alleviate the ambiguities of T 1 variation caused either by the T 1 relaxometry inaccuracies or the underlying pathological conditions of the patients.

estimates' accuracy for four different approaches: NLSE with the new TSI sequence (blue), NLSE with the original TSI sequence (cyan), DESPOT1 (green), and NLSE with the SPGR sequence (black).
The accuracy is measured in terms of the relative bias %, calculated as jðMean ofT 1 À True T 1 Þ=True T 1 j. This experiment uses a nominal T 1 value of 1500 ms. Simulated SNR levels are calibrated equivalently as 5 SNR 60 for the TSI sequences and 40 SNR 480 for the SPGR sequence. The new TSI sequence achieves the lowest overall relative bias and therefore highest accuracy among the four approaches. Comparing T 1 sensitivity and orthogonality As described in the Theory section, the CRB provides a more complete model for the factors controlling T 1 estimate precision. There are three factors contributing to CRB improvement: SNR M 0 /σ, sensitivity jj @h @T 1 jj, and orthogonality sin ϕ. Improving SNR often requires increasing the B 0 field strength, employing lower noise receiver coils, or increasing the number of signals averaged, all of which increase hardware costs or scan time. Therefore, redesigning pulse sequences to increase the measurement sensitivity jj @h @T 1 jj and orthogonality sin ϕ can improve the T 1 estimate precision without requiring improved hardware or additional scan time. Here we evaluate and compare the sensitivity and orthogonality terms separately for each pulse sequence and investigate how the different terms contribute to the CRB as T 1 varies. Fig 5 compares the T 1 sensitivity and orthogonality for the new TSI sequence, the original TSI sequence, and the SPGR sequence as a function of T 1 over the range of 700-5000 ms. The top panel shows the new TSI sequence exhibits the highest sensitivity among the three sequences, especially for the T 1 range of 700-2000 ms corresponding to normal brain tissues. This high sensitivity of the new TSI sequence explains its low mean T 1 estimation error shown in Fig 3. For the original TSI sequence, the sensitivity decreases sharply for T 1 within 700-1500 ms, which implies this sequence is relatively unstable for T 1 relaxometry within the WM and GM regions. For all tested T 1 values, the SPGR sequence has the lowest T 1 sensitivity. The bottom panel in Fig 5 shows the TSI sequences have greater orthogonality than the SPGR sequence. Specifically, within the tested T 1 range, the new TSI sequence has an average orthogonality value of 0.88, compared with the original TSI sequence of 0.92 and the SPGR sequence of 0.57. Comparing the sensitivity, orthogonality, and the equivalent TSI-SPGR SNR factor, we find that the optimized TSI sequence owes its improved T 1 estimation capability to its high T 1 sensitivity and orthogonality of the TSI-family sequences. In contrast, the DESPOT1 sequence owes much of its T 1 estimation capability to its inherently high SNR, largely due to very short TR intervals. The dramatic improved sensitivity and modest improved orthogonality explain the superior performance of the new TSI sequence over DESPOT1 in spite of a conceding factor of 8 in SNR.

Sequence design approach
Prior to any validation from phantom or in vivo experiments, the Monte Carlo simulation results show that the new TSI sequence achieves improved T 1 precision and accuracy over the popular SPGR approaches. This improvement in precision is prominent for a broad range of brain T 1 values, not only the ones corresponding to normal brain parenchyma, but also to more advanced lesioned tissues and the CSF region. Both T 1 precision and accuracy show desirable stability to T 1 variation and varying SNR levels. The relative mean T 1 estimation errors using NLSE for both the new TSI sequence and the SPGR sequence confirm that the CRB is a reliable approach to predict the performance of a T 1 relaxometry approach.
This work uses a qualitatively new approach in the MRI pulse sequence design, which exploits the Cramér-Rao Bound (CRB) on the variance (stability) achievable by any estimator [28]. This bound-based design allows us to predict the performance of an MRI relaxometry approach based on the pulse sequence parameters, and then optimize the performance of that sequence by adjusting these parameters. We optimized the stability of a T 1 image by finding the pulse times and flip angles that produce the smallest variance for all T 1 s within a range of interest. This differs from prior pulse sequence design approaches that either presumed nominal tissue T 1 values a priori [14] or focused on improving the signal dynamic range [10,13]  The norm of T 1 sensitivity jj @h @T 1 jj is calculated as the Euclidean norm of the derivative of the signal weighting vector with respect to T 1 . Increasing the norm of the sensitivity will increase the impact of small changes in T 1 on the overall signal weighting vector h. sin ϕ is the orthogonality term defined in Fig 1. The more orthogonal these vectors are, the easier it becomes to ascribe changes in the observed signal to M 0 or T 1 unambiguously. The top panel shows the new TSI sequence has the best sensitivity among the three sequences and SPGR has very poor sensitivity for T 1 estimation. The bottom panel shows the TSI-family sequences have greater orthogonality (above 0.7) than the SPGR sequence (equal to 0.58) for the tested T 1 range.
doi:10.1371/journal.pone.0172573.g005 without assessing the implications for the T 1 map stability. Another advantage of the proposed pulse design approach is considering a range of nominal T 1 values in the pulse optimization process. For prior approaches [10,14], the optimized pulse times and flip angles have strong dependence on the assumed tissue T 1 values. Although T 1 measurements for brain tissues are widely available [32], the inherent T 1 variation within the target tissue degrades the performance of the sequences designed assuming a specific T 1 value for each tissue type. In contrast, the min-max approach provides the best worst case performance, resulting in robust stability for T 1 throughout the region of interest and not just estimates in the neighborhoods of the nominal T 1 values.
From the pulse sequence design perspective, the proposed bound-based framework was derived assuming a general signal model under additive white Gaussian noise. This framework is illustrated with a case study for the TSI sequence to show the effectiveness of the CRB-based framework in optimizing sequence parameters for T 1 relaxometry. However, the proposed framework is readily extendible to other MRI relaxometry sequence designs and sequence parameter optimizations. For example, the proposed framework can be applied to optimize the inversion time and number of data measurements when using the inversion recovery (IR) sequence for T 1 relaxometry. Another example is using the proposed framework to optimize the number of excitation pulses and their inter-pulse times and flip angles when using the Look-Locker (LL) sequence for fast T 1 relaxometry. Alternatively, the proposed framework is also extendible for designing new MRI sequences and optimizing their parameters to simultaneously extract T 1 and T 2 information from the signal measurements. One such example is the DESPOT approach that jointly uses the SPGR and steady state free precession sequences (SSFP) to collect a series of steady state images over a range of flip angles for fast T 1 and T 2 mappings [10]. The proposed framework can be adapted to optimize the flip angles and pulse repetition times for the SPGR and SSFP sequences for optimal precision of joint T 1 /T 2 estimation. We believe our pulse design framework is promising for producing sequences with improved stability, precision, and accuracy over a wide range of T 1 /T 2 values corresponding to normal tissue and advanced pathologies as well.

Practical considerations
A dominant practical error impacting the T 1 estimation performance of a pulse sequence is the pulse flip angle perturbations. In practice, the flip angle perturbations are mainly caused by patient-induced B1 inhomogeneities due to distortions of the radio-frequency field generated by the transmit coils. Although the use of B1-insensitive adiabatic pulses can result in accurate inversions, the B1 inhomogeneities could affect the imaging pulses with the flip angles perturbations up to ±20% of their nominal values [14]. Observing how the relative error ffiffiffiffiffiffiffiffiffi CRB p =T 1 degrades due to the flip angle perturbations quantifies how B1 field inhomogeneities affect T 1 estimate precision. Fig 6 demonstrates the degradation of the relative errors due to flip angle perturbations up to ±20% of their nominal values for the new TSI sequence at SNR = 125 (top panel), the original TSI sequence at SNR = 125 (middle panel), and the SPGR sequence at SNR = 1000 (bottom panel). The top panel shows the new TSI sequence is robust to the B1 field inhomogeneities for a wide T 1 range within 700-5000 ms. The flip angle perturbations alter the relative error curve, but not dramatically. Specifically, the relative errors fall within ±0.1% of the unperturbed nominal values for the normal brain tissue range (T 1 700-2000 ms). For T 1 up to 5000 ms, the relative error still falls within ±0.4% of the unperturbed value of 6.5%. Compared against the new TSI sequence, the flip angle perturbations affect the original TSI sequence more severely. As shown in the middle panel, the relative errors fall within ±0.9% of the unperturbed nominal values for the normal The relative error is calculated as the square root of the CRB on T 1 estimates normalized by the true T 1 value. All curves are generated keeping the flip angle of the inversion pulse as 180˚and simultaneously varying the flip angles of brain tissue range (T 1 700-2000 ms). The relative error is most sensitive to flip angle perturbations for T 1 of 1300-2100 ms, where most brain GM resides. Compared against the TSIfamily sequences, the bottom panel shows that the flip angle perturbations affect the SPGR sequence more severely for all T 1 s within 700-5000 ms. This result agrees with the earlier findings on the SPGR-based techniques' sensitivity to B1 variations [10]. Specifically, for the normal brain tissue range (700-2000 ms), the relative errors fall within ±0.4% of the unperturbed nominal values. As T 1 increases, the relative error degrades more severely (off the chart) and falls within ±0.65% of the unperturbed value of 10.5% for T 1 = 5000 ms. The relative robustness to B1 inhomogeneities and stability to T 1 variation make the new TSI sequence appropriate for studying the effect of diseases such as MS, where the presence of lesions at different degrees of severity may lead to T 1 variations in the order of thousands of milliseconds, rendering measurements by other techniques unreliable.
Regarding the underlying MR physics for different imaging sequences, we assume that the same physical object is scanned by the same MRI scanner under the same physical conditions. The B 0 field strength, the voxel volume size, and the number of data measurements are assumed the same between TSI and DESPOT1. The simulations adjust the distinct SNR levels between TSI and DESPOT1 for different echo times during T Ã 2 decay and different receiver bandwidths. However, we ignore several MR physical factors which might affect the performance of the optimized pulse sequence in practice. For example, TSI usually requires about three additional repetitions of the five-pulse sequence before the start of data acquisition to allow magnetization to reach equilibrium condition [14]. Similarly, for the DESPOT1, we have not considered either the time required to reach steady state for each acquisition or the wait time between different sub-sequences in the actual acquisitions. Therefore, a natural next step is to confirm the improved performance of the new TSI sequence in actual MRI experiments. As a first validation, our next effort will focus on phantom experiments with known T 1 values to demonstrate the reliability of the proposed technique.

Conclusion
This paper designed a new TSI T 1 relaxometry MRI sequence that achieves improved precision and accuracy over a broad range of T 1 values, covering both healthy and lesioned brain tissues. The new sequence demonstrates robustness to B1 inhomogeneity and stability over the wide range of T 1 variations encountered in neuro-degenerative diseases under clinically feasible scan time. This suggests that the improved TSI sequence may be helpful in the study of neurodegenerative diseases such as multiple sclerosis. The improved performance of the new sequence illustrates the value of the min-max design strategy minimizing the Cramér-Rao bound on T 1 estimation. The geometric interpretation of the Cramér-Rao bound developed here illuminates three factors controlling relaxometry performance: improving the SNR (or signal dynamic range), increasing the signal's sensitivity to T 1 variation, and increasing the orthogonality between the signal vector and the sensitivity vector. The second and third approaches can be implemented on existing systems to improve the T 1 relaxometry performance without requiring hardware improvements to magnets or measurement coils. The improved relaxometry performance predicted by the Cramér-Rao bound and demonstrated in Monte-Carlo simulations suggests that the proposed benefits of the design strategy may prove portable to other relaxometry sequences.