Class: Psych 204a
Tutorial: Basic MR
Duration: 90 minutes
Authors: Originally written in Matlab by Brian Wandell in 2004
Checked and updated:
Copyright: Stanford University
This tutorial explains basic principles of magnetic resonance signals. As the first tutorial in the series, it also gives you an opportunity to examine basic Python commands. (For a short tutorial on the Jupyter/IPython notebook interface, have a look at this video: https://www.youtube.com/watch?v=lmoNmY-cmSI.)
In this tutorial, principles are illustrated for signals derived from a bulk substance, such as a beaker of water. Various terms used to describe the MR principles, including spins, parallel and anti-parallel, Boltzman distribution are introduced.
Also, tissue properties such as T1 (longitudinal) and T2 (spin-spin) interactions, are explained and the way in which signal contrast depends on these parameters is explained.
The next tutorial, MRImaging, takes up the topic of how to make images that measure these tissue properties in a non-uniform volume, such as a head.
References to help with this tutorial:
1st v. 2nd edition of text: The course text is Huettel et al, 2nd edition. The first and second editions are quite similar, especially in the earlier chapters. References followed by "ii" are 2nd edition, "i" first edition. Hence p. 51i, 60ii means p 51 in the first edition, p 60 in the second. And so on.
OK, let's get started with the tutorial! Each cell that has a gray box with 'In [ ]:' in front of it contains code. To run the code, click on the cell (it will be highlighted with a green outline when you do) and then hit the Enter key while holding shift. Some later cells depend on values set in earlier cells, so be sure to execute them in order. In the first cell of code just below this paragraph, we set some general parameters for python.
In [ ]:
%pylab inline
import matplotlib as mpl
mpl.rcParams["figure.figsize"] = (8, 6)
mpl.rcParams["axes.grid"] = True
from IPython.display import display
from ipywidgets import interact,FloatSlider
When an atomic species with net spin is placed in a steady magnetic field, it precesses at a frequency that is characteristic of that species and that varies with the strength of the magnetic field. We can compute the precession frequency from a following simple formula. Note the units. (see chapter 3 and ca. p. 51i, 60ii for a discussion of precession).
In [ ]:
B0 = 1.5 # Magnetic field strength (Tesla)
g = 42.58e6 # Gyromagnetic constant for hyrdogen (Hz / Tesla)
v = g * B0 # The resonant frequence of hydrogen, also called its Larmor frequency
What are the units of the Larmor frequency, $v$, for hydrogen as expressed above?
In [ ]:
# Python string interpolation is accomplished through the % operator
print "The resonant frequency of spins in hydrogen is %0.4f (MHz) at %.2f Tesla" % (v / (10 ** 6), B0)
The gyromagnetic constant of sodium is $11.27 \times 10^6 Hz/Tesla$. Compute the Larmor frequency of sodium in a 3T magnet.
Ordinarily, the resonant frequencies are expressed in units of MegaHertz (millions of hertz)
In [ ]:
# Compute your answer here
The energy in precessing spin is proportional to its resonant frequency, $v$. The constant of proportionality between energy and frequency is a famous constant called Planck's constant. (Extra credit: Find out something interesting about Max Planck).
Hence, the amount of energy in a hydrogen molecular in this magnetic field is
In [ ]:
h = 6.626e-34 # Planck's constant (Joules-seconds)
E = h * v
print "E = %.4g" % E
What are the units of E?
At steady state, say when the subject enters the magnet and prior to any measurements, the dipoles align parallel or anti-parallel to the magnetic field. (The textbook illustrates this with a nice metaphor of gravity in Figure 3.7ii, 3,6i.)
The energy difference between these two states is proportional to the mean magnetic field. This is called the Zeeman effect (see Figure 3.8i, 3.9ii in the book).
The proportion of dipoles aligned parallel (low energy) and anti-parallel (high energy) to the main field is described by the Boltzmann distribution. The formula that determines the fraction of dipoles in the low and high energy states is in Hornak (Chapter 2).
In [ ]:
k = 1.3805e-23 # Boltzmann's constant, J/Kelvin
T = 300 # Degrees Kelvin at room temperature
dE = h * g * B0 # Transition energy
ratio_high_to_low = exp(-dE / (k * T)) # Boltzmann's formula, Hornak, Chapter 3; Huettel, p. 76ii
print "High-low ratio: %.6f" % ratio_high_to_low
At low temperatures (near absolute zero), very few of the dipoles enter the high (anti-parallel) energy state.
No, this is not important for us. But I find the numbers interesting and can see why people might want to work on low temperature physics for a while.
In [ ]:
T = logspace(-3, 2.5, 50)
r = exp(-dE / (k*T))
plot(T, r)
semilogx()
xlabel('Temperature (K)')
ylabel('Ratio of high/low energy state dipoles')
ylim([0, 1.1]);
Where is human body temperature on the graph? Given human body temperature, what is the ratio of high/low energy?
Would it matter if we kept the room cooler?
The tendency of dipoles to align (either parallel or anti-parallel) with the B0 magnetic field imposes an opportunity to probe tissue properties using magnetic resonance. MR signals are derived by perturbing the ordered spins (excitation) and then measuring the emissions (reception) as the dipoles return to their low energy state. The way the dipoles return to their low energy state provides information about the local tissue.
Summing the effect of the many many dipoles within a voxel, we obtain a measure of the dipoles within the voxel called the net magnetization. This is represented by a single vector (see the many examples in Hornak). Most reasoning about the MR signal is based on understanding models of the net magnetization vector and how it recovers after being perturbed by the radio frequency pulses in the presence of changing magnetic fields.
The MR measurements that we obtain describe the net magnetization in the direction perpendicular to the main axis (B0 field). This direction is illustrated in the following figures.
First, here is a 3D plot that shows the net magnetization as a red circle in the steady-state. The black lines show the three axes.
In [ ]:
# This function just saves us some typing below
# It's not relevant to the conceptual material
# But it does show a nice trick on the last line
def axisplot(ax):
"""Convenience function to plot axis lines in a 3D plot"""
ax.plot([-1, 1], [0, 0], [0, 0], "k:")
ax.plot([0, 0], [-1, 1], [0, 0], "k:")
ax.plot([0, 0], [0, 0], [-1, 1], "k:")
for axis in ["x", "y", "z"]:
getattr(ax, "set_%slabel" % axis)(axis)
In [ ]:
from mpl_toolkits.mplot3d import Axes3D
f = figure()
ax = f.add_subplot(111, projection='3d', aspect="equal")
ax.plot([0], [0], [1], "ro")
axisplot(ax)
az0 = 332.5
el0 = 30
ax.view_init(el0, az0)
The size of the MR signal we measure depends on how far the magnetization deviates from the main z-axis. We can see this more easily by looking at the figure straight from the top. In initial situation, the point is at $(0, 0)$ in the $(x, y)$ plane.
In [ ]:
az1 = 0
el1 = 90
ax.view_init(el1, az1)
display(f)
Notice that from this view, looking down the z-axis, we see only the x- and y-axes. The net magnetization is at $(0,0)$ so we will not see any signal.
Suppose we excite the tissue and place the net magnetization along the x-axis.
In [ ]:
f = figure()
ax = f.add_subplot(111, projection='3d', aspect="equal")
ax.plot([1], [0], [0], "go")
axisplot(ax)
ax.view_init(el0, az0)
When we look from the top now, we see that there is a large magnetization component. The green point is removed from $(0,0)$ and falls along the x-axis.
In [ ]:
ax.view_init(el1, az1)
display(f)
When we change the net magnetization from the steady-state position (red circle) to the excited position (green circle), it is like introducing a 90 deg rotation in the magnetization direction. This is usually called the flip angle. This is one of the parameters that you select when doing MR imaging.
In [ ]:
ax = subplot(111, projection='3d', aspect="equal")
ax.plot([0], [0], [1], "ro")
ax.plot([1], [0], [0], "go")
axisplot(ax)
ax.view_init(el0, az0)
As the net magnetization returns to the steady-state position, the distance from the origin decreases, reducing the signal. This is illustrated by the collection of points plotted here that show the net magnetization rotating back towards the z-axis.
In [ ]:
theta = ((pi / 2) * (linspace(0, 10, 11) / 10)).reshape(11, 1)
x = cos(theta)
y = zeros_like(x)
z = sin(theta)
f = figure()
ax = f.add_subplot(111, projection='3d', aspect="equal")
clut = linspace(1, 0, 11)
for i, (x_i, y_i, z_i) in enumerate(zip(x, y, z)):
ax.plot(x_i, y_i, z_i, "o", color=cm.RdYlGn(clut[i]))
axisplot(ax)
ax.view_init(el0, az0)
When viewed from the top, you can see that the green points head back towards the origin.
In [ ]:
ax.view_init(el1, az1)
display(f)
After receiving an excitation pulse, this component of the MR signal decays gradually over time. The spin-lattice decay has been measured and, in general, it follows an exponential decay rate.
Specifically, here is the rate of recovery of the T1 magnetization for hydrogen molecules in gray matter.
In [ ]:
T1_gray = 0.88 # Time constant units of S
t_T1 = arange(0.02, 6, 0.02) # Time in seconds
Mo = 1 # Set the net magnetization in the steady state to 1 and ignore.
This is the exponential rate of recovery of the T1 magnetization, after it has been set to zero by flipping the net magnetization 90 degrees.
In [ ]:
MzG_T1 = Mo * (1 - exp(-t_T1 / T1_gray))
Plotted is a graph we have the magnetization of gray matter as
In [ ]:
plot(t_T1, MzG_T1)
xlabel('Time (s)')
ylabel('Longitudinal magnetization (T1)');
The decay rates for various brain tissues (summarized by the parameter T1 above) differs both with the material and with the level of the B0 field. The value T1 = 0.88 seconds above is typical of gray matter at 1.5T.
The T1 value for white matter is slightly smaller. Comparing the two we see that white matter recovers slightly faster (has a smaller T1):
In [ ]:
T1_white = 0.64;
MzW_T1 = Mo * (1 - exp(-t_T1 / T1_white))
plot(t_T1, MzG_T1, 'black', label="Gray")
plot(t_T1, MzW_T1, 'gray', linestyle="--", label="White")
xlabel('Time (s)')
ylabel('Longitudinal magnetization (T1)')
legend(loc="best");
Notice that the time to recover all the way back to Mo is fairly long, on the order of 3-4 sec. This is a limitation in the speed of acquisition for T1 imaging.
The difference in the T1 component of the signal from gray matter and white matter changes over time. This difference is plotted in the next graph.
In [ ]:
plot(t_T1, abs(MzW_T1 - MzG_T1))
xlabel('Time (s)')
ylabel('Magnetization difference');
If you are seeking to make a measurement that optimizes the signal to noise ratio between these two materials, at what time would you measure the recovery of the T1 signal?
Look up the T1 value of cerebro-spinal fluid (CSF). Plot the T1 recovery of CSF. At what time you would measure to maximize the white/CSF contrast.
We can visualize this difference as follows. Suppose that we have two beakers, adjacent to one another, containing materials with different T1 values. Suppose we make a pair of images in which the intensity of each image is set to the T1 value over time. What would the T1 images look like?
The beakers start with the same, Mo, magnetization.
In [ ]:
beaker1 = Mo * ones((32, 32))
beaker2 = Mo * ones((32, 32))
They will have different T1 relaxation values
In [ ]:
T1 = (0.64, 0.88) # White, gray T1
Here is a movie showing the relative intensities of the images over a 4 sec period, measured every 100 ms. You can play the movie by moving the slider to change the time point displayed.
In [ ]:
f = figure()
beakers = [beaker1, beaker2]
def draw_beakers(time):
for i, beaker in enumerate(beakers):
subplot(1, 2, i + 1)
img = beaker * (1 - exp(-time / T1[i]))
imshow(img, vmin=0, vmax=1, cmap="gray")
grid(False)
xticks([])
yticks([])
ylim([-5, 32])
title("Beaker %d" % (i + 1))
text(8, -3, 'Time: %.2f sec' % time, fontdict=dict(size=12))
interact(draw_beakers, time=FloatSlider(min=0, max=4.0, value=0));
As you can see, if we make a picture of the net magnetization around 0.6-1.0 sec during the decay, there will be a good contrast difference between the gray and white matter. Measured earlier or later, the picture will have less contrast.
As a preview of further work, later in the course, we should add just a little noise to the measurements. After all, all measurements have some noise. Let's look again.
In [ ]:
def draw_beakers(time, noise_level):
for i, beaker in enumerate(beakers):
subplot(1, 2, i + 1)
img = beaker * (1 - exp(-time / T1[i])) + randn(*beaker.shape)*noise_level
img = abs(img)
imshow(img, vmin=0, vmax=1, cmap="gray")
grid(False)
xticks([])
yticks([])
ylim([-5, 32])
title("Beaker %d" % (i + 1))
text(8, -3, 'Time: %.2f sec' % time, fontdict=dict(size=12))
interact(draw_beakers, time=FloatSlider(min=0, max=4.0, value=0), noise_level=FloatSlider(min=0, max=1.0, value=0.1));
There is a second physical mechanism, in addition to the spin-lattice measurement, that influences the MR signal. This second mechanism is called spin-spin interaction (transverse relaxation). This signaling mechanism is particularly important for functional magnetic resonance imaging and the BOLD signal.
In describing the T1 signal, we treated the MR signal as a single unified vector. In the example above, we explored what happens to the net magnetization of the MR signal when the vector is rotated 90 deg into the x-y plane.
But we omitted any discussion the fact that the dipoles are assumed to be continuously precessing around the main axis together, in unison. Perhaps in a perfectly homogeneous environment, these rotating dipoles would precess at the Larmor frequency in perfect synchronization and we could treat the single large vector as we have.
But in practice, the dipoles within a single voxel of a functional image each experience slightly different magnetic field environments. Consequently, they each have their own individual Larmor frequencies, proportional to the magnetic field that they experience. An important second mechanism of MR is a consequence of the fact that the individual dipoles each have their own local magnetic field and the synchrony soon dissipate.
Suppose we have a large sample of dipoles that are spinning together in perfect phase. We can specify their orientation as an angle in this plane, theta. Let's assume they all share a common angle
In [ ]:
n_samples = 10000
theta = zeros(n_samples)
The position of the spins in the $(x, y)$ plane will be
In [ ]:
spins = [cos(theta), sin(theta)]
And they will all fall at the same position
In [ ]:
subplot(111, aspect="equal")
plot(spins[0], spins[1], "o")
xlim(-2, 2)
ylim(-2, 2)
xlabel("x")
ylabel("y");
The total magnetization, summed across all the dipoles, is the vector length of the sum of these spins
In [ ]:
avg_pos = sum(spins, axis=1) / n_samples
net_mag = sqrt(sum(square(avg_pos)))
print "Net magnetization: %.4f" % net_mag
Now, suppose that spins are precessing at slightly different rates. So after a few moments in time they do not fall at exactly the same angle. We can express this by creating a new vector theta that has some variability in it.
In [ ]:
theta = rand(n_samples) * 0.5 # Uniform random number generator
Here is the distribution of the angles
In [ ]:
hist(theta)
xlabel('Angle')
ylabel('Number of spins')
xlim(0, 2 * pi)
xticks([0, pi / 2, pi, 3 * pi / 2, 2 * pi],
["0", r"$\frac{\pi}{2}$", "$\pi$", r"$\frac{3\pi}{2}$", "$2\pi$"],
size=14);
Going through the same process we can make a plot of the spin positions
In [ ]:
spins = [cos(theta), sin(theta)]
subplot(111, aspect="equal")
plot(spins[0], spins[1], 'o')
xlim(-2, 2)
ylim(-2, 2)
xlabel('x'),
ylabel('y');
Now, you can see that the net magnetization is somewhat smaller
In [ ]:
avg_pos = sum(spins, axis=1) / n_samples
net_mag = sqrt(sum(square(avg_pos)))
print "Net magnetization: %.4f" % net_mag
If the spins grow even further out of phase, say they are spread over a full $\pi$ radians (180 degrees), then we have a dramatic reduction in the net magnetization
In [ ]:
theta = rand(n_samples) * pi
hist(theta)
xlabel('Angle')
ylabel('Number of spins')
xlim(0, 2 * pi)
xticks([0, pi / 2, pi, 3 * pi / 2, 2 * pi],
["0", r"$\frac{\pi}{2}$", "$\pi$", r"$\frac{3\pi}{2}$", "$2\pi$"],
size=14);
In [ ]:
spins = [cos(theta), sin(theta)]
subplot(111, aspect="equal")
plot(spins[0], spins[1], 'o')
xlim(-2, 2)
ylim(-2, 2)
xlabel('x')
ylabel('y')
avg_pos = sum(spins, axis=1) / n_samples
net_mag = sqrt(sum(square(avg_pos)))
print "Net magnetization: %.4f" % net_mag
In a typical experiment, in which the spins are in an inhomogeneous environment, the spins spread out more and more with time, and the transverse magnetization declines. The loss of signal from this spin-dephasing mechanism follows an exponential time constant.
In [ ]:
t_T2 = arange(0.01, 0.3, 0.01) # Time in secs
T2_white = 0.08 # T2 for white matter
T2_gray = 0.11 # T2 for gray matter
MzG_T2 = Mo * exp(-t_T2 / T2_gray)
MzW_T2 = Mo * exp(-t_T2 / T2_white)
plot(t_T2, MzG_T2, 'k', label="Gray")
plot(t_T2, MzW_T2, 'gray', linestyle="--", label="White")
xlabel('Time (s)')
ylabel('Transverse magnetization (T2)')
legend();
Experimental measurements of spin-spin decay shows that it occurs at a much faster rate than spin-lattice. Comparison of T1 (spin-lattice) and T2 (spin-spin) decay constants at various B0 field strengths are:
| T1 | T2 | |||||
| Field | 1.5T | 3.0T | 4.0T | 1.5T | 3.0T | 4.0T |
| White | 0.64 | 0.86 | 1.04 | 0.08 | 0.08 | 0.05 |
| Gray | 0.88 | 1.20 | 1.40 | 0.08 | 0.11 | 0.05 |
Source: Jezzard and Clare, Chapter 3 in the Oxford fMRI Book
Also, notice that the peak difference occurs a very short time compared to the T1 difference.
In [ ]:
plot(t_T1, abs(MzG_T1 - MzW_T1), label="T1")
plot(t_T2, abs(MzG_T2 - MzW_T2), label="T2")
xlabel('Time (s)')
ylabel('Transverse magnetization difference')
legend();
Suppose you wanted to measure brain structure, and you were particularly interested in gray/white differences. Based on the T1 and T2 curves we have drawn, would you choose to distinguish these two tissues using T1 or T2?
The T2 difference is so short that measurement of the T2 signal would be a problem were it not for the invention of an important measurement method called 'Spin Echo'. We will begin the next tutorial by explaining that idea.
In fact, there are many essential imaging ideas that we have not yet explored. For example, how can we measure T1 and T2 separately? Perhaps most importantly, how can we form an image? These are the topics we will take up in the next tutorial.