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

!pip install cmrseq --index-url https://gitlab.ethz.ch/api/v4/projects/30300/packages/pypi/simple

import sys
sys.path.append("../")
sys.path.insert(0, "../../")
sys.path.insert(0, "../../../cmrseq")

import os
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
os.environ['CUDA_VISIBLE_DEVICES'] = "3"

import tensorflow as tf
gpu = tf.config.get_visible_devices("GPU")
if gpu:
    tf.config.experimental.set_memory_growth(gpu[0], True)

# 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 cmrseq
import cmrsim
import cmrsim.utils.particle_properties as part_factory

Define CPMG sequence with hard-pulses

The sequence consist of a initial 90 degree and consecutive 180 degree pulses and it also includes contious readout blocks between the RF hardpulses to observe the evolution of magnetization

%matplotlib inline
# Define MR-system specifications
system_specs = cmrseq.SystemSpec(max_grad=Quantity(80, "mT/m"),
                                 max_slew=Quantity(200., "mT/m/ms"),
                                 rf_peak_power=Quantity(300, "uT"),
                                 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.HardRFPulse(system_specs, flip_angle=Quantity(90, "degrees"),
                                       duration=Quantity(0.01, "ms"), name="90")
pulse180 = cmrseq.bausteine.HardRFPulse(system_specs, flip_angle=Quantity(180, "degrees"),
                                        duration=Quantity(0.01, "ms"), name="180")

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([pulse90], system_specs=system_specs)
seq.extend([delay, adc_block90, delay,
            pulse180, delay, adc_block180, delay,
            pulse180, delay, adc_block180, delay,
            pulse180, delay, adc_block180, delay,
            end_delay], copy=True)

plt.close("all")
f = cmrseq.plotting.plot_sequence(seq, axes="single", format_axes=True, add_legend=True, legend_position="lower left");

Extending Sequence: 100%|██████████| 16/16 [00:00<00:00, 799.82it/s]

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)

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1711371318.306910 3560283 device_compiler.h:186] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

Results

As visualized in the plot below, the intermittent 180 degree refocussing pulses, restore the transversal magnetization decay due to T2 star. The green curve corresponds to pure FID decay including T2 star and the orange line illustrates the pure T2 decay according to the set material properties.

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.)'))