Display content:

In [1]:

```
% **Blog post code introduction**
%
% Congrats on activating the "All cells" option in this interactive blog post =D
%
% Below, several new HTML blocks have appears prior to the figures, displaying the Octave/MATLAB code that was used to generate the figures in this blog post.
%
% If you want to reproduce the data on your own local computer, you simply need to have qMRLab installed in your Octave/MATLAB path and run the "startup.m" file, as is shown below.
%
% If you want to get under the hood and modify the code right now, you can do so in the Jupyter Notebook of this blog post hosted on MyBinder. The link to it is in the introduction above.
```

In [2]:

Dictionary-based MRI techniques capable of generating T_{1} maps are increasing in popularity, due to their growing availability on clinical scanners, rapid scan times, and fast post-processing computation time, thus making quantitative T_{1} mapping accessible for clinical applications. Generally speaking, dictionary-based quantitative MRI techniques use numerical dictionaries—databases of pre-calculated signal values simulated for a wide range of tissue and protocol combinations—during the image reconstruction or post-processing stages. Popular examples of dictionary-based techniques that have been applied to T_{1} mapping are MR Fingerprinting (MRF) (Ma et al. 2013), certain flavours of compressed sensing (CS) (Doneva et al. 2010; Li et al. 2012), and Magnetization Prepared 2 Rapid Acquisition Gradient Echoes (MP2RAGE) (Marques et al. 2010). Dictionary-based techniques can usually be classified into one of two categories: techniques that use information redundancy from parametric data to assist in accelerated imaging (e.g. CS, MRF), or those that use dictionaries to estimate quantitative maps using the MR images after reconstruction. Because MP2RAGE is a technique implemented primarily for T_{1} mapping, and it is becoming increasingly available as a standard pulse sequence on many MRI systems, the remainder of this section will focus solely on this technique. However, many concepts discussed are shared by other dictionary-based techniques.

MP2RAGE is an extension of the conventional MPRAGE pulse sequence widely used in clinical studies (Haase et al. 1989; Mugler & Brookeman 1990). A simplified version of the MP2RAGE pulse sequence is shown in Figure 1. MP2RAGE can be seen as a hybrid between the inversion recovery and VFA pulse sequences: a 180° inversion pulse is used to prepare the magnetization for T_{1} sensitivity at the beginning of each TRMP2RAGE, and then two images are acquired at different inversion times using gradient recalled echo (GRE) imaging blocks with low flip angles and short repetition times (TR). During a given GRE imaging block, each excitation pulse is followed by a constant in-plane (“y”) phase encode weighting (varied for each TRMP2RAGE), but with different 3D (“z”) phase encoding gradients (varied at each TR). The center of k-space for the 3D phase encoding direction is acquired at the TI for each GRE imaging block. The main motivation for developing the MP2RAGE pulse sequence was to provide a metric similar to MPRAGE, but with self-bias correction of the static (B_{0}) and receive (B_{1}^{-}) magnetic fields, and a first order correction of the transmit magnetic field (B_{1}^{+}). However, because two images at different TIs are acquired (unlike MPRAGE, which only acquires data at a single TI), information about the T_{1} values can also be inferred, thus making it possible to generate quantitative T_{1} maps using this data.

Prior to considering the full signal equations, we will first introduce the equation for the MP2RAGE parameter (*S*_{MP2RAGE}) that is calculated in addition to the T_{1} map. For complex data (magnitude and phase, or real and imaginary), the MP2RAGE signal (*S*_{MP2RAGE}) is calculated from the images acquired at two TIs (*S*_{GRE,TI1} and *S*_{GRE,TI2}) using the following expression (Marques et al. 2010):

This value is bounded between [-0.5, 0.5], and helps reduce some B_{0} inhomogeneity effects using the phase data. For real data, or magnitude data with polarity restoration, this metric is instead calculated as:

Because MP2RAGE is a hybrid of pulse sequences used for inversion recovery and VFA, the resulting signal equations are more complex. Typically, a steady state is not achieved during the short train of GRE imaging blocks, so the signal at the center of k-space for each readout (which defines the contrast weighting) will depend on the number of phase-encoding steps. For simplicity, the equations presented here assume that the 3D phase-encoding dimension is fully sampled (no partial Fourier or parallel imaging acceleration). For this case (see appendix of (Marques et al. 2010) for derivation details), the signal equations are:

where B_{1}^{-} is the receive field sensitivity, “eff” is the adiabatic inversion pulse efficiency, ER = exp(-TR/T_{1}), EA = exp(-TA/T_{1}), EB = exp(-TB/T_{1}), EC = exp(-TC/T_{1}). The variables TA, TB, and TC are the three different delay times (TA: time between inversion pulse and beginning of the GRE_{1} block, TB: time between the end of GRE_{1} and beginning of GRE_{2}, TC: time between the end of GRE_{2} and the end of the TR). If no k-space acceleration is used (e.g. no partial Fourier or parallel imaging acceleration), then these values are TA = TI_{1} - (n/2)TR, TB = TI_{2} - TI_{1} + (n/2)TR, and TC = TR_{MP2RAGE} - TI_{2} + (n/2)TR, where n is the number of voxels acquired in the 3D phase encode direction varied within each GRE block. The value m_{z,ss} is the steady-state longitudinal magnetization prior to the inversion pulse, and is given by:

From Eqs. (3–6), it is evident that the MP2RAGE parameter *S*_{MP2RAGE} (Eqs. (1.11) and (1.12)) cancels out the effects of receive field sensitivity, T_{2}^{*}, and M_{0}. The signal sensitivity related to the transmit field (B_{1}^{+}), hidden in Eqs. (3–6) within the flip angle values *θ _{1}* and

_{1}), and then interpolation is done within this dictionary of values to estimate the T_{1} value that matches the observed signal. This approach results in rapid post-processing times because the dictionaries can be simulated/generated prior to scanning and interpolating between these values is much faster than most fitting algorithms. This means that the quantitative image can be produced and displayed directly on the MRI scanner console rather than needing to be fitted offline.

In [3]:

```
%% MATLAB/OCTAVE CODE
% Code used to generate the data required for Figure 5 of the blog post
clear all
cd ../qMRLab
startup
MP2RAGE.B0=7; % in Tesla
MP2RAGE.TR=6; % MP2RAGE TR in seconds
MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout
MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding
MP2RAGE.NZslices=[35 72];% Slices Per Slab * [PartialFourierInSlice-0.5 0.5]
MP2RAGE.FlipDegrees=[4 5];% Flip angle of the two readouts in degrees
invEFF=0.99;
B1_vector=0.005:0.05:1.9;
T1_vector=0.5:0.05:5.2;
npoints=40;
%% creates a lookup table of MP2RAGE intensities as a function of B1 and T1
k=0;
for b1val=B1_vector
k=k+1;
[Intensity T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,b1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF);
MP2RAGEmatrix(k,:)=interp1(T1vector,Intensity,T1_vector);
end;
%% make the matrix MP2RAGEMatrix into T1_matrix(B1,ratio)
MP2RAGE_vector=linspace(-0.5,0.5,npoints);
k=0;
for b1val=B1_vector
k=k+1;
try
T1matrix(k,:)=interp1(MP2RAGEmatrix(k,:),T1_vector,MP2RAGE_vector,'pchirp');
catch
temp=MP2RAGEmatrix(k,:);
temp(isnan(temp))=linspace(-0.5,-1,length(find(isnan(temp)==1)));
temp=interp1(temp,T1_vector,MP2RAGE_vector);
T1matrix(k,:)=temp;
end
end
```

In [4]:

In [5]:

_{1} maps with good accuracy and precision using dictionary-based interpolation methods, it is important that the signal curves are unique for each parameter value. MP2RAGE can produce good T_{1} maps by using a dictionary with only dimensions (T_{1}, *S*_{MP2RAGE}), since *S*_{MP2RAGE} is unique for each T_{1} value for a given protocol (Marques et al. 2010). However, as was noted above, *S*_{MP2RAGE} is also sensitive to B_{1} because of *θ _{1}* and

In [6]:

```
%% MATLAB/OCTAVE CODE
% Download variable flip angle brain MRI data for Figure 7 of the blog post
cmd = ['curl -L -o mp2rage.zip https://osf.io/8x2c9/download/'];
[STATUS,MESSAGE] = unix(cmd);
unzip('mp2rage.zip');
```

In [7]:

```
%% MATLAB/OCTAVE CODE
% Code used to generate the data required for Figure 5 of the blog post
clear all
cd ../qMRLab
startup
MP2RAGE.B0=7; % in Tesla
MP2RAGE.TR=6; % MP2RAGE TR in seconds
MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout
MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding
MP2RAGE.NZslices=[35 72];% Slices Per Slab * [PartialFourierInSlice-0.5 0.5]
MP2RAGE.FlipDegrees=[4 5];% Flip angle of the two readouts in degrees
MP2RAGE.filenameUNI='MP2RAGE_UNI.nii'; % file with UNI
MP2RAGE.filenameINV1='MP2RAGE_INV1.nii'; % file with UNI
MP2RAGE.filenameINV2='MP2RAGE_INV2.nii'; % file with INV2
% load the MP2RAGE data - it can be either the SIEMENS one scaled from
% 0 4095 or the standard -0.5 to 0.5
MP2RAGEimg=load_untouch_nii(MP2RAGE.filenameUNI, [], [], [], [], [], []);
MP2RAGEINV1img=load_untouch_nii(MP2RAGE.filenameINV1, [], [], [], [], [], []);
MP2RAGEINV2img=load_untouch_nii(MP2RAGE.filenameINV2, [], [], [], [], [], []);
[T1map , M0map , R1map]=T1M0estimateMP2RAGE(MP2RAGEimg,MP2RAGEINV2img,MP2RAGE,0.96);
% Code used to re-orient the images to make pretty figures, and to assign variables with the axis lengths.
T1_map = imrotate(squeeze(T1map.img(:,130,:)),180);
T1_map(T1map.img>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
S_INV1 = imrotate(squeeze(MP2RAGEINV1img.img(:,130,:)),180);
S_INV2 = imrotate(squeeze(MP2RAGEINV2img.img(:,130,:)),180);
B1map = imrotate(-0.5+1/4095*double(squeeze(MP2RAGEimg.img(:,130,:))),180);
```

In [8]:

In [9]: