Pulseq-CEST
The following live script demonstrates how to generate an pulseq-file for an APT-weighted CEST experiment with Pulseq-CEST.
Installation
The Pulseq-CEST MATLAB Code can be cloned from git and is set up using the Installation script
% if you do not have git installed you can use the commented lines to unzip
% the code directly from GitHub:
% unzip("https://github.com/kherz/pulseq-cest/archive/master.zip");
% movefile('pulseq-cest-master', 'pulseq-cest');
system('git clone https://github.com/kherz/pulseq-cest');
cd pulseq-cest;
install_pulseqcest;
cd ..
The pulseq-cest-library repository was cloned during the installation:
dir('pulseq*')
The pulseq-cest folder contains the core funtionality and all the functions to simulate CEST protocols, whereas the pulseq-cest-library is a database for different protocols.
For this tutorial, we have a detailled look into the generation of an APT-weighted protocol from the pulseq-cest-library.
cd pulseq-cest-library\seq-library\APTw_3T_001_2uT_36SincGauss_DC90_2s_braintumor\
Each protol in the library comes with a MATLAB and python script to generate the pulseq-file, as well as a short description and of course the pulseq-file itself. Let's have a detailled look into the MATLAB generation file.
The seq_def struct contains everything that gets later stored as a definition in the pulseq-file. This can contain anything you would need later for post-processing etc.
seqid = 'APTw_3T_001_2uT_36SincGauss_DC90_2s_braintumor';
seq_defs.n_pulses = 36 ; % number of pulses
seq_defs.tp = 50e-3 ; % pulse duration [s]
seq_defs.td = 5e-3 ; % interpulse delay [s]
seq_defs.Trec = 3.5 ; % recovery time [s]
seq_defs.Trec_M0 = 3.5 ; % recovery time before M0 [s]
seq_defs.M0_offset = -1560 ; % m0 offset [ppm]
seq_defs.DCsat = (seq_defs.tp)/(seq_defs.tp+seq_defs.td); % duty cycle
seq_defs.offsets_ppm = [seq_defs.M0_offset -4:0.25:4]; % offset vector [ppm]
seq_defs.num_meas = numel(seq_defs.offsets_ppm) ; % number of repetition
seq_defs.Tsat = seq_defs.n_pulses*(seq_defs.tp+seq_defs.td) - ...
seq_defs.td ; % saturation time [s]
seq_defs.B0 = 3 ; % B0 [T]
seq_defs.seq_id_string = seqid ; % unique seq id
You can find them in the existing pulseq-file in the [DEFINITIONS]:
fid = fopen('APTw_3T_001_2uT_36SincGauss_DC90_2s_braintumor.seq');
while ~strcmp(fgetl(fid), '[DEFINITIONS]')
continue;
end
def = '[DEFINITIONS]';
while ~isempty(def)
disp(def);
def = fgetl(fid);
end
fclose(fid);
After the definitions, we prepare the variables we need to fill the pulseq-file.
Pulseq takes care of scanner limits, such as gradient slew rates, rf times, etc. We provide a function to set some low-demand limits.
lims = getScannerLimits()
In the next step we create the saturation pulses. Pulseq comes with some basic pulse shahes that can be further adapted with additional parameters. For this example we use a Sinc pulse with an adjusted time-bandwidth product and apodization.
B1pa = 1.78; % mean sat pulse b1 [uT]
gyroRatio_hz = 42.5764; % for H [Hz/uT]
gyroRatio_rad = gyroRatio_hz*2*pi; % [rad/uT]
fa_sat = B1pa*gyroRatio_rad*seq_defs.tp; % flip angle of sat pulse
satPulse = mr.makeSincPulse(fa_sat, 'Duration', seq_defs.tp, 'system', lims,'timeBwProduct', 2,'apodization', 0.15);
Pulseq uses a samle rate of 1 us for RF pulses, which can result in large files and long simulation and plotting times. One workaround is to downsample pulse and upsample them again with a nearest neighbour interpolation. This way, we keep the sample rate, but the internally used run-length-encoding algorithm is more effective. Since pulse objects are ususlly simulates with less samples, we do not loose any specificity here.
% resample pulse for reduced file size and io time
nPulseSamples = 1000;
satPulse = resamplePulseForRLE(satPulse, nPulseSamples);
It is also possible to calculate the different power equivalents from the saturation pulse and the duty cycle of the entire preparation. We use that to store the B1 continous wave power equivalent in the pulseq-file:
[B1cwpe,B1cwae,B1cwae_pure,alpha]= calculatePowerEquivalents(satPulse,seq_defs.tp,seq_defs.td,1,gyroRatio_hz);
seq_defs.B1cwpe = B1cwpe;
Now we have all the objects we need and can fill the pulseq-file by looping through the offsets we want to simulate / measure.
% init sequence with system limits
seq = SequenceSBB(lims);
% pulseq uses offsets in Hz
offsets_Hz = seq_defs.offsets_ppm*gyroRatio_hz*seq_defs.B0;
% loop through offsets
for currentOffset = offsets_Hz
% if the current offset is an M0 scan we use the M0 recovery time
if currentOffset == seq_defs.M0_offset*gyroRatio_hz*seq_defs.B0
if seq_defs.Trec_M0 > 0
seq.addBlock(mr.makeDelay(seq_defs.Trec_M0));
end
else % if no M0 scan use Trec
if seq_defs.Trec > 0
seq.addBlock(mr.makeDelay(seq_defs.Trec)); % recovery time
end
end
% set frequency offset of the pulse
satPulse.freqOffset = currentOffset;
% take care of the accumulated phase during the saturation
accumPhase=0;
% loop through pulses
for np = 1:seq_defs.n_pulses
% set accumulated pahse from previous rf pulse
satPulse.phaseOffset = mod(accumPhase,2*pi);
% add satURATION PULSE pulse
seq.addBlock(satPulse)
% calculate the phase for the next rf pulse
accumPhase = mod(accumPhase + currentOffset*2*pi*(numel(find(abs(satPulse.signal)>0))*1e-6),2*pi);
% delay between pulses
if np < seq_defs.n_pulses
seq.addBlock(mr.makeDelay(seq_defs.td)); % add delay
end
end
% add the standard spoiling gradients
seq.addSpoilerGradients();
% add the readout trigger event
seq.addPseudoADCBlock();
end
As a last step we add the definitions to the pulseq-file that we defined in the beginning and write the file.
%% set definitions
def_fields = fieldnames(seq_defs);
for n_id = 1:numel(def_fields)
seq.setDefinition(def_fields{n_id}, seq_defs.(def_fields{n_id}));
end
% write the pulseq-file with author info
seq_filename = 'Demo.seq';
author = 'John Doe';
seq.write(seq_filename, author);
Lets have a look at the first saturation phase:
seq.plotSaturationPhase();
This file could now be used at the scanner with the Pulseq-CEST hybrid sequence. But it is of course also possible to use it in our simulation. For this purpose we provide simulation parameters in the pulseq-cest-library. Lets have a look at the content of the simulation parameters file for gray matter at 3T.
param_filename = '../../sim-library/WM_3T_001_bmsim.yaml';
sim_parameters = yaml.ReadYaml(param_filename)
sim_parameters.water_pool
sim_parameters.mt_pool
sim_parameters.cest_pool
sim_parameters.cest_pool.amide
As you can see, all relevant parameters for simulation are stored in the human-readbale yaml-file. We can use this in combination with our pulseq-file directly in the simulation.
M_z = simulate_pulseqcest(seq_filename,param_filename);
The simulation returns the Z-magnetization which we can now plot together with the MTRasym.
plotSimulationResults(M_z, seq_defs.offsets_ppm, seq_defs.M0_offset);
