CoCalc -- Bloch-simulation-with-exchange.ipynb

Download

IPython notebook Bloch-simulation-with-exchange.ipynb

Bloch Simulation

Bloch-simulation-with-exchange.ipynb

⁶² views

Kernel: Unknown Kernel

Crash course in ODEs with scipy for Bloch simulations

Let's say we have $\frac{dy}{dt} = \dot{y} = a \cdot y$

Clearly the solution is $y=y_0 \,e^{at} .$

In [3]:

import numpy
from scipy.integrate import odeint, ode

a=-1/.05
# in this trivial example the derivative does not explicitly depend on t
def ydot(y, t):
    return a*y

def ydot_dopri(t, y): # note the order of parameters is different!
    return a*y

# we need a vector of times where we want to evaluate the solution
delta_t = 0.002
t_output = delta_t * numpy.arange(100)

The simple odeint interface is easy to handle. Just give it a function handle for the derivative, a starting point, and some points where you want answers. Done.

In [4]:

y0 = 1.0
y_result = odeint(ydot, y0, t_output)
y_result = y_result[:,0] # let's make this result a 1D vector
y_analytics = y0 * numpy.exp(a*t_output)

More setup required for the generic ode interface but in exchange for a choice of integrator (Runge Kutta 4-5 in this case) comes more foot work.

In [5]:

r = ode(ydot_dopri)
r = r.set_integrator('dopri5')
r = r.set_initial_value(y0, t=0)
#r = r.set_f_params
y_result_dopri = numpy.empty_like(t_output)+numpy.nan
y_result_dopri[0]=y0
idx = 1
for tf in t_output[1:]:
    y_result_dopri[idx] = r.integrate(tf)
    idx +=1

Let's make sure the ODE result and the analytical solution agree.

In [12]:

import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] =(10,5)
pylab.subplot(221)
pylab.plot(t_output, y_result, label='ode solution')
pylab.plot(t_output, y_analytics, label='analytical soln')
pylab.legend()
pylab.subplot(222)
pylab.plot(t_output, (y_analytics-y_result)/(y_analytics))

pylab.subplot(223)
pylab.plot(t_output, y_result_dopri, label='ode(dopri) solution')
pylab.plot(t_output, y_analytics, label='analytical soln')
pylab.legend()
pylab.subplot(224)
pylab.plot(t_output, (y_analytics-y_result_dopri)/(y_analytics))
pylab.title('orders of magnitude better')

Out[12]:

<matplotlib.text.Text at 0x7fdd8585dad0>

DOPRI looks better but the real question of course is to find out whether this is due to denser points since you have to compare Cox with Granny Smith not Valencias. Short answer, dopri is more efficient.

Onwards to higher orders and applications of Bloch equations

Let's try this with some real stuff. For Bloch equations we have to move on from 1D. The ODE is: $\frac{d\vec{M}}{dt}=\mathbf{A}\cdot \vec{M} - \vec{B}$ or more specifically:

ParseError: KaTeX parse error: Undefined control sequence: \eqalignno at position 1: \̲e̲q̲a̲l̲i̲g̲n̲n̲o̲{ &{dM_{x}(t)\…

where the B1 field of the RF pulse is along x ("phase=0"); so rotation of $\vec{M}$ only happens around x-axis. The frequency of the RF pulse is $\omega_{RF}$ and the magnetization is being observed in the rotating frame which rotates with $\omega_{RF}$ which can be off-resonance from the Larmor frequency $\omega_0=\gamma B_0$ by $\Delta\omega=\omega_0-\omega_{RF}$ . This causes an apparent rotation of the magnetization $\vec{M^{\prime}}$ (short $\vec{M}$ since we drop the $\prime$ ) with that frequency $\Delta\omega$ .

With the above in mind, let's define a $\frac{d\mathbf{M}}{dt}$ that can be used for integration:

In [15]:

gamma = 2*numpy.pi*42.6e6 # rad/(s T)
B1 = 1e-6 # Tesla
omega1 = gamma * B1
T2=0.10 # seconds
T1=numpy.infty # seconds
def dMdt(t, M_vec, T1=T1, T2=T2, M0=1, domega=0, omega1=omega1):
    # on-resonance pulse: domega=0
    # B1=1uT -> omega1= 2pi * 42.6e6 rad/(s T) * 1uT = 267 rad/s  
    # inhomogeneous term in ODE
    
#    B = numpy.array([M_vec[0]/T2, M_vec[1]/T2, (M_vec[2]-M0)/T1])
    B = numpy.array([0, 0, -M0/T1])

    A = numpy.array([[-1./T2 ,  +domega , 0],
                     [-domega, -1./T2       , +omega1],
                     [0,   -omega1  , -1./T1]])
    return numpy.dot(A,M_vec) - B

In [17]:

dt = 0.002
n_pts = 100
max_t = n_pts*dt

r = ode(dMdt)
r = r.set_integrator('dopri5')
#r = r.set_integrator('dop853')#, atol=1e-10, rtol=1e-10)
Mstart = numpy.array([0,0,1])
#Mstart=array([0,1,0])
r = r.set_initial_value(Mstart, t=0)
#r = r.set_f_params
t_output = numpy.empty(n_pts)
M_result = numpy.empty((n_pts,3))
t_output[0]=0
M_result[0,:]=Mstart
idx = 1
while r.successful() and r.t < max_t and idx<n_pts:
    M_result[idx,:] = r.integrate(r.t+dt)
    t_output[idx]=r.t
    idx += 1

It might be nice to check the results by using some analytical expressions for this. Some helpful discussions can be found on Brian Hargreave's website on this. He has a tutorial for this: http://mrsrl.stanford.edu/~brian/bloch/

For the near analytical solution we need to define some operators for what would constitute an RF-pulse as well as one for evolution under free precession:

In [18]:

import numpy

def xrot(phi):
    return numpy.array([[1,0,0],
                        [0, numpy.cos(phi), -numpy.sin(phi)],
                        [0, numpy.sin(phi), numpy.cos(phi)]])

def yrot(phi):
    return numpy.array([[numpy.cos(phi), 0, numpy.sin(phi)],
                        [0,1,0],
                        [-numpy.sin(phi), 0, numpy.cos(phi)]])

def zrot(phi):
    return numpy.array([[numpy.cos(phi), -numpy.sin(phi), 0],
                        [numpy.sin(phi), numpy.cos(phi), 0],
                        [0,0,1]])

def freeprecess(dt, T1=numpy.inf, T2=numpy.inf, Mo=1.0, domega=0):
    ''' return the A matrix and B vector for the dM/dt magnetization evolution '''
    phi = domega*dt	 # Resonant precession, radians.
    E1 = numpy.exp(-dt/T1)
    E2 = numpy.exp(-dt/T2)

    B = numpy.array([0, 0, Mo*(1 - E1)])

    A = numpy.array([[E2,0 , 0],
                     [0, E2, 0],
                     [0, 0 , E1]])
    return numpy.dot(A, zrot(phi)),B

Before we compare the 'analytical' solution to the ODE one, let's see whether some straightforward free precession with T1 and T2 relaxation gives us what we expect: An FID.

In [21]:

T1_fp = .600 
T2_fp = .170
df_fp = 10.*2*numpy.pi
dt_fp = 0.001
n_pts_fp = 1000
A_fp, B_fp = freeprecess(dt_fp, T1=T1_fp, T2=T2_fp, domega=df_fp)
print A_fp, B_fp

M_result_fp = numpy.empty((n_pts_fp,3))
# we are starting with Mx=1 and My=Mz=0
# this is the situation **after** a hard 90 pulse
M_result_fp[0,:] = numpy.array([1,0,0])

# let's evolve the magnetization step by step
for i in xrange(n_pts_fp-1):
    M_result_fp[i+1,:] = numpy.dot(A_fp,M_result_fp[i,:])+B_fp

Out[21]:

[[ 0.99217322 -0.06242225  0.        ]
 [ 0.06242225  0.99217322  0.        ]
 [ 0.          0.          0.99833472]] [ 0.          0.          0.00166528]

In [24]:

pylab.subplot(121)
pylab.plot(M_result_fp[:,0],M_result_fp[:,1])
pylab.title('signal in transverse plane')
pylab.subplot(122)
t = dt_fp*numpy.arange(n_pts_fp)
pylab.plot(t,M_result_fp[:,0],label='Mx')
pylab.plot(t,M_result_fp[:,1],label='My')
pylab.plot(t,M_result_fp[:,2],label='Mz')
pylab.plot(t, numpy.exp(-t/T2_fp)) 
pylab.title('FID')
pylab.legend()

Out[24]:

<matplotlib.legend.Legend at 0x7fdd85b40b50>

Now, let's see whether the evolution using our ODE further up corresponds to what we find the direct way. This is a simple sanity check.

In [25]:

idx = 1
pulse = xrot(dt*omega1)
A,B = freeprecess(dt, T2=T2)
M_result_analytical = numpy.empty((n_pts,3))
M_result_analytical[0,:]=Mstart
print A,B
for idx in xrange(1,n_pts):
    M_result_analytical[idx,:] = numpy.dot(A, numpy.dot(M_result_analytical[idx-1,:],pulse)) + B
    idx += 1

Out[25]:

[[ 0.98019867  0.          0.        ]
 [ 0.          0.98019867  0.        ]
 [ 0.          0.          1.        ]] [ 0.  0.  0.]

In [29]:

pylab.subplot(221)
pylab.plot(t_output, M_result[:,1], label='integrated My')
pylab.plot(t_output, M_result_analytical[:,1], label='analytical My')
pylab.plot(t_output, numpy.exp(-t_output/T2/2), label='T2 envelope')
pylab.legend()
pylab.subplot(222)
pylab.plot(t_output, M_result[:,1] - M_result_analytical[:,1], label='difference')
pylab.legend()
pylab.subplot(223)
pylab.plot(t_output, M_result[:,2], label='integrated Mz')
pylab.plot(t_output, M_result_analytical[:,2], label='analytical Mz')
pylab.plot(t_output, numpy.exp(-t_output/T2/2), label='T2 envelope')
pylab.legend()
pylab.subplot(224)
pylab.plot(t_output, M_result[:,2] - M_result_analytical[:,2], label='difference')
pylab.legend()

Out[29]:

<matplotlib.legend.Legend at 0x7fdd85631e50>

Bonus points: Why is the envelope an exponential that decays with 2*T2?

Next steps

two-pool version
exchange terms
extension to composite evolutions: saturation pulses with exchange, hard pulse readout (FLASH), and relaxation...
larger number of pools
could this even be used for fitting?!

In [ ]:

Crash course in ODEs with scipy for Bloch simulations

Onwards to higher orders and applications of Bloch equations

Next steps

Product

Resources

Company