import marimo
__generated_with = "0.23.2"
app = marimo.App(
width="medium",
app_title="MR Physics: From Protons to Brain Images",
)
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
# MR Physics: From Protons to Brain Images
*Written by Luke Chang & Ben Graul*
How do MRI scanners work? And what signal are we actually measuring?
Understanding the physics underlying MRI is essential for interpreting
neuroimaging results, yet many people who use MRI don't fully grasp
how it works -- largely because MR physics can be unintuitive.
This course primarily focuses on Blood Oxygenated Level Dependent (BOLD) fMRI signals. Gaining a deep understanding of the MR physics and physiological basis for the BOLD fMRI signal is beyond the scope of this course and we refer the interested reader to the excellent [Huettel, Song, & McCarthy (2004) Functional magnetic resonance imaging textbook](https://www.amazon.com/Functional-Magnetic-Resonance-Imaging-Huettel/dp/0878936270/ref=pd_sbs_14_1/144-9493364-1935804?_encoding=UTF8&pd_rd_i=0878936270&pd_rd_r=ac61b1df-17bf-47c5-8db5-25dfa36bcd16&pd_rd_w=J61zv&pd_rd_wg=d1O2i&pf_rd_p=703f3758-d945-4136-8df6-a43d19d750d1&pf_rd_r=PCEXDFT3TQQ4JW7FD8HF&psc=1&refRID=PCEXDFT3TQQ4JW7FD8HF) for a more in depth conceptual and quantitative overview.
The goal of this tutorial is to build your intuition from the ground
up, starting with things you already know (magnets and compasses) and
working toward functional brain imaging. Along the way, you'll interact
with simulations that let you **see** and **feel** how changing physical
parameters affects the signal.
> **Note on simplifications:** We'll use *classical* physics explanations
> throughout, following the approach advocated by
> [Lars G. Hanson](https://www.drcmr.dk/bloch). While not quantum-mechanically
> complete, this classical picture provides strong intuitions that will
> serve you well. Where the full story differs, we'll flag it in
> optional "Deep Dive" sections.
This notebook covers three major topics:
1. **Magnetism & Resonance** — magnetic fields, proton alignment, precession, RF excitation, and the FID signal
2. **Signal & Contrast** — T₁/T₂ relaxation, the Bloch equations, tissue contrast, and pulse sequences
3. **Imaging & fMRI** — gradients, spatial encoding, k-space, the BOLD signal, and EPI
""")
return
@app.cell(hide_code=True)
def _():
import marimo as mo
import numpy as np
import plotly.graph_objects as go
from pathlib import Path
from plotly.subplots import make_subplots
from dartbrains_tools.mr_simulations import (
GAMMA_H, GAMMA, TISSUE_PROPERTIES,
rotation_x, rotation_y, rotation_z,
apply_rf_pulse, apply_relaxation, simulate_bloch,
fid_signal, compute_spectrum,
t1_recovery, t2_decay,
spin_echo_signal, gradient_echo_signal,
hrf, image_to_kspace, kspace_to_image, mask_kspace,
plot_magnetization_3d, plot_signal_timeline,
plot_contrast_bars, plot_pulse_sequence,
plot_kspace_and_image,
)
from dartbrains_tools.mr_widgets import (
CompassWidget, NetMagnetizationWidget, PrecessionWidget,
SpinEnsembleWidget, EncodingWidget, KSpaceWidget, ConvolutionWidget,
)
# Resolve IMG_DIR relative to book.yml so image paths work on every
# build host (locally + GitHub Actions runners). Hardcoded absolute
# paths get baked into the rendered HTML and 404 on the deployed site.
_ROOT = next(
p for p in (Path.cwd(), *Path.cwd().resolve().parents)
if (p / "book.yml").exists()
)
IMG_DIR = _ROOT / "images" / "signal_generation"
return (
CompassWidget,
ConvolutionWidget,
EncodingWidget,
GAMMA,
GAMMA_H,
IMG_DIR,
KSpaceWidget,
NetMagnetizationWidget,
PrecessionWidget,
SpinEnsembleWidget,
TISSUE_PROPERTIES,
compute_spectrum,
fid_signal,
go,
hrf,
make_subplots,
mo,
np,
plot_contrast_bars,
spin_echo_signal,
t1_recovery,
t2_decay,
)
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
Let's first start with a short video on the basics of MR physics by Martin Lindquist.
---
# Part 1: Magnetism & Resonance
---
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.Html("""
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
## 1. What is a Magnetic Field?
You've used a compass before -- a tiny magnetized needle that aligns
with the Earth's magnetic field to point north. This simple device
illustrates the core idea behind MRI: **magnetic objects tend to align
with an applied magnetic field.**
The Earth's magnetic field is weak -- only about 25-65 *micro*Teslas
(µT). An MRI scanner produces a field that is roughly **100,000 times
stronger**: typically 1.5 or 3 Tesla (T), with research scanners
going up to 7T or beyond.
Let's start with a simulation inspired by the
[DRCMR Compass MR Simulator](https://www.drcmr.dk/CompassMR/).
Use the slider below to change the strength of the external magnetic
field (\(B_0\)) and watch how a compass needle responds.
""")
return
@app.cell(hide_code=True)
def _(mo):
b0_compass_slider = mo.ui.slider(
start=0.0, stop=5.0, step=0.1, value=3.0,
label="B\u2080 field strength (mT)",
full_width=True,
)
return (b0_compass_slider,)
@app.cell(hide_code=True)
def _(CompassWidget, b0_compass_slider, mo):
_widget = CompassWidget(b0=float(b0_compass_slider.value))
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
_wrapped,
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
**Key observations:**
- With no field (\(B_0 = 0\)), the needle just sits wherever it is -- no preferred direction.
- As you increase \(B_0\), the needle oscillates back toward alignment **faster**.
- The **frequency** of oscillation increases with field strength.
- The oscillation is detected by a nearby coil as a **signal** that decays over time.
These same principles apply inside an MRI scanner, except instead of compass
needles, we're working with hydrogen nuclei -- tiny magnets inside your body.
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 2. From Compass to Proton
We don't have compass needles inside the body, but do have a large number of **hydrogen nuclei** (protons).
Hydrogen is the most abundant element in the body (it's in every water
molecule), and each hydrogen nucleus acts like a tiny magnet because of
a quantum property called **spin**.
In the absence of an external magnetic field, these tiny magnets point
in random directions, and their effects cancel out -- there's no net
magnetization. But when we place them in a strong \(B_0\) field, a
slight majority align *with* the field rather than *against* it. This
tiny surplus creates a measurable **net magnetization** vector, called
\(M_0\), that points along \(B_0\).
Use the sliders below to see how random spins create a net magnetization
when a field is applied.
""")
return
@app.cell(hide_code=True)
def _(mo):
n_protons_slider = mo.ui.slider(
start=10, stop=500, step=10, value=100,
label="Number of protons",
)
b0_on_toggle = mo.ui.switch(label="B\u2080 field ON", value=False)
return b0_on_toggle, n_protons_slider
@app.cell(hide_code=True)
def _(NetMagnetizationWidget, b0_on_toggle, mo, n_protons_slider):
_widget = NetMagnetizationWidget(
n_protons=int(n_protons_slider.value),
b0_on=bool(b0_on_toggle.value),
)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
mo.hstack([n_protons_slider, b0_on_toggle], justify="start", gap=2),
_wrapped,
mo.callout(
mo.md(
"With the field **OFF**, spins point in random directions and jitter freely -- "
"the net magnetization (red arrow) is near zero. Toggle B\u2080 **ON** to watch "
"the spins gradually align, and a net magnetization emerges along z. "
"**Drag to rotate** the 3D view."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 3. Precession & the Larmor Frequency
When a spinning top is tilted, it doesn't just fall over -- it
**precesses**, tracing a circle as it wobbles around the vertical axis.
Protons in a magnetic field do the same thing: their spin axes precess
around the direction of \(B_0\).
The rate of precession is governed by the **Larmor equation**:
$$\omega_0 = \gamma \cdot B_0$$
where:
- \(\omega_0\) is the Larmor (precession) frequency
- \(\gamma\) is the **gyromagnetic ratio** -- a constant unique to each nucleus
- \(B_0\) is the magnetic field strength
For hydrogen protons, \(\gamma = 42.576\) MHz/T. This means at 3T,
protons precess at about **127.7 MHz** -- in the radiofrequency (RF)
range!
Adjust the field strength below and see how the precession frequency changes.
""")
return
@app.cell(hide_code=True)
def _(mo):
b0_larmor_slider = mo.ui.slider(
start=0.5, stop=7.0, step=0.1, value=3.0,
label="B\u2080 (Tesla)",
full_width=True,
)
return (b0_larmor_slider,)
@app.cell(hide_code=True)
def _(GAMMA_H, PrecessionWidget, b0_larmor_slider, mo):
_b0 = b0_larmor_slider.value
_larmor_freq = GAMMA_H * _b0
_widget = PrecessionWidget(b0=_b0, flip_angle=30.0, show_relaxation=False)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
b0_larmor_slider,
_wrapped,
])
return
@app.cell(hide_code=True)
def _(GAMMA, GAMMA_H, b0_larmor_slider, mo):
_b0 = b0_larmor_slider.value
_larmor_freq = GAMMA_H * _b0
mo.vstack([
mo.callout(
mo.md(
f"At **B\u2080 = {_b0:.1f} T**, the Larmor frequency for hydrogen is "
f"**{_larmor_freq:.1f} MHz** ({_larmor_freq/1e3:.4f} GHz). "
f"This is in the **radiofrequency** range -- the same part of the "
f"electromagnetic spectrum used by FM radio!\n\n"
"**Drag to rotate** the 3D view. The signal traces on the right show "
"|Mxy| (red) and Mz (teal) in real time."
),
kind="success",
),
mo.md("**Larmor frequencies of different nuclei at this field strength:**"),
mo.md(
"| Nucleus | \u03b3 (MHz/T) | Larmor freq at "
+ f"{_b0:.1f}T |\n|---------|-----------|----------|\n"
+ "\n".join(
f"| {name} | {ratio:.3f} | {ratio * _b0:.2f} MHz |"
for name, ratio in GAMMA.items()
)
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 4. RF Excitation & Resonance
In an MRI machine, $B_0$ is the strong, static magnetic field that runs along the bore of the scanner (the **longitudinal plane** or z-axis). It aligns hydrogen protons in your body roughly parallel to the field, creating a net magnetization $M_0$.
""")
return
@app.cell(hide_code=True)
def _(IMG_DIR, mo):
mo.image(str(IMG_DIR / "b0.png"))
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
Unfortunately, we can't directly measure magnetization
along z -- our receiver coils are only sensitive to magnetization
in the **transverse (xy) plane**.
To get a signal, we need to **tip** the magnetization away from z.
We do this by applying a second, much weaker magnetic field called
\(B_1\), oriented perpendicular to \(B_0\), and oscillating at the
**Larmor frequency**. This is the radiofrequency (RF) pulse.
Think of pushing a child on a swing: if you push at the swing's
natural frequency, each push adds energy and the swing goes higher.
Push at the wrong frequency, and nothing much happens. This is
**resonance** -- and it's the "R" in MRI.
The angle by which the magnetization is tipped is called the **flip
angle** (\(\alpha\)). A 90° pulse tips \(M\) entirely into the
xy-plane. A 180° pulse inverts it.
Try adjusting the flip angle and the B1 frequency below.
""")
return
@app.cell(hide_code=True)
def _(mo):
flip_angle_slider = mo.ui.slider(
start=0, stop=180, step=5, value=90,
label="Flip angle (degrees)",
full_width=True,
)
return (flip_angle_slider,)
@app.cell(hide_code=True)
def _(PrecessionWidget, flip_angle_slider, mo):
_flip = flip_angle_slider.value
_widget = PrecessionWidget(b0=3.0, flip_angle=float(_flip), show_relaxation=False)
_wrapped = mo.ui.anywidget(_widget)
_mxy = abs(round(float(__import__('math').sin(__import__('math').radians(_flip))), 2))
_mz = round(float(__import__('math').cos(__import__('math').radians(_flip))), 2)
mo.vstack([
flip_angle_slider,
_wrapped,
mo.md(
f"Transverse magnetization |Mxy| = **{_mxy:.2f}** "
f"(this is our detectable signal strength)\n\n"
f"Longitudinal magnetization Mz = **{_mz:.2f}**"
),
mo.callout(
mo.md(
"**Try this:** Set the flip angle to 90\u00b0 and watch the vector tip "
"fully into the xy-plane -- maximum signal! At 180\u00b0, the vector "
"inverts completely. The signal traces on the right show how |Mxy| and "
"Mz change. **Drag to rotate** the 3D view."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 5. The Signal We Measure: Free Induction Decay
After the RF pulse tips the magnetization into the transverse plane,
the precessing \(M_{xy}\) component induces a voltage in the receiver
coil -- just like the oscillating compass needle produced a signal
in Section 1. This signal is called the **Free Induction Decay (FID)**.
The FID has two key features:
1. It **oscillates** at the Larmor frequency (the precessing magnetization)
2. It **decays** over time as the transverse magnetization dephases
(characterized by the time constant \(T_2^*\))
We can extract the frequency content of the FID using a **Fourier
Transform** (FFT) -- the same mathematical tool introduced in the
Dartbrains Signal Processing chapter. The FFT reveals a peak at the
Larmor frequency, confirming that the signal came from precessing protons.
""")
return
@app.cell(hide_code=True)
def _(mo):
fid_flip_slider = mo.ui.slider(
start=5, stop=90, step=5, value=90,
label="Flip angle (\u00b0)",
)
fid_t2star_slider = mo.ui.slider(
start=10, stop=200, step=10, value=50,
label="T\u2082* (ms)",
)
fid_show_fft = mo.ui.switch(label="Show frequency spectrum (FFT)", value=False)
return fid_flip_slider, fid_show_fft, fid_t2star_slider
@app.cell(hide_code=True)
def _(
compute_spectrum,
fid_flip_slider,
fid_show_fft,
fid_signal,
fid_t2star_slider,
go,
mo,
np,
):
_flip = fid_flip_slider.value
_t2star = fid_t2star_slider.value
_show_fft = fid_show_fft.value
_dt = 0.1
_t = np.arange(0, 500, _dt)
_f0 = 0.05
_sig = fid_signal(_t, _t2star, f0=_f0, flip_angle_deg=_flip)
_fig = go.Figure()
_fig.add_trace(go.Scatter(
x=_t, y=np.real(_sig), mode="lines",
line=dict(color="#636EFA", width=2), name="FID (real part)",
))
_envelope = np.sin(np.radians(_flip)) * np.exp(-_t / _t2star)
_fig.add_trace(go.Scatter(
x=_t, y=_envelope, mode="lines",
line=dict(color="red", width=1.5, dash="dash"), name="T\u2082* envelope",
))
_fig.add_trace(go.Scatter(
x=_t, y=-_envelope, mode="lines",
line=dict(color="red", width=1.5, dash="dash"), showlegend=False,
))
_fig.update_layout(
title=f"Free Induction Decay (flip={_flip}\u00b0, T\u2082*={_t2star} ms)",
xaxis_title="Time (ms)", yaxis_title="Signal (a.u.)",
yaxis=dict(range=[-1.1, 1.1]),
width=700, height=300, margin=dict(l=60, r=20, t=40, b=40),
)
_elements = [
mo.hstack([fid_flip_slider, fid_t2star_slider, fid_show_fft], justify="start", gap=2),
_fig,
]
if _show_fft:
_freqs, _mag = compute_spectrum(_sig, _dt)
_fft_fig = go.Figure()
_fft_fig.add_trace(go.Scatter(
x=_freqs, y=_mag, mode="lines",
line=dict(color="#00CC96", width=2), name="Magnitude spectrum",
))
_fft_fig.update_layout(
title="Frequency Spectrum (FFT of FID)",
xaxis_title="Frequency (kHz)", yaxis_title="Magnitude",
yaxis=dict(range=[0, 0.55]),
width=700, height=300, margin=dict(l=60, r=20, t=40, b=40),
)
_elements.append(_fft_fig)
_elements.append(
mo.md(
f"The FFT shows a peak at **{_f0*1000:.0f} Hz** -- the precession "
f"frequency of our protons. The width of the peak is inversely "
r"related to \(T_2^*\): shorter \(T_2^*\) means faster decay, "
f"which means a broader spectral peak."
)
)
_elements.append(
mo.callout(
mo.md(
"**Key takeaways:**\n"
"- The FID amplitude depends on the flip angle (maximum at 90\u00b0)\n"
"- The decay rate is determined by T\u2082* (we'll explore this next)\n"
"- The Fourier transform reveals the precession frequency\n"
"- **This is the fundamental signal that MRI is built upon!**"
),
kind="success",
)
)
mo.vstack(_elements)
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
# Part 2: Signal & Contrast
We've learned how an RF pulse tips the net magnetization into the
transverse plane, producing a Free Induction Decay signal. But
**why does the signal decay?** And how do MRI scanners produce images
where gray matter, white matter, and CSF all look different?
The answers lie in **relaxation** -- the processes by which the
magnetization returns to equilibrium after excitation. Different
tissues relax at different rates, and MRI pulse sequences exploit
these differences to create **contrast**.
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
## 6. T1 Relaxation
After a 90° RF pulse tips all the magnetization into the xy-plane,
the longitudinal component \(M_z\) is zero. Over time, \(M_z\)
**recovers** back to its equilibrium value \(M_0\) through a process
called **T1 relaxation** (also known as spin-lattice relaxation).
This recovery follows an exponential curve:
$$M_z(t) = M_0 \left(1 - e^{-t/T_1}\right)$$
The **T1 time constant** is the time it takes for \(M_z\) to recover
to about 63% of \(M_0\). Crucially, **different tissues have different
T1 values**:
- **Fat** has a short T1 (fast recovery) because fat molecules are
large and tumble slowly, efficiently transferring energy to the lattice
- **CSF** has a long T1 (slow recovery) because water molecules are
small and tumble quickly, making energy transfer less efficient
- **Gray and white matter** fall in between
""")
return
@app.cell(hide_code=True)
def _(mo):
t1_slider = mo.ui.slider(
start=100, stop=5000, step=100, value=1300,
label="T₁ (ms)",
full_width=True,
)
t1_field_select = mo.ui.dropdown(
options={"1.5T": "1.5T", "3T": "3T"},
value="3T",
label="Field strength",
)
t1_show_tissues = mo.ui.switch(label="Show tissue curves", value=True)
return t1_field_select, t1_show_tissues, t1_slider
@app.cell(hide_code=True)
def _(
TISSUE_PROPERTIES,
go,
mo,
np,
t1_field_select,
t1_recovery,
t1_show_tissues,
t1_slider,
):
_t1_custom = t1_slider.value
_field = t1_field_select.value
_show_tissues = t1_show_tissues.value
_t = np.linspace(0, 8000, 500)
_fig = go.Figure()
# Custom T1 curve
_recovery = t1_recovery(_t, _t1_custom)
_fig.add_trace(go.Scatter(
x=_t, y=_recovery, mode="lines",
line=dict(color="#636EFA", width=3),
name=f"Custom T₁ = {_t1_custom} ms",
))
# 63% line
_fig.add_hline(y=0.63, line_dash="dot", line_color="gray",
annotation_text="63% recovery (t = T₁)")
if _show_tissues:
_tissue_colors = {
"Gray Matter": "#808080", "White Matter": "#C8A882",
"CSF": "#4169E1", "Fat": "#FFD700", "Muscle": "#CD5C5C"
}
for _tissue, _props in TISSUE_PROPERTIES[_field].items():
_r = t1_recovery(_t, _props["T1"])
_fig.add_trace(go.Scatter(
x=_t, y=_r, mode="lines",
line=dict(color=_tissue_colors.get(_tissue, "gray"), width=2, dash="dash"),
name=f"{_tissue} (T₁={_props['T1']} ms)",
))
_fig.update_layout(
title=f"T₁ Recovery at {_field}",
xaxis_title="Time after excitation (ms)",
yaxis_title="Mz / M₀",
yaxis=dict(range=[0, 1.05]),
width=750, height=400,
margin=dict(l=60, r=20, t=40, b=40),
)
mo.vstack([
mo.hstack([t1_slider, t1_field_select, t1_show_tissues], justify="start", gap=2),
_fig,
mo.callout(
mo.md(
"Notice how **fat** (short T₁) recovers quickly while **CSF** "
"(long T₁) takes much longer. Gray and white matter are in between "
"but have *different* T₁ values -- this difference is what creates "
"**T₁-weighted contrast** in MRI images."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.accordion(
{
"Deep Dive: The physics of T₁ relaxation": mo.md(
r"""
T₁ relaxation describes the return of longitudinal magnetization
to equilibrium. The recovery follows the Bloch equation for \(M_z\):
$$\frac{dM_z}{dt} = \frac{M_0 - M_z}{T_1}$$
with the solution (after a 90° excitation):
$$M_z(t) = M_0 \left(1 - e^{-t/T_1}\right)$$
The physical mechanism involves energy exchange between the
excited spin system and the surrounding molecular lattice
(hence "spin-lattice" relaxation). The efficiency depends on
molecular tumbling rates:
- **Optimal T₁ relaxation** occurs when the molecular tumbling
frequency matches the Larmor frequency
- **Small, fast-tumbling molecules** (like free water) are
inefficient → long T₁
- **Large, slow-tumbling molecules** (like fat) are more
efficient → short T₁
- T₁ **increases with field strength** because the Larmor
frequency moves further from typical tumbling frequencies
"""
)
}
)
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 7. T2 Relaxation
While \(M_z\) is recovering along the z-axis, the transverse
magnetization \(M_{xy}\) is **decaying** in the xy-plane. This happens
because individual proton spins gradually lose phase coherence --
they precess at slightly different frequencies due to interactions
with neighboring spins. This is **T2 relaxation** (spin-spin relaxation):
$$M_{xy}(t) = M_{xy}(0) \cdot e^{-t/T_2}$$
In practice, additional dephasing from \(B_0\) field inhomogeneities
makes the signal decay even faster, characterized by \(T_2^*\):
$$\frac{1}{T_2^*} = \frac{1}{T_2} + \frac{1}{T_2'}$$
where \(T_2'\) represents dephasing from field inhomogeneities.
**T2 is always shorter than or equal to T1** -- the transverse signal
always decays before the longitudinal signal fully recovers.
""")
return
@app.cell(hide_code=True)
def _(mo):
t2_slider = mo.ui.slider(
start=10, stop=300, step=10, value=80,
label="T₂ (ms)",
full_width=True,
)
t2prime_slider = mo.ui.slider(
start=10, stop=500, step=10, value=100,
label="T₂' (ms) — field inhomogeneity",
full_width=True,
)
t2_show_tissues = mo.ui.switch(label="Show tissue curves", value=True)
t2_field_select = mo.ui.dropdown(
options={"1.5T": "1.5T", "3T": "3T"}, value="3T", label="Field strength",
)
return t2_field_select, t2_show_tissues, t2_slider, t2prime_slider
@app.cell(hide_code=True)
def _(
TISSUE_PROPERTIES,
go,
mo,
np,
t2_decay,
t2_field_select,
t2_show_tissues,
t2_slider,
t2prime_slider,
):
_t2_custom = t2_slider.value
_t2prime = t2prime_slider.value
_field = t2_field_select.value
_t2star = 1.0 / (1.0 / _t2_custom + 1.0 / _t2prime)
_t = np.linspace(0, 500, 500)
_fig = go.Figure()
# T2 decay
_fig.add_trace(go.Scatter(
x=_t, y=t2_decay(_t, _t2_custom), mode="lines",
line=dict(color="#636EFA", width=3),
name=f"T₂ = {_t2_custom} ms",
))
# T2* decay
_fig.add_trace(go.Scatter(
x=_t, y=t2_decay(_t, _t2star), mode="lines",
line=dict(color="#EF553B", width=3, dash="dash"),
name=f"T₂* = {_t2star:.0f} ms",
))
# 37% line (1/e)
_fig.add_hline(y=0.37, line_dash="dot", line_color="gray",
annotation_text="37% remaining (t = T₂)")
if t2_show_tissues.value:
_tissue_colors = {
"Gray Matter": "#808080", "White Matter": "#C8A882",
"CSF": "#4169E1", "Fat": "#FFD700", "Muscle": "#CD5C5C"
}
for _tissue, _props in TISSUE_PROPERTIES[_field].items():
_d = t2_decay(_t, _props["T2"])
_fig.add_trace(go.Scatter(
x=_t, y=_d, mode="lines",
line=dict(color=_tissue_colors.get(_tissue, "gray"), width=2, dash="dot"),
name=f"{_tissue} (T₂={_props['T2']} ms)",
))
_fig.update_layout(
title=f"T₂ and T₂* Decay at {_field}",
xaxis_title="Time after excitation (ms)",
yaxis_title="Mxy / Mxy(0)",
yaxis=dict(range=[0, 1.05]),
width=750, height=400,
margin=dict(l=60, r=20, t=40, b=40),
)
mo.vstack([
mo.vstack([
mo.hstack([t2_slider, t2prime_slider], gap=2),
mo.hstack([t2_field_select, t2_show_tissues], justify="start", gap=2),
]),
_fig,
mo.md(
f"T₂* = **{_t2star:.0f} ms** (always shorter than T₂ = {_t2_custom} ms due "
f"to field inhomogeneities)"
),
mo.callout(
mo.md(
"**T₂ vs T₂*:** T₂ is an intrinsic tissue property (spin-spin interactions). "
"T₂* includes *additional* dephasing from magnetic field non-uniformities. "
"Spin echo sequences can recover the T₂' component (getting back to 'true' T₂), "
"but gradient echo sequences are sensitive to T₂*. This distinction will "
"become important for fMRI later in this notebook!"
),
kind="warn",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 8. T1 and T2 Together: The Full Picture
After an RF excitation pulse, **both relaxation processes happen
simultaneously**: \(M_z\) recovers (T1) while \(M_{xy}\) decays (T2).
The combined motion of the magnetization vector traces a **spiral**
path as it precesses, dephases in the transverse plane, and recovers
along the longitudinal axis.
This is the complete Bloch equation picture. The visualization below
shows the 3D trajectory of the magnetization vector after a 90° pulse.
""")
return
@app.cell(hide_code=True)
def _(mo):
bloch_t1_slider = mo.ui.slider(
start=200, stop=4000, step=100, value=1300,
label="T\u2081 (ms)",
)
bloch_t2_slider = mo.ui.slider(
start=10, stop=300, step=10, value=80,
label="T\u2082 (ms)",
)
bloch_b0_slider = mo.ui.slider(
start=0.5, stop=7.0, step=0.5, value=3.0,
label="B\u2080 (T)",
)
return bloch_b0_slider, bloch_t1_slider, bloch_t2_slider
@app.cell(hide_code=True)
def _(PrecessionWidget, bloch_b0_slider, bloch_t1_slider, bloch_t2_slider, mo):
_t1 = bloch_t1_slider.value
_t2 = bloch_t2_slider.value
_b0 = bloch_b0_slider.value
_widget = PrecessionWidget(
b0=float(_b0), flip_angle=90.0,
t1=float(_t1), t2=float(_t2),
show_relaxation=True,
)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
mo.hstack([bloch_t1_slider, bloch_t2_slider, bloch_b0_slider], justify="start", gap=2),
_wrapped,
mo.callout(
mo.md(
"**The spiral tells the whole story:** The magnetization simultaneously "
"precesses (circles in xy), dephases (|Mxy| decays via T\u2082), and "
"recovers (Mz returns to equilibrium via T\u2081). Try making T\u2082 very short -- the "
"signal collapses quickly. Make T\u2081 very short -- Mz snaps back fast.\n\n"
"The signal traces on the right show |Mxy| (red) and Mz (teal) in real time. "
"**Drag to rotate** the 3D view."
),
kind="success",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.accordion(
{
"Deep Dive: The complete Bloch equations": mo.md(
r"""
The Bloch equations describe the time evolution of the
magnetization vector \(\vec{M} = (M_x, M_y, M_z)\) in a
magnetic field:
$$\frac{dM_x}{dt} = \gamma (\vec{M} \times \vec{B})_x - \frac{M_x}{T_2}$$
$$\frac{dM_y}{dt} = \gamma (\vec{M} \times \vec{B})_y - \frac{M_y}{T_2}$$
$$\frac{dM_z}{dt} = \gamma (\vec{M} \times \vec{B})_z - \frac{M_z - M_0}{T_1}$$
The cross-product term describes precession, while the decay
terms describe relaxation. In matrix form, for a single time
step \(\Delta t\) with only \(B_0\) along z:
1. **Precession**: Rotate \(M\) about z by angle \(\Delta\phi = \gamma B_0 \Delta t\)
2. **T2 decay**: Multiply \(M_x\) and \(M_y\) by \(e^{-\Delta t / T_2}\)
3. **T1 recovery**: Update \(M_z \rightarrow M_z \cdot e^{-\Delta t / T_1} + M_0 (1 - e^{-\Delta t / T_1})\)
This is exactly how our simulation module implements the Bloch
equations -- using rotation matrices for precession and
exponential factors for relaxation.
"""
)
}
)
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 9. Tissue Contrast: How TE and TR Create Different Images
This is where MRI gets its remarkable soft-tissue contrast. By choosing
when we **repeat** the excitation pulse (\(TR\) = Repetition Time) and
when we **read** the signal (\(TE\) = Echo Time), we can make the image
sensitive to different tissue properties:
| Weighting | TR | TE | Sensitive to |
|-----------|----|----|-------------|
| **T1-weighted** | Short (~500ms) | Short (~10ms) | T1 differences (anatomy) |
| **T2-weighted** | Long (~3000ms) | Long (~80ms) | T2 differences (pathology) |
| **PD-weighted** | Long (~3000ms) | Short (~10ms) | Proton density |
The spin echo signal equation combines both relaxation effects:
$$S = PD \cdot (1 - e^{-TR/T_1}) \cdot e^{-TE/T_2}$$
The first term (\(1 - e^{-TR/T_1}\)) is T1-weighting: short TR means
tissues haven't fully recovered, so T1 differences matter. The second
term (\(e^{-TE/T_2}\)) is T2-weighting: long TE means more T2 decay
has occurred, so T2 differences matter.
**This is the most important interactive plot in this notebook.** Use the
sliders below to explore how TE and TR change the contrast between tissues.
""")
return
@app.cell(hide_code=True)
def _(mo):
te_slider = mo.ui.slider(
start=5, stop=200, step=5, value=10,
label="TE (ms)",
full_width=True,
)
tr_slider = mo.ui.slider(
start=100, stop=6000, step=100, value=500,
label="TR (ms)",
full_width=True,
)
contrast_field = mo.ui.dropdown(
options={"1.5T": "1.5T", "3T": "3T"}, value="3T", label="Field strength",
)
return contrast_field, te_slider, tr_slider
@app.cell(hide_code=True)
def _(
TISSUE_PROPERTIES,
contrast_field,
go,
mo,
np,
plot_contrast_bars,
spin_echo_signal,
t1_recovery,
t2_decay,
te_slider,
tr_slider,
):
_te = te_slider.value
_tr = tr_slider.value
_field = contrast_field.value
_tissues = TISSUE_PROPERTIES[_field]
# Calculate signal for each tissue
_tissue_names = list(_tissues.keys())
_signals = []
for _name in _tissue_names:
_p = _tissues[_name]
_s = spin_echo_signal(_te, _tr, _p["T1"], _p["T2"], _p["PD"])
_signals.append(_s)
# Determine weighting type
if _tr < 800 and _te < 30:
_weighting = "T₁-weighted"
elif _tr > 2000 and _te > 60:
_weighting = "T₂-weighted"
elif _tr > 2000 and _te < 30:
_weighting = "PD-weighted"
else:
_weighting = "Mixed weighting"
# Bar chart (compact)
_bar_fig = plot_contrast_bars(_tissue_names, _signals,
title=f"{_weighting}")
_bar_fig.update_layout(width=300, height=300, margin=dict(l=40, r=5, t=30, b=60),
xaxis_tickangle=-35, xaxis_tickfont=dict(size=9))
# T1 recovery and T2 decay as separate compact plots
_t_t1 = np.linspace(0, 6000, 500)
_t_t2 = np.linspace(0, 300, 500)
_tissue_colors = {
"Gray Matter": "#808080", "White Matter": "#C8A882",
"CSF": "#4169E1", "Fat": "#FFD700", "Muscle": "#CD5C5C"
}
# T1 recovery plot
_t1_fig = go.Figure()
for _name in _tissue_names:
_p = _tissues[_name]
_color = _tissue_colors.get(_name, "gray")
_t1_fig.add_trace(go.Scatter(
x=_t_t1, y=t1_recovery(_t_t1, _p["T1"], _p["PD"]),
mode="lines", line=dict(color=_color, width=2), name=_name,
))
_t1_fig.add_vline(x=_tr, line_dash="dash", line_color="red",
annotation_text=f"TR={_tr}")
_t1_fig.update_layout(
title="T\u2081 Recovery", width=350, height=300,
margin=dict(l=40, r=5, t=30, b=40),
xaxis_title="Time (ms)", yaxis_title="Signal",
yaxis=dict(range=[0, 1.05]),
legend=dict(font=dict(size=8), x=0.55, y=0.35),
)
# T2 decay plot
_t2_fig = go.Figure()
for _name in _tissue_names:
_p = _tissues[_name]
_color = _tissue_colors.get(_name, "gray")
_t2_fig.add_trace(go.Scatter(
x=_t_t2, y=_p["PD"] * t2_decay(_t_t2, _p["T2"]),
mode="lines", line=dict(color=_color, width=2),
name=_name, showlegend=False,
))
_t2_fig.add_vline(x=_te, line_dash="dash", line_color="red",
annotation_text=f"TE={_te}")
_t2_fig.update_layout(
title="T\u2082 Decay", width=350, height=300,
margin=dict(l=40, r=5, t=30, b=40),
xaxis_title="Time (ms)", yaxis_title="Signal",
yaxis=dict(range=[0, 1.05]),
)
mo.vstack([
mo.hstack([te_slider, tr_slider, contrast_field], justify="start", gap=2),
mo.hstack([_bar_fig, _t1_fig, _t2_fig], justify="center"),
mo.callout(
mo.md(
f"**Current weighting: {_weighting}** \u2014 "
"Try: TR=500, TE=10 (T\u2081w) | TR=4000, TE=80 (T\u2082w) | TR=4000, TE=10 (PD)\n\n"
"The red dashed lines show where TR and TE sample the recovery/decay curves. "
"The bar chart shows the resulting signal for each tissue."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 10. Spin Echo & Gradient Echo
We've been talking about the signal at time TE after excitation. But
how do we actually *collect* the signal at a specific echo time? There
are two fundamental approaches:
### Spin Echo (SE)
A **180° refocusing pulse** is applied at time TE/2 after the initial
90° excitation. This reverses the dephasing caused by static field
inhomogeneities, causing the spins to **rephase** and form an echo
at time TE. The spin echo signal depends on **T2** (not T2*), because
the refocusing pulse undoes the T2' dephasing.
""")
return
@app.cell(hide_code=True)
def _(IMG_DIR, mo):
mo.image(str(IMG_DIR / "spin_echo_pulse_sequence.svg"))
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
The diagram shows one full spin echo cycle.
A **spin echo pulse sequence** begins with a 90° RF pulse that tips the net magnetization from the longitudinal axis into the transverse plane. Immediately after excitation, spins begin to dephase due to local field inhomogeneities, producing a decaying free induction decay (FID).
A 180° refocusing pulse applied at time TE/2 reverses this dephasing, causing spins to rephase and form an echo at time **TE**.
The sequence repeats every **TR** (indicated by the faded 90° pulse at the right edge).
Three spatial encoding gradients operate alongside the RF pulses:
- **Gss (slice select)** is applied during both RF pulses to restrict excitation to a single imaging slice
- **Gpe (phase encode)** applies a brief gradient of varying amplitude on each repetition, stepping through k-space line by line
- **Gro (readout/frequency encode)** is switched on during the echo to spatially encode signal along the remaining axis.
Together, these three gradients provide the spatial information needed to reconstruct a 2D image from the acquired signal.
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
### Gradient Echo (GRE)
Instead of a 180° pulse, a **gradient reversal** is used to form the
echo. This is faster but does NOT refocus static field inhomogeneities,
so the signal depends on **T2*** (not T2). Gradient echoes are the
basis of most fMRI sequences (because fMRI needs T2* sensitivity!).
""")
return
@app.cell(hide_code=True)
def _(IMG_DIR, mo):
mo.image(str(IMG_DIR / "gradient_echo_pulse_sequence.svg"))
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
A **gradient echo pulse sequence** begins with a low flip angle (α) RF pulse — typically much less than 90° — that tips a fraction of the longitudinal magnetization into the transverse plane.
Unlike the spin echo, there is no 180° refocusing pulse.
Instead, the echo is formed entirely by gradient manipulation: a negative dephasing lobe on the **readout gradient (Gro)** deliberately dephases the spins, and then a positive lobe of opposite polarity rephases them to produce the gradient echo at time TE.
Because the 180° pulse is absent, static field inhomogeneities are not corrected, so the signal decays with T2* rather than T2.
The **slice select gradient (Gss)** is applied during the α pulse to restrict excitation to a single slice, and the **phase encode gradient (Gpe)** steps through different amplitudes on each repetition to fill k-space.
The combination of a small flip angle and no refocusing pulse allows very short TR and TE, making gradient echo sequences the basis for fast imaging methods such as FLASH, SPGR, and MPRAGE.
Here is an interactive figure showing this process over time for each sequence type.
""")
return
@app.cell(hide_code=True)
def _(mo):
echo_type = mo.ui.radio(
options={"Spin Echo (90\u00b0-180\u00b0)": "spin_echo", "Gradient Echo": "gradient_echo"},
value="Spin Echo (90\u00b0-180\u00b0)",
label="Sequence type",
)
echo_speed = mo.ui.slider(
start=0.3, stop=3.0, step=0.1, value=1.0,
label="Animation speed",
)
return echo_speed, echo_type
@app.cell(hide_code=True)
def _(SpinEnsembleWidget, echo_speed, echo_type, mo):
_seq_type = echo_type.value
_speed = echo_speed.value
_widget = SpinEnsembleWidget(sequence_type=_seq_type, speed=_speed)
_wrapped = mo.ui.anywidget(_widget)
if _seq_type == "spin_echo":
_seq_desc = (
"**Spin Echo**: 90\u00b0 \u2192 wait TE/2 \u2192 180\u00b0 (refocus) \u2192 wait TE/2 \u2192 Echo\n\n"
"The 180\u00b0 pulse reverses all static dephasing. The echo signal "
"reflects **true T\u2082** decay only. Used for anatomical imaging (T\u2081w, T\u2082w)."
)
else:
_seq_desc = (
"**Gradient Echo**: \u03b1\u00b0 \u2192 gradient dephase \u2192 gradient rephase \u2192 Echo\n\n"
"No refocusing pulse, so static field inhomogeneities are NOT reversed. "
"The echo signal reflects **T\u2082*** decay. Faster than spin echo. "
"Used for **fMRI** (BOLD signal depends on T\u2082*!)."
)
mo.vstack([
mo.hstack([echo_type, echo_speed], justify="start", gap=2),
_wrapped,
mo.md(_seq_desc),
mo.callout(
mo.md(
"**Why does this matter for fMRI?** The BOLD signal (which we'll "
"cover next) depends on local magnetic field changes caused "
"by deoxyhemoglobin. These are T\u2082* effects -- they would be refocused "
"away by a spin echo! That's why fMRI uses **gradient echo** sequences.\n\n"
"Watch the signal trace on the right: in **spin echo**, the signal fully "
"recovers at the echo. In **gradient echo**, residual dephasing remains."
),
kind="warn",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
# Part 3: Imaging & fMRI
We now know how to generate a signal and how different tissues produce
different signals. But there's a critical problem: **where did the
signal come from?**
If the entire volume inside the scanner experiences the same \(B_0\)
field, all protons precess at the same frequency, and we can't tell
whether the signal is coming from your frontal cortex or your cerebellum.
We need **spatial encoding** -- a way to tag different locations with
different frequencies or phases so we can reconstruct an image.
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.Html("""
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 11. The Localization Problem: Why We Need Gradients
In a perfectly uniform \(B_0\) field, every proton in the body
precesses at the same Larmor frequency. Our receiver coil picks up
the *sum* of all these signals, with no way to tell where each
contribution came from.
The solution: **gradient coils**. These are additional electromagnetic
coils that make the magnetic field vary *linearly* across space. If
we add a gradient along the x-axis:
$$B(x) = B_0 + G_x \cdot x$$
then the Larmor frequency also varies with position:
$$f(x) = \gamma \cdot (B_0 + G_x \cdot x)$$
Now protons at different x-positions precess at different frequencies
-- and we can separate them using the Fourier transform!
""")
return
@app.cell(hide_code=True)
def _(mo):
grad_strength_slider = mo.ui.slider(
start=0.0, stop=40.0, step=1.0, value=0.0,
label="Gradient strength Gx (mT/m)",
full_width=True,
)
return (grad_strength_slider,)
@app.cell(hide_code=True)
def _(GAMMA_H, go, grad_strength_slider, make_subplots, mo, np):
_gx = grad_strength_slider.value # mT/m
_b0 = 3.0 # Tesla
_positions = np.linspace(-0.15, 0.15, 200) # meters (30cm FOV)
# B field at each position
_b_field = _b0 + (_gx / 1000) * _positions # Convert mT/m to T/m
_freq = GAMMA_H * _b_field # MHz
_fig = make_subplots(rows=1, cols=2,
subplot_titles=["Magnetic Field B(x)", "Larmor Frequency f(x)"])
# B field
_fig.add_trace(go.Scatter(
x=_positions * 100, y=_b_field,
mode="lines", line=dict(color="#636EFA", width=3),
name="B(x)",
), row=1, col=1)
# Frequency
_fig.add_trace(go.Scatter(
x=_positions * 100, y=_freq,
mode="lines", line=dict(color="#EF553B", width=3),
name="f(x)",
), row=1, col=2)
_fig.update_xaxes(title_text="Position (cm)", row=1, col=1)
_fig.update_xaxes(title_text="Position (cm)", row=1, col=2)
_fig.update_yaxes(title_text="B (Tesla)", row=1, col=1)
_fig.update_yaxes(title_text="Frequency (MHz)", row=1, col=2)
_freq_range = _freq[-1] - _freq[0]
_fig.update_layout(
width=800, height=350, margin=dict(l=60, r=20, t=40, b=40),
yaxis=dict(range=[2.99, 3.01]),
yaxis2=dict(range=[127.3, 128.2]),
)
mo.vstack([
mo.hstack([grad_strength_slider], justify="start", gap=2),
_fig,
mo.md(
f"**Gradient: {_gx:.0f} mT/m** → Frequency spread across 30cm FOV: "
f"**{abs(_freq_range * 1000):.1f} kHz** "
f"({'uniform field, no spatial encoding!' if _gx == 0 else 'spatial encoding active'})"
),
mo.callout(
mo.md(
"With **no gradient** (Gx=0), the field and frequency are the same "
"everywhere -- no spatial information. As you increase the gradient, "
"frequency varies across space, and the Fourier transform of the signal "
"directly maps to spatial position!"
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 12. Slice Selection
MRI typically images one **slice** at a time. To select a specific
slice, we apply a gradient along the z-axis (head-to-foot) during
the RF pulse. The RF pulse has a limited bandwidth -- it only excites
protons within a narrow range of frequencies. Since the z-gradient
makes frequency vary with position, only protons in one slice
are on-resonance and get excited.
**Slice thickness** is controlled by:
$$\Delta z = \frac{\Delta f}{\gamma \cdot G_z}$$
where \(\Delta f\) is the RF pulse bandwidth and \(G_z\) is the
slice-select gradient strength. Stronger gradient = thinner slice.
Wider RF bandwidth = thicker slice.
""")
return
@app.cell(hide_code=True)
def _(mo):
ss_grad_slider = mo.ui.slider(
start=5, stop=40, step=1, value=20,
label="Slice-select gradient Gz (mT/m)",
)
ss_bw_slider = mo.ui.slider(
start=500, stop=5000, step=100, value=2000,
label="RF bandwidth (Hz)",
)
return ss_bw_slider, ss_grad_slider
@app.cell(hide_code=True)
def _(GAMMA_H, go, mo, np, ss_bw_slider, ss_grad_slider):
_gz = ss_grad_slider.value / 1000 # T/m
_bw = ss_bw_slider.value # Hz
_b0 = 3.0
# Slice thickness in mm
_slice_thickness = (_bw / (GAMMA_H * 1e6 * _gz)) * 1000 # mm
# Positions along z
_z_pos = np.linspace(-100, 100, 500) # mm
_freq_hz = GAMMA_H * 1e6 * (_b0 + _gz * _z_pos / 1000) - GAMMA_H * 1e6 * _b0 # Hz offset from center
# RF excitation profile (sinc-like bandwidth)
_center_freq = 0 # Hz (excite at isocenter)
_excited = np.abs(_freq_hz - _center_freq) < (_bw / 2)
_fig = go.Figure()
# Frequency vs position
_fig.add_trace(go.Scatter(
x=_z_pos, y=_freq_hz / 1000,
mode="lines", line=dict(color="#636EFA", width=2),
name="Frequency offset",
))
# Excited region
_fig.add_trace(go.Scatter(
x=_z_pos[_excited], y=_freq_hz[_excited] / 1000,
mode="lines", line=dict(color="red", width=4),
name="Excited slice",
))
# RF bandwidth band
_fig.add_hrect(y0=-_bw / 2000, y1=_bw / 2000,
fillcolor="rgba(255, 0, 0, 0.1)", line_width=0,
annotation_text="RF bandwidth")
_fig.update_layout(
title=f"Slice Selection: thickness = {_slice_thickness:.1f} mm",
xaxis_title="Position along z (mm)",
yaxis_title="Frequency offset (kHz)",
yaxis=dict(range=[-5, 5]),
width=700, height=400,
margin=dict(l=60, r=20, t=40, b=40),
)
mo.vstack([
mo.hstack([ss_grad_slider, ss_bw_slider], justify="start", gap=2),
_fig,
mo.md(
f"With Gz = **{_gz * 1000:.0f} mT/m** and RF bandwidth = **{_bw} Hz**, "
f"the selected slice is **{_slice_thickness:.1f} mm** thick.\n\n"
"The red segment shows the positions where protons are within the RF bandwidth "
"and will be excited. All other protons are off-resonance and unaffected."
),
mo.callout(
mo.md(
"**Try this:** Increase the gradient strength -- the slice gets thinner "
"(better resolution but less signal). Increase the RF bandwidth -- the "
"slice gets thicker (more signal but lower resolution). This is a "
"fundamental tradeoff in MRI."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 13. Frequency & Phase Encoding
After selecting a slice, we still need to encode position within
that 2D slice. MRI uses two clever tricks:
**Frequency encoding (readout direction, x):** During signal readout,
a gradient along x makes protons at different x-positions precess at
different frequencies. The Fourier transform of the readout signal
directly gives us the spatial profile along x.
**Phase encoding (y-direction):** Before readout, a brief gradient
pulse along y gives protons at different y-positions different *phases*.
This is repeated many times with different gradient strengths to fill
out the y-dimension. Each repetition fills one line of k-space.
The combination gives each voxel a unique (frequency, phase) signature
that can be decoded by a 2D Fourier transform.
""")
return
@app.cell(hide_code=True)
def _(EncodingWidget, mo):
_widget = EncodingWidget(speed=1.0)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
_wrapped,
mo.md(
"Watch how the spin arrows in each voxel respond to gradients:\n"
"- **No gradients**: all spins precess at the same rate\n"
"- **Frequency encoding (Gx)**: columns spin at different rates\n"
"- **Phase encoding (Gy)**: rows accumulate different phase offsets\n"
"- **Both**: each voxel gets a unique (frequency, phase) signature\n\n"
"The 2D Fourier transform decodes these signatures back into an image."
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.accordion(
{
"Deep Dive: The Fourier encoding equations": mo.md(
r"""
The MRI signal at time \(t\) during readout, for a given
phase-encode step with gradient area \(A_y\), is:
$$S(t) = \int\int \rho(x, y) \cdot e^{-i 2\pi (\gamma G_x x t + \gamma A_y y)} \, dx \, dy$$
where \(\rho(x, y)\) is the spin density (our image).
If we define spatial frequencies:
- \(k_x = \gamma G_x t\) (varies during readout)
- \(k_y = \gamma A_y\) (set by phase-encode gradient)
then the signal equation becomes:
$$S(k_x, k_y) = \int\int \rho(x, y) \cdot e^{-i 2\pi (k_x x + k_y y)} \, dx \, dy$$
This is exactly the **2D Fourier transform** of the image!
Therefore, the image is simply the inverse Fourier transform
of the acquired data:
$$\rho(x, y) = \mathcal{F}^{-1}\{S(k_x, k_y)\}$$
The space of \((k_x, k_y)\) is called **k-space**.
"""
)
}
)
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 14. K-Space: Where MRI Data Lives
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.Html("""
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
The raw data from an MRI scan doesn't look like an image -- it lives
in **k-space**, the spatial frequency domain. K-space and the image
are related by the 2D Fourier transform:
$$\text{Image} = \text{FFT}^{-1}(\text{K-space})$$
Different regions of k-space encode different features:
- **Center of k-space** → overall contrast, brightness, low-frequency
shapes
- **Edges of k-space** → fine details, sharp edges, high spatial
frequency
This means you don't need ALL of k-space to get a useful image.
Keeping only the center gives you a blurry but recognizable image.
Keeping only the edges gives you an edge map with no contrast.
Try masking different regions below to build intuition!
""")
return
@app.cell(hide_code=True)
def _(mo):
kspace_mask_type = mo.ui.dropdown(
options={
"Progressive fill (animated)": "progressive",
"Center only (low frequencies)": "center",
"Periphery only (high frequencies)": "periphery",
"Every 4th line (undersampled)": "undersampled",
},
value="Progressive fill (animated)",
label="K-space mode",
)
kspace_radius = mo.ui.slider(
start=0.05, stop=0.5, step=0.05, value=0.2,
label="Mask radius",
)
kspace_speed = mo.ui.slider(
start=0.5, stop=5.0, step=0.5, value=2.0,
label="Fill speed",
)
return kspace_mask_type, kspace_radius, kspace_speed
@app.cell(hide_code=True)
def _(KSpaceWidget, kspace_mask_type, kspace_radius, kspace_speed, mo):
_widget = KSpaceWidget(
mask_type=kspace_mask_type.value,
radius_fraction=float(kspace_radius.value),
speed=float(kspace_speed.value),
)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
mo.hstack([kspace_mask_type, kspace_radius, kspace_speed], justify="start", gap=2),
_wrapped,
mo.md(
"**Progressive fill** shows how an MRI scanner acquires k-space line by line, "
"with the image emerging gradually. Try other modes:\n"
"- **Center only**: Blurry but recognizable -- contrast lives in the center\n"
"- **Periphery only**: Only edges visible -- detail lives at the edges\n"
"- **Undersampled**: Aliasing artifacts from skipping lines"
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 15. BOLD & the Hemodynamic Response
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.Html("""
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
Now we arrive at **functional MRI (fMRI)** -- using MRI to measure
brain *activity*. But MRI doesn't measure neural activity directly.
Instead, it measures a downstream consequence: changes in **blood
oxygenation**.
The key insight: **deoxyhemoglobin is paramagnetic** (slightly magnetic),
while **oxyhemoglobin is diamagnetic** (not magnetic). When a brain
region becomes active:
1. Neurons fire and consume oxygen locally
2. The vascular system responds by increasing blood flow to that region
3. Blood flow *overshoots* metabolic demand -- more oxygen arrives
than is consumed
4. The result: less deoxyhemoglobin locally
5. Less deoxyhemoglobin → less local field distortion → **increased
T₂* signal**
This is the **Blood Oxygen Level Dependent (BOLD)** signal. The
temporal shape of the BOLD response to a brief neural event is called
the **Hemodynamic Response Function (HRF)**.
""")
return
@app.cell(hide_code=True)
def _(mo):
hrf_peak_slider = mo.ui.slider(
start=3.0, stop=9.0, step=0.5, value=6.0,
label="Peak time (s)",
)
hrf_undershoot_slider = mo.ui.slider(
start=10.0, stop=25.0, step=1.0, value=16.0,
label="Undershoot time (s)",
)
hrf_ratio_slider = mo.ui.slider(
start=0.0, stop=0.5, step=0.05, value=0.167,
label="Undershoot ratio",
)
return hrf_peak_slider, hrf_ratio_slider, hrf_undershoot_slider
@app.cell(hide_code=True)
def _(
go,
hrf,
hrf_peak_slider,
hrf_ratio_slider,
hrf_undershoot_slider,
mo,
np,
):
_peak = hrf_peak_slider.value
_undershoot = hrf_undershoot_slider.value
_ratio = hrf_ratio_slider.value
_t = np.linspace(0, 30, 300)
_h = hrf(_t, peak_time=_peak, undershoot_time=_undershoot,
undershoot_ratio=_ratio)
_fig = go.Figure()
_fig.add_trace(go.Scatter(
x=_t, y=_h, mode="lines",
line=dict(color="#636EFA", width=3), name="HRF",
))
_fig.add_hline(y=0, line_dash="dot", line_color="gray")
# Annotate key features
_peak_idx = np.argmax(_h)
_fig.add_annotation(x=_t[_peak_idx], y=_h[_peak_idx],
text=f"Peak (~{_t[_peak_idx]:.1f}s)",
arrowhead=2, ay=-40)
if _ratio > 0:
_undershoot_idx = np.argmin(_h[_peak_idx:]) + _peak_idx
if _h[_undershoot_idx] < 0:
_fig.add_annotation(x=_t[_undershoot_idx], y=_h[_undershoot_idx],
text=f"Undershoot (~{_t[_undershoot_idx]:.1f}s)",
arrowhead=2, ay=40)
_fig.update_layout(
title="Hemodynamic Response Function (HRF)",
xaxis_title="Time after stimulus (s)",
yaxis_title="BOLD signal change (a.u.)",
yaxis=dict(range=[-0.3, 1.1]),
width=700, height=350,
margin=dict(l=60, r=20, t=40, b=40),
)
mo.vstack([
mo.hstack([hrf_peak_slider, hrf_undershoot_slider, hrf_ratio_slider], justify="start", gap=2),
_fig,
mo.md(
"The canonical HRF has several notable features:\n"
f"- **Initial dip** (small, often not detectable at typical fMRI resolution)\n"
f"- **Peak** at ~{_peak:.0f}s after stimulus\n"
f"- **Post-stimulus undershoot** at ~{_undershoot:.0f}s\n"
"- Returns to baseline after ~25-30s\n\n"
"This sluggish hemodynamic response means fMRI has poor **temporal resolution** "
"(~seconds) compared to EEG (milliseconds), even though the neural events "
"happen on the millisecond timescale."
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
### Convolution: From Events to Predicted BOLD Signal
In an fMRI experiment, we present multiple stimuli over time. The
predicted BOLD signal is the **convolution** of the stimulus timing
with the HRF:
$$\text{Predicted BOLD}(t) = \text{stimulus}(t) * \text{HRF}(t)$$
This is the foundation of the **General Linear Model (GLM)** used
in fMRI analysis (covered in detail in later Dartbrains chapters).
Try adjusting the stimulus timing below to see how the predicted BOLD
signal changes.
""")
return
@app.cell(hide_code=True)
def _(mo):
stim_pattern = mo.ui.dropdown(
options={
"Single event": "single",
"3 events (spaced)": "spaced",
"Block design (10s on/off)": "block",
"Rapid event-related": "rapid",
},
value="Single event",
label="Stimulus pattern",
)
conv_speed = mo.ui.slider(
start=0.5, stop=4.0, step=0.5, value=1.5,
label="Animation speed",
)
return conv_speed, stim_pattern
@app.cell(hide_code=True)
def _(ConvolutionWidget, conv_speed, mo, stim_pattern):
_widget = ConvolutionWidget(
pattern=stim_pattern.value,
speed=float(conv_speed.value),
)
_wrapped = mo.ui.anywidget(_widget)
mo.vstack([
mo.hstack([stim_pattern, conv_speed], justify="start", gap=2),
_wrapped,
mo.callout(
mo.md(
"Watch how each stimulus event spawns its own HRF response (purple ghost curves), "
"and the predicted BOLD signal (blue) is their sum. In block designs, the BOLD "
"builds up during ON periods. In rapid event-related designs, individual HRFs "
"overlap -- this is **linear superposition**."
),
kind="info",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 16. fMRI Pulse Sequences: Echo Planar Imaging (EPI)
Standard MRI acquires one line of k-space per TR -- far too slow for
fMRI, which needs to image the whole brain every 1-2 seconds.
The solution is **Echo Planar Imaging (EPI)**: after a single RF
excitation, rapidly oscillating gradients traverse *all* of k-space
in about 50-100 milliseconds. This is the workhorse sequence for fMRI.
**Typical fMRI parameters:**
| Parameter | Structural (T₁w MPRAGE) | Functional (GRE-EPI) |
|-----------|------------------------|---------------------|
| Weighting | T₁ | T₂* |
| TR | ~2000 ms | 500-2000 ms |
| TE | ~3 ms | ~30 ms |
| Flip angle | ~9° | ~70-90° |
| Resolution | ~1 mm³ | ~2-3 mm³ |
| Volumes | 1 (whole brain) | 100s-1000s |
| Duration | ~5 min | 5-60 min |
""")
return
@app.cell(hide_code=True)
def _(mo):
epi_tr_slider = mo.ui.slider(
start=500, stop=3000, step=100, value=1000,
label="TR (ms)",
)
epi_te_slider = mo.ui.slider(
start=15, stop=50, step=5, value=30,
label="TE (ms)",
)
epi_matrix_slider = mo.ui.slider(
start=32, stop=128, step=16, value=64,
label="Matrix size",
)
return epi_matrix_slider, epi_te_slider, epi_tr_slider
@app.cell(hide_code=True)
def _(epi_matrix_slider, epi_te_slider, epi_tr_slider, go, mo):
_tr = epi_tr_slider.value
_te = epi_te_slider.value
_matrix = epi_matrix_slider.value
# Calculate derived parameters
_fov = 220 # mm (typical brain FOV)
_voxel_size = _fov / _matrix
_readout_time = _matrix * 0.5 # ~0.5ms per readout line (simplified)
_n_slices = int(_tr / (_readout_time + 5)) # rough estimate
# Tradeoff visualization
_metrics = {
"Spatial Resolution": f"{_voxel_size:.1f} mm",
"Temporal Resolution": f"{_tr/1000:.1f} s",
"Slices per TR": f"~{_n_slices}",
"BOLD sensitivity": f"{'Good' if 25 <= _te <= 35 else 'Suboptimal'} (TE={_te}ms)",
"Readout duration": f"~{_readout_time:.0f} ms",
}
# K-space trajectory: line density reflects matrix size
_fig = go.Figure()
_n_lines = _matrix
_kx_all = []
_ky_all = []
_colors = []
for _line in range(_n_lines):
_y = _line / _n_lines - 0.5
if _line % 2 == 0:
_kx_all.extend([-0.5, 0.5, None])
else:
_kx_all.extend([0.5, -0.5, None])
_ky_all.extend([_y, _y, None])
_fig.add_trace(go.Scatter(
x=_kx_all, y=_ky_all, mode="lines",
line=dict(color="#636EFA", width=max(0.5, 3 - _n_lines / 40)),
name=f"{_matrix} PE lines",
))
_fig.update_layout(
title=f"EPI Trajectory ({_matrix} lines)",
xaxis_title="kx", yaxis_title="ky",
xaxis=dict(range=[-0.6, 0.6]),
yaxis=dict(range=[-0.6, 0.6], scaleanchor="x"),
width=350, height=350,
margin=dict(l=50, r=10, t=35, b=40),
)
# Parameters table
_table_md = "| Parameter | Value |\n|-----------|-------|\n"
for _k, _v in _metrics.items():
_table_md += f"| {_k} | {_v} |\n"
mo.vstack([
mo.hstack([epi_tr_slider, epi_te_slider, epi_matrix_slider], justify="start", gap=2),
mo.hstack([_fig, mo.md(_table_md)], justify="center"),
mo.callout(
mo.md(
"**Key tradeoffs in fMRI:**\n"
"- **Larger matrix** → better spatial resolution but slower readout, "
"more distortion\n"
"- **Shorter TR** → better temporal resolution but fewer slices\n"
"- **TE ~30ms at 3T** → optimal BOLD sensitivity (matches T₂* of "
"gray matter ~40-50ms)\n"
"- Modern **multiband/simultaneous multi-slice** techniques can "
"excite multiple slices at once, allowing sub-second whole-brain "
"TR with good resolution"
),
kind="warn",
),
])
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## 17. Putting It All Together: From Proton to Activation Map
Let's trace the complete chain from physics to functional brain imaging:
```
Strong magnetic field (B₀)
│
▼
Protons align → net magnetization (M₀)
│
▼
RF pulse at Larmor frequency → tips M into transverse plane
│
▼
Precessing Mxy induces signal in coil (FID)
│
▼
Gradients encode spatial position
│
▼
Raw signal fills k-space
│
▼
2D Inverse FFT → image for each slice
│
▼
Repeat with T₂*-weighted GRE-EPI every ~1s
│
▼
BOLD signal changes reflect neural activity (via hemodynamics)
│
▼
Statistical analysis (GLM) → activation maps
```
**What we gain:**
- Non-invasive measurement of brain activity
- Whole-brain coverage
- Reasonable spatial resolution (~2-3mm)
- No ionizing radiation (unlike PET/CT)
**What we lose / must be aware of:**
- Temporal resolution is limited (~seconds, not milliseconds)
- BOLD is an *indirect* measure of neural activity
- Signal can be affected by head motion, physiological noise,
susceptibility artifacts (especially near sinuses and ear canals)
- The hemodynamic response varies across brain regions and individuals
These limitations -- and how to address them -- are covered in the
subsequent Dartbrains chapters on preprocessing, the GLM, and
statistical analysis.
""")
return
@app.cell(hide_code=True)
def _(mo):
mo.md(r"""
---
## Summary
| Topic | Key Concepts |
|-------|-------------|
| **Magnetism & Resonance** | B₀ alignment, precession, Larmor equation, RF excitation, FID |
| **Signal & Contrast** | T₁/T₂ relaxation, Bloch equations, TE/TR contrast, spin echo vs gradient echo |
| **Imaging & fMRI** | Gradients, spatial encoding, k-space, BOLD signal, HRF, EPI |
These concepts form the foundation for everything else in neuroimaging.
Understanding *why* the signal looks the way it does -- and *what can go
wrong* -- is essential for designing good experiments, preprocessing data
correctly, and interpreting results with appropriate caution.
**Continue your learning:** The subsequent Dartbrains chapters cover
signal processing, preprocessing with fMRIPrep, the General Linear Model,
group analysis, and multivariate pattern analysis.
""")
return
@app.cell
def _():
return
if __name__ == "__main__":
app.run()