Author Archives: Steven Whitaker

Mastering MRI Bloch Simulations with BlochSim.jl in Julia

By: Steven Whitaker

Re-posted from: https://glcs.hashnode.dev/mastering-mri-bloch-simulations-with-blochsimjl-in-julia

Julia is a relatively new programming language that excels especially in scientific computing.

In addition to the functionality provided by the language itself, Julia boasts a large and growing collection of packages that can be seamlessly installed and loaded and that provide even more functionality.

In this post, we will learn how to use BlochSim.jl, a Julia package that provides functionality for running MRI Bloch simulations.

Note that this post assumes a basic understanding of the Bloch equations and MRI in general. See our previous post for a very brief overview of the Bloch equations (and a walkthrough of how to write your own Bloch simulation code in Julia).

Installing Julia

You can head over to Julia’s website to download the language on your own computer, or you can try out Julia in your web browser first.

Installing and Loading BlochSim.jl

We need to install BlochSim.jl. We can do so by running the following from the Julia prompt (also known as the REPL):

julia> using Pkg; Pkg.add("BlochSim")

Then we can load BlochSim.jl into our current Julia session.

julia> using BlochSim

Great, BlochSim.jl is installed and loaded, so now we can start learning how to use it!

The Bloch equations manipulate magnetization vectors, so first we will learn how BlochSim.jl represents a magnetization vector.

Spin: Representing a Magnetization Vector

BlochSim.jl uses Spins to represent magnetization vectors. More accurately, a Spin represents an isochromat, i.e., a magnetization vector with associated values for M0, T1, T2, and off-resonance frequency.

Thus, to construct a Spin, we need to pass those values to the constructor.

julia> (M0, T1, T2, df) = (1, 1000, 80, 100)(1, 1000, 80, 100)julia> spin = Spin(M0, T1, T2, df)Spin{Float64}: M = Magnetization(0.0, 0.0, 1.0) M0 = 1.0 T1 = 1000.0 ms T2 = 80.0 ms f = 100.0 Hz pos = Position(0.0, 0.0, 0.0) cm

Some notes:

  • spin.M is the magnetization vector associated with the isochromat. This vector is what will be manipulated by the Bloch equations.

  • spin.M by default starts in thermal equilibrium.

  • We can grab the x-component of the magnetization with spin.M.x or getproperty(spin.M, :x) (similarly for y and z).

  • T1 and T2 are given in units of ms, while off-resonance uses Hz.

  • The spin optionally can be given a position (useful if simulating B0 gradients).

Now that we have a spin, we can manipulate it with the Bloch equations.

Running a Bloch Simulation

We will simulate a 1 ms excitation pulse followed by 9 ms of free precession with a time step of 10 s.

julia> (t_ex, t_fp, dt) = (1, 9, 0.01)(1, 9, 0.01)julia> t_total = t_ex + t_fp10

We will want to plot the magnetization over time, so we will need to store the magnetization vector at each time step of the simulation. We can create a Vector of Magnetizations with enough pre-allocated storage for the simulation.

julia> nsteps = round(Int, t_total / dt)1000julia> M = Vector{Magnetization{Float64}}(undef, nsteps + 1)1001-element Vector{Magnetization{Float64}}: #undef  #undef

RF: Representing an RF Excitation Pulse

BlochSim.jl represents excitation pulses with the RF object. The constructor takes in the RF waveform (in units of Gauss) and the time step (in ms) to use (and optionally a phase offset and/or a B0 gradient).

For the excitation pulse, we will use a sinc pulse that gives a flip angle of 90. First we will specify the shape.

julia> nrf = round(Int, t_ex / dt)100julia> waveform = sinc.(range(-3, 3, nrf));

(Note the dot (.) when calling sinc. This calls sinc, a scalar function, on each element of the input array. See our blog post on broadcasting for more information.)

And now we need to scale our waveform to get the desired flip angle.

julia> flip = /21.5707963267948966julia> waveform .*= flip / sum(waveform) / (GAMMA * dt / 1000);

Now that we have our waveform we can construct an RF object.

julia> rf_full = RF(waveform, dt)RF{Float64,Gradient{Int64}}:  = [0.0, 0.001830277817402092  0.001830277817402092, 0.0] rad  = [0.0, 0.0  0.0, 0.0] rad t = 0.01 ms  (initial) = 0.0 rad  (current) = 0.0 rad grad = Gradient(0, 0, 0) G/cm

Using this RF object would give us the magnetization vector at the end of excitation, but wouldn’t allow us to plot the magnetization during excitation. So instead, we will pretend we have nrf separate excitations.

julia> rf = [RF([x], dt) for x in waveform]100-element Vector{RF{Float64, Gradient{Int64}}}: RF([0.0], [0.0], 0.01, 0.0, Gradient(0, 0, 0)) RF([0.001830277817402092], [0.0], 0.01, 0.0, Gradient(0, 0, 0))  RF([0.0], [0.0], 0.01, 0.0, Gradient(0, 0, 0))

(Note that above we used a comprehension.)

excite and freeprecess

BlochSim.jl provides two main Bloch simulation functions: excite for simulating excitation and freeprecess for simulating free precession. They each return a matrix A and a vector B such that A * M + B applies the dynamics to a magnetization vector M.

We will use these functions in our simulation.

julia> M[1] = copy(spin.M);julia> for i = 1:nsteps           if i <= nrf               (A, B) = excite(spin, rf[i])           else               (A, B) = freeprecess(spin, dt)           end           M[i+1] = A * M[i] + B       end

(Note that if we used M[1] = spin.M above, i.e., without copying, any modifications to spin.M would also change M[1], and vice versa, because they both would refer to the same computer memory location.)

And with that, we have successfully run a Bloch simulation!

Now we want to visualize the results.

Simulation Results

To plot the results, we will use the Plots.jl package, installed via

julia> using Pkg; Pkg.add("Plots")

Our plot will have two subplots. The first will show the RF waveform, and the second will show the x-, y-, and z-components of the magnetization vector over time.

julia> using Plotsjulia> t = (0:nsteps) .* dt0.0:0.01:10.0julia> rf_plot = [0; waveform; zeros(nsteps - nrf)];julia> prf = plot(t, rf_plot; label = "RF", xlabel = "t (ms)", ylabel = "Amplitude (G)");julia> Mx = getproperty.(M, :x); My = getproperty.(M, :y); Mz = getproperty.(M, :z);julia> pM = plot(t, Mx; label = "Mx", xlabel = "t (ms)", ylabel = "Magnetization");julia> plot!(pM, t, My; label = "My");julia> plot!(pM, t, Mz; label = "Mz");julia> plot(prf, pM; layout = (2, 1))

And there you have it! We can see the excitation occur followed by off-resonance precession. If we simulated free precession for longer we would also see T1 recovery and T2 decay.

Running a Bloch Simulation Without Plotting

But what if you don’t want to plot because you only care about the magnetization at the end of the simulation?

BlochSim.jl was made with this use case in mind.

Let’s repeat the Bloch simulation above, but this time we want only the magnetization at the end of the simulation.

We will reuse the following variables:

  • spin (Spin object with the initial magnetization vector and tissue parameters)

  • t_fp (free precession duration)

  • rf_full (RF object of duration 1 ms with a time step of 10 s and a flip angle of 90)

Since we don’t need to store intermediate magnetization vectors, we can operate directly on spin.M, reusing it instead of creating new Magnetization objects. To do this, we will use the functions excite! and freeprecess!.

Note the ! at the end of the above function names; this is the Julia convention for indicating a function mutates, or modifies, one or more of its inputs (typically the first input).

The functions excite! and freeprecess! take the same inputs as their non-mutating counterparts. But instead of returning a matrix A and a vector B, they modify spin.M directly.

Also, since we aren’t storing intermediate magnetization vectors, instead of looping for nsteps steps, we call excite! and freeprecess! just one time each.

julia> excite!(spin, rf_full)julia> freeprecess!(spin, t_fp)

Now spin.M is the magnetization at the end of the Bloch simulation.

julia> spin.MMagnetization vector with eltype Float64: Mx = 0.8506279504983056 My = 0.25506720555905055 Mz = 0.015720687775520957

And if we want the complex-valued MRI signal, we just use the signal function.

julia> signal(spin)0.8506279504983056 + 0.25506720555905055im

We can also verify that the two approaches give the same result (up to floating point rounding errors) at the end of the simulation.

julia> M[end]  spin.M # To type , type \approx<tab>true

So now we can run Bloch simulations in two ways:

  1. where we keep track of the intermediate magnetization vectors for plotting purposes, and

  2. where we call excite! and freeprecess! one time each to compute the magnetization only at the end of the simulation.

Summary

In this post, we have seen how to use BlochSim.jl to run Bloch simulations. In particular, we learned how to construct a Spin object and an RF object, and we manipulated magnetization vectors using excite/excite! and freeprecession/freeprecession!.

Additional Notes about BlochSim.jl

  • excite! and freeprecess! could have been used in the plotting version of the Bloch simulation. Then, instead of M[i+1] = A * M[i] + B, we would need M[i+1] = copy(spin.M).

  • InstantaneousRF can be used in place of RF if one wants to assume the RF pulse is very fast compared to T1/T2/off-resonance effects. (This is a typical assumption.)

  • For advanced MRI users/researchers: BlochSim.jl provides a SpinMC object for simulating multiple compartments using the Bloch-McConnell equations.

Additional Links

Simulating MRI Physics with the Bloch Equations

By: Steven Whitaker

Re-posted from: https://blog.glcs.io/simulating-mri-physics-with-the-bloch-equations

MRI is an important imaging modality with impressive diagnostic capabilities.

The equations that specify the physics of how MRI works are called the Bloch equations.

People who study MRI use the Bloch equations to simulate MRI scans. These simulations are called Bloch simulations.

In this post, we will learn how to simulate MRI physics in the Julia programming language, a free and open source programming language that excels especially in scientific computing.

Note that this post assumes a basic understanding of the Bloch equations and MRI in general.

Installing Julia

You can head over to Julia’s website to download the language on your own computer, or you can try out Julia in your web browser first.

Notation

Here is the mathematical notation used in this post.

  • \( \mathbf{M} = [M_x, M_y, M_z] \) : Magnetization vector

  • \( M_{xy} = M_x + i M_y \) : Transverse magnetization ( \( i = \sqrt{-1} \) )

  • \( \mathbf{M}_0 = [0, 0, M_0] \) : Equilibrium magnetization vector

  • \( M_0 \) : Equilibrium magnetization

  • \( T_1 \) : T1 time constant

  • \( T_2 \) : T2 time constant

  • \( \Delta\omega \) : Off-resonance frequency

  • \( \omega_{1,x} \) : Rotational frequency due to the excitation RF field (along x)

  • \( \omega_{1,y} \) : Rotational frequency due to the excitation RF field (along y)

Bloch Equations

MRI scans manipulate magnetization vectors, and the Bloch equations explain the physics of how that happens. The system of ordinary differential equations (ODE) that make up the Bloch equations is as follows.

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} -\frac{1}{T_2} & \Delta\omega & -\omega_{1,y}(t) \\ -\Delta\omega & -\frac{1}{T_2} & \omega_{1,x}(t) \\ \omega_{1,y}(t) & -\omega_{1,x}(t) & -\frac{1}{T_1} \end{bmatrix} \mathbf{M}(t) + \begin{bmatrix} 0 \\ 0 \\ \frac{M_0}{T_1} \end{bmatrix}\]

Normally in MRI, however, we make important assumptions that lead to closed-form solutions.

Excitation

One assumption is that the RF excitation pulses are fast (instantaneous) relative to all the other dynamics. This leads to

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & \omega_{1,x}(t) \\ 0 & -\omega_{1,x}(t) & 0 \end{bmatrix} \mathbf{M}(t).\]

(For simplicity, we assume the excitation pulse is aligned with the x-axis.)

The solution to this ODE is rotation, i.e.,

\[\mathbf{M}(t_0^{+}) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\alpha) & \sin(\alpha) \\ 0 & -\sin(\alpha) & \cos(\alpha) \end{bmatrix} \mathbf{M}(t_0^{-}),\]

where \( \alpha \) is the flip angle.

Julia Code: Excitation

To simulate excitation in Julia, we need a magnetization vector and a rotation matrix. Let’s say we start with \( \mathbf{M}(t_0^{-}) = [0, 0, 1] \) , i.e., magnetization along the z-axis. In Julia, we create a variable M and assign it the value of [0, 0, 1]. Here’s how it looks in the Julia prompt (REPL):

julia> M = [0, 0, 1]
3-element Vector{Int64}:
 0
 0
 1

Next we need the rotation matrix. First, we specify the flip angle, and then we make the rotation matrix.

julia> alpha =  / 2 # To type , type \pi<tab>
1.5707963267948966

julia> R = [0 0 0; 0 cos(alpha) sin(alpha); 0 -sin(alpha) cos(alpha)]
3x3 Matrix{Float64}:
 0.0   0.0          0.0
 0.0   6.12323e-17  1.0
 0.0  -1.0          6.12323e-17

(Note that the result of cos(alpha) is 6.12323e-17, not 0.0, because of floating point rounding errors, i.e., a computer can’t exactly represent \( \frac{\pi}{2} \) .)

Finally, the excitation can be simulated by multiplying the magnetization by the rotation matrix.

julia> M_ex = R * M
3-element Vector{Float64}:
 0.0
 1.0
 6.123233995736766e-17

We started with magnetization along the z-axis, applied a 90 excitation pulse along the x-axis, and ended up with magnetization along the y-axis as expected.

We have officially simulated MRI physics!

But excitation is just part of the physics, so let’s move on to free precession.

Free Precession

Free precession is what occurs when the RF excitation pulse is turned off. In this case, the Bloch equations become

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} -\frac{1}{T_2} & \Delta\omega & 0 \\ -\Delta\omega & -\frac{1}{T_2} & 0 \\ 0 & 0 & -\frac{1}{T_1} \end{bmatrix} \mathbf{M}(t) + \begin{bmatrix} 0 \\ 0 \\ \frac{M_0}{T_1} \end{bmatrix}.\]

The solution to this ODE is free precession, which includes the effects of off-resonance precession, T1 recovery, and T2 decay.

\[\begin{align} \mathbf{M}(t – t_0) &= \mathbf{M}_0 + \mathbf{E}(t – t_0) \mathbf{F}(t – t_0) (\mathbf{M}(t_0) – \mathbf{M}_0) \\ \mathbf{E}(t) &= \begin{bmatrix} e^{-t/T_2} & 0 & 0 \\ 0 & e^{-t/T_2} & 0 \\ 0 & 0 & e^{-t/T_1} \end{bmatrix} \\ \mathbf{F}(t) &= \begin{bmatrix} \cos(\Delta\omega t) & \sin(\Delta\omega t) & 0 \\ -\sin(\Delta\omega t) & \cos(\Delta\omega t) & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{align}\]

Julia Code: Free Precession

To simulate free precession in Julia, we need a magnetization vector, and we need to specify some constants. Let’s say we start with \( \mathbf{M}(t_0) = [0, 1, 0] \) , i.e., magnetization along the y-axis, and let’s use the following constant values:

  • \( M_0 = 1 \)

  • \( T_1 = 1 \) second

  • \( T_2 = 0.1 \) seconds

  • \( \Delta\omega = 25\pi \) rad/s

  • \( t – t_0 = 0.01 \) seconds

We can create M and variables for the constants and equilibrium magnetization.

julia> M = [0, 1, 0]
3-element Vector{Int64}:
 0
 1
 0

julia> (M0, T1, T2, dw, dt) = (1, 1, 0.1, 25, 0.01)
(1, 1, 0.1, 78.53981633974483, 0.01)

julia> M0_vec = [0, 0, M0]
3-element Vector{Int64}:
 0
 0
 1

Next we can create the \( \mathbf{E} \) and \( \mathbf{F} \) matrices from above.

julia> E = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 exp(-dt/T1)]
3x3 Matrix{Float64}:
 0.904837  0.0       0.0
 0.0       0.904837  0.0
 0.0       0.0       0.99005

julia> phi = dw * dt; F = [cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]
3x3 Matrix{Float64}:
  0.707107  0.707107  0.0
 -0.707107  0.707107  0.0
  0.0       0.0       1.0

Finally, the free precession can be simulated.

julia> M_fp = M0_vec + E * F * (M - M0_vec)
3-element Vector{Float64}:
 0.6398166741645539
 0.639816674164554
 0.009950166250831893

Awesome, now we can simulate excitation and free precession, two major building blocks of MRI physics!

We can use these building blocks to perform a more complete simulation of MRI physics.

A More Complete Simulation

Let’s simulate how a magnetization vector evolves over time.

We will simulate a 1 ms excitation pulse followed by 9 ms of free precession with a time step of 10 s.

julia> (t_ex, t_fp, dt) = (0.001, 0.009, 0.00001)
(0.001, 0.009, 1.0e-5)

julia> t_total = t_ex + t_fp
0.009999999999999998

In this case, however, we don’t want to assume the excitation is instantaneous. Instead, we will split the RF pulse into segments of width \( \Delta t \) . For each segment, we will

  • simulate free precession for a duration of \( \Delta t / 2 \) ,

  • simulate an instantaneous RF pulse whose flip angle is determined by the current value of the RF waveform, and

  • again simulate free precession for a duration of \( \Delta t / 2 \) .

Thus we can still use our building blocks while simulating a more accurate excitation.

We will want to plot the magnetization over time, so we will need to store the magnetization vector at each time step of the simulation. We can create a Vector of Vectors with enough pre-allocated storage for the simulation.

julia> nsteps = round(Int, t_total / dt)
1000

julia> M = Vector{Vector{Float64}}(undef, nsteps + 1)
1001-element Vector{Vector{Float64}}:
 #undef
 
 #undef

We also need to specify the RF pulse. We will make a sinc pulse that gives a flip angle of 90.

julia> nrf = round(Int, t_ex / dt)
100

julia> rf = sinc.(range(-3, 3, nrf));

Above, range(-3, 3, nrf) creates nrf equally spaced points from -3 to 3, inclusive. And sinc is a scalar function, so we broadcast it, or apply it elementwise, by adding the dot (.) after the function name. The flip angle is not correct, so we need to scale rf.

julia> rf .*= /2 / sum(rf);

Finally, let’s start with magnetization in equilibrium and specify the other constants.

julia> (M0, T1, T2, dw) = (1, 1, 0.1, 25)
(1, 1, 0.1, 78.53981633974483)

julia> M[1] = [0, 0, M0]
3-element Vector{Int64}:
 0
 0
 1

To aid with the simulation, we will make a function to do excitation and another to do free precession.

julia> function excite(M, alpha, M0, T1, T2, dw, dt)

           M_next = freeprecess(M, M0, T1, T2, dw, dt/2)
           R = [0 0 0; 0 cos(alpha) sin(alpha); 0 -sin(alpha) cos(alpha)]
           M_next = R * M_next
           M_next = freeprecess(M_next, M0, T1, T2, dw, dt/2)

           return M_next

       end
excite (generic function with 1 method)

julia> function freeprecess(M, M0, T1, T2, dw, dt)

           M0_vec = [0, 0, M0]
           E = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 exp(-dt/T1)]
           phi = dw * dt
           F = [cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]
           M_next = M0_vec + E * F * (M - M0_vec)

           return M_next

       end
freeprecess (generic function with 1 method)

Now we can run the simulation.

julia> for i = 1:nsteps
           if i <= nrf
               M[i+1] = excite(M[i], rf[i], M0, T1, T2, dw, dt)
           else
               M[i+1] = freeprecess(M[i], M0, T1, T2, dw, dt)
           end
       end

And with that, we have successfully simulated a realistic excitation pulse followed by free precession.

Now we want to visualize the results.

Simulation Results

To plot the results, we will use the Plots.jl package. Our plot will have two subplots. The first will show the RF waveform, and the second will show \( M_x \) , \( M_y \) , and \( M_z \) .

julia> using Plots # Load Plots.jl package

julia> t = (0:nsteps) .* (dt * 1000)
0.0:0.01:10.0

julia> rf = [0; rf; zeros(nsteps - nrf)];

julia> prf = plot(t, rf; label = "RF", xlabel = "t (ms)", ylabel = "Flip Angle (rad)");

julia> Mx = getindex.(M, 1); My = getindex.(M, 2); Mz = getindex.(M, 3);

julia> pM = plot(t, Mx; label = "Mx", xlabel = "t (ms)", ylabel = "Magnetization");

julia> plot!(pM, t, My; label = "My");

julia> plot!(pM, t, Mz; label = "Mz");

julia> plot(prf, pM; layout = (2, 1))

We can also plot the magnitude and phase of the transverse component of the magnetization (i.e., \( M_{xy} = M_x + i M_y \) ).

julia> Mxy = complex.(Mx, My);

julia> pmag = plot(t, abs.(Mxy); label = "|Mxy|", xlabel = "t (ms)", ylabel = "Magnitude", color = 1);

# To type , type \angle<tab>
julia> pphase = plot(t, angle.(Mxy); label = "Mxy", xlabel = "t (ms)", ylabel = "Phase (rad)", color = 2);

julia> plot(pmag, pphase; layout = (2, 1))

Finally, if we want to see the magnetization recover to (almost) equilibrium, we can increase t_fp (the free precession duration) to, say, 3 * T1. To reduce the number of time steps, it is advisable to use a larger time step for free precession than for excitation (the results below used a 1 ms time step for free precession).

And now we can easily see the characteristic T2 decay, off-resonance precession, and T1 recovery of the magnetization vector, indicating we have done a successful Bloch simulation!

Summary

In this post, we have reviewed the Bloch equations and seen how to implement excitation and free precession in Julia. We have also seen how to run a Bloch simulation that tracks a magnetization vector over time, starting with an excitation pulse and continuing with free precession.

Additional Links

Simulating MRI Physics with the Bloch Equations

By: Steven Whitaker

Re-posted from: https://glcs.hashnode.dev/simulating-mri-physics-with-the-bloch-equations

MRI is an important imaging modality with impressive diagnostic capabilities.

The equations that specify the physics of how MRI works are called the Bloch equations.

People who study MRI use the Bloch equations to simulate MRI scans. These simulations are called Bloch simulations.

In this post, we will learn how to simulate MRI physics in the Julia programming language, a free and open source programming language that excels especially in scientific computing.

Note that this post assumes a basic understanding of the Bloch equations and MRI in general.

Installing Julia

You can head over to Julia’s website to download the language on your own computer, or you can try out Julia in your web browser first.

Notation

Here is the mathematical notation used in this post.

  • \( \mathbf{M} = [M_x, M_y, M_z] \) : Magnetization vector

  • \( M_{xy} = M_x + i M_y \) : Transverse magnetization ( \( i = \sqrt{-1} \) )

  • \( \mathbf{M}_0 = [0, 0, M_0] \) : Equilibrium magnetization vector

  • \( M_0 \) : Equilibrium magnetization

  • \( T_1 \) : T1 time constant

  • \( T_2 \) : T2 time constant

  • \( \Delta\omega \) : Off-resonance frequency

  • \( \omega_{1,x} \) : Rotational frequency due to the excitation RF field (along x)

  • \( \omega_{1,y} \) : Rotational frequency due to the excitation RF field (along y)

Bloch Equations

MRI scans manipulate magnetization vectors, and the Bloch equations explain the physics of how that happens. The system of ordinary differential equations (ODE) that make up the Bloch equations is as follows.

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} -\frac{1}{T_2} & \Delta\omega & -\omega_{1,y}(t) \\ -\Delta\omega & -\frac{1}{T_2} & \omega_{1,x}(t) \\ \omega_{1,y}(t) & -\omega_{1,x}(t) & -\frac{1}{T_1} \end{bmatrix} \mathbf{M}(t) + \begin{bmatrix} 0 \\ 0 \\ \frac{M_0}{T_1} \end{bmatrix}\]

Normally in MRI, however, we make important assumptions that lead to closed-form solutions.

Excitation

One assumption is that the RF excitation pulses are fast (instantaneous) relative to all the other dynamics. This leads to

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & \omega_{1,x}(t) \\ 0 & -\omega_{1,x}(t) & 0 \end{bmatrix} \mathbf{M}(t).\]

(For simplicity, we assume the excitation pulse is aligned with the x-axis.)

The solution to this ODE is rotation, i.e.,

\[\mathbf{M}(t_0^{+}) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\alpha) & \sin(\alpha) \\ 0 & -\sin(\alpha) & \cos(\alpha) \end{bmatrix} \mathbf{M}(t_0^{-}),\]

where \( \alpha \) is the flip angle.

Julia Code: Excitation

To simulate excitation in Julia, we need a magnetization vector and a rotation matrix. Let’s say we start with \( \mathbf{M}(t_0^{-}) = [0, 0, 1] \) , i.e., magnetization along the z-axis. In Julia, we create a variable M and assign it the value of [0, 0, 1]. Here’s how it looks in the Julia prompt (REPL):

julia> M = [0, 0, 1]3-element Vector{Int64}: 0 0 1

Next we need the rotation matrix. First, we specify the flip angle, and then we make the rotation matrix.

julia> alpha =  / 2 # To type , type \pi<tab>1.5707963267948966julia> R = [0 0 0; 0 cos(alpha) sin(alpha); 0 -sin(alpha) cos(alpha)]3x3 Matrix{Float64}: 0.0   0.0          0.0 0.0   6.12323e-17  1.0 0.0  -1.0          6.12323e-17

(Note that the result of cos(alpha) is 6.12323e-17, not 0.0, because of floating point rounding errors, i.e., a computer can’t exactly represent \( \frac{\pi}{2} \) .)

Finally, the excitation can be simulated by multiplying the magnetization by the rotation matrix.

julia> M_ex = R * M3-element Vector{Float64}: 0.0 1.0 6.123233995736766e-17

We started with magnetization along the z-axis, applied a 90 excitation pulse along the x-axis, and ended up with magnetization along the y-axis as expected.

We have officially simulated MRI physics!

But excitation is just part of the physics, so let’s move on to free precession.

Free Precession

Free precession is what occurs when the RF excitation pulse is turned off. In this case, the Bloch equations become

\[\frac{d}{dt} \mathbf{M}(t) = \begin{bmatrix} -\frac{1}{T_2} & \Delta\omega & 0 \\ -\Delta\omega & -\frac{1}{T_2} & 0 \\ 0 & 0 & -\frac{1}{T_1} \end{bmatrix} \mathbf{M}(t) + \begin{bmatrix} 0 \\ 0 \\ \frac{M_0}{T_1} \end{bmatrix}.\]

The solution to this ODE is free precession, which includes the effects of off-resonance precession, T1 recovery, and T2 decay.

\[\begin{align} \mathbf{M}(t – t_0) &= \mathbf{M}_0 + \mathbf{E}(t – t_0) \mathbf{F}(t – t_0) (\mathbf{M}(t_0) – \mathbf{M}_0) \\ \mathbf{E}(t) &= \begin{bmatrix} e^{-t/T_2} & 0 & 0 \\ 0 & e^{-t/T_2} & 0 \\ 0 & 0 & e^{-t/T_1} \end{bmatrix} \\ \mathbf{F}(t) &= \begin{bmatrix} \cos(\Delta\omega t) & \sin(\Delta\omega t) & 0 \\ -\sin(\Delta\omega t) & \cos(\Delta\omega t) & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{align}\]

Julia Code: Free Precession

To simulate free precession in Julia, we need a magnetization vector, and we need to specify some constants. Let’s say we start with \( \mathbf{M}(t_0) = [0, 1, 0] \) , i.e., magnetization along the y-axis, and let’s use the following constant values:

  • \( M_0 = 1 \)

  • \( T_1 = 1 \) second

  • \( T_2 = 0.1 \) seconds

  • \( \Delta\omega = 25\pi \) rad/s

  • \( t – t_0 = 0.01 \) seconds

We can create M and variables for the constants and equilibrium magnetization.

julia> M = [0, 1, 0]3-element Vector{Int64}: 0 1 0julia> (M0, T1, T2, dw, dt) = (1, 1, 0.1, 25, 0.01)(1, 1, 0.1, 78.53981633974483, 0.01)julia> M0_vec = [0, 0, M0]3-element Vector{Int64}: 0 0 1

Next we can create the \( \mathbf{E} \) and \( \mathbf{F} \) matrices from above.

julia> E = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 exp(-dt/T1)]3x3 Matrix{Float64}: 0.904837  0.0       0.0 0.0       0.904837  0.0 0.0       0.0       0.99005julia> phi = dw * dt; F = [cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]3x3 Matrix{Float64}:  0.707107  0.707107  0.0 -0.707107  0.707107  0.0  0.0       0.0       1.0

Finally, the free precession can be simulated.

julia> M_fp = M0_vec + E * F * (M - M0_vec)3-element Vector{Float64}: 0.6398166741645539 0.639816674164554 0.009950166250831893

Awesome, now we can simulate excitation and free precession, two major building blocks of MRI physics!

We can use these building blocks to perform a more complete simulation of MRI physics.

A More Complete Simulation

Let’s simulate how a magnetization vector evolves over time.

We will simulate a 1 ms excitation pulse followed by 9 ms of free precession with a time step of 10 s.

julia> (t_ex, t_fp, dt) = (0.001, 0.009, 0.00001)(0.001, 0.009, 1.0e-5)julia> t_total = t_ex + t_fp0.009999999999999998

In this case, however, we don’t want to assume the excitation is instantaneous. Instead, we will split the RF pulse into segments of width \( \Delta t \) . For each segment, we will

  • simulate free precession for a duration of \( \Delta t / 2 \) ,

  • simulate an instantaneous RF pulse whose flip angle is determined by the current value of the RF waveform, and

  • again simulate free precession for a duration of \( \Delta t / 2 \) .

Thus we can still use our building blocks while simulating a more accurate excitation.

We will want to plot the magnetization over time, so we will need to store the magnetization vector at each time step of the simulation. We can create a Vector of Vectors with enough pre-allocated storage for the simulation.

julia> nsteps = round(Int, t_total / dt)1000julia> M = Vector{Vector{Float64}}(undef, nsteps + 1)1001-element Vector{Vector{Float64}}: #undef  #undef

We also need to specify the RF pulse. We will make a sinc pulse that gives a flip angle of 90.

julia> nrf = round(Int, t_ex / dt)100julia> rf = sinc.(range(-3, 3, nrf));

Above, range(-3, 3, nrf) creates nrf equally spaced points from -3 to 3, inclusive. And sinc is a scalar function, so we broadcast it, or apply it elementwise, by adding the dot (.) after the function name. The flip angle is not correct, so we need to scale rf.

julia> rf .*= /2 / sum(rf);

Finally, let’s start with magnetization in equilibrium and specify the other constants.

julia> (M0, T1, T2, dw) = (1, 1, 0.1, 25)(1, 1, 0.1, 78.53981633974483)julia> M[1] = [0, 0, M0]3-element Vector{Int64}: 0 0 1

To aid with the simulation, we will make a function to do excitation and another to do free precession.

julia> function excite(M, alpha, M0, T1, T2, dw, dt)           M_next = freeprecess(M, M0, T1, T2, dw, dt/2)           R = [0 0 0; 0 cos(alpha) sin(alpha); 0 -sin(alpha) cos(alpha)]           M_next = R * M_next           M_next = freeprecess(M_next, M0, T1, T2, dw, dt/2)           return M_next       endexcite (generic function with 1 method)julia> function freeprecess(M, M0, T1, T2, dw, dt)           M0_vec = [0, 0, M0]           E = [exp(-dt/T2) 0 0; 0 exp(-dt/T2) 0; 0 0 exp(-dt/T1)]           phi = dw * dt           F = [cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]           M_next = M0_vec + E * F * (M - M0_vec)           return M_next       endfreeprecess (generic function with 1 method)

Now we can run the simulation.

julia> for i = 1:nsteps           if i <= nrf               M[i+1] = excite(M[i], rf[i], M0, T1, T2, dw, dt)           else               M[i+1] = freeprecess(M[i], M0, T1, T2, dw, dt)           end       end

And with that, we have successfully simulated a realistic excitation pulse followed by free precession.

Now we want to visualize the results.

Simulation Results

To plot the results, we will use the Plots.jl package. Our plot will have two subplots. The first will show the RF waveform, and the second will show \( M_x \) , \( M_y \) , and \( M_z \) .

julia> using Plots # Load Plots.jl packagejulia> t = (0:nsteps) .* (dt * 1000)0.0:0.01:10.0julia> rf = [0; rf; zeros(nsteps - nrf)];julia> prf = plot(t, rf; label = "RF", xlabel = "t (ms)", ylabel = "Flip Angle (rad)");julia> Mx = getindex.(M, 1); My = getindex.(M, 2); Mz = getindex.(M, 3);julia> pM = plot(t, Mx; label = "Mx", xlabel = "t (ms)", ylabel = "Magnetization");julia> plot!(pM, t, My; label = "My");julia> plot!(pM, t, Mz; label = "Mz");julia> plot(prf, pM; layout = (2, 1))

We can also plot the magnitude and phase of the transverse component of the magnetization (i.e., \( M_{xy} = M_x + i M_y \) ).

julia> Mxy = complex.(Mx, My);julia> pmag = plot(t, abs.(Mxy); label = "|Mxy|", xlabel = "t (ms)", ylabel = "Magnitude", color = 1);# To type , type \angle<tab>julia> pphase = plot(t, angle.(Mxy); label = "Mxy", xlabel = "t (ms)", ylabel = "Phase (rad)", color = 2);julia> plot(pmag, pphase; layout = (2, 1))

Finally, if we want to see the magnetization recover to (almost) equilibrium, we can increase t_fp (the free precession duration) to, say, 3 * T1. To reduce the number of time steps, it is advisable to use a larger time step for free precession than for excitation (the results below used a 1 ms time step for free precession).

And now we can easily see the characteristic T2 decay, off-resonance precession, and T1 recovery of the magnetization vector, indicating we have done a successful Bloch simulation!

Summary

In this post, we have reviewed the Bloch equations and seen how to implement excitation and free precession in Julia. We have also seen how to run a Bloch simulation that tracks a magnetization vector over time, starting with an excitation pulse and continuing with free precession.

Additional Links