<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Diffusion
Class
Step1:  Self-diffusion of water
The self-diffusion coefficient of water (in micrometers<sup>2</sup>/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
Step2:  The self-diffusion of water at body temperature and standard pressure, in micrometers<sup>2</sup>/millimeter, is
Step3:  Now we'll plot D for a biologically meaningful range of temperatures
Step4:  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?
Step5:  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 micrometers<sup>2</sup>/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.
Step6:  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.
Step7:  Question 2
a. What is the average position change of each particle in the X dimension? In Y? (Hint
Step8:  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.
Step9:  Question 3
a. What are the units of the ADC?
b. What are the units of the PDD?
Step10:  Now lets show the particle starting positions a little line segment showing where each moved to.
Step11:  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.
Step12:  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
Step13:  Question 4
a. What is the Larmor frequency of Hydrogen spins at 3T? 
b. What is it at 7T?
Step14:  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
Step15:  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
Step16:  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)
Step17:  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?
Step18:  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.
<ASSISTANT_TASK:>
Python Code:
%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
Explanation: 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.
End of explanation
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
Explanation: Self-diffusion of water
The self-diffusion coefficient of water (in micrometers<sup>2</sup>/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).
End of explanation
D = selfDiffusionOfWater(37)
print("%f micrometers^2/millisecond" % D)
Explanation: The self-diffusion of water at body temperature and standard pressure, in micrometers<sup>2</sup>/millimeter, is:
End of explanation
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)
Explanation: Now we'll plot D for a biologically meaningful range of temperatures
End of explanation
# compute your answer here
Explanation: 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?
End of explanation
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)
Explanation: 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 micrometers<sup>2</sup>/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.
End of explanation
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()
Explanation: 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.
End of explanation
# compute your answer here
Explanation: 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.)
End of explanation
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))
Explanation: 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.
End of explanation
# compute your answer here
Explanation: Question 3
a. What are the units of the ADC?
b. What are the units of the PDD?
End of explanation
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')
Explanation: Now lets show the particle starting positions a little line segment showing where each moved to.
End of explanation
# compute your answer here
Explanation: 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.
End of explanation
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
Explanation: 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:
End of explanation
# compute your answer here
Explanation: Question 4
a. What is the Larmor frequency of Hydrogen spins at 3T? 
b. What is it at 7T?
End of explanation
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)
Explanation: 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
End of explanation
spinFreq = (sx * gx + sy * gy) * gyromagneticRatio + B0 * gyromagneticRatio
print(spinFreq)
Explanation: 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 + Tesla*Hz/Tesla leaves us with MHz
End of explanation
relativeSpinFreq = (sx * gx + sy * gy) * gyromagneticRatio
print(relativeSpinFreq)
Explanation: 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):
End of explanation
# compute your answer here
Explanation: 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?
End of explanation
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()
Explanation: 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.
End of explanation
<END_TASK>