Benefits 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_{0}$, $kM_{0}$]), 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 2.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 T₁ mapping Boudreau et al., 2023.

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. Equations Eq. 2.1 and Eq. 2.4), and up to 5 T₁ for long TR approximated models (Equation 2.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.

% 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