Class: Psych 204a
Tutorial: Diffusion
Author: Dougherty
Date: 2001.10.31
Duration: 90 minutes
Copyright: Stanford University, Robert F. Dougherty
Translated to Python by Bob Dougherty, 11/2012 and Grace Tang 10/13
The purpose of this tutorial is to illustrate the nature of the data acquired in a diffusion-weighted imaging scan. The computational methods available for interpreting these data are also introduced.
First, we'll set up the python environment and define some utility functions that will be used below.
In [ ]:
%pylab inline
rcParams["figure.figsize"] = (8, 6)
rcParams["axes.grid"] = True
from IPython.display import display, clear_output
from mpl_toolkits.axes_grid1 import make_axes_locatable
from time import sleep
from __future__ import division
def cart2pol(x, y):
theta = arctan2(y, x)
r = sqrt(x ** 2 + y ** 2)
return theta, r
def pol2cart(theta, r):
x = r * cos(theta)
y = r * sin(theta)
return x, y
The self-diffusion coefficient of water (in micrometers2/millisecond) is dependent on the temperature and pressure. Several groups have derived quantitative descriptions of the relationship between temperature, pressure, and the diffusion coefficient of water. Here we use the formula presented in:
Let's start by defining a function that implements the Krynicki formula. The the default value for the pressure parameter will be set to the standard atmospheric pressure at sea level: 101.325 kilo Pascals (kPa).
In [ ]:
def selfDiffusionOfWater(T, P=101.325):
# Implements the Krynicki formula; returns the self-diffusion of water (in micrometers^2/millisec)
# given the temperature T (in Centigrade) and the pressure P (in kPa).
d = 12.5 * exp(P * -5.22 * 1e-6) * sqrt(T + 273.15) * exp(-925 * exp(P * -2.6 * 1e-6)/(T + 273.15 - (95 + P * 2.61 * 1e-4)))
return d
The self-diffusion of water at body temperature and standard pressure, in micrometers2/millimeter, is:
In [ ]:
D = selfDiffusionOfWater(37)
print("%f micrometers^2/millisecond" % D)
Now we'll plot D for a biologically meaningful range of temperatures
In [ ]:
T = arange(25,41)
D = selfDiffusionOfWater(T)
figure()
plot(T, D, 'k')
xlabel('Temperature (Centigrade)', fontsize=14)
ylabel('Self-diffusion ($\mu$m$^2$/ms)', fontsize=14)
plot([37,37], [2,3.4], 'r-')
text(37, 3.45, 'Body Temperature', ha='center', color='r', fontsize=12)
a. The average atmospheric pressure in Denver Colorado is about 84 kPa. How different (in percent) is the self-diffusion coefficient of water at body temperature in Denver relative to that in Palo Alto, which is about at sea level?
b. Suppose you are running a fever of 40 deg Centigrade. Compared to someone without a fever, how much higher (in percent) is the water diffusion coefficient in your body?
In [ ]:
# compute your answer here
Here we simulate the brownian motion in a small chunk of tissue. First, we define some parameters, including the size of the simulated voxel (in micrometers), the time step (in milliseconds), and the Apparent Diffusion Coefficient (ADC) to simulate (in micrometers2/millisecond). Our tissue model will include simple barriers that will roughly approximate the effects of cell membranes, which are relatively impermeable to the free diffusion of water. So we will also define the barrier spacing (in micrometers). Finally, we'll specify the number of particles and time-steps to run.
In [ ]:
voxel_size = 50.0 # micrometers
ADC = 2.0 # micrometers^2/millisecond)
barrier_spacing = 10.0 # micrometers (set this to 0 for no barriers)
num_particles = 500
def draw_particles(ax, title, xy, particle_color, voxel_size, barrier_spacing):
ax.set_xlabel('X position $(\mu m)$')
ax.set_ylabel('Y position $(\mu m)$')
ax.axis('equal')
ax.set_title(title)
ax.set_xlim(-voxel_size/2, voxel_size/2)
ax.set_ylim(-voxel_size/2, voxel_size/2)
if barrier_spacing > 0:
compartments = unique(np.round(xy[1,:] / barrier_spacing))
for c in range(compartments.size):
ax.hlines(compartments[c]*barrier_spacing, -voxel_size/2, voxel_size/2, linewidth=4, colors=[.7, .7, .8], linestyles='solid')
particles = []
for p in range(xy.shape[1]):
particles.append(Circle(xy[:,p], 0.3, color=particle_color[p]))
ax.add_artist(particles[p])
return particles
# Place some particles randomly distributed in the volume
xy = random.rand(2, num_particles) * voxel_size - voxel_size/2.
start_xy = xy
particle_color = [((xy[0,p] + voxel_size/2) / voxel_size, (xy[1,p] + voxel_size/2) / voxel_size, .5) for p in range(num_particles)]
# draw 'em
fig,ax = subplots(1, 1, figsize=(6, 6))
particles = draw_particles(ax, 'initial particle positions', xy, particle_color, voxel_size, barrier_spacing)
In this loop, we update all the particle positions at each time step. The diffusion equation tells us that the final position of a particle moving in Brownian motion can be described by a Gaussian distribution with a standard deviation of sqrt(ADC*timeStep). So we update the current particle position by drawing numbers from a Gaussian with this standard deviation.
In [ ]:
import time, sys
from IPython.core.display import clear_output
time_step = 0.1 # milliseconds
nTimeSteps = 100
fig,ax = subplots(1, 1, figsize=(6, 6))
total_time = 0
for t_i in range(nTimeSteps):
dxy = np.random.standard_normal(xy.shape) * sqrt(2 * ADC * time_step)
new_xy = xy + dxy
if barrier_spacing>0:
curCompartment = np.round(xy[1,:]/barrier_spacing)
newCompartment = np.round(new_xy[1,:]/barrier_spacing)
for p in range(newCompartment.size):
if newCompartment[p]!=curCompartment[p]:
# approximate particles reflecting off the impermeable barrier
new_xy[1,p] = xy[1,p]
xy = new_xy
title = 'elapsed time: %5.2f ms' % (time_step * t_i)
particles = draw_particles(ax, title, xy, particle_color, voxel_size, barrier_spacing)
clear_output(wait=True)
display(fig,ax)
ax.cla()
total_time += time_step
close()
In [ ]:
# compute your answer here
By comparing the particle ending positions (in xy) with their starting positions (in start_xy), we can compute the diffusion tensor. This is essentially just a 2-d Gaussian fit to the position differences, computed using the covariance function (cov). We also need to normalize the positions by the total time that we diffused.
The eigensystem of the diffusion tensor (computed using 'eig') describes an isoprobability ellipse through the data points.
In [ ]:
Dt = cov(start_xy - xy) / (2 * total_time)
[val,vec] = eig(Dt)
estimatedADC = val / total_time
principalDiffusionDirection = vec[0]
print('ADC = ' + str(estimatedADC))
print('Principal Diffusion Direction (PDD) = ' + str(principalDiffusionDirection))
In [ ]:
# compute your answer here
Now lets show the particle starting positions a little line segment showing where each moved to.
In [ ]:
fig,ax = subplots(1, 1, figsize=(6, 6))
start_p = draw_particles(ax, 'initial particle positions', start_xy, particle_color, voxel_size, barrier_spacing)
for p in range(num_particles):
ax.plot((start_xy[0,p], xy[0,p]), (start_xy[1,p], xy[1,p]), linewidth=1, color=[.5, .5, .5], linestyle='solid')
a. Run the simulation with and without barriers by adjusting the 'barrierSpacing' variable. How does the diffusion tensor change?
b. Adjust the barrier spacing. What effect does this have on the princpal diffusion direction? On the estimatedADC values?
c. With barriers in place, reduce the number of time steps (nTimeSteps). How does this affect the estimated ADC values? Explore the interaction between the barrier spacing and the number of time steps.
In [ ]:
# compute your answer here
We'll simulate an image that represents a vial of water in a 3 Tesla magnetic field. The image intensity at each point will represent the local magnetic field strength, expressed as the Larmor frequency difference between that region and the frequency at 3 T.
First let's define some parameters, such as the simulated filed strength (B0) and the gyromagnetic ratio for Hydrogen:
In [ ]:
B0 = 3.0 # Magnetic field strength (Tesla)
gyromagneticRatio = 42.58e+6 # Gyromagnetic constant for hydrogen (Hz / Tesla)
# The Larmor frequency (in Hz) of Hydrogen spins in this magnet is:
spinFreq = gyromagneticRatio * B0
In [ ]:
# compute your answer here
We start by defining the size of our simulated voxel, in micrometers. For the diffusion gradients, we'll start with the Y gradient turned off and the X gradient set to 5e-8 Tesla per micrometer (that's 50 mT/m, a typical gradient strengh for clinical MR scanners). We'll also quantize space into 100 regions and use meshgrid to lay these out into 2d arrays that be used to compute a 100x100 image. Finally, we compute the relative field strength at each spatial location across our simulated voxel. Gradient strengths are symmetric about the center. To understand the following expression, work through the units: (micrometers T/um + micrometers T/um) leaves us with T. We scale this by 1e6 to express the resulting image in micro-Teslas
In [ ]:
voxelSize = 100.0 # micrometers
gx = 5e-8 # Tesla / micrometer
gy = 0.0 # Tesla / micrometer
def draw_spins(ax, title, field_image, im_unit, sx, sy, px, py):
# a function to draw spin-packets
# draw the relative magnetic field map image
im = ax.imshow(field_image, extent=im_unit+im_unit, cmap=matplotlib.cm.bone)
ax.set_ylabel('Y position (micrometers)')
ax.set_xlabel('X position (micrometers)')
ax.set_title(title)
# Place some spin packets in there:
ax.scatter(x=sx+px, y=sy+py, color='r', s=3)
ax.scatter(x=sx, y=sy, color='c', s=3)
ax.set_xlim(im_unit)
ax.set_ylim(im_unit)
# add a colorbar
divider = make_axes_locatable(ax)
cax = divider.append_axes("bottom", size="7%", pad=0.5)
cbl = fig.colorbar(im, cax=cax, orientation='horizontal')
cbl.set_label('Relative field strength (micro Tesla)')
im_unit = (-voxelSize/2, voxelSize/2)
x = linspace(im_unit[0], im_unit[1], 100)
y = linspace(im_unit[0], im_unit[1], 100)
[X,Y] = meshgrid(x,y)
# micrometers * Tesla/micrometer * 1e6 = microTesla
relative_field_image = (X*gx + Y*gy) * 1e6
locations = linspace(-voxelSize/2+voxelSize/10, voxelSize/2-voxelSize/10, 5)
sx,sy = meshgrid(locations, locations);
sx = sx.flatten()
sy = sy.flatten()
# set the phase/magnitude to be zero
px = zeros(sx.shape)
py = zeros(sy.shape)
fig,ax = subplots(1, 1)
draw_spins(ax, 'Spin packets at rest in a gradient', relative_field_image, im_unit, sx, sy, px, py)
Calculate the relative spin frequency at each spin location. Our gradient strengths are expressed as T/cm and are symmetric about the center of the voxelSize. To understand the following expression, work through the units: (centimeters T/cm + centimeters T/cm) Hz/Tesla + TeslaHz/Tesla leaves us with MHz
In [ ]:
spinFreq = (sx * gx + sy * gy) * gyromagneticRatio + B0 * gyromagneticRatio
print(spinFreq)
Note that including the B0 field in this equation simply adds an offset to the spin frequency. For most purposes, we usually only care about the spin frequency relative to the B0 frequency (i.e., the rotating frame of reference), so we can leave that last term off and compute relative frequencies (in Hz):
In [ ]:
relativeSpinFreq = (sx * gx + sy * gy) * gyromagneticRatio
print(relativeSpinFreq)
In [ ]:
# compute your answer here
When we first apply an RF pulse, all the spins will all precess in phase. If they are all experienceing the same magnetic field, they will remain in phase. However, if some spins experience a different local field, they will become out of phase with the others. Let's show this with a movie, where the phase will be represented with color. Our timestep is 1 millisecond.
In [ ]:
fig,ax = subplots(1, 1)
timeStep = .001
# Initialize the transverse magnitization to reflect a 90 deg RF pulse.
# The scale here is arbitrary and is selected simply to make a nice plot.
Mxy0 = 5
# Set the T2 value of the spins (in seconds)
T2 = 0.07
curPhase = zeros(sx.size)
t = 0.
nTimeSteps = 50
for ti in range(nTimeSteps):
# Update the current phase based on the spin's precession rate, which is a function
# of its local magnetic field.
curPhase = curPhase + 2*pi*gyromagneticRatio * (sx*gx+sy*gy) * timeStep
# Do a 180 flip at the TE/2:
if ti==round(nTimeSteps/2.):
curPhase = -curPhase
# The transverse magnitization magnitude decays with the T2:
curMagnitude = Mxy0 * exp(-t/T2)
px = sin(curPhase) * curMagnitude
py = cos(curPhase) * curMagnitude
# Summarize the total (relative) MR signal for this iteration
S = sqrt(sum(px**2 + py**2)) / sx.size
title = 'elapsed time: %5.1f/%5.1f ms' % (t*1000., timeStep*(nTimeSteps-1)*1000)
draw_spins(ax, title, relative_field_image, im_unit, sx, sy, px, py)
clear_output(wait=True)
display(fig,ax)
ax.cla()
t = t+timeStep
close()
Developer notes: