qmt_spgr: quantitative Magnetizatoion Transfer (qMT) using Spoiled Gradient Echo

Contents

• 1. Print qmt_spgr information
• 2. Setting model parameters
• 2.a. Create qmt_spgr object
• 2.b. Set protocol and options
• 2.b.1 Set protocol the CLI way
• 2.b.2 Set protocol and options the GUI way
• 3. Fit MRI data
• 3.a. Load input data
• 3.b. Execute fitting process
• 3.c. Display FitResults
• 3.d. Save fit results
• 3.e. Re-use or share fit configuration files
• 4. Simulations
• 4.a. Single Voxel Curve
• 4.b. Sensitivity Analysis
• 5. Notes
• 5.a. Notes specific to qmt_spgr
• 5.a.1 BIDS-like convention
• 5.b. Generic notes
• 5.b.1. Batch friendly option and protocol conventions
• 5.b.2 Parallelization
• 6. Citations

% This m-file has been automatically generated using qMRgenBatch(qmt_spgr)
% for publishing documentation.
% Command Line Interface (CLI) is well-suited for automatization
% purposes and Octave.
%
% Please execute this m-file section by section to get familiar with batch
% processing for qmt_spgr on CLI.
%
% Demo files are downloaded into qmt_spgr_data folder.
%
% Written by: Agah Karakuzu, 2017
% ==============================================================================

1. Print qmt_spgr information

qMRinfo('qmt_spgr');

 qmt_spgr:  quantitative Magnetizatoion Transfer (qMT) using Spoiled Gradient Echo
Pulse Sequence Diagram

Assumptions:
FILL

Inputs:
MTdata              Magnetization Transfer data
R1map               Longitudinal relaxation rate (VFA RECOMMENDED Boudreau 2017 MRM).
(B1map)             Normalized transmit excitation field map (B1+). B1+ is defined
as a  normalized multiplicative factor such that:
FA_actual = B1+ * FA_nominal. (OPTIONAL). (=1 if not provided)
(B0map)             B0 field map, used for offset correction  (OPTIONAL) (=0Hz if not provided)
(Mask)              Binary mask to accelerate the fitting (OPTIONAL)

Outputs:
F                   Ratio of number of restricted pool to free pool, defined
as F = M0r/M0f = kf/kr.
kr                  Exchange rate from the free to the restricted pool
(note that kf and kr are related to one another via the
definition of F. Changing the value of kf will change kr
accordingly, and vice versa).
R1f                 Longitudinal relaxation rate of the free pool
(R1f = 1/T1f).
R1r                 Longitudinal relaxation rate of the restricted pool
(R1r = 1/T1r).
T2f                 Tranverse relaxation time of the free pool (T2f = 1/R2f).
T2r                 Tranverse relaxation time of the restricted pool (T2r = 1/R2r).
(kf)                Exchange rate from the restricted to the free pool.
(resnorm)           Fitting residual.

Protocol:
MTdata              Array [Nb of volumes x 2]
Angle             MT pulses angles (degree)
Offset            Offset frequencies (Hz)

TimingTable         Vector [5x1]
Tmt               Duration of the MT pulses (s)
Ts                Free precession delay between the MT and excitation pulses (s)
Tp                Duration of the excitation pulse (s)
Tr                Free precession delay after the excitation pulse, before
the next MT pulse (s)
TR                Repetition time of the whole sequence (TR = Tmt + Ts + Tp + Tr)


Options:
MT Pulse
Shape                 Shape of the MT pulse.
Available shapes are:
- hard
- gaussian
- gausshann (gaussian pulse with Hanning window)
- sinc
- sinchann (sinc pulse with Hanning window)
- singauss (sinc pulse with gaussian window)
- fermi
Sinc TBW              Time-bandwidth product for the sinc MT pulses
(applicable to sinc, sincgauss, sinchann MT
pulses).
Bandwidth             Bandwidth of the gaussian MT pulse (applicable
to gaussian, gausshann and sincgauss MT pulses).
Fermi transition (a)  slope 'a' (related to the transition width)
of the Fermi pulse (applicable to fermi MT
pulse).
Assuming pulse duration at 60 dB (from the Bernstein handbook)
and t0 = 10a,
slope = Tmt/33.81;
# of MT pulses        Number of pulses used to achieve steady-state
before a readout is made.
Fitting constraints
Use R1map to         By checking this box, you tell the fitting
constrain R1f          algorithm to check for an observed R1map and use
its value to constrain R1f. Checking this box
will automatically set the R1f fix box to true
in the Fit parameters table.
Fix R1r = R1f        By checking this box, you tell the fitting
algorithm to fix R1r equal to R1f. Checking this
box will automatically set the R1r fix box to
true in the Fit parameters table.
Fix R1f*T2f          By checking this box, you tell the fitting
algorithm to compute T2f from R1f value. R1f*T2f
value is set in the next box.
R1f*T2f =            Value of R1f*T2f (no units)

Model                  Model you want to use for fitting.
Available models are:
- SledPikeRP (Sled & Pike rectangular pulse),
- SledPikeCW (Sled & Pike continuous wave),
- Yarkykh (Yarnykh & Yuan)
- Ramani
Note: Sled & Pike models will show different
options than Yarnykh or Ramani.
Lineshape              The absorption lineshape of the restricted pool.
Available lineshapes are:
- Gaussian
- Lorentzian
- SuperLorentzian
Read pulse alpha       Flip angle of the excitation pulse.
Compute SfTable        By checking this box, you compute a new SfTable

Command line usage:
qMRusage(qmt_spgr

Author: Ian Gagnon, 2017

References:
Please cite the following if you use this module:
Sled, J.G., Pike, G.B., 2000. Quantitative interpretation of magnetization transfer in spoiled gradient echo MRI sequences. J. Magn. Reson. 145, 24?36.
In addition to citing the package:
Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F.,
Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab:
Quantitative MRI analysis, under one umbrella doi: 10.21105/joss.02343

Reference page in Doc Center
doc qmt_spgr

2. Setting model parameters

2.a. Create qmt_spgr object

Model = qmt_spgr;

2.b. Set protocol and options

Protocol: MRI acquisition parameters that are accounted for by the respective model.

For example: TE, TR, FA FieldStrength. The assigned protocol values are subjected to a sanity check to ensure that they are in agreement with the data attributes.

Options: Fitting preferences that are left at user's discretion.

For example: linear fit, exponential fit, drop first echo.

2.b.1 Set protocol the CLI way

If you are using Octave, or would like to serialize your operations any without GUI involvement, you can assign protocol directly in CLI:

Angle = [142.0000; 426.0000; 142.0000; 426.0000; 142.0000; 426.0000; 142.0000; 426.0000; 142.0000; 426.0000];
% Angle is a vector of [10X1]
Offset = [443.0000; 443.0000; 1088.0000; 1088.0000; 2732.0000; 2732.0000; 6862.0000; 6862.0000; 17235.0000; 17235.0000];
% Offset is a vector of [10X1]
Model.Prot.MTdata.Mat = [ Angle Offset];

Tmt  = 0.0102;
Ts  = 0.003;
Tp  = 0.0018;
Tr  = 0.01;
TR  = 0.025;
Model.Prot.TimingTable.Mat = [ Tmt  Ts  Tp  Tr  TR ];

  See the generic notes section below for further information.

2.b.2 Set protocol and options the GUI way

The following command opens a panel to set protocol and options (if GUI is available to the user):

Model = Custom_OptionsGUI(Model);

If available, you need to close this panel for the remaining of the script to proceed.

  Using this panel, you can save qMRLab protocol files that can be used in both interfaces. See the generic notes section below for details.

3. Fit MRI data

3.a. Load input data

This section shows how you can load data into a(n) qmt_spgr object.

• At the CLI level, qMRLab accepts structs containing (double) data in the fields named in accordance with a qMRLab model.

  See the generic notes section below for BIDS compatible wrappers and scalable
      qMRLab workflows.

%          |- qmt_spgr object needs 5 data input(s) to be assigned:
%          |-   MTdata
%          |-   R1map
%          |-   B1map
%          |-   B0map
%          |-   Mask

data = struct();

% MTdata.mat contains [88  128    1   10] data.
load('qmt_spgr_data/MTdata.mat');
% R1map.mat contains [88  128] data.
load('qmt_spgr_data/R1map.mat');
% B1map.mat contains [88  128] data.
load('qmt_spgr_data/B1map.mat');
% B0map.mat contains [88  128] data.
load('qmt_spgr_data/B0map.mat');
% Mask.mat contains [88  128] data.
load('qmt_spgr_data/Mask.mat');
data.MTdata= double(MTdata);
data.R1map= double(R1map);
data.B1map= double(B1map);
data.B0map= double(B0map);
data.Mask= double(Mask);

3.b. Execute fitting process

This section will fit the loaded data.

FitResults = FitData(data,Model,0);

Visit the generic notes section below for instructions to accelerate fitting by
      parallelization using ParFitData.

3.c. Display FitResults

You can display the current outputs by:

qMRshowOutput(FitResults,data,Model);

A representative fit curve will be plotted if available.

To render images in this page, we will load the fit results that had been saved before. You can skip the following code block;

% Load FitResults that comes with the example dataset.
FitResults_old = load('FitResults/FitResults.mat');
qMRshowOutput(FitResults_old,data,Model);

3.d. Save fit results

Outputs can be saved as *.nii.(gz) if NIfTI inputs are available:

% Generic function call to save nifti outputs
FitResultsSave_nii(FitResults, 'reference/nifti/file.nii.(gz)');

If not, FitResults.mat file can be saved. This file contains all the outputs as workspace variables:

% Generic function call to save FitResults.mat
FitResultsSave_mat(FitResults);

  FitResults.mat files can be loaded to qMRLab GUI for visualization and ROI
       analyses.

The section below will be dynamically generated in accordance with the example data format (mat or nii). You can substitute FitResults_old with FitResults if you executed the fitting using example dataset for this model in section 3.b..

FitResultsSave_mat(FitResults_old);

3.e. Re-use or share fit configuration files

qMRLab's fit configuration files (qmt_spgr_Demo.qmrlab.mat) store all the options and protocol in relation to the used model and the release version.

  *.qmrlab.mat files can be easily shared with collaborators to allow them fit their own
      data or run simulations using identical option and protocol configurations.

Model.saveObj('my_qmt_spgr_config.qmrlab.mat');

4. Simulations

4.a. Single Voxel Curve

Simulates single voxel curves

1. Analytically generate synthetic MRI data
2. Add rician noise
3. Fit and plot the respective curve

      x = struct;
x.F = 0.16;
x.kr = 30;
x.R1f = 1;
x.R1r = 1;
x.T2f = 0.03;
x.T2r = 1.3e-05;
% Set simulation options
Opt.SNR = 50;
Opt.Method = 'Analytical equation';
Opt.ResetMz = false;
% run simulation
figure('Name','Single Voxel Curve Simulation');
FitResult = Model.Sim_Single_Voxel_Curve(x,Opt);

4.b. Sensitivity Analysis

Simulates sensitivity to fitted parameters

1. Iterate fitting parameters from lower (lb) to upper (ub) bound
2. Run Sim_Single_Voxel_Curve for Nofruns times
3. Compute the mean and std across runs

            F             kr            R1f           R1r           T2f           T2r

      OptTable.st = [0.16          30            1             1             0.03          1.3e-05]; % nominal values
OptTable.fx = [0             1             1             1             1             1]; %vary F...
OptTable.lb = [0.0001        0.0001        0.05          0.05          0.003         3e-06]; %...from 0.0001
OptTable.ub = [0.5           1e+02         5             5             0.5           5e-05]; %...to 0.5
% Set simulation options
Opt.SNR = 50;
Opt.Method = 'Analytical equation';
Opt.ResetMz = false;
Opt.Nofrun = 5;
% run simulation
SimResults = Model.Sim_Sensitivity_Analysis(OptTable,Opt);
figure('Name','Sensitivity Analysis');
SimVaryPlot(SimResults, 'F' ,'F' );

5. Notes

5.a. Notes specific to qmt_spgr

5.a.1 BIDS-like convention

qMT-SPGR has not been defined in BIDS yet. However, following convention is suggested:

|== sub-01/
|~~~~~~ anat/
|---------- sub-01_flip-1_mt-1_QMTSPGR.json
|---------- sub-01_flip-1_mt-1_QMTSPGR.nii
|---------- sub-01_flip-2_mt-1_QMTSPGR.json
|---------- sub-01_flip-2_mt-1_QMTSPGR.nii
|---------- sub-01_flip-1_mt-2_QMTSPGR.json
|---------- sub-01_flip-1_mt-2_QMTSPGR.nii
|---------- sub-01_flip-2_mt-2_QMTSPGR.json
|---------- sub-01_flip-2_mt-2_QMTSPGR.nii
|---------- .
|---------- .
|---------- sub-01_flip-1_mt-N_QMTSPGR.json
|---------- sub-01_flip-1_mt-N_QMTSPGR.nii
|---------- sub-01_flip-2_mt-N_QMTSPGR.json
|---------- sub-01_flip-2_mt-N_QMTSPGR.nii
|
|== derivatives/
|~~~~~~ qMRLab/
|---------- dataset_description.json
|~~~~~~~~~~ sub-01/anat/
|-------------- sub-01_R1map.nii.gz (pre-exists)
|-------------- sub-01_R1map.json
|-------------- sub-01_Fmap.nii.gz
|-------------- sub-01_Fmap.json
|-------------- sub-01_kRmap.nii.gz
|-------------- sub-01_kRmap.json
|-------------- sub-01_R1Fmap.nii.gz
|-------------- sub-01_R1Fmap.json
|-------------- sub-01_R1Rmap.nii.gz
|-------------- sub-01_R1Rmap.json
|-------------- sub-01_T2Fmap.nii.gz
|-------------- sub-01_T2Fmap.json
|-------------- sub-01_T2Rmap.nii.gz
|-------------- sub-01_T2Rmap.json
|~~~~~~~~~~ sub-01/fmap/
|-------------- sub-01_TB1map.nii.gz
|-------------- sub-01_TB1map.json

Please visit sequence-specific metadata fields for applicable json keys.

5.b. Generic notes

5.b.1. Batch friendly option and protocol conventions

If you would like to load a desired set of options / protocols programatically, you can use *.qmrlab.mat files. To save a configuration from the protocol panel of qmt_spgr, first open the respective panel by running the following command in your MATLAB command window (MATLAB only):

Custom_OptionsGUI(qmt_spgr);

In this panel, you can arrange available options and protocols according to your needs, then click the save button to save my_qmt_spgr.qmrlab.mat file. This file can be later loaded into a qmt_spgr object in batch by:

Model = qmt_spgr;
Model = Model.loadObj('my_qmt_spgr.qmrlab.mat');

  Model.loadObj('my_qmt_spgr.qmrlab.mat') call won't update the fields in the Model object, unless the output is assigned to the object as shown above. This compromise on convenience is to retain Octave CLI compatibility.

If you don't have MATLAB, hence cannot access the GUI, two alternatives are available to populate options:

1. Use qmrlab/mcrgui:latest Docker image to access GUI. The instructions are available here.
2. Set options and protocols in CLI:

• List available option fields using tab completion in Octave's command prompt (or window)

Model = qmt_spgr;
Model.option. % click the tab button on your keyboard and list the available fields.

• Assign the desired field. For example, for a mono_t2 object:

Model = mono_t2;
Model.options.DropFirstEcho = true;
Model.options.OffsetTerm = false;

Some option fields may be mutually exclusive or interdependent. Such cases are handled by the GUI options panel; however, not exposed to the CLI. Therefore, manual CLI options assignments may be challenging for some involved methods such as qmt_spgr or qsm_sb. If above options are not working for you and you cannot infer how to set options solely in batch, please feel free to open an issue in qMRLab and request the protocol file you need.

Similarly, in CLI, you can inspect and assign the protocols:

Model = qmt_spgr;
Model.Prot. % click the tab button on your keyboard and list the available fields.

Each protocol field has two subfields of Format and Mat. The first one is a cell indicating the name of the protocol parameter (such as EchoTime (ms)) and the latter one contains the respective values (such as 30 x 1 double array containing EchoTimes).

The default Mat protocol values are set according to the example datasets served via OSF.

5.b.2 Parallelization

Beginning from release 2.5.0, you can accelerate fitting for the voxelwise models using parallelization.

Available in MATLAB only. Requires parallel processing toolbox.

In CLI, you can perform parallel fitting by:

parpool();
FitResults = ParFitData(data,Model);

If a parpool exists, the ParFitData will use it. If not, a new pool will be created using the local profile. By default, ParFitData saves outputs automatically every 5 minutes. You can disable this feature by:

FitResults = ParFitData(data, Model, 'AutosaveEnabled', false);

Alternatively, you can change the autosave interval (min 1 min) by:

FitResults = ParFitData(data,Model,'AutoSaveInterval',10);

If something went wrong during the fitting (e.g. your computer had to be restarted), you can recover the autosaved data by:

FitResults = ParFitData(data,Model,'RecoverDirectory','/ParFitTempResults_*/folder/from/the/previous/session');

GUI users will be prompted a question about whether they would like to use parallelization after clicking the Fit Data button, if the conditions are met. When called from GUI, ParFitData will be run with default options:

• Save temporary results every 5 minutes or whenever a chunk has finished processing
• Split data into chunks with a granularity factor of 3
• Do not remove temporary fit results upon completion

For further information:

help ParFitData

The default parallelization options can be changed in the preferences.json file located at the root qMRLab directory.

6. Citations

qMRLab JOSS article

Karakuzu A., Boudreau M., Duval T.,Boshkovski T., Leppert I.R., Cabana J.F., Gagnon I., Beliveau P., Pike G.B., Cohen-Adad J., Stikov N. (2020), qMRLab: Quantitative MRI analysis, under one umbrella 10.21105/joss.02343

Reference article for qmt_spgr

Sled, J.G., & Pike, G.B. (2000). Quantitative interpretation of magnetization transfer in spoiled gradient echo MRI sequences. Journal of Magnetic Resonance, 145(1), 24-36. 10.1006/jmre.2000.2059

---

Quantitative MRI, under one umbrella.

| qMRPullseq | qMRFlow | Interactive Tutorials |

NeuroPoly Lab, Montreal, Canada

---

Published with MATLAB® R2018b