Simulation

Simulation of CEST is based on the Bloch-McConnell equation. You will first need to define which species you want to simulate. Each species is a instance of the Species class.

A single saturation block-pulse can be simulated as follows, with a set of frequency offsets to obtain a Z-spectrum. This examples is based on the two-pools Bloch-McConnell equations, using the formalism of Numerical solution of the Bloch equations provides insights into the optimum design of PARACEST agents for MRI (Woessner et al. Magnetic Resonance in Medicine, 53(4), 2005).

import cest
import matplotlib.pyplot
import numpy

# Experimental conditions: gyromagnetic ratio (rad/s/T), main magnetic field (T and rad/s)
gamma = 267522187.44 # rad/s/T
B0 = 3 # T
w0 = gamma * B0 # rad/s

# Block saturation pulse: amplitude and duration
B1 = 12e-6 # T
w1 = gamma * B1 # rad/m
tau = 10 # s

# Species and their exchange rate (B -> A)
species_a = cest.Species(1, 0.2, 0, 1)
species_b = cest.Species(0.1, 0.1, 160, 0.0003636)
Cb = 500 # Hz

# Simulate between -100 and +200 ppms
offsets = numpy.linspace(-100, +200, 401)
magnetization = numpy.array([
    cest.bm(
        species_a, species_b, Cb, w0, offset, w1, tau,
        [0, 0, species_a.M0, 0, 0, species_b.M0])
    for offset in offsets])

# Plot the results
matplotlib.pyplot.plot(offsets, magnetization[:,2]/species_a.M0, lw=1)
matplotlib.pyplot.xlabel("$\Delta\omega_{RF}$ (ppm)")
matplotlib.pyplot.ylim(0, 1)
matplotlib.pyplot.ylabel("$M_z/M_0$")
matplotlib.pyplot.gca().invert_xaxis()

Shaped pulses can be defined using the pre-defined shapes or by adding your own. Each shape is normalized and discretized, so it needs to be scaled:

# Discretization step
step = 10e-3 # s

# Shaped pulse. The RF pulses are normalized: multiply by ω₁⋅N to get a
# saturation equivalent to the block pulse defined above.
steps = int(round(tau/step)) # Unitless
pulse = cest.pulses.gaussian(steps) * w1 * steps # rad/s

The shaped pulse is then simulated using the same function as above:

magnetization = numpy.array([
    cest.bm(
        species_a, species_b, Cb, w0, offset, pulse, step,
        [0, 0, species_a.M0, 0, 0, species_b.M0])
    for offset in offsets])

# Plot the results
matplotlib.pyplot.plot(offsets, magnetization[:,2]/species_a.M0, lw=1)