Diffusion

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

Self-diffusion of water

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:

  • Krynicki, Green & Sawyer (1978) Pressure and temperature dependence of self-diffusion in water. Faraday Discuss. Chem. Soc., 66, 199 - 208.
  • Mills (1973) Self-Diffusion in Normal and Heavy Water. JPhysChem 77(5), pg. 685 - 688.
  • Also see http://www.lsbu.ac.uk/water/explan5.html.

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)

Question 1

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

Brownian motion

Set up the diffusion simulation

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)

Run the diffusion simulation

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()

Question 2

a. What is the average position change of each particle in the X dimension? In Y? (Hint: start_xy contains the starting positions.)

b. What is the average distance that each particle moved? (Hint: compute the Euclidean distance that each moved.)


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

Question 3

a. What are the units of the ADC?

b. What are the units of the PDD?


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

Question 4

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

The effect of diffusion on the MR signal

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

Question 4

a. What is the Larmor frequency of Hydrogen spins at 3T?

b. What is it at 7T?


In [ ]:
# compute your answer here

Simulate spins in an MR experiment

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)

Question 5

a. Speculate on why the relative spin frequency is the most important value to calculate here.

b. Do you think the B0 field strength will play a role in the calculation of the diffusion tensor?


In [ ]:
# compute your answer here

Display the relative frequencies (in Hz)

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: