Quantitative T1 MRI - Quantitative T1 MRI

Abstract¶

This NeuroLibre Reproducible preprint is an interactive tutorial on quantitative T1 mapping MRI. It is an interactive version of two subsections of the chapter “Quantitative T1 and T1r Mapping” in the book Quantitative Magnetic Resonance Imaging Boudreau et al., 2020.

\awesomebox[red]{2pt}{\faExclamationCircle}{red}{\textbf{NOTE}}

The following section in this document repeats the narrative content exactly as found in the corresponding NeuroLibre Reproducible Preprint (NRP). The content was automatically incorporated into this PDF using the NeuroLibre publication workflow Karakuzu et al., 2022 to credit the referenced resources. The submitting author of the preprint has verified and approved the inclusion of this section through a GitHub pull request made to the source repository from which this document was built. Please note that the figures and tables have been excluded from this (static) document. To interactively explore such outputs and re-generate them, please visit the corresponding NRP. For more information on integrated research objects (e.g., NRPs) that bundle narrative and executable content for reproducible and transparent publications, please refer to DuPre et al. (2022). NeuroLibre is sponsored by the Canadian Open Neuroscience Platform (CONP) Harding et al., 2023.

This Jupyter Book is a series of interactive tutorials about quantitative T1 mapping, powered by qMRLab Karakuzu et al., 2020.

Widely considered the gold standard for T1 mapping, the inversion recovery technique estimates T1 values by fitting the signal recovery curve acquired at different delays after an inversion pulse (180°). In a typical inversion recovery experiment (Figure 1), the magnetization at thermal equilibrium is inverted using a 180° RF pulse. After the longitudinal magnetization recovers through spin-lattice relaxation for predetermined delay (“inversion time”, TI), a 90° excitation pulse is applied, followed by a readout imaging sequence (typically a spin-echo or gradient-echo readout) to create a snapshot of the longitudinal magnetization state at that TI.

Inversion recovery was first developed for NMR in the 1940s Hahn, 1949Drain, 1949, and the first T1 map was acquired using a saturation-recovery technique (90° as a preparation pulse instead of 180°) by Pykett & Mansfield, 1978. Some distinct advantages of inversion recovery are its large dynamic range of signal change and an insensitivity to pulse sequence parameter imperfections Stikov et al., 2015. Despite its proven robustness at measuring T1, inversion recovery is scarcely used in practice, because conventional implementations require repetition times (TRs) on the order of 2 to 5 T1 Steen et al., 1994, making it challenging to acquire whole-organ T1 maps in a clinically feasible time. Nonetheless, it is continuously used as a reference measurement during the development of new techniques, or when comparing different T1 mapping techniques, and several variations of the inversion recovery technique have been developed, making it practical for some applications Messroghli et al., 2004Piechnik et al., 2010.

M\_{z}(TI) = M\_0 \frac{1-\text{cos}(\theta\_{180})e^{- \frac{TR}{T\_1}} -[1-\text{cos}(\theta\_{180})]e^{- \frac{TI}{T\_1}}}{1 - \text{cos}(\theta\_{180}) \text{cos}(\theta\_{90}) e^{- \frac{TR}{T\_1}}}

(1)

where Mz is the longitudinal magnetization prior to the θ90 pulse. If the in-phase real signal is desired, it can be calculated by multiplying Eq. 1 by ksin(θ90)e-TE/T2, where k is a constant. This general equation can be simplified by grouping together the constants for each measurements regardless of their values (i.e. at each TI, same TE and θ90 are used) and assuming an ideal inversion pulse:

M\_z(TI) = C(1-2e^{- \frac{TI}{T\_1}} + e^{- \frac{TR}{T\_1}})

(2)

M\_z(TI) = C(1-2e^{- \frac{TI}{T\_1}})

(3)

S(TI) = a + be^{- \frac{TI}{T\_1}}

(4)

where a and b are complex values. If magnitude-only data is available, a 3-parameter model can be sufficient by taking the absolute value of Eq. 4. While the RD-NLS algorithms are too complex to be presented here (the reader is referred to the paper, (Barral et al. 2010)), the code for these algorithms was released open-source along with the original publication, and is also available as a qMRLab T1 mapping model. One important thing to note about Eq. 4 is that it is general – no assumption is made about TR – and is thus as robust as Eq. 1 as long as all pulse sequence parameters other than TI are kept constant between each measurement. Figure 4 compares simulated data (Eq. 1) using a range of TRs (1.5T1 to 5T1) fitted using either RD-NLS & Eq. 4 or a Levenberg-Marquardt fit of Eq. 2.

Figure 5 displays an example brain dataset from an inversion recovery experiment, along with the T1 map fitted using the RD-NLS technique.

Despite a widely acknowledged robustness for measuring accurate T1 maps, inversion recovery is not often used in studies. An important drawback of this technique is the need for long TR values, generally on the order of a few T1 for general models (e.g. Eqs. 1 and 4), and up to 5T1 for long TR approximated models (Eq. 3). It takes about to 10-25 minutes to acquire a single-slice T1 map using the inversion recovery technique, as only one TI is acquired per TR (2-5 s) and conventional cartesian gradient readout imaging acquires one phase encode line per excitation (for a total of ~100-200 phase encode lines). The long acquisition time makes it challenging to acquire whole-organ T1 maps in clinically feasible protocol times. Nonetheless, it is useful as a reference measurement for comparisons against other T1 mapping methods, or to acquire a single-slice T1 map of a tissue to get T1 estimates for optimization of other pulse sequences.

0.1Other Saturation-Recovery T1 Mapping techniques¶

1Variable Flip Angle T1 Mapping¶

Variable flip angle (VFA) T1 mapping Christensen et al., 1974Fram et al., 1987Gupta, 1977, also known as Driven Equilibrium Single Pulse Observation of T1 (DESPOT1) Homer & Beevers, 1985Deoni et al., 2003, is a rapid quantitative T1 measurement technique that is widely used to acquire 3D T1 maps (e.g. whole-brain) in a clinically feasible time. VFA estimates T1 values by acquiring multiple spoiled gradient echo acquisitions, each with different excitation flip angles (θn for n = 1, 2, .., N and θi ≠ θj). The steady-state signal of this pulse sequence (Figure 1) uses very short TRs (on the order of magnitude of 10 ms) and is very sensitive to T1 for a wide range of flip angles.

Figure 1. Simplified pulse sequence diagram of a variable flip angle (VFA) pulse sequence with a gradient echo readout. TR: repetition time, θn: excitation flip angle for the nth measurement, IMG: image acquisition (k-space readout), SPOIL: spoiler gradient.

\tag{1} M\_{z}(\theta\_n) = M\_0 \frac{1-e^{- \frac{TR}{T\_1}}}{1-\text{cos}(\theta\_n) e^{- \frac{TR}{T\_1}}} \text{sin}(\theta\_n)

(5)

where Mz is the longitudinal magnetization, M0 is the magnetization at thermal equilibrium, TR is the pulse sequence repetition time (Figure 1), and θn is the excitation flip angle. The Mz curves of different T1 values for a range of θn and TR values are shown in Figure 2.

From Figure 2, it is clearly seen that the flip angle at which the steady-state signal is maximized is dependent on the T1 and TR values. This flip angle is a well known quantity, called the Ernst angle Ernst & Anderson, 1966, which can be solved analytically from Equation 1 using properties of calculus:

\tag{2} \theta\_{Ernst} = \text{acos}(e^{- \frac{TR}{T\_1}})

(6)

Figure 4. Signal curves estimated using Bloch simulations for three categories of signal spoiling: (1) ideal spoiling (blue), gradient & RF Spoiling (orange), and no spoiling (green). Simulations details: TR = 25 ms, T1 = 900 ms, Te = 100 ms, TE = 5 ms, 100 spins. For the ideal spoiling case, the transverse magnetization is set to zero at the end of each TR. For the gradient & RF spoiling case, each spin is rotated by different increments of phase (2𝜋 / # of spins) to simulate complete decoherence from gradient spoiling, and the RF phase of the excitation pulse is ɸn = ɸn-1 + nɸ0 = ½ ɸ0(n2 + n + 2) Bernstein et al., 2004 with ɸ0 = 117° Zur et al., 1991 after each TR.

\tag{3} \frac{S\_n}{ \text{sin}(\theta\_n)} = e^{- \frac{TR}{T\_1}} \frac{S\_n}{ \text{tan}(\theta\_n)} + C (1-e^{- \frac{TR}{T\_1}})

(7)

As the third term does not change between measurements (it is constant for each θn), it can be grouped into the constant for a simpler representation:

\tag{4} \frac{S\_n}{ \text{sin}(\theta\_n)} = e^{- \frac{TR}{T\_1}} \frac{S\_n}{ \text{tan}(\theta\_n)} + C

(8)

With this rearranged form of Equation 1, T1 can be simply estimated from the slope of a linear regression calculated from Sn/sin(θn) and Sn/tan(θn) values:

\tag{5} T\_1 = - \frac{TR}{ \text{ln}(slope)}

(9)

\tag{6} \theta\_{corrected} = B\_1 \theta\_{nominal}

(10)

B1 in this context is normalized, meaning that it is unitless and has a value of 1 in voxels where the RF field has the expected amplitude (i.e. where the nominal flip angle is the actual flip angle). Figure 6 displays fitted VFA T1 values from a Monte Carlo dataset simulated using biased flip angle values, and fitted without/with B1 correction.

Figure 6. Mean and standard deviations of fitted VFA T1 values for a set of Monte Carlo simulations (SNR = 100, N = 1000), simulated using a wide range of biased flip angles and fitted without (blue) or with (red) B1 correction. Simulation parameters: TR = 25 ms, T1 = 900 ms, θnominal = 6° and 32° (optimized values for this TR/T1 combination). Notice how even after B1 correction, fitted T1 values at B1 values far from the nominal case (B1 = 1) exhibit larger variance, as the actual flip angles of the simulated signal deviate from the optimal values for this TR/T1 (Deoni et al. 2003).

Figure 7 displays an example VFA dataset and a B1 map in a healthy brain, along with the T1 map estimated using a linear fit (Equations 4 and 5).

Figure 7. Example variable flip angle dataset and B1 map of a healthy adult brain (left). The relevant VFA protocol parameters used were: TR = 15 ms, θnominal = 3° and 20°. The T1 map (right) was fitted using a linear regression (Equations 4 and 5).

It has been well reported in recent years that the accuracy of VFA T1 estimates is very sensitive to pulse sequence implementations Baudrexel et al., 2017Lutti & Weiskopf, 2013Stikov et al., 2015, and as such is less robust than the gold standard inversion recovery technique. In particular, the signal bias resulting from insufficient spoiling can result in inaccurate T1 estimates of up to 30% relative to inversion recovery estimated values Stikov et al., 2015. VFA T1 map accuracy and precision is also strongly dependent on the quality of the measured B1 map Lee et al., 2017, which can vary substantially between implementations Boudreau et al., 2017. Modern rapid B1 mapping pulse sequences are not as widely available as VFA, resulting in some groups attempting alternative ways of removing the bias from the T1 maps like generating an artificial B1 map through the use of image processing techniques Liberman et al., 2013 or omitting B1 correction altogether Yuan et al., 2012. The latter is not recommended, because most MRI scanners have default pulse sequences that, with careful protocol settings, can provide B1 maps of sufficient quality very rapidly Boudreau et al., 2017Samson et al., 2006Wang et al., 2005.

Despite some drawbacks, VFA is still one of the most widely used T1 mapping methods in research. Its rapid acquisition time, rapid image processing time, and widespread availability makes it a great candidate for use in other quantitative imaging acquisition protocols like quantitative magnetization transfer imaging Cercignani et al., 2005Yarnykh, 2002 and dynamic contrast enhanced imaging Li et al., 2018Sung et al., 2013.

Keywords:MRIqMRIQuantitativeT1Relaxometry¶

2About¶

This Jupyter Book is a series of interactive tutorials about quantitative T₁ mapping, powered by qMRLab Karakuzu et al., 2020.

3Inversion Recovery T₁ Mapping¶

Widely considered the gold standard for T₁ mapping, the inversion recovery technique estimates T₁ values by fitting the signal recovery curve acquired at different delays after an inversion pulse (180°). In a typical inversion recovery experiment (Figure 1), the magnetization at thermal equilibrium is inverted using a 180° RF pulse. After the longitudinal magnetization recovers through spin-lattice relaxation for predetermined delay (“inversion time”, TI), a 90° excitation pulse is applied, followed by a readout imaging sequence (typically a spin-echo or gradient-echo readout) to create a snapshot of the longitudinal magnetization state at that TI.

Inversion recovery was first developed for NMR in the 1940s Hahn, 1949Drain, 1949, and the first T₁ map was acquired using a saturation-recovery technique (90° as a preparation pulse instead of 180°) by Pykett & Mansfield, 1978. Some distinct advantages of inversion recovery are its large dynamic range of signal change and an insensitivity to pulse sequence parameter imperfections Stikov et al., 2015. Despite its proven robustness at measuring T₁, inversion recovery is scarcely used in practice, because conventional implementations require repetition times (TRs) on the order of 2 to 5 T₁ Steen et al., 1994, making it challenging to acquire whole-organ T₁ maps in a clinically feasible time. Nonetheless, it is continuously used as a reference measurement during the development of new techniques, or when comparing different T₁ mapping techniques, and several variations of the inversion recovery technique have been developed, making it practical for some applications Messroghli et al., 2004Piechnik et al., 2010.

Figure 1. Pulse sequence of an inversion recovery experiment.

4Signal Modelling¶

The steady-state longitudinal magnetization of an inversion recovery experiment can be derived from the Bloch equations for the pulse sequence {θ₁₈₀ – TI – θ₉₀ – (TR-TI)}, and is given by:

M_{z}(TI) = M_0 \frac{1-\text{cos}(\theta_{180})e^{- \frac{TR}{T_1}} -[1-\text{cos}(\theta_{180})]e^{- \frac{TI}{T_1}}}{1 - \text{cos}(\theta_{180}) \text{cos}(\theta_{90}) e^{- \frac{TR}{T_1}}}

(11)

where M_z is the longitudinal magnetization prior to the θ₉₀ pulse. If the in-phase real signal is desired, it can be calculated by multiplying Eq. 1 by ksin(θ₉₀)e^{-TE/T ₂}, where k is a constant. This general equation can be simplified by grouping together the constants for each measurements regardless of their values (i.e. at each TI, same TE and θ₉₀ are used) and assuming an ideal inversion pulse:

M_z(TI) = C(1-2e^{- \frac{TI}{T_1}} + e^{- \frac{TR}{T_1}})

(12)

where the first three terms and the denominator of Eq. 1 have been grouped together into the constant C. If the experiment is designed such that TR is long enough to allow for full relaxation of the magnetization (TR > 5T₁), we can do an additional approximation by dropping the last term in Eq. 2:

M_z(TI) = C(1-2e^{- \frac{TI}{T_1}})

(13)

The simplicity of the signal model described by Eq. 3, both in its equation and experimental implementation, has made it the most widely used equation to describe the signal evolution in an inversion recovery T₁ mapping experiment. The magnetization curves are plotted in Figure 2 for approximate T₁ values of three different tissues in the brain. Note that in many practical implementations, magnitude-only images are acquired, so the signal measured would be proportional to the absolute value of Eq. 3.

Source:

Signal Modelling

Figure 1:Inversion recovery curves (Eq. 2) for three different T₁ values, approximating the main types of tissue in the brain.

Practically, Eq. 1 is the better choice for simulating the signal of an inversion recovery experiment, as the TRs are often chosen to be greater than 5T₁ of the tissue-of-interest, which rarely coincides with the longest T₁ present (e.g. TR may be sufficiently long for white matter, but not for CSF which could also be present in the volume). Equation 3 also assumes ideal inversion pulses, which is rarely the case due to slice profile effects. Figure 3 displays the inversion recovery signal magnitude (complete relaxation normalized to 1) of an experiment with TR = 5 s and T₁ values ranging between 250 ms to 5 s, calculated using both equations.

Source:

Signal Modelling

Figure 2:Signal recovery curves simulated using Eq. 3 (solid) and Eq. 1 (dotted) with a TR = 5 s for T₁ values ranging between 0.25 to 5 s.

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 2.


% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in seconds
% All flip angles are in degrees

params.TR = 5.0;
params.TI = linspace(0.001, params.TR, 1000);
            
params.TE = 0.004;
params.T2 = 0.040;
            
params.EXC_FA = 90;  % Excitation flip angle
params.INV_FA = 180; % Inversion flip angle

params.signalConstant = 1;

%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the
% gradient echo readout inversion recovery experiment The option
% '4' is a flag that selects the long TR approximation of the 
% analytical solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run:
% `help inversion_recovery.analytical_solution`

% White matter
params.T1 = 0.900; % in seconds

signal_WM = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);

% Grey matter
params.T1 = 1.500;  % in seconds
signal_GM = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);

% CSF
params.T1 = 4.000;  % in seconds
signal_CSF = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 3.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in seconds
% All flip angles are in degrees

params.TR = 5.0;
params.TI = linspace(0.001, params.TR, 1000);
            
params.TE = 0.004;
params.T2 = 0.040;
            
params.EXC_FA = 90;  % Excitation flip angle
params.INV_FA = 180; % Inversion flip angle

params.signalConstant = 1;

T1_range = 0.25:0.25:5; % in seconds

%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the
% gradient echo readout inversion recovery experiment. The option
% '1' is a flag that selects full analytical solution equation 
% (no approximation), Eq. 1 of the blog post. The option '4' is a
% flag that selects the long TR approximation of the analytical 
% solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run:
% `help inversion_recovery.analytical_solution`

for ii = 1:length(T1_range)
    params.T1 = T1_range(ii);
    
    signal_T1_Eq1{ii} = inversion_recovery.analytical_solution(params, 'GRE-IR', 1);

    signal_T1_Eq3{ii} = inversion_recovery.analytical_solution(params, 'GRE-IR', 4);
end

5Data Fitting¶

Several factors impact the choice of the inversion recovery fitting algorithm. If only magnitude images are available, then a polarity-inversion is often implemented to restore the non-exponential magnitude curves (Figure 3) into the exponential form (Figure 2). This process is sensitive to noise due to the Rician noise creating a non-zero level at the signal null. If phase data is also available, then a phase term must be added to the fitting equation Barral et al., 2010. Equation 3 must only be used to fit data for the long TR regime (TR > 5T₁), which in practice is rarely satisfied for all tissues in subjects.

Early implementations of inversion recovery fitting algorithms were designed around the computational power available at the time. These included the “null method” Pykett et al., 1983, assuming that each T₁ value has unique zero-crossings (see Figure 2), and linear fitting of a rearranged version of Eq. 3 on a semi-log plot Fukushima, 1981. Nowadays, a non-linear least-squares fitting algorithm (e.g. Levenberg-Marquardt) is more appropriate, and can be applied to either approximate or general forms of the signal model (Eq. 3 or Eq. 1). More recent work Barral et al., 2010 demonstrated that T₁ maps can also be fitted much faster (up to 75 times compared to Levenberg-Marquardt) to fit Eq. 1 – without a precision penalty – by using a reduced-dimension non-linear least squares (RD-NLS) algorithm. It was demonstrated that the following simplified 5-parameter equation can be sufficient for accurate T₁ mapping:

S(TI) = a + be^{- \frac{TI}{T_1}}

(14)

where a and b are complex values. If magnitude-only data is available, a 3-parameter model can be sufficient by taking the absolute value of Eq. 4. While the RD-NLS algorithms are too complex to be presented here (the reader is referred to the paper, (Barral et al. 2010)), the code for these algorithms was released open-source along with the original publication, and is also available as a qMRLab T₁ mapping model. One important thing to note about Eq. 4 is that it is general – no assumption is made about TR – and is thus as robust as Eq. 1 as long as all pulse sequence parameters other than TI are kept constant between each measurement. Figure 4 compares simulated data (Eq. 1) using a range of TRs (1.5T₁ to 5T₁) fitted using either RD-NLS & Eq. 4 or a Levenberg-Marquardt fit of Eq. 2.

Source:

Data Fitting

Figure 3:Fitting comparison of simulated data (blue markers) with T₁ = 1 s and TR = 1.5 to 5 s, using fitted using RD-NLS & Eq. 4 (green) and Levenberg-Marquardt & Eq. 2 (orange, long TR approximation).

Figure 5 displays an example brain dataset from an inversion recovery experiment, along with the T₁ map fitted using the RD-NLS technique.

Source:

Data Fitting

Figure 4:Example inversion recovery dataset of a healthy adult brain (left). Inversion times used to acquire this magnitude image dataset were 30 ms, 530 ms, 1030 ms, and 1530 ms, and the TR used was 1550 ms. The T₁ map (right) was fitted using a RD-NLS algorithm.

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 4.


% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

params.TI = 50:50:1500;
TR_range = 1500:50:5000;

params.EXC_FA = 90;
params.INV_FA = 180;

params.T1 = 1000;

%% Calculate signals
%
% The option 'GRE-IR' selects the analytical equations for the gradient echo readout inversion recovery experiment
% The option '1' is a flag that selects full analytical solution equation (no approximation), Eq. 1 of the blog post.
%
% To see all the options available, run `help inversion_recovery.analytical_solution`

for ii = 1:length(TR_range)
    params.TR = TR_range(ii);
    Mz_analytical(ii,:) = inversion_recovery.analytical_solution(params, 'GRE-IR', 1);
end

%% Fit data using Levenberg-Marquardt with the long TR approximation equation
%
% The option '4' is a flag that selects the long TR approximation of the analytical solution (TR>>T1), Eq. 3 of the blog post.
%
% To see all the options available, run `help inversion_recovery.fit_lm`

for ii=1:length(TR_range)
    fitOutput_lm{ii} = inversion_recovery.fit_lm(Mz_analytical(ii,:), params, 4);
    T1_lm(ii) = fitOutput_lm{ii}.T1;
end

%% Fit data using the RDLS method (Barral), Eq. 4 of the blog post.
%

% Create a qMRLab inversion recovery model object and load protocol values
irObj = inversion_recovery();
irObj.Prot.IRData.Mat = params.TI';

for ii=1:length(TR_range)

    data.IRData = Mz_analytical(ii,:);

    fitOutput_barral{ii} = irObj.fit(data);

    T1_barral(ii) = fitOutput_barral{ii}.T1;

end

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 5.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup

clear all

%% MATLAB/OCTAVE CODE

% Load data into environment, and rotate mask to be aligned with IR data
load('data/ir_dataset/IRData.mat');
load('data/ir_dataset/IRMask.mat');

IRData = data;
Mask = imrotate(Mask,180);
clear data

% Format qMRLab inversion_recovery model parameters, and load them into the Model object
Model = inversion_recovery; 
TI = [30; 530; 1030; 1530];
Model.Prot.IRData.Mat = [TI];

% Format data structure so that they may be fit by the model
data = struct();
data.IRData= double(IRData);
data.Mask= double(Mask);

FitResults = FitData(data,Model,0); % The '0' flag is so that no wait bar is shown.

% Code used to re-orient the images to make pretty figures, and to assign variables with the axis lengths.

T1_map = imrotate(FitResults.T1.*Mask,-90);
xAxis = [0:size(T1_map,2)-1];
yAxis = [0:size(T1_map,1)-1];

% Raw MRI data at different TI values
TI_0030 = imrotate(squeeze(IRData(:,:,:,1).*Mask),-90);
TI_0530 = imrotate(squeeze(IRData(:,:,:,2).*Mask),-90);
TI_1030 = imrotate(squeeze(IRData(:,:,:,3).*Mask),-90);
TI_1530 = imrotate(squeeze(IRData(:,:,:,4).*Mask),-90);

References

6Benefits and Pitfalls¶

The conventional inversion recovery experiment is considered the gold standard T₁ mapping technique for several reasons:

A typical protocol has a long TR value and a sufficient number of inversion times for stable fitting (typically 5 or more) covering the range [0, TR].
It offers a wide dynamic range of signals (up to [-kM₀, kM₀]), allowing a number of inversion times where high SNR is available to sample the signal recovery curve Fukushima, 1981.
T₁ maps produced by inversion recovery are largely insensitive to inaccuracies in excitation flip angles and imperfect spoiling Stikov et al., 2015, as all parameters except TI are constant for each measurement and only a single acquisition is performed (at TI) during each TR.

One important protocol design consideration is to avoid acquiring at inversion times where the signal for T₁ values of the tissue-of-interest is nulled, as the magnitude images at this TI time will be dominated by Rician noise which can negatively impact the fit under low SNR circumstances (Figure 6). Inversion recovery can also often be acquired using commonly available standard pulse sequences available on most MRI scanners by setting up a customized acquisition protocol, and does not require any additional calibration measurements. For an example, please visit the interactive preprint of the ISMRM Reproducible Research Group 2020 Challenge on inversion recovery T1 mapping Boudreau et al., 2023.

Source:

Benefits and Pitfalls

Figure 5:Monte Carlo simulations (mean and standard deviation (STD), blue markers) and fitted T₁ values (mean and STD, red and green respectively) generated for a T₁ value of 900 ms and 5 TI values linearly spaced across the TR (ranging from 1 to 5 s). A bump in T₁ STD occurs near TR = 3000 ms, which coincides with the TR where the second TI is located near a null point for this T₁ value.

Despite a widely acknowledged robustness for measuring accurate T₁ maps, inversion recovery is not often used in studies. An important drawback of this technique is the need for long TR values, generally on the order of a few T₁ for general models (e.g. Eqs. 1 and 4), and up to 5T₁ for long TR approximated models (Eq. 3). It takes about to 10-25 minutes to acquire a single-slice T₁ map using the inversion recovery technique, as only one TI is acquired per TR (2-5 s) and conventional cartesian gradient readout imaging acquires one phase encode line per excitation (for a total of ~100-200 phase encode lines). The long acquisition time makes it challenging to acquire whole-organ T₁ maps in clinically feasible protocol times. Nonetheless, it is useful as a reference measurement for comparisons against other T₁ mapping methods, or to acquire a single-slice T₁ map of a tissue to get T₁ estimates for optimization of other pulse sequences.

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 6.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in seconds
% All flip angles are in degrees

TR_range = 1000:100:5000; % in seconds

x = struct;
x.T1 = 900; % in seconds

Opt.SNR = 25;
Opt.M0 = 1;
Opt.FAexcite = 90; % Excitation flip angle
Opt.FAinv = 180;   % Inversion flip angle

%% Monte Carlo data simulation
% Simulate noisy signal data 1,000 time, fit the data, then calculate the means and standard deviations of the data and fitted T1
% Data is calculated by calculating the a and b values of Eq. 4 from the full analytical equations (Eq. 1)

Model = inversion_recovery; 

for ii = 1:length(TR_range)
    Opt.TR = TR_range(ii);
    Opt.T1 = x.T1;
    TI_lowres(ii,:) = linspace(0.05, Opt.TR, 6)';
    Model.Prot.IRData.Mat = [TI_lowres(ii,:)];
    [ra,rb] = Model.ComputeRaRb(x,Opt);
    x.rb = rb;
    x.ra = ra;
    for jj = 1:1000
        [FitResult{ii,jj}, noisyData{ii,jj}] = Model.Sim_Single_Voxel_Curve(x,Opt,0); 
        fittedT1(ii,jj) = FitResult{ii,jj}.T1;
        noisyData_array(ii,jj,:) = noisyData{ii,jj}.IRData;
    end
        
    for kk=1:length(TI_lowres(ii,:))
        data_mean(ii,kk) = mean(noisyData_array(ii,:,kk));
        data_std(ii,kk) = std(noisyData_array(ii,:,kk));
    end
    
    T1_mean(ii) = mean(fittedT1(ii,:));
    T1_std(ii) = std(fittedT1(ii,:));
end

%% Calculate the noiseless data at a higher TI resolution to plot the ideal signal curve.
%

for ii = 1:length(TR_range)
    TI_highres(ii,:) = linspace(0.05, TR_range(ii), 500);
    Model.Prot.IRData.Mat = [TI_highres(ii,:)];
    Opt.TR = TR_range(ii);
    Opt.T1 = x.T1;
    [ra,rb] = Model.ComputeRaRb(x,Opt);
    x.rb = rb;
    x.ra = ra;

    data_noiseless(ii,:) = Model.equation(x);
end

References

7Other Saturation-Recovery T₁ Mapping techniques¶

Several variations of the inversion recovery pulse sequence were developed to overcome challenges like those specified above. Amongst them, the Look-Locker technique Look & Locker, 1970 stands out as one of the most widely used in practice. Instead of a single 90° acquisition per TR, a periodic train of small excitation pulses θ are applied after the inversion pulse, {θ₁₈₀ – 𝛕 – θ – 𝛕 – θ – ...}, where 𝛕 = TR/n and n is the number of sampling acquisitions. This pulse sequence samples the inversion time relaxation curve much more efficiently than conventional inversion recovery, but at a cost of lower SNR. However, because the magnetization state of each TI measurement depends on the previous series of θ excitation, it has higher sensitivity to B₁-inhomogeneities and imperfect spoiling compared to inversion recovery Gai et al., 2013Stikov et al., 2015. Nonetheless, Look-Locker is widely used for rapid T₁ mapping applications, and variants like MOLLI (Modified Look-Locker Inversion recovery) and ShMOLLI (Shortened MOLLI) are widely used for cardiac T₁ mapping Messroghli et al., 2004Piechnik et al., 2010.

Another inversion recovery variant that’s worth mentioning is saturation recovery, in which the inversion pulse is replaced with a saturation pulse: {θ₉₀ – TI – θ₉₀}. This technique was used to acquire the very first T₁ map Pykett & Mansfield, 1978. Unlike inversion recovery, this pulse sequence does not need a long TR to recover to its initial condition; every θ₉₀ pulse resets the longitudinal magnetization to the same initial state. However, to properly sample the recovery curve, TIs still need to reach the order of ~T₁, the dynamic range of signal potential is cut in half ([0, M₀]), and the short TIs (which have the fastest acquisition times) have the lowest SNRs.

8Variable Flip Angle T₁ Mapping¶

Variable flip angle (VFA) T₁ mapping Christensen et al., 1974Fram et al., 1987Gupta, 1977, also known as Driven Equilibrium Single Pulse Observation of T₁ (DESPOT1) Homer & Beevers, 1985Deoni et al., 2003, is a rapid quantitative T₁ measurement technique that is widely used to acquire 3D T₁ maps (e.g. whole-brain) in a clinically feasible time. VFA estimates T₁ values by acquiring multiple spoiled gradient echo acquisitions, each with different excitation flip angles (θ _n for n = 1, 2, .., N and θ _i ≠ θ _j). The steady-state signal of this pulse sequence (Figure 1) uses very short TRs (on the order of magnitude of 10 ms) and is very sensitive to T₁ for a wide range of flip angles.

VFA is a technique that originates from the NMR field, and was adopted because of its time efficiency and the ability to acquire accurate T₁ values simultaneously for a wide range of values Christensen et al., 1974Gupta, 1977. For imaging applications, VFA also benefits from an increase in SNR because it can be acquired using a 3D acquisition instead of multislice, which also helps to reduce slice profile effects. One important drawback of VFA for T₁ mapping is that the signal is very sensitive to inaccuracies in the flip angle value, thus impacting the T₁ estimates. In practice, the nominal flip angle (i.e. the value set at the scanner) is different than the actual flip angle experienced by the spins (e.g. at 3.0 T, variations of up to ±30%), an issue that increases with field strength. VFA typically requires the acquisition of another quantitative map, the transmit RF amplitude (B₁⁺, or B₁ for short), to calibrate the nominal flip angle to its actual value because of B₁ inhomogeneities that occur in most loaded MRI coils Sled & Pike, 1998. The need to acquire an additional B₁ map reduces the time savings offered by VFA over saturation-recovery techniques, and inaccuracies/imprecisions of the B₁ map are also propagated into the VFA T₁ map Boudreau et al., 2017Lee et al., 2017.

Figure 1. Simplified pulse sequence diagram of a variable flip angle (VFA) pulse sequence with a gradient echo readout. TR: repetition time, θ_n: excitation flip angle for the nth measurement, IMG: image acquisition (k-space readout), SPOIL: spoiler gradient.

9Signal Modelling¶

The steady-state longitudinal magnetization of an ideal variable flip angle experiment can be analytically solved from the Bloch equations for the spoiled gradient echo pulse sequence {θ _n–TR}:

\tag{1} M_{z}(\theta_n) = M_0 \frac{1-e^{- \frac{TR}{T_1}}}{1-\text{cos}(\theta_n) e^{- \frac{TR}{T_1}}} \text{sin}(\theta_n)

(15)

where M _z is the longitudinal magnetization,

₀ is the magnetization at thermal equilibrium, TR is the pulse sequence repetition time (Figure 1), and θ _n is the excitation flip angle. The M _z curves of different T₁ values for a range of θ _n and TR values are shown in Figure 2.

Source:

Signal Modelling

Figure 6:Variable flip angle technique signal curves (Eq. 1) for three different T₁ values, approximating the main types of tissue in the brain at 3T.

From Figure 2, it is clearly seen that the flip angle at which the steady-state signal is maximized is dependent on the T₁ and TR values. This flip angle is a well known quantity, called the Ernst angle Ernst & Anderson, 1966, which can be solved analytically from Equation 1 using properties of calculus:

\tag{2} \theta_{Ernst} = \text{acos}(e^{- \frac{TR}{T_1}})

(16)

The closed-form solution (Equation 1) makes several assumptions which in practice may not always hold true if care is not taken. Mainly, it is assumed that the longitudinal magnetization has reached a steady state after a large number of TRs, and that the transverse magnetization is perfectly spoiled at the end of each TR. Bloch simulations – a numerical approach at solving the Bloch equations for a set of spins at each time point – provide a more realistic estimate of the signal if the number of repetition times is small (i.e. a steady-state is not achieved). As can be seen from Figure 3, the number of repetitions required to reach a steady state not only depends on T₁, but also on the flip angle; flip angles near the Ernst angle need more TRs to reach a steady state. Preparation pulses or an outward-in k-space acquisition pattern are typically sufficient to reach a steady state by the time that the center of k-space is acquired, which is where most of the image contrast resides.

Source:

Signal Modelling

Figure 7:Signal curves simulated using Bloch simulations (orange) for a number of repetitions ranging from 1 to 150, plotted against the ideal case (Equation 1 – blue). Simulation details: TR = 25 ms, T₁ = 900 ms, 100 spins. Ideal spoiling was used for this set of Bloch simulations (transverse magnetization was set to 0 at the end of each TR).

Sufficient spoiling is likely the most challenging parameter to control for in a VFA experiment. A combination of both gradient spoiling and RF phase spoiling Bernstein et al., 2004Zur et al., 1991 are typically recommended (Figure 4). It has also been shown that the use of very strong gradients, introduces diffusion effects (not considered in Figure 4), further improving the spoiling efficacy in the VFA pulse sequence Yarnykh, 2010.

Source:

Signal Modelling

Figure 8:Signal curves estimated using Bloch simulations for three categories of signal spoiling: (1) ideal spoiling (blue), gradient & RF Spoiling (orange), and no spoiling (green). Simulations details: TR = 25 ms, T₁ = 900 ms, T_e = 100 ms, TE = 5 ms, 100 spins. For the ideal spoiling case, the transverse magnetization is set to zero at the end of each TR. For the gradient & RF spoiling case, each spin is rotated by different increments of phase (2𝜋 / # of spins) to simulate complete decoherence from gradient spoiling, and the RF phase of the excitation pulse is ɸ_n = ɸ_n-1 + nɸ₀ = ½ ɸ₀(n² + n + 2) Bernstein et al., 2004 with ɸ₀ = 117° Zur et al., 1991 after each TR.

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 1.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

TR_range = 5:5:200;

params.EXC_FA = 1:90;

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

for ii = 1:length(TR_range)
    params.TR = TR_range(ii);
    
    % White matter
    params.T1 = 900; % in milliseconds

    signal_WM(ii,:) = vfa_t1.analytical_solution(params);

    % Grey matter
    params.T1 = 1500;  % in milliseconds
    signal_GM(ii,:) = vfa_t1.analytical_solution(params);

    % CSF
    params.T1 = 4000;  % in milliseconds
    signal_CSF(ii,:) = vfa_t1.analytical_solution(params);
end

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 3.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

% White matter
params.T1 = 900; % in milliseconds
params.T2 = 10000;
params.TR = 25;
params.TE = 5;
params.EXC_FA = 1:90;
Nex_range = 1:1:150;

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

for ii = 1:length(Nex_range)
    params.Nex = Nex_range(ii);
    
    signal_analytical(ii,:) = vfa_t1.analytical_solution(params);

    [~, complex_signal] = vfa_t1.bloch_sim(params);
    signal_blochsim(ii,:) = abs(complex(complex_signal));
end

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 4.


% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

% White matter
params.T1 = 900; % in milliseconds
params.T2 = 100;
params.TR = 25;
params.TE = 5;
params.EXC_FA = 1:90;
Nex_range = [1:9, 10:10:100];

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

for ii = 1:length(Nex_range)
    params.Nex = Nex_range(ii);
    
    params.crushFlag = 1;
    
    [~, complex_signal] = vfa_t1.bloch_sim(params);
    signal_ideal_spoil(ii,:) = abs(complex_signal);
    
    
    params.inc = 117;
    params.partialDephasing = 1;
    params.partialDephasingFlag = 1;
    params.crushFlag = 0;
    
    [~, complex_signal] = vfa_t1.bloch_sim(params);
    signal_optimal_crush_and_rf_spoil(ii,:) = abs(complex_signal);
    
    params.inc = 0;
    params.partialDephasing = 0;

    [~, complex_signal] = vfa_t1.bloch_sim(params);
    signal_no_gradient_and_rf_spoil(ii,:) = abs(complex_signal);
end

References

10Data Fitting¶

At first glance, one could be tempted to fit VFA data using a non-linear least squares fitting algorithm such as Levenberg-Marquardt with Eq. 1, which typically only has two free fitting variables (T₁ and

₀). Although this is a valid way of estimating T₁ from VFA data, it is rarely done in practice because a simple refactoring of Equation 1 allows T₁ values to be estimated with a linear least square fitting algorithm, which substantially reduces the processing time. Without any approximations, Equation 1 can be rearranged into the form y = mx+b Gupta, 1977:

\tag{3} \frac{S_n}{ \text{sin}(\theta_n)} = e^{- \frac{TR}{T_1}} \frac{S_n}{ \text{tan}(\theta_n)} + C (1-e^{- \frac{TR}{T_1}})

(17)

As the third term does not change between measurements (it is constant for each θ _n), it can be grouped into the constant for a simpler representation:

\tag{4} \frac{S_n}{ \text{sin}(\theta_n)} = e^{- \frac{TR}{T_1}} \frac{S_n}{ \text{tan}(\theta_n)} + C

(18)

With this rearranged form of Equation 1, T₁ can be simply estimated from the slope of a linear regression calculated from S _n/sin(θ _n) and S _n/tan(θ _n) values:

\tag{5} T_1 = - \frac{TR}{ \text{ln}(slope)}

(19)

If data were acquired using only two flip angles – a very common VFA acquisition protocol – then the slope can be calculated using the elementary slope equation. Figure 5 displays both Equation 1 and 4 plotted for a noisy dataset.

Source:

Data Fitting

Figure 9:Mean and standard deviation of the VFA signal plotted using the nonlinear form (Equation 1 – blue) and linear form (Equation 4 – red). Monte Carlo simulation details: SNR = 25, N = 1000. VFA simulation details: TR = 25 ms, T₁ = 900 ms.

There are two important imaging protocol design considerations that should be taken into account when planning to use VFA: (1) how many and which flip angles to use to acquire VFA data, and (2) correcting inaccurate flip angles due to transmit RF field inhomogeneity. Most VFA experiments use the minimum number of required flip angles (two) to minimize acquisition time. For this case, it has been shown that the flip angle choice resulting in the best precision for VFA T₁ estimates for a sample with a single T₁ value (i.e. single tissue) are the two flip angles that result in 71% of the maximum possible steady-state signal (i.e. at the Ernst angle) Deoni et al., 2003Schabel & Morrell, 2008.

Time allowing, additional flip angles are often acquired at higher values and in between the two above, because greater signal differences between tissue T₁ values are present there (e.g. Figure 2). Also, for more than two flip angles, Equations 1 and 4 do not have the same noise weighting for each fitting point, which may bias linear least-square T₁ estimates at lower SNRs. Thus, it has been recommended that low SNR data should be fitted with either Equation 1 using non-linear least-squares (slower fitting) or with a weighted linear least-squares form of Equation 4 Chang et al., 2008.

Accurate knowledge of the flip angle values is very important to produce accurate T₁ maps. Because of how the RF field interacts with matter Sled & Pike, 1998, the excitation RF field (B₁⁺, or B₁ for short) of a loaded RF coil results in spatial variations in intensity/amplitude, unless RF shimming is available to counteract this effect (not common at clinical field strengths). For quantitative measurements like VFA which are sensitive to this parameter, the flip angle can be corrected (voxelwise) relative to the nominal value by multiplying it with a scaling factor (B₁) from a B₁ map that is acquired during the same session:

\tag{6} \theta_{corrected} = B_1 \theta_{nominal}

(20)

B₁ in this context is normalized, meaning that it is unitless and has a value of 1 in voxels where the RF field has the expected amplitude (i.e. where the nominal flip angle is the actual flip angle). Figure 6 displays fitted VFA T₁ values from a Monte Carlo dataset simulated using biased flip angle values, and fitted without/with B₁ correction.

Source:

Data Fitting

Figure 10:Mean and standard deviations of fitted VFA T₁ values for a set of Monte Carlo simulations (SNR = 100, N = 1000), simulated using a wide range of biased flip angles and fitted without (blue) or with (red) B₁ correction. Simulation parameters: TR = 25 ms, T₁ = 900 ms, θ_nominal = 6° and 32° (optimized values for this TR/T₁ combination). Notice how even after B₁ correction, fitted T₁ values at B₁ values far from the nominal case (B₁ = 1) exhibit larger variance, as the actual flip angles of the simulated signal deviate from the optimal values for this TR/T₁ (Deoni et al. 2003).

Figure 7 displays an example VFA dataset and a B₁ map in a healthy brain, along with the T₁ map estimated using a linear fit (Equations 4 and 5).

Source:

Data Fitting

Figure 11:Example variable flip angle dataset and B₁ map of a healthy adult brain (left). The relevant VFA protocol parameters used were: TR = 15 ms, θ_nominal = 3° and 20°. The T₁ map (right) was fitted using a linear regression (Equations 4 and 5).

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 5.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

params.EXC_FA = [1:4,5:5:90];

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

params.TR = 0.025;
params.EXC_FA = [2:9,10:5:90];

% White matter
x.M0 = 1;
x.T1 = 0.900; % in milliseconds

Model = vfa_t1; 

Opt.SNR = 25;
Opt.TR = params.TR;
Opt.T1 = x.T1;

clear Model.Prot.VFAData.Mat(:,1) 
Model.Prot.VFAData.Mat = zeros(length(params.EXC_FA),2);
Model.Prot.VFAData.Mat(:,1) = params.EXC_FA';
Model.Prot.VFAData.Mat(:,2) = Opt.TR;

for jj = 1:1000
    [FitResult{jj}, noisyData{jj}] = Model.Sim_Single_Voxel_Curve(x,Opt,0); 
    fittedT1(jj) = FitResult{jj}.T1;
    noisyData_array(jj,:) = noisyData{jj}.VFAData;
    noisyData_array_div_sin(jj,:) = noisyData_array(jj,:) ./ sind(Model.Prot.VFAData.Mat(:,1))';
    noisyData_array_div_tan(jj,:) = noisyData_array(jj,:) ./ tand(Model.Prot.VFAData.Mat(:,1))';
end
        
for kk=1:length(params.EXC_FA)
    data_mean(kk) = mean(noisyData_array(:,kk));
    data_std(kk) = std(noisyData_array(:,kk));
    
    data_mean_div_sin(kk) = mean(noisyData_array_div_sin(:,kk));
    data_std_div_sin(kk) = std(noisyData_array_div_sin(:,kk));
    
    data_mean_div_tan(kk) = mean(noisyData_array_div_tan(:,kk));
    data_std_div_tan(kk) = std(noisyData_array_div_tan(:,kk));
end

%% Setup parameters
% All times are in milliseconds
% All flip angles are in degrees

params_highres.EXC_FA = [2:1:90];

%% Calculate signals
%
% To see all the options available, run `help vfa_t1.analytical_solution`

params_highres.TR = params.TR * 1000; % in milliseconds
    
% White matter
params_highres.T1 = x.T1*1000; % in milliseconds

signal_WM = vfa_t1.analytical_solution(params_highres);
signal_WM_div_sin = signal_WM ./ sind(params_highres.EXC_FA);
signal_WM_div_tan = signal_WM ./ tand(params_highres.EXC_FA);

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 6.

% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% Setup parameters
% All times are in seconds
% All flip angles are in degrees

params.TR = 0.025; % in seconds

% White matter
params.T1 = 0.900; % in seconds

% Calculate optimal flip angles for a two flip angle VFA experiment (for this T1 and TR)
% The will be the nominal flip angles (the flip angles assumed by the "user", before a 
% "realistic"B1 bias is applied)

nominal_EXC_FA = vfa_t1.find_two_optimal_flip_angles(params); % in degrees
disp('Nominal flip angles:')
disp(nominal_EXC_FA)

% Range of B1 values biasing the excitation flip angle away from their nominal values
B1Range = 0.1:0.1:2;

x.M0 = 1;
x.T1 = params.T1; % in seconds

Model = vfa_t1; 

Opt.SNR = 100;
Opt.TR = params.TR;
Opt.T1 = x.T1;

% Monte Carlo signal simulations
for ii = 1:1000
    for jj = 1:length(B1Range)
        B1 = B1Range(jj);
        actual_EXC_FA = B1 * nominal_EXC_FA;
 
        params.EXC_FA = actual_EXC_FA;

        clear Model.Prot.VFAData.Mat(:,1)
        Model.Prot.VFAData.Mat = zeros(length(params.EXC_FA),2);
        Model.Prot.VFAData.Mat(:,1) = params.EXC_FA';
        Model.Prot.VFAData.Mat(:,2) = Opt.TR;

        [FitResult{ii,jj}, noisyData{ii,jj}] = Model.Sim_Single_Voxel_Curve(x,Opt,0); 
        noisyData_array(ii,jj,:) = noisyData{ii,jj}.VFAData;
    end
end
%
Model = vfa_t1; 
    
FlipAngle = nominal_EXC_FA';
TR = params.TR .* ones(size(FlipAngle));

Model.Prot.VFAData.Mat = [FlipAngle TR];

data.VFAData(:,:,1,1) = noisyData_array(:,:,1);
data.VFAData(:,:,1,2) = noisyData_array(:,:,2);
data.Mask = repmat(ones(size(B1Range)),[size(noisyData_array,1),1]);

data.B1map = repmat(ones(size(B1Range)),[size(noisyData_array,1),1]);
FitResults_noB1Correction = FitData(data,Model,0);

data.B1map = repmat(B1Range,[size(noisyData_array,1),1]);
FitResults_withB1Correction = FitData(data,Model,0);

%%
%

mean_T1_noB1Correction = mean(FitResults_noB1Correction.T1);
mean_T1_withB1Correction = mean(FitResults_withB1Correction.T1);
std_T1_noB1Correction = std(FitResults_noB1Correction.T1);
std_T1_withB1Correction = std(FitResults_withB1Correction.T1);

Click here to view the qMRLab (MATLAB/Octave) code that generated Figure 7.


% Verbosity level 0 overrides the disp function and supresses warnings.
% Once executed, they cannot be restored in this session
% (kernel needs to be restarted or a new notebook opened.)
VERBOSITY_LEVEL = 0;

if VERBOSITY_LEVEL==0
    % This hack was used to supress outputs from external tools
    % in the Jupyter Book.
    function disp(x)
    end
    warning('off','all')
end

try
    cd qMRLab
catch
    try
        cd ../../../qMRLab
    catch
        cd ../qMRLab
    end
end

startup
clear all

%% MATLAB/OCTAVE CODE

% Load data into environment, and rotate mask to be aligned with IR data
load('data/vfa_dataset/VFAData.mat');
load('data/vfa_dataset/B1map.mat');
load('data/vfa_dataset/Mask.mat');

% Format qMRLab vfa_t1 model parameters, and load them into the Model object
Model = vfa_t1; 
FlipAngle = [    3;     20];
TR        = [0.015; 0.0150];

Model.Prot.VFAData.Mat = [FlipAngle, TR];

% Format data structure so that they may be fit by the model
data = struct();
data.VFAData= double(VFAData);
data.B1map= double(B1map);
data.Mask= double(Mask);

FitResults = FitData(data,Model,0); % The '0' flag is so that no wait bar is shown.

T1_map = imrotate(FitResults.T1.*Mask,-90);
T1_map(T1_map>5)=0;
T1_map = T1_map*1000; % Convert to ms

xAxis = [0:size(T1_map,2)-1];
yAxis = [0:size(T1_map,1)-1];

% Raw MRI data at different TI values
FA_03 = imrotate(squeeze(VFAData(:,:,:,1).*Mask),-90);
FA_20 = imrotate(squeeze(VFAData(:,:,:,2).*Mask),-90);
B1map = imrotate(squeeze(B1map.*Mask),-90);

References

11Benefits and Pitfalls¶

It has been well reported in recent years that the accuracy of VFA T₁ estimates is very sensitive to pulse sequence implementations Baudrexel et al., 2017Lutti & Weiskopf, 2013Stikov et al., 2015, and as such is less robust than the gold standard inversion recovery technique. In particular, the signal bias resulting from insufficient spoiling can result in inaccurate T₁ estimates of up to 30% relative to inversion recovery estimated values Stikov et al., 2015. VFA T₁ map accuracy and precision is also strongly dependent on the quality of the measured B₁ map Lee et al., 2017, which can vary substantially between implementations Boudreau et al., 2017. Modern rapid B₁ mapping pulse sequences are not as widely available as VFA, resulting in some groups attempting alternative ways of removing the bias from the T₁ maps like generating an artificial B₁ map through the use of image processing techniques Liberman et al., 2013 or omitting B₁ correction altogether Yuan et al., 2012. The latter is not recommended, because most MRI scanners have default pulse sequences that, with careful protocol settings, can provide B₁ maps of sufficient quality very rapidly Boudreau et al., 2017Samson et al., 2006Wang et al., 2005.

Despite some drawbacks, VFA is still one of the most widely used T₁ mapping methods in research. Its rapid acquisition time, rapid image processing time, and widespread availability makes it a great candidate for use in other quantitative imaging acquisition protocols like quantitative magnetization transfer imaging Cercignani et al., 2005Yarnykh, 2002 and dynamic contrast enhanced imaging Li et al., 2018Sung et al., 2013.

Acknowledgments¶

This work was supported by the Canadian Open Neuroscience Platform (CONP), the Quebec Bio-Imaging Network (QBIN), and the Montreal Heart Institute Foundation (MHI).

References¶

Boudreau, M., Keenan, K. E., & Stikov, N. (2020). Chapter 2 - Quantitative T1 and T1ρ Mapping. In N. Seiberlich, V. Gulani, F. Calamante, A. Campbell-Washburn, M. Doneva, H. H. Hu, & S. Sourbron (Eds.), Quantitative Magnetic Resonance Imaging (Vol. 1, pp. 19–45). Academic Press. 10.1016/B978-0-12-817057-1.00004-4
Boudreau, M., Keenan, K. E., & Stikov, N. (2023). Quantitative T1 MRI. 10.55458/neurolibre.00019
Karakuzu, A., DuPre, E., Tetrel, L., Bermudez, P., Boudreau, M., Chin, M., Poline, J.-B., Das, S., Bellec, P., & Stikov, N. (2022). NeuroLibre : A preprint server for full-fledged reproducible neuroscience. OSF Preprints. 10.31219/osf.io/h89js
DuPre, E., Holdgraf, C., Karakuzu, A., Tetrel, L., Bellec, P., Stikov, N., & Poline, J.-B. (2022). Beyond advertising: New infrastructures for publishing integrated research objects. PLOS Computational Biology, 18(1), e1009651. 10.1371/journal.pcbi.1009651
Harding, R. J., Bermudez, P., Bernier, A., Beauvais, M., Bellec, P., Hill, S., Karakuzu, A., Knoppers, B. M., Pavlidis, P., Poline, J.-B., Roskams, J., Stikov, N., Stone, J., Strother, S., Consortium, C., & Evans, A. C. (2023). The Canadian Open Neuroscience Platform—An open science framework for the neuroscience community. PLOS Computational Biology, 19(7), 1–14. 10.1371/journal.pcbi.1011230
Karakuzu, A., Boudreau, M., Duval, T., Boshkovski, T., Leppert, I., Cabana, J.-F. ̧cois, Gagnon, I., Beliveau, P., Pike, G., Cohen-Adad, J., & Stikov, N. (2020). qMRLab: Quantitative MRI analysis, under one umbrella. J. Open Source Softw., 5(53), 2343. 10.21105/joss.02343
Hahn, E. L. (1949). An accurate nuclear magnetic resonance method for measuring spin-lattice relaxation times. Physical Review, 76(1), 145. 10.1103/PhysRev.76.145
Drain, L. (1949). A direct method of measuring nuclear spin-lattice relaxation times. Proceedings of the Physical Society. Section A, 62(5), 301. 10.1088/0370-1298/62/5/306
Pykett, I., & Mansfield, P. (1978). A line scan image study of a tumorous rat leg by NMR. Physics in Medicine & Biology, 23(5), 961. 10.1088/0031-9155/23/5/012
Stikov, N., Boudreau, M., Levesque, I. R., Tardif, C. L., Barral, J. K., & Pike, G. B. (2015). On the accuracy of T1 mapping: Searching for common ground. Magnetic Resonance in Medicine, 73(2), 514–522. 10.1002/mrm.25135
Steen, R. G., Gronemeyer, S. A., Kingsley, P. B., Reddick, W. E., Langston, J. S., & Taylor, J. S. (1994). Precise and accurate measurement of proton T1 in human brain in vivo: validation and preliminary clinical application. Journal of Magnetic Resonance Imaging, 4(5), 681–691. 10.1002/jmri.1880040511
Messroghli, D. R., Radjenovic, A., Kozerke, S., Higgins, D. M., Sivananthan, M. U., & Ridgway, J. P. (2004). Modified Look-Locker inversion recovery (MOLLI) for high-resolution T1 mapping of the heart. Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine, 52(1), 141–146. 10.1002/mrm.20110
Piechnik, S. K., Ferreira, V. M., Dall’Armellina, E., Cochlin, L. E., Greiser, A., Neubauer, S., & Robson, M. D. (2010). Shortened Modified Look-Locker Inversion recovery (ShMOLLI) for clinical myocardial T1-mapping at 1.5 and 3 T within a 9 heartbeat breathhold. Journal of Cardiovascular Magnetic Resonance, 12(1), 1–11. 10.1186/1532-429X-12-69
Christensen, K. A., Grant, D. M., Schulman, E. M., & Walling, C. (1974). Optimal determination of relaxation times of Fourier transform nuclear magnetic resonance. Determination of spin-lattice relaxation times in chemically polarized species. The Journal of Physical Chemistry, 78(19), 1971–1977. 10.1021/j100612a022
Fram, E. K., Herfkens, R. J., Johnson, G. A., Glover, G. H., Karis, J. P., Shimakawa, A., Perkins, T. G., & Pelc, N. J. (1987). Rapid calculation of T1 using variable flip angle gradient refocused imaging. Magnetic Resonance Imaging, 5(3), 201–208. 10.1016/0730-725x(87)90021-x

Abstract¶

0.1Other Saturation-Recovery T1 Mapping techniques¶

1Variable Flip Angle T1 Mapping¶

2About¶

3Inversion Recovery T1 Mapping¶

4Signal Modelling¶

5Data Fitting¶

8Variable Flip Angle T1 Mapping¶

9Signal Modelling¶

Acknowledgments¶

3Inversion Recovery T₁ Mapping¶

8Variable Flip Angle T₁ Mapping¶