Diffusion-induced Signal Attenuation
The purpose of this tutorial is to showcase the simulation of diffusion-related effects. For this, we are going to define a Path motion to simulate the Brownian motion of spins. This is not the most efficient way of simulating diffusion, but it is a good way to understand the phenomenon. In particular, we will going to simulate isotropic diffusion, characterized by the Apparent Diffusion Coefficient (ADC).
Creating a phantom with isotropic diffusion
First we will define a Phantom without motion containing $10,000$ spins. The spins will have the same relaxation times $T_1 = 1000\,\mathrm{ms}$ and $T_2 = 100\,\mathrm{ms}$, and will be placed at the origin.
Nspins = 10_000
obj = Phantom(;
x = zeros(Nspins),
T1 = ones(Nspins) * 1000e-3,
T2 = ones(Nspins) * 100e-3,
)Phantom{Float64}
name: String "spins"
x: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
y: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
z: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
ρ: Array{Float64}((10000,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
T1: Array{Float64}((10000,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 … 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
T2: Array{Float64}((10000,)) [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 … 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1]
T2s: Array{Float64}((10000,)) [1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6 … 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6, 1.0e6]
Δw: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Dλ1: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Dλ2: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Dθ: Array{Float64}((10000,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
motion: NoMotion{Float64} NoMotion{Float64}()
Now we will define the Brownian motion of the spins using the Path motion definition. The motion will be defined by the displacements in the x, y, and z directions (dx, dy, and dz) of the spins. The displacements will be generated by a random walk with mean square displacement
\[\mathbb{E}\left[x^{2}\right]=2D Δt,\]
where $D$ is the diffusion coefficient and $Δt$ is time step duration.
D = 2e-9 # Diffusion Coefficient of water in m^2/s
T = 100e-3 # Duration of the motion
Nt = 100 # Number of time steps
Δt = T / (Nt - 1) # Time sep
Δr = sqrt(2 * D * Δt) # √ Mean square displacement2.0100756305184243e-6Random walk is defined as the cumulative sum of random displacements:
rng = MersenneTwister(1234) # Setting up the random seed
dx = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)
dy = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)
dz = cumsum([zeros(Nspins) Δr .* randn(rng, Nspins, Nt - 1)]; dims=2)10000×100 Matrix{Float64}:
0.0 -1.07467e-6 -9.19486e-7 2.83151e-6 … 1.22958e-5 7.97302e-6
0.0 -2.9865e-6 -1.47145e-6 -2.51876e-6 -4.83153e-5 -4.94091e-5
0.0 2.6588e-7 1.01477e-6 1.02467e-6 -9.73165e-7 -8.43845e-7
0.0 -2.0266e-6 -1.96299e-6 -1.89612e-6 -1.18977e-5 -1.08615e-5
0.0 2.72498e-6 2.0258e-6 2.61778e-6 7.2308e-6 4.26138e-6
0.0 1.18205e-6 9.84498e-7 -1.05401e-6 … 5.74571e-6 9.05569e-6
0.0 -1.85799e-6 7.36487e-8 -2.84953e-6 -1.7914e-5 -1.55996e-5
0.0 7.9729e-7 -1.70114e-6 -4.85136e-6 -2.25009e-6 -6.24905e-7
0.0 -3.81476e-9 9.19291e-7 -1.93225e-7 -1.43333e-5 -1.58779e-5
0.0 -2.96106e-6 -5.7953e-6 -5.33949e-6 -7.37129e-6 -4.45957e-6
⋮ ⋱
0.0 3.06237e-6 1.37898e-6 1.89968e-6 -2.73319e-5 -2.3879e-5
0.0 -7.94671e-7 2.85117e-7 -1.64305e-6 -2.02476e-5 -2.13101e-5
0.0 2.14378e-6 4.02033e-6 5.64978e-6 6.98896e-6 9.42077e-6
0.0 4.24531e-7 2.36026e-6 9.46319e-7 1.2153e-9 2.72995e-6
0.0 3.18803e-8 1.94519e-6 7.07642e-8 … 1.61758e-5 1.46001e-5
0.0 1.83955e-6 3.97057e-6 3.21251e-6 -6.45348e-6 -6.45838e-6
0.0 1.41076e-7 -3.17962e-7 -4.17676e-6 -4.63959e-5 -4.81074e-5
0.0 5.90846e-7 -3.96405e-7 2.7393e-6 3.83404e-8 -1.78516e-7
0.0 -2.96346e-6 -2.00388e-6 1.73694e-6 -7.63724e-6 -4.87867e-6Including the random_walk into the Phantom definition:
random_walk = Path(dx, dy, dz, TimeRange(0.0, T))
obj.motion = MotionList(random_walk)
p1 = plot_phantom_map(obj, :T1; time_samples=Nt÷4, height=450)"../assets/6-displacements.html"The plot shows the random walk of spins due to diffusion, also known as Brownian motion. This motion was named after Robert Brown, who first described the phenomenon in 1827 while looking at pollen suspended in water under a microscope.
Pulse Gradient Spin Echo (PGSE) sequence
The classical sequence used to measure diffusion is the pulse gradient spin echo (PGSE) introduced by Stejskal and Tanner in 1965. This sequence is characterized by the use of two diffusion gradients, placed right before and right after the inversion RF pulse. The duration of each gradient is defined by the δ parameter and the distance between the beginning of both gradients is described by the Δ parameter. In this tutorial square-shaped gradients will be used.
First, we generate the RF pulses:
sys = Scanner()
durRF = 1e-3
B1 = (π / 2) / (2π * γ * durRF)
rf90 = PulseDesigner.RF_hard(B1, durRF, sys)
rf180 = (0.0 + 2im) * rf90Sequence[ τ = 1.0 ms | blocks: 1 | ADC: 0 | GR: 0 | RF: 1 | DEF: 0 ]Now we generate the gradients:
G = 30e-3 # Gradient amplitude
δ = 30e-3 # Duration of the gradient
Δ = durRF + δ # Time between the two gradients
gx_diff = Grad(G, δ)Grad(30.0 mT, 30.0 ms)Finally, we generate the ADC:
adc_dwell_time = 1e-6
adc = ADC(1, adc_dwell_time, durRF/2 - adc_dwell_time/2) # ADCs with N=1 are positioned at the centerADC(1, 1.0e-6, 0.0004995, 0.0, 0.0)Obtaining the PGSE sequence:
seq = Sequence()
seq += rf90
seq += gx_diff
seq += rf180
seq += gx_diff
seq += adc
p2 = plot_seq(seq; show_adc=true) # Plotting the sequence"../assets/6-pgse_sequence.html"For the isotropic diffusion, the signal attenuation is given by the Stejskal-Tanner formula:
\[E = \exp\left(-b D\right)\]
where (b) is the b-value, defined as:
\[b = (\gamma G \delta)^{2} \cdot \left(\Delta - \delta/3\right)\]
where $\gamma$ is the gyromagnetic ratio, $G$ is the gradient amplitude, $\delta$ is the gradient duration, and $\Delta$ is the time between the two gradients.
function bvalue(seq)
block, axis = 2, 1 # Gx from second block
G = seq.GR[axis, block].A
δ = seq.GR[axis, block].T
Δ = dur(seq[2:3]) # Because there are no gaps
b = (2π * γ * G * δ)^2 * (Δ - δ/3)
return b * 1e-6
endbvalue (generic function with 1 method)Diffusion Weighted Imaging (DWI)
For DWI, multiple b-values are acquired to determine the tissue's ADC. For this, we will scale the gradient's amplitude of the previous sequence to obtain a desired b-value. We will store the sequences in a vector seqs and simulate the signal for each one of them.
seqs = Sequence[] # Vector of sequences
bvals = [0, 250, 500, 1000, 1500, 2000] # b-values in s/mm^2
for bval_target in bvals
gradient_scaling = sqrt(bval_target / bvalue(seq))
seq_b = gradient_scaling * seq
push!(seqs, seq_b)
endTo simulate, we will broadcast the simulate function over the sequences and store the signals in a vector Sb. The Ref's are used to avoid broadcasting the obj and sys arguments (they will remain constant for all seqs).
sim_params = KomaMRICore.default_sim_params()
sim_params["return_type"] = "mat"
sim_params["Δt"] = Δt # Set max. grad. time step to fit diffusion time step
signals = simulate.(Ref(obj), seqs, Ref(sys); sim_params) # simulate broadcasted over seqs
Sb = [sb[1] for sb in signals] # Reshaping the simulated signals
bvals_si = bvals .* 1e6 # Convert b-values from s/mm^2 to s/m^2
E_simulated = abs.(Sb) ./ abs.(Sb[1])
E_theoretical = exp.(-bvals_si .* D)"../assets/6-pgse_signal_attenuation.html"The plot shows the signal attenuation as a function of the b-value. The simulated signal attenuation matches the theoretical curve, showing the expected exponential decay with the b-value.
This page was generated using Literate.jl.