Friday, June 28, 2019

The Innovators

The Innovators: How a Group of Hackers, Geniuses, and Geeks Created the Digital Revolution, by Walter Isaacson, superimposed on Intermediate Physics for Medicine and Biology.
The Innovators: How a Group
of Hackers, Geniuses, and Geeks
Created the Digital Revolution
,
by Walter Isaacson.
Computers play a large and growing role in biomedical research. Intermediate Physics for Medicine and Biology contains four computer programs that illustrate the computations behind topics such as the Hodgkin-Huxley model and computed tomography. To learn more about the history of computers, I read The Innovators: How a Group of Hackers, Geniuses, and Geeks Created the Digital Revolution, by Walter Isaacson. The introduction begins
The computer and the Internet are among the most important inventions of our era, but few people know who created them. They were not conjured up in a garret or garage by solo inventors suitable to be singled out on magazine covers or put into a pantheon with Edison, Bell, and Morse. Instead, most of the innovations of the digital age were done collaboratively. There were a lot of fascinating people involved, some ingenious and a few even geniuses. This is the story of these pioneers, hackers, inventors, and entrepreneurs—who they were, how their minds worked, and what made them so creative. It’s also a narrative of how they collaborated and why their ability to work as teams made them even more creative.
Mr. Ray Dennis, my FORTRAN teacher at Shawnee Mission South High School.
Mr. Ray Dennis,
my FORTRAN teacher at
Shawnee Mission South High School.
From the 1977-1978 school yearbook.
Transistors, microchips and computers were invented before I was born, but I witnessed the rise of personal computers and the creation of the Internet in the 70s, 80s. and 90s. My first experience with computers was during my senior year of high school (1978), when I took a programming class based on the computer language FORTRAN. It was the most useful class I ever took; I programmed in FORTRAN almost every day for decades. At that time, many computers used punch cards to input programs. I remember typing up a stack of cards during class and giving them to my teacher, Ray Dennis, who ran them on an off-site computer after school. The next day he would return our cards along with our output on folded, perforated paper. Mistype one comma and you lost a day. I learned to program with care.

The next year I attended the University of Kansas, where they were phasing out punch cards and replacing them with time-sharing terminals. I thought it was the greatest advance ever. During my freshmen year, my roommate was an engineering graduate student from California who owned an Apple II, and he let me use it to play a primitive Star Trek video game. The Physics Department purchased a Vax mainframe computer that I used for undergraduate research analyzing light scattering data.

At Vanderbilt, my advisor John Wikswo made sure each of his graduate students had their own computer (an IBM clone), which I thought was amazing. I also got my first email account. During my grad school years I followed the rivalry between Bill Gates and Steve Jobs. After starting as a PC guy,  I moved to the National Institutes of Health in 1988 and switched over to an Apple Macintosh, and have been a loyal Mac man ever since. I remember going to a lecture at NIH about something called the World Wide Web, and thinking that it sounded silly. Since then, I’ve seen the rise of Yahoo!, Google, and Wikipedia, and I write these posts ever week using blogger. All these developments and more are described in The Innovators

I’m a fan of Walter Isaacson. He’s written several “genius biographies (Steve Jobs, Leonardo da Vinci, Albert Einstein, Ben Franklin), but my favorite of his books is The Wise Men, the story of six influential leaders—Averell Harriman, Dean Acheson, George Kennan, Robert Lovett, John McCloy, and Charles Bohlen—during the cold war. The Innovators reminds me of The Wise Men in that it focuses on the interactions among a group, rather than the contributions of an individual.

For readers of IPMB who use computers for biomedical research, I recommend The Innovators. It provides insight into how the digital world was invented, and the role of collaboration in science and technology. Enjoy!

Friday, June 21, 2019

More About Implantable Microcoils for Intracortical Magnetic Stimulation

Three years ago, I wrote a post in this blog discussing an article by Seung Woo Lee and his coworkers (“Implantable Microcoils for Intracortical Magnetic Stimulation,” Science Advances, 2:e1600889, 2016). They claimed to have performed magnetic stimulation of nerves by passing a 40 mA, 3 kHz current through a single-turn microcoil with a size less than a millimeter. I claimed that the electric field induced in the surrounding tissue by such a coil would be much smaller than Lee et al. predicted. In their Figure 2 they calculated that a 1 mA current induced an electric field on the order of 2 V/m. I calculated an electric field about a million times smaller, and concluded “their results are too strange to believe and too important to ignore.”

I didn’t ignore them. Recently a graduate student here at Oakland University, Mohammed Alzahrani, and I tested the hypothesis that excitation using microcoils is caused by capacitive coupling rather than magnetic stimulation. The picture below shows our model. The current at the left end of the microcoil passes through the capacitance of the insulation and enters the surrounding tissue. It then flows through the tissue, possibly exciting neurons along its path, until reentering the wire through the capacitance near the right end.


Does this model look familiar? It’s similar to the cable model for a nerve axon (for more about the cable model, see Section 6.11 of Intermediate Physics for Medicine and Biology). The wire in our model is analogous to the intracellular space of the axon in the traditional cable model, and the insulation surrounding the wire is analogous to the cell membrane. Our model’s even simpler than the traditional cable model because the conductance of the insulation is so low that it can be taken as zero; the only way for current to leave the wire is through the capacitance. This model is not new; it was derived in the 19th century to describe current through the transatlantic telegraph cable.

Our goal was to calculate the electric field assuming capacitive coupling, to see whether it’s larger or smaller than what you’d expect from magnetic stimulation. We concluded
In summary, we predict an electric field in the tissue due to capacitive coupling of about 4 mV/m for a current of 1 mA and 3 kHz. The electric field produced by magnetic stimulation would be thousands of times less, on the order of 0.002 mV/m. Therefore, capacitive coupling should be the dominant mechanism for stimulation with a microcoil.
We haven’t published our results yet, but you can download a preprint on my ResearchGate page (doi 10.13140/RG.2.2.19222.40006). I’m curious what you think.

What’s the moral to this story? As I wrote at the end of my previous post, experiments “need to be consistent with the fundamental physical laws outlined in Intermediate Physics for Medicine and Biology.”

Friday, June 14, 2019

Isotopes to Worry About: Cesium-137, Iodine-131, and Strontium-90

The radiation hazard symbol.
Three hazardous isotopes released by a nuclear explosion or accident are cesium-137, iodine-131, and strontium-90.

Cesium-137

Isotopes with short half lives often have disappeared within days of a nuclear accident, after they emit lots of radiation. Isotopes with long half lives may linger for millennia, but aren’t very radioactive. Isotopes with half lives of a few decades are “just right” for being an environmental hazard; their lifetimes are short enough that they release a lot of radioactivity, but are long enough that they cause decades of danger. Cesium-137 (137Cs) has a half-life of 30 years. It is the main source of remaining radiation at the site of the Chernobyl accident.

Cesium-137 is volatile, meaning it evaporates at high temperatures, allowing it to mix with the air and spread with the wind. In addition, cesium is in the same column of the periodic table as sodium and potassium, so it forms water soluble salts that distribute throughout the body. It beta decays (0.51 MeV) to a meta-stable state of barium-137, which then decays by emission of a gamma ray (0.66 MeV).

Accidental uptake of caesium-137 can be treated with Prussian blue, which binds to it and reduces its biological half-life from 70 to 30 days.

Iodine-131

Iodine-131 (131I) has a half-life of eight days, so it is dangerous for only a few weeks after a nuclear explosion or accident. However, radioactive iodine is concentrated in the thyroid gland, causing thyroid cancer. Iodine-131 undergoes beta decay (0.61 MeV) to xenon-131, which then emits a 0.36 MeV gamma ray. It’s so potent a radioisotope that it’s used for cancer therapy (see Section 17.11 of Intermediate Physics for Medicine and Biology).

After a nuclear accident, people take potassium iodide pills to flood the thyroid gland with non-radioactive iodine, suppressing the uptake of the radioactive isotope.

Strontium-90

Like cesium-137, strontium-90 (90Sr) has a half life of about 30 years. It undergoes beta decay (0.55 MeV) to yttrium-90, which in turn beta decays (2.3 MeV) to stable zirconium-90. Strontium-90 is in the same column of the periodic table as calcium, so it is taken up by bones (it’s a “bone seeker”) and therefore has a long biological half-life (about 18 years).

Dirty Bombs

Want something to worry about? Consider this radiological nightmare: a terrorist dirty bomb consisting of equal parts cesium-137, iodine-131, and strontium-90. Good luck sleeping tonight.

Friday, June 7, 2019

Slice Selection During Magnetic Resonance Imaging

In Chapter 18 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I discuss magnetic resonance imaging. The first step in most MRI pulse sequences is to excite the spins in a slice of the sample. Three magnetic fields are required: A static field B0 in the z direction that aligns the spins, a magnetic field gradient G in the z direction that selects a slice with thickness Δz, and a radiofrequency field in the x direction Bx(t) that excites the spins.

The three magnetic fields applied during slice selection for magnetic resonance imaging: A static field, a gradient field, and a radiofrequency field.

The RF signal is called a π/2 pulse because it rotates the spins through an angle of 90 degrees (π/2 radians), from the z direction to the x-y plane. The Larmor frequency (the angular frequency at which the spins precess at the center of the slice, z = 0) is ω0 = γB0, where γ is the gyromagnetic ratio. Figure 18.23a of IPMB shows the RF signal as a function of time t: an oscillation at the Larmor frequency modulated by the sinc function.

Figure 18.23a of Intermediate Physics for Medicine and Biology, showing the radiofrequency signal used as a π/2 pulse during slice selection.
Figure 18.23a of IPMB, showing the radiofrequency signal
used as a π/2 pulse during slice selection.

Figure 18.24 illustrates part of the pulse sequence, showing when the RF pulse (Bx, top) and the gradient (G, bottom) are applied.

From Figure 18.24 of Intermediate Physics for Medicine and Biology, showing the part of the pulse sequence responsible for slice selection.
From Figure 18.24 of IPMB, showing the part of the
pulse sequence responsible for slice selection.

Did you ever wonder how the spins behave during this part of the pulse sequence? Do the spins flip all at once near the center of the π/2 pulse, or gradually throughout its duration? Why is there a second negative lobe of the gradient after the π/2 pulse is complete? To answer these questions, let’s analyze slice selection using the Bloch equations—the differential equations governing how spins precess in a magnetic field. Our analysis is similar to that in Section 18.5 of IPMB

If we ignore relaxation, the Bloch equations are

The Bloch equations describing the magnetization during slice selection in magnetic resonance imaging.

where Mx, My, and Mz are the components of the magnetization (a measure of the average spin orientation).

Next, let’s transform to a coordinate system (Mx', My', Mz') rotating at the Larmor frequency. I will leave the details of the transformation to you (they are similar to those in Section 18.5.2 of IPMB). We obtain

The Bloch equations in a rotating coordinate system, describing the magnetization during slice selection in magnetic resonance imaging.

In general the Larmor frequency γB0 is much higher than the modulation frequency γGΔz/2 (thousands of oscillations occur within the central lobe of the sinc function rather than the dozen shown in Fig. 18.23a), so we average over many cycles (as in Section 18.5.3 of IPMB). The Bloch equations simplify to

The Bloch equations in a rotating coordinate system, averaged over time, when describing the magnetization during slice selection in magnetic resonance imaging.

Before we can solve this system of three coupled ordinary differential equations, we need to specify the initial conditions. At t = −∞, Mx' = My' = 0 and Mz' = M0, where M0 is the equilibrium value of the magnetization.

At this point, we define dimensionless variables

Dimensionless variables for time, position, and magnetization, used during slice selection for magnetic resonance imaging.

where M = (Mx, My, Mz). The differential equations become

The Bloch equations in a rotating coordinate system, averaged over time, in terms of dimensionless variables, when describing the magnetization during slice selection in magnetic resonance imaging.

where z' ranges from -1 to 1, and the initial conditions at t' = −∞ are M'x' = M'y' = 0 and M'z' = 1.

I like this elegant system of equations, except all the primes are annoying. For the rest of this post I’ll drop the primes. Your job is to remember that these are dimensionless coordinates in a rotating coordinate system.

The Bloch equations in a rotating coordinate system, averaged over time, in terms of dimensionless variables, when describing the magnetization during slice selection in magnetic resonance imaging. The annoying primes are removed.

I wanted to solve these equations analytically, but couldn’t. I tried weird guesses involving the sine integral and sometimes I thought I was close, but none of my attempts worked. If you, dear reader, find an analytical solution, send it to me. You will be my hero forever.

Having no analytical solution, I wrote a simple Matlab program to solve the equations numerically.
dt=pi/1000;
T=10*pi;
nt=2*T/dt;
dz=0.05;
nz=2/dz+1;
for j=1:nz
    mx(1,j)=0;
    my(1,j)=0;
    mz(1,j)=1;
    z(j)=-1+(j-1)*dz;
    end
t(1)=-T;
for i=1:nt-1
    t(i+1)=t(i)+dt;
    for j=1:nz
        mx(i+1,j)=mx(i,j) + dt*z(j)*my(i,j);
        my(i+1,j)=my(i,j) + dt*(0.5*sin(t(i))/t(i)*mz(i,j) - z(j)*mx(i,j));
        mz(i+1,j)=mz(i,j) - dt*0.5*sin(t(i))/t(i)*my(i,j);
        end
    end
A few notes:
  • I didn’t use a fancy numerical technique, just Euler’s method. To ensure accuracy, I used a tiny time step: Δt = π/1000. If you are not familiar with solving differential equations numerically, see Section 6.14 in IPMB.
  • The sinc function, sin(t)/t, ranges from −∞ to +∞. I had to truncate it, so I applied the RF pulse from −10π to +10π. The middle, tall phase of the sinc function is centered at t = 0.
  • I actually used the software Octave, which is a free version of Matlab. I recommend it.
My results are shown below: Mx(t) (blue), My(t) (red), and Mz(t) (yellow) for different values of z, plotted for t from −10π to +10π.

The magnetization Mx, My, and Mz as functions of time and at several positions within the slice, during slice selection in magnetic resonance imaging.
More notes:
  • Mx is an odd function of z, whereas My and Mz are even. 
  • Mz changes most rapidly during the central phase of the sinc function. 
  • At z = 0, Mz rotates to My and then sits there constant. Remember, we are in the rotating coordinate system; in the laboratory coordinate system the spins precess at the Larmor frequency.
  • As z increases, the magnetization rotates in the x-y plane for times greater than zero. The larger |z|, the higher the frequency; the local precession frequency increasingly differs from the frequency of the rotating coordinate system. 
  • I don’t show z = ±1. Why not? The RF pulse was designed to contain frequencies corresponding to z between -1 and 1, so z = 1 is just on the edge of the slice. What happens outside the slice? The figure below shows that the RF pulse is ineffective at exciting spins at z = 1.25 and 1.5. At exactly z = 1 the signal is just weird.
The magnetization Mx, My, and Mz as functions of time and at several positions outside the slice, during slice selection in magnetic resonance imaging.

In an MRI experiment we don’t measure the magnetization at each value of z individually, but instead detect a signal from the entire slice. Therefore, I averaged the magnetization over z. Mx is odd, so its average is zero. The averages of My and Mz are shown below.

The magnetization Mx, My, and Mz as functions of time, averaged over the slice, during slice selection in magnetic resonance imaging.

The behavior of My is interesting. After its central peak it decays back to zero because at all the different values of z the spins precess at different frequencies and interfer with each other (they dephase). When averaged over the entire slice, My is nearly zero except briefly around t = 0.

We want My to get big and stay big. Therefore, we need to rephase the spins after they dephase. The second lobe of the gradient field does this. All the plots so far are for times from −10π to +10π, and the gradient G is on the entire time. Now, let’s extend the calculation so that for times between 10π and 20π the gradient is present but with opposite sign (like in Fig. 18.24) and the RF pulse is off. After time 20π, the gradient also turns off and the only magnetic field present is B0. The figure below shows the behavior of the spins.

The magnetization Mx, My, and Mz as functions of time, averaged over the slice, during slice selection in magnetic resonance imaging. This calculation includes the negative lobe of the gradient magnetic field that produces the echo.

We did it! When averaged over the entire slice, My (the red trace) eventually increases, reaches one, and then stays one. It remains constant for times greater than 20π because we ignore relaxation. The bottom line is that we found a way to excite the spins in the slice so they are all aligned in the y direction (in the rotating coordinate system), and we are ready to apply other fancy pulse sequences to extract information and form an image.

How did this miracle happen? The key insight is that the gradient produces an echo, much like in echo planar imaging. My decays after time zero because the spins diphase, but changing the sign of the gradient rephases the spins, so they are all lined up again when we turn the gradient off at t = 20π. We produce an echo without using a π pulse; it’s a gradient echo (if you didn’t understand that last sentence, study Figures 18.31 and 18.32 in IPMB).

Until I did this calculation, I didn’t realize that the second phase of the gradient is so crucial. I thought it corrected for a small phase shift and was only needed for high accuracy. No! The second phase produces the echo. Without it, you get no signal at all. And you have to turn it off at just the right time, or the spins will dephase again and the echo will vanish. The second lobe of the gradient is essential; the whole process fails without it.

I learned a lot from analyzing slice selection; I hope you did too.