T2star decay with CMPG sequence

In this notebook, the functionality of the T2star module is demonstrated by simulating a CMPG/Turbo-spin-echo sequence without spatial encoding, but the signal is acquired on every simulation step.

Imports

import os
os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'

# Load TF and check for GPUs
import tensorflow as tf
gpus = tf.config.list_physical_devices('GPU')
for gpu in gpus:
    tf.config.experimental.set_memory_growth(gpu, True)
logical_gpus = tf.config.list_logical_devices('GPU')
print(len(gpus), "Physical GPUs,", len(logical_gpus), "Logical GPUs")


# 3rd Party dependencies
from pint import Quantity
import pyvista
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
%matplotlib inline

# Project library cmrsim
import sys
sys.path.append("../../")
sys.path.append("../../../cmrseq")
import cmrseq
import cmrsim
import cmrsim.utils.particle_properties as part_factory

1 Physical GPUs, 1 Logical GPUs

Define CPMG sequence with hard-pulses

# Define MR-system specifications
system_specs = cmrseq.SystemSpec(max_grad=Quantity(80, "mT/m"), max_slew=Quantity(200., "mT/m/ms"),
                                 grad_raster_time=Quantity(0.01, "ms"),
                                 rf_raster_time=Quantity(0.01, "ms"),
                                 adc_raster_time=Quantity(0.01,"ms"))

# Set up relevant sequence blocks

delay = cmrseq.bausteine.Delay(system_specs=system_specs, duration = Quantity(0.3,'ms'))
end_delay = cmrseq.bausteine.Delay(system_specs=system_specs, duration = Quantity(10,'ms'))

pulse90 = cmrseq.bausteine.SincRFPulse(system_specs=system_specs, flip_angle=Quantity(90,'degrees'),
                                      duration=Quantity(0.01,'ms'),
                                      time_bandwidth_product=0,
                                      frequency_offset=Quantity(0,'Hz'),
                                      center=0.5, delay=Quantity(0, "ms"),
                                      apodization=0., name="90")

pulse180 = cmrseq.bausteine.SincRFPulse(system_specs=system_specs, flip_angle=Quantity(180,'degrees'),
                                      duration=Quantity(0.01,'ms'),
                                      time_bandwidth_product=0,
                                      frequency_offset=Quantity(0,'Hz'),
                                      center=0.5, delay=Quantity(0, "ms"),
                                      apodization=0., name="90")

adc_block90 = cmrseq.bausteine.SymmetricADC(system_specs=system_specs, num_samples=100,
                                            duration=Quantity(25,'ms'),
                                            delay=Quantity(0,'ms'))

adc_block180 = cmrseq.bausteine.SymmetricADC(system_specs=system_specs, num_samples=100,
                                            duration=Quantity(50,'ms'),
                                            delay=Quantity(0,'ms'))


# Create sequence with 90-180-180-180...
# Delays are required around each RF object to prevent sampling during RF
seq = cmrseq.Sequence([delay],system_specs=system_specs)
seq.append(pulse90)
seq.append(delay)
seq.append(adc_block90)
seq.append(delay)

seq.append(pulse180)
seq.append(delay)
seq.append(adc_block180)
seq.append(delay)

seq.append(pulse180)
seq.append(delay)
seq.append(adc_block180)
seq.append(delay)

seq.append(pulse180)
seq.append(delay)
seq.append(adc_block180)
seq.append(delay)
seq.append(end_delay)

plt.close("all")
f, a = plt.subplots(1, 1,figsize=(16, 5))
cmrseq.plotting.plot_sequence(seq, axes=[a,a,a,a], format_axes=False, add_legend=True)

Set up Bloch operator

# Define T2* module
T2s = cmrsim.bloch.submodules.T2starDistributionModel()
submodules = [T2s]

# Define bloch simulator module
time, rf_grid, grad_grid, adc_on_grid = [np.stack(v) for v in cmrseq.utils.grid_sequence_list([seq], force_uniform_grid=False)]

bloch_module = cmrsim.bloch.GeneralBlochOperator(name="FID", gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
                                                 time_grid=time[0],
                                                 gradient_waveforms=grad_grid,
                                                 rf_waveforms=rf_grid,
                                                 adc_events=adc_on_grid,
                                                 device="GPU:0", submodules=submodules)

Set up particles and relaxation properties

n_particles = 100000
T1 = 10000
T2 = 300
T2p = 30

properties = {"M0": part_factory.uniform(default_value=1.)(n_particles),
              "T1": part_factory.uniform(default_value=T1)(n_particles),
              "T2": part_factory.uniform(default_value=T2)(n_particles),
              "omega_t2s": part_factory.t2star_lorentzian(T2p = T2p,prctile_cutoff=0.002)(n_particles),
              "magnetization": part_factory.norm_magnetization()(n_particles),
              "initial_position": np.zeros([n_particles,3]).astype(np.float32)}

# Run simulation
m_init = properties.pop("magnetization")
initial_position = properties.pop("initial_position")
bloch_module.reset()

Run Bloch-simulation

_,_ = bloch_module(initial_position=initial_position, magnetization=m_init, **properties)

Plot Signal decay

sig= np.abs(bloch_module.time_signal_acc[0].numpy())
sig = sig/sig[0]

t = time[adc_on_grid[:,:,0]==1]

T2s = 1/(1/T2+1/T2p)

sig = sig*np.exp(-t[0]/T2s)

et2 = np.exp(-t/T2)
et2s = np.exp(-t/T2s)

fig, axes = plt.subplots(1, 1, figsize=(14, 5))
axes.plot(t, sig, label='Simulated CPMG')
axes.plot(t, et2, label='Expected T2')
axes.plot(t, et2s, label='Expected T2s')
axes.legend()
axes.set_xlabel('Time (ms)'), axes.set_ylabel('Normalized Signal (a.u.)')

(Text(0.5, 0, 'Time (ms)'), Text(0, 0.5, 'Normalized Signal (a.u.)'))