Discrete (arbitrary gradient moment) EPG#
In the discrete EPG model (C++, Python), the gradient moments may vary across time intervals. In this model, the orders of the model are stored in bins of user-specified width (hence the term “discrete”), expressed in rad/m. The interface is otherwise similar to that of regular EPG.
The following code sample show the full simulation of the DW-DESS sequence, with its two read-out modules and its diffusion module. Since the read-out gradient moment is largely smaller than the diffusion gradient moment, its simulation using regular EPG would require to subdivide every time interval to a common duration: with the parameters of the following simulation, even if we neglect the pre-phasing, each repetition would generate over 20 new states with regular EPG, while at most 12 will be generated with discrete EPG.
The discrete EPG model, as well as other models, includes a threshold member (C++, Python) which controls how states with a low population are discarded. This yields better peformance, at the expanse of a slight loss in precision. As an example, using a threshold of \(10^{-6}\) in the DESS simulation reduces the computation time from 23 ms to 15 ms (i.e. by a factor 1.5) while conserving a good precision, as shown in the figure below.
import time
import numpy
import sycomore
from sycomore.units import *
species = sycomore.Species(1000*ms, 100*ms, 1700*um**2/s)
# Sequence parameters
flip_angle = 20*deg
TE = 2*ms, 8*ms
TR = 10*ms
slice_thickness = 1*mm
tau_readout = 1*ms
q = 70/cm
# Motion to k-space extremity and its associated gradient amplitude
k_max = 0.5 * 2*numpy.pi/slice_thickness
G_readout = k_max / sycomore.gamma / (tau_readout/2)
# Diffusion gradient
tau_diffusion = (TE[1]-tau_readout/2) - (TE[0]+tau_readout/2)
G_diffusion = q/(sycomore.gamma_bar*tau_diffusion)
models = [sycomore.epg.Discrete(species), sycomore.epg.Discrete(species)]
models[1].threshold = 1e-6
repetitions = int((4*species.T1/TR))
S_plus = numpy.zeros((len(models), repetitions), dtype=complex)
S_minus = numpy.zeros((len(models), repetitions), dtype=complex)
for index, model in enumerate(models):
begin = time.time()
for repetition in range(repetitions):
# RF-pulse and idle until the first read-out
model.apply_pulse(flip_angle)
model.apply_time_interval(TE[0]-tau_readout)
# First echo
model.apply_time_interval(tau_readout/2, -G_readout)
model.apply_time_interval(tau_readout/2, G_readout)
S_plus[index, repetition] = model.echo
model.apply_time_interval(tau_readout/2, G_readout)
# Diffusion gradient between the two echoes
model.apply_time_interval(tau_diffusion, G_diffusion)
# Second echo
model.apply_time_interval(tau_readout/2, G_readout)
S_minus[index, repetition] = model.echo
model.apply_time_interval(tau_readout/2, G_readout)
model.apply_time_interval(tau_readout/2, -G_readout)
# Idle until the end of the TR
model.apply_time_interval(TR-TE[1]-tau_readout)
# Make sure the sequence timing is correct
if repetition == 0:
assert((model.elapsed-TR)/TR < 1e-6)
end = time.time()
print(
"Threshold:", model.threshold, len(model), "orders",
1e3*(end-begin), "ms")
#include <chrono>
#include <cmath>
#include <iostream>
#include <sycomore/epg/Discrete.h>
#include <sycomore/Species.h>
#include <sycomore/sycomore.h>
#include <sycomore/units.h>
using Clock = std::chrono::high_resolution_clock;
using Duration = std::chrono::duration<double>;
int main()
{
using namespace sycomore::units;
sycomore::Species species(1000*ms, 100*ms, 1700*std::pow(um, 2)/s);
// Sequence parameters
auto flip_angle = 20*deg;
sycomore::Vector2Q TE{2*ms, 8*ms};
auto TR = 10*ms, slice_thickness = 1*mm, tau_readout = 1*ms, q = 70/cm;
// Motion to k-space extremity and its associated gradient amplitude
auto k_max = 0.5 * 2*M_PI/slice_thickness;
auto G_readout = k_max / sycomore::gamma / (tau_readout/2);
// Diffusion gradient
auto tau_diffusion = (TE[1]-tau_readout/2) - (TE[0]+tau_readout/2);
auto G_diffusion = q/(sycomore::gamma_bar*tau_diffusion);
std::vector<sycomore::epg::Discrete>models{
sycomore::epg::Discrete(species), sycomore::epg::Discrete(species)};
models[1].threshold = 1e-6;
std::size_t repetitions = std::lround(4*species.T1()/TR);
sycomore::TensorC<2> S_plus({models.size(), repetitions}, 0);
sycomore::TensorC<2> S_minus({models.size(), repetitions}, 0);
for(std::size_t index=0; index!=models.size(); ++index)
{
auto & model = models[index];
auto const begin = Clock::now();
for(std::size_t repetition=0; repetition!=repetitions; ++repetition)
{
// RF-pulse and idle until the first read-out
model.apply_pulse(flip_angle);
model.apply_time_interval(TE[0]-tau_readout);
// First echo
model.apply_time_interval(tau_readout/2, -G_readout);
model.apply_time_interval(tau_readout/2, G_readout);
S_plus(index, repetition) = model.echo();
model.apply_time_interval(tau_readout/2, G_readout);
// Diffusion gradient between the two echoes
model.apply_time_interval(tau_diffusion, G_diffusion);
// Second echo
model.apply_time_interval(tau_readout/2, G_readout);
S_minus(index, repetition) = model.echo();
model.apply_time_interval(tau_readout/2, G_readout);
model.apply_time_interval(tau_readout/2, -G_readout);
// Idle until the end of the TR
model.apply_time_interval(TR-TE[1]-tau_readout);
// Make sure the sequence timing is correct
if(repetition == 0)
{
assert((model.elapsed()-TR)/TR < 1e-6);
}
}
auto const end = Clock::now();
std::cout
<< "Threshold:" << model.threshold << ", "
<< model.size() << " orders, "
<< 1e3*Duration(end-begin).count() << " ms\n";
}
return 0;
}
Note
The default bin width of \(\delta_k=1\ \mathrm{rad/m}\) corresponds to a gradient area of \(\delta_k/\gamma \approx 3\ \mathrm{\mu T/m\, ms}\). The dephasing orders are stored as 64-bits, signed integers; the default bin width hence allows to handle dephasing orders in the range \(\left[-9\cdot 10^{18}\ \mathrm{rad/m}, +9\cdot 10^{18}\ \mathrm{rad/m} \right]\), i.e. \(\left[-3\cdot 10^{16}\ \mathrm{mT/m\, ms}, +3\cdot 10^{16}\ \mathrm{mT/m\, ms} \right]\)
Simulation of DW-DESS with discrete EPG, without threshold (plain lines) and with threshold (dots)#