GENsim - Matrix Method Diffusion MRI Simulator

Introduction

GENsim is a Matlab toolkit designed to simulate diffusion MRI signal for a large variety of pulse sequences in diffusion substrates which include restriction. It provides the means to understand the diffusion patterns in synthetic tissue models (white matter, cancer, etc) for various diffusion sequences.

This tutorial describes the use of GENsim Matlab toolkit for simulating diffusion MRI data in various substrates.

The novelty of GENSim is that it allows the user to create a flexible diffusion gradient waveform. It assumes a basic understanding of diffusion MRI pulse sequences, including more advanced sequences such as square and trapezoidal oscillating gradients. The diffusion signal is calculated using the matrix formalism detailed in Drobnjak et al. 2010 and Drobnjak et al. 2011 . The existing parametrizations for square and trapezoidal oscillating gradients follow the description in Ianus et al. 2012 . The diffusion substrates available in GENSim are the basic compartments presented in Panagiotaki et al. 2011, as well as some two- and three-compartment models of white matter.

Terms and Conditions

GENSim is distributed under the Artistic License 2.0. The full text of the license can be found here . The Owner draws your attention to the fact that the Software has been developed for and is intended for use in a research environment only. No endorsement can be given for other use including, but not limited to, use in a clinical environment.

If you are using GENSim in your research, please cite Drobnjak et al. 2010 , Drobnjak et al. 2011 and Ianus et al. 2012

Download

Coming soon

Step-by-step getting started tutorial

Include GENSim toolbox in the Matlab search path:
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=matlab_path button=1 :)

Assuming that you have unpacked the GENSim toolbox in the directory /usr/local/GENSim, then run

addpath(genpath('/usr/local/GENSim'))

Take a minute to understand the format of the gradient waveform necessary for the simulation:
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=GEN_sequence button=1 :)

The information about the diffusion protocol is passed to the simulator as the fields of a Matab structure, which we conveniently call "protocol" in the examples. To generate the diffusion signal for a generalized waveform, the structure protocol should have the following fields: protocol.pulseseq = 'GEN' - the name of the sequence, protocol.G - the gradient waveform, protocol.tau - time interval between two points on the gradient waveform and optionally protocol.smalldel and protocol.delta which are the pulse duration and time interval between the onset of the first and second pulse. It is not mandatory to specify the last two parameters, but they save computation time if the diffusion sequence has two gradient intervals. Additionally the user can specify if the second gradient is repeated protocol.mirror = 0 - default value or reflected protocol.mirror = 1 .

The gradient waveform protocol.G is a M x 3K matrix which stores the values of the diffusion gradient at each time point. Each row represents one measurement and contains the values of the gradients in x, y and z direction. M is the total number of measurements and K is the number of gradient points in each direction for one measurement:

G1x(1) G1y(1) G1z(1) G1x(2) G1y(2) G1z(2) ... G1x(K) G1y(K) G1z(K)

G2x(1) G2y(1) G2z(1) G2x(2) G2y(2) G2z(2) ... G2x(K) G2y(K) G2z(K)

...

GMx(1) GMy(1) GMz(1) GMx(2) GMy(2) GMz(2) ... GMx(K) GMy(K) GMz(K)

protocol.G includes the complete diffusion gradient and the user should take into account any realistic situations such as the duration of the rf pulses, crusher gradients, etc.

protocol.G should satisfy the echo condition, that the integral of the gradient over time should be 0.

Take a minute to understand the diffusion substrates available in the simulation:
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=model button=1 :)

GENSym uses the 3D extension of the Matrix Method (first introduced by Callaghan JMR 95) presented in Drobnjak et al 2011. For the diffusion substrates, we follow the naming scheme presented in Panagiotaki et al. 2011. The following substrates are available in GENSim:

Basic compartments with Gaussian diffusion:

Ball (isotropic free diffusion)

Zeppelin (anisotropic, cylindrically symmetric diffusion tensor)

Tensor (full diffusion tensor)

Basic compartments with restricted diffusion:

Stick

AstroSticks (isotropically oriented sticks)

Cylinder

AstroCylinders (isotropically oriented cylinders)

Sphere

Dot (fully restricted)

Multi-compartment models:

ZeppelinCylinder

TortZeppelinCylinder (same as ZeppelinCylinder, but with tortuosity constraint on volume fraction)

TortZeppelinCylinderBall

Other substrates can be easily implemented by the user once familiarized with the code.

All the information related to the diffusion substrate is stored in a Matlab structure, which we call “model” in the examples. For the purpose of signal generation, “model” has only two fields: model.name – the name of the model (listed above) and model.params - the values of the model parameters in S.I. (m, s, etc.) of each model. See documentation for details.

Create an example gradient waveform for square oscillating gradients
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=example_SWOGSE button=1 :)

We provide the means of generating the gradient waveform for several parametrized sequences by calling the wave_form function. The available sequences are: pulse gradient spin echo sequence PGSE, square oscillating gradient - SWOGSE, trapezoidal oscillating gradients - TWOGSE, square oscillating gradients with different parameters in x,y and z directions - SWOGSE_3D. See documentation for details.

Here we show an example for SWOGSE:

* Setup an initial protocol which includes the necessary information to compute the discrete waveform for a SWOGSE sequence

 % add some combinations of parameters:
delta = 0.015:0.005:0.04; % duration in s
smalldel = delta - 0.005;  % duration in s
Nosc = 1:5;  % number of lobes in the oscillating gradient. A full period has 2 lobes
G = 0.08; % gradient strength in T/m;
tau = 1E-4; % time interval for waveform discretization
it = 0;
protocol_init.pulseseq = 'SWOGSE';
for i = 1:length(delta)
    for j = 1:length(Nosc)
        it = it+1;
        protocol_init.smalldel(it) = smalldel(i);
        protocol_init.delta(it) = delta(i);
        protocol_init.omega(it) = Nosc(j).*pi./smalldel(i);
        protocol_init.G(it) = G;
        protocol_init.grad_dirs(it,:) = [1 0 0]; % gradient in x direction
    end
end
protocol_init.tau = tau;

* Create the GEN protocol

protocolGEN.pulseseq = 'GEN';
protocolGEN.G = wave_form(protocol_init);
protocolGEN.tau = tau;
% include smalldel and delta as they make computation slightly faster
protocolGEN.delta = protocol_init.delta;
protocolGEN.smalldel = protocol_init.smalldel;

View the example waveform for square oscillating gradients
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=view_SWOGSE button=1 :)

The previous diffusion protocol has 30 measurements with varying timing parameters and number of oscillations. To visualize the gradient waveform for the first measurement, run the following code:

figure();
Gx = protocolGEN.G(1,1:3:end);
plot((0:tau:(length(Gx)-1)*tau)*1E3,Gx*1000,'LineWidth',2)
xlim([0,(protocolGEN.smalldel(1)+protocolGEN.delta(1)+10*tau)*1E3])
ylim([min(Gx)*1200 max(Gx)*1200])
xlabel('time (ms)');
ylabel('Gx (mT/m)');

Synthesize the diffusion signal for the previous example
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=run_SWOGSE button=1 :)

In order to synthesize diffusion signal, we need to choose a diffusion substrate.

Chose a white-matter model:

model.name = 'ZeppelinCylinder';
di = 1.7E-9; % intrinsic diffusivity
dh = 1.2E-9; % hindered diffusivity
rad = 4E-6; % cylinder radius
% angles in spherical coordinates describing the cylinder orientation; 
theta = 0; % angle from z axis
phi = 0; % azimuthal angle
ficvf = 0.7; % intracellular volume fraction
model.params = [ficvf di dh rad theta phi];

Add the matrices and other constants required by the Matrix Method:

protocolGEN = MMConstants(model,protocolGEN);

Compute the diffusion signal for the given protocol and model:

signal = SynthMeas(model,protocolGEN);

Plot the synthesized diffusion signal
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=plot_signal button=1 :)

To visualize the diffusion signal computed in the previous step as a function of diffusion time, run the following code:

signal_matrix = reshape(signal,length(Nosc),length(delta));
figure();
colors = jet(length(Nosc));
hold on;
for i  = 1:length(Nosc)  
    legend_name{i} = ['Nosc = ' num2str(Nosc(i))];
    plot(delta*1000,signal_matrix(i,:),'LineWidth',2,'Color',colors(i,:));
end
xlabel('\Delta (ms)','FontSize',16);
ylabel('Diffusion Signal','FontSize',16);
set(gca,'FontSize',16);
title('Diffusion signal as a function of \Delta for various number of oscillations, G = 80mT/m')
legend(legend_name,'FontSize',16)

A true advantage of GENSim is that it can compute the diffusion signal for any sequence, not only parametrized ones. An example with random values
(:toggle init=hide hide="Hide the detail" show="Show the detail" div=rand_protocol button=1 :)

Here you will see how to create a valid gradient waveform from a vector which specifies the value of the gradient at each point in time. The direction of the gradient is given by the field grad_dirs, similarly to the parametrized versions.

We can use the function wave_form with a protocol which has pulseseq = 'GEN' and a field GENGx which is the vector with the gradient values. To understand the resulting waveform, think of a PGSE where the constant gradient was replaced by GENGx.

Set up the initial protocol:

delta = 0.04; % duration in s

smalldel = delta - 0.005; % duration in s G = 0.08; % maximum gradient strength in T/m; protocol_rand.pulseseq = 'GEN'; protocol_rand.grad_dirs = [1/sqrt(2) 1/sqrt(2) 0 ]; % gradient along x and y protocol_rand.smalldel = smalldel; protocol_rand.delta = delta; % protocol_rand consists of a random waveform with 40 points which will be repeated after the 180 pulse protocol_rand.GENGx = rand(1,40)*G; protocol_rand.tau = protocol_rand.smalldel(1)./length(protocol_rand.GENGx);% calculate the time interval based on the number of points

Create the protocol needed for the diffusion simulation

protocolGEN.pulseseq = 'GEN';
protocolGEN.G = wave_form(protocol_rand);
protocolGEN.tau = tau;
% include smalldel and delta as they make computation slightly faster
protocolGEN.delta = protocol_rand.delta;
protocolGEN.smalldel = protocol_rand.smalldel;

View the protocol

figure();
plot(protocol_rand.GENGx)
figure();
Gx = protocolGEN.G(1,1:3:end);
plot((0:protocolGEN.tau:(length(Gx)-1)*protocolGEN.tau)*1E3,Gx*1000,'LineWidth',2)
xlabel('time (ms)');
ylabel('Gx (mT/m)');
xlim([0,(protocolGEN.smalldel(1)+protocolGEN.delta(1)+2*protocolGEN.tau)*1E3])

Choose the diffusion substrate (same as before)

model.name = 'ZeppelinCylinder';
di = 1.7E-9; % intrinsic diffusivity
dh = 1.2E-9; % hindered diffusivity
rad = 4E-6; % cylinder radius
% angles in spherical coordinates describing the cylinder orientation; 
theta = 0; % angle from z axis
phi = 0; % azimuthal angle
ficvf = 0.7; % intracellular volume fraction
model.params = [ficvf di dh rad theta phi];

Add the matrices and other constants required by the Matrix Method:

protocolGEN = MMConstants(model,protocolGEN);

Compute the diffusion signal for the given protocol and model:

signal = SynthMeas(model,protocolGEN);

Create your own diffusion protocol.
To help you start, the same commands are in GENSim/example/RunGEN.m