Author Archives: Steven Whitaker

MRI Scan Simulation Made Easy with BlochSim.jl

By: Steven Whitaker

Re-posted from: https://blog.glcs.io/mri-scan-simulation-made-easy-with-blochsimjl

MRI scans manipulate magnetization vectors to produce an observable signal.

Often it can be useful to simulate what signal an MRI scan will generate. Example scenarios include simulating MRI signals for educational purposes, for developing a new MRI scan, or for fitting model parameters to an acquired signal.

In this post, we will learn how to use the Julia programming language, a free and open source programming language that excels especially in scientific computing, to simulate the signal generated by an MRI scan. Specifically, we will use the BlochSim.jl package to simulate a balanced steady-state free precession (bSSFP) scan.

Note that this post assumes a basic understanding of MRI pulse sequences and Bloch simulations. See our previous posts for an overview of the Bloch equations and an overview of BlochSim.jl.

Prerequisites

To run the code in this post, both Julia and BlochSim.jl need to be installed. See our previous post for installation instructions.

bSSFP

First, we need to review bSSFP so we know how to simulate it.

Here is a pulse sequence diagram of bSSFP (source).

For this post, we are not simulating signal localization, so we will ignore the imaging gradients.

Thus, for each TR we will simulate

  1. excitation,

  2. free precession until the echo time, and

  3. free precession until the end of the TR.

This, by itself, would be essentially the same simulation as done in our previous post. However, two aspects of bSSFP make simulating it a bit more interesting.

  1. bSSFP, as the name states, is a steady-state scan. This means we don’t simulate just one excitation followed by free precession to the echo time. Instead, we need to simulate multiple TRs until we get to steady-state. We can do this either by manually simulating multiple TRs or by mathematically computing the steady-state magnetization. More on this later.

  2. bSSFP utilizes phase-cycling, where the phase of the RF pulse is incremented after each TR.

Therefore, our goal is to write a Julia function (that we will call bssfp) that computes the steady-state signal given by a bSSFP scan that uses phase-cycling.

Throughout this post, we will be editing a text file called bssfp.jl that will load BlochSim.jl and define our bssfp function.

Function Outline

We will begin with an outline of what our bssfp function needs to do.

# bssfp.jl

using BlochSim

# Step 0: Specify the inputs to the function.
function bssfp(inputs)

    # Step 1: Create a `Spin` object.
    # Step 2: Create the RF pulse.
    # Step 3: Compute the steady-state signal.

end

First, we need to specify what inputs bssfp needs. We need values for M0, T1, T2, and off-resonance frequency \( \Delta f \) to create a Spin object, a flip angle \( \alpha \) to create the RF pulse, and the TR, TE, and phase-cycling factor \( \Delta\phi \) to use. We will also add an input that specifies how to compute the steady-state magnetization.

# bssfp.jl

using BlochSim

function bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)

    # Step 1: Create a `Spin` object.
    # Step 2: Create the RF pulse.
    # Step 3: Compute the steady-state signal.

end

Here we see that Julia allows using Greek letters in variable names, which helps the code look more like the math, enhancing readability. (To type, e.g., , type \alpha and then press TAB, and similarly for other Greek letters.)

We also made manual_ss an optional argument. If the user does not specify manual_ss, it will default to false.

Now we can focus on writing code for each step outlined above.

Step 1: Create a Spin Object

As discussed in our previous post, a Spin object is created as follows.

spin = Spin(M0, T1, T2, f)

That’s it!

So let’s update bssfp.

# bssfp.jl

using BlochSim

function bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)

    # Step 1: Create a `Spin` object.
    spin = Spin(M0, T1, T2, f)

    # Step 2: Create the RF pulse.
    # Step 3: Compute the steady-state signal.

end

Step 2: Create the RF Pulse

In this simulation, we will assume the excitation pulse is very short relative to relaxation and off-resonance effects. In this case, we will use an InstantaneousRF object for the RF pulse that will simulate excitation by just rotating the magnetization vector.

rf = InstantaneousRF()

Let’s add that to bssfp.

# bssfp.jl

using BlochSim

function bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)

    # Step 1: Create a `Spin` object.
    spin = Spin(M0, T1, T2, f)

    # Step 2: Create the RF pulse.
    rf = InstantaneousRF()

    # Step 3: Compute the steady-state signal.

end

Step 3: Compute the Steady-State Signal

Now we are ready to tackle the core of the simulation. As mentioned earlier, there are two ways to go about computing the steady-state signal:

  1. we can manually simulate multiple TRs, or

  2. we can mathematically compute the steady-state magnetization.

We will discuss how to do both.

Version 1: Simulate Multiple TRs

The idea behind simulating multiple TRs is to simulate for long enough so that the magnetization reaches a steady-state. A good rule of thumb is to simulate for a few T1’s worth of time. Thus, the number of TRs we will simulate is 3T1 / TR, rounded up so we get an integer value.

nTR = ceil(Int, 3 * T1 / TR)

(Of course, the longer we simulate the closer to the true steady-state we will get.)

Then for each TR, we have excitation followed by free precession. But we also need to remember to increment the phase of the RF pulse. That means we will need to create a new InstantaneousRF object for each excitation. In this case, we will pass in the flip angle as we did previously, but we will also pass in the phase of the RF pulse (which will be the phase of the previous pulse plus the phase-cycling increment).

for i = 1:nTR
    excite!(spin, rf)
    freeprecess!(spin, TR)
    rf = InstantaneousRF(, rf. + )
end

Now we have the magnetization at the end of the TR, but we want the signal at the echo time. So we need to excite one more time, and then free precess until the echo time.

excite!(spin, rf)
freeprecess!(spin, TE)

Then the steady-state signal can be obtained. One detail to remember, however, is that when doing phase-cycling, the phase of the receiver coil is matched to the phase of the RF pulse. To simulate this, we will multiply our signal by a phase factor.

sig = signal(spin) * cis(-rf.)

(cis() is an efficient way to compute \( e^{i \theta} \) .)

And there we have our steady-state signal!

Let’s update bssfp.

# bssfp.jl

using BlochSim

function bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)

    # Step 1: Create a `Spin` object.
    spin = Spin(M0, T1, T2, f)

    # Step 2: Create the RF pulse.
    rf = InstantaneousRF()

    # Step 3: Compute the steady-state signal.
    if manual_ss
        # Version 1: Simulate multiple TRs.
        nTR = ceil(Int, 3 * T1 / TR)
        for i = 1:nTR
            excite!(spin, rf)
            freeprecess!(spin, TR)
            rf = InstantaneousRF(, rf. + )
        end
        excite!(spin, rf)
        freeprecess!(spin, TE)
        sig = signal(spin) * cis(-rf.)
    else
        # Version 2: Mathematically compute the steady-state.
        # This version is used by default if `manual_ss` is not provided.
    end

    return sig

end

Now let’s see how to compute the steady-state signal a different way.

Version 2: Mathematically Compute the Steady-State

In steady-state, a magnetization vector immediately after excitation (call it \( \mathbf{M}^{+} \) ) is the same as the one immediately after the previous excitation (call it \( \mathbf{M} \) ). Let \( \mathbf{A} \) and \( \mathbf{B} \) be such that \( \mathbf{A} \mathbf{M} + \mathbf{B} \) applies free precession to the magnetization vector, and let \( \mathbf{R} \) be the rotation matrix that represents the instantaneous excitation. Then in steady-state we have

\[
\mathbf{M}^{+}=\mathbf{R}(\mathbf{A}\mathbf{M}+\mathbf{B})=\mathbf{M}.
\]

Solving for \( \mathbf{M} \) gives us

\[
\mathbf{M}=(\mathbf{I}\mathbf{R}\mathbf{A})^{1}\mathbf{R}\mathbf{B},
\]

where \( \mathbf{I} \) is the identity matrix.

This derivation assumes a constant RF phase, so we need to account for phase-cycling somehow. Phase-cycling in bSSFP shifts the off-resonance profile, so it will suffice to use a Spin object with a suitably chosen off-resonance frequency to simulate phase-cycling.

Let’s use this derivation in our code.

spin_pc = Spin(M0, T1, T2, f - ( / 2 / (TR / 1000)))
(R,) = excite(spin_pc, rf)
(A, B) = freeprecess(spin_pc, TR)
M = (I - R * A) \ (R * B)

Here, I comes from the LinearAlgebra standard library module, so we will need to make sure to load it.

Notice we didn’t modify spin, so we need to make sure to copy the computed steady-state magnetization to spin.

copyto!(spin.M, M)

The above gives us the steady-state magnetization immediately after excitation, so now we need to simulate free precession until the echo time.

freeprecess!(spin, TE)

Then we can get the steady-state signal.

sig = signal(spin)

Updating bssfp gives us the following.

# bssfp.jl

using BlochSim, LinearAlgebra

function bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)

    # Step 1: Create a `Spin` object.
    spin = Spin(M0, T1, T2, f)

    # Step 2: Create the RF pulse.
    rf = InstantaneousRF()

    # Step 3: Compute the steady-state signal.
    if manual_ss
        # Version 1: Simulate multiple TRs.
        nTR = ceil(Int, 3 * T1 / TR)
        for i = 1:nTR
            excite!(spin, rf)
            freeprecess!(spin, TR)
            rf = InstantaneousRF(, rf. + )
        end
        excite!(spin, rf)
        freeprecess!(spin, TE)
        sig = signal(spin) * cis(-rf.)
    else
        # Version 2: Mathematically compute the steady-state.
        # This version is used by default if `manual_ss` is not provided.
        spin_pc = Spin(M0, T1, T2, f - ( / 2 / (TR / 1000)))
        (R,) = excite(spin_pc, rf)
        (A, B) = freeprecess(spin_pc, TR)
        M = (I - R * A) \ (R * B)
        copyto!(spin.M, M)
        freeprecess!(spin, TE)
        sig = signal(spin)
    end

    return sig

end

Now bssfp is complete and ready to use!

Using bssfp

To use bssfp, we need to run the bssfp.jl file that loads BlochSim.jl and defines bssfp.

julia> include("bssfp.jl")
bssfp (generic function with 2 methods)

First, let’s make sure the two methods for computing the steady-state magnetization give approximately the same result.

julia> (M0, T1, T2, f, , TR, TE, ) = (1, 1000, 80, 0, /6, 5, 2.5, /2)
(1, 1000, 80, 0, 0.5235987755982988, 5, 2.5, 1.5707963267948966)

julia> sig_manual = bssfp(M0, T1, T2, f, , TR, TE, , true)
0.09889476203144602 + 0.09290409266118088im

julia> sig_exact = bssfp(M0, T1, T2, f, , TR, TE, )
0.09880282817239999 + 0.09281666742806782im

julia> isapprox(sig_exact, sig_manual; rtol = 0.001)
true

The results are less than 0.1% different from each other, indicating the manually simulated steady-state does approach the mathematically computed steady-state, as expected.

Now let’s use bssfp to plot the characteristic bSSFP off-resonance profile. First we compute the steady-state signal for various off-resonance values.

julia> (M0, T1, T2, , TR, TE, , N) = (1, 1000, 80, /6, 5, 2.5, , 801)
(1, 1000, 80, 0.5235987755982988, 5, 2.5, , 801)

julia> f = range(-2 / (TR / 1000), 2 / (TR / 1000), N)
-400.0:1.0:400.0

julia> sig = bssfp.(M0, T1, T2, f, , TR, TE, );

Note the dot (.) when calling bssfp. This essentially is shorthand for

sig = zeros(ComplexF64, length(f))
for i = 1:length(f)
    sig[i] = bssfp(M0, T1, T2, f[i], , TR, TE, )
end

(See our blog post on broadcasting for more information.)

And now we can plot the results.

julia> using Plots

julia> pmag = plot(f, abs.(sig); label = "", title = "bSSFP Off-Resonance Profile", ylabel = "Magnitude", color = 1);

julia> pphase = plot(f, angle.(sig); label = "", xlabel = "Off-Resonance (Hz)", ylabel = "Phase (rad)", yticks = ([-, 0, ], ["-", "0", ""]), color = 2);

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

Summary

In this post, we have learned how to use BlochSim.jl to simulate the steady-state signal obtained from a bSSFP scan. We used the bssfp function we wrote to plot the characteristic bSSFP off-resonance profile to demonstrate that our code produces the expected results.

Additional Links

MRI Scan Simulation Made Easy with BlochSim.jl

By: Steven Whitaker

Re-posted from: https://glcs.hashnode.dev/mri-scan-simulation-made-easy-with-blochsimjl

MRI scans manipulate magnetization vectors to produce an observable signal.

Often it can be useful to simulate what signal an MRI scan will generate. Example scenarios include simulating MRI signals for educational purposes, for developing a new MRI scan, or for fitting model parameters to an acquired signal.

In this post, we will learn how to use the Julia programming language, a free and open source programming language that excels especially in scientific computing, to simulate the signal generated by an MRI scan. Specifically, we will use the BlochSim.jl package to simulate a balanced steady-state free precession (bSSFP) scan.

Note that this post assumes a basic understanding of MRI pulse sequences and Bloch simulations. See our previous posts for an overview of the Bloch equations and an overview of BlochSim.jl.

Prerequisites

To run the code in this post, both Julia and BlochSim.jl need to be installed. See our previous post for installation instructions.

bSSFP

First, we need to review bSSFP so we know how to simulate it.

Here is a pulse sequence diagram of bSSFP (source).

For this post, we are not simulating signal localization, so we will ignore the imaging gradients.

Thus, for each TR we will simulate

  1. excitation,

  2. free precession until the echo time, and

  3. free precession until the end of the TR.

This, by itself, would be essentially the same simulation as done in our previous post. However, two aspects of bSSFP make simulating it a bit more interesting.

  1. bSSFP, as the name states, is a steady-state scan. This means we don’t simulate just one excitation followed by free precession to the echo time. Instead, we need to simulate multiple TRs until we get to steady-state. We can do this either by manually simulating multiple TRs or by mathematically computing the steady-state magnetization. More on this later.

  2. bSSFP utilizes phase-cycling, where the phase of the RF pulse is incremented after each TR.

Therefore, our goal is to write a Julia function (that we will call bssfp) that computes the steady-state signal given by a bSSFP scan that uses phase-cycling.

Throughout this post, we will be editing a text file called bssfp.jl that will load BlochSim.jl and define our bssfp function.

Function Outline

We will begin with an outline of what our bssfp function needs to do.

# bssfp.jlusing BlochSim# Step 0: Specify the inputs to the function.function bssfp(inputs)    # Step 1: Create a `Spin` object.    # Step 2: Create the RF pulse.    # Step 3: Compute the steady-state signal.end

First, we need to specify what inputs bssfp needs. We need values for M0, T1, T2, and off-resonance frequency \( \Delta f \) to create a Spin object, a flip angle \( \alpha \) to create the RF pulse, and the TR, TE, and phase-cycling factor \( \Delta\phi \) to use. We will also add an input that specifies how to compute the steady-state magnetization.

# bssfp.jlusing BlochSimfunction bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)    # Step 1: Create a `Spin` object.    # Step 2: Create the RF pulse.    # Step 3: Compute the steady-state signal.end

Here we see that Julia allows using Greek letters in variable names, which helps the code look more like the math, enhancing readability. (To type, e.g., , type \alpha and then press TAB, and similarly for other Greek letters.)

We also made manual_ss an optional argument. If the user does not specify manual_ss, it will default to false.

Now we can focus on writing code for each step outlined above.

Step 1: Create a Spin Object

As discussed in our previous post, a Spin object is created as follows.

spin = Spin(M0, T1, T2, f)

That’s it!

So let’s update bssfp.

# bssfp.jlusing BlochSimfunction bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)    # Step 1: Create a `Spin` object.    spin = Spin(M0, T1, T2, f)    # Step 2: Create the RF pulse.    # Step 3: Compute the steady-state signal.end

Step 2: Create the RF Pulse

In this simulation, we will assume the excitation pulse is very short relative to relaxation and off-resonance effects. In this case, we will use an InstantaneousRF object for the RF pulse that will simulate excitation by just rotating the magnetization vector.

rf = InstantaneousRF()

Let’s add that to bssfp.

# bssfp.jlusing BlochSimfunction bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)    # Step 1: Create a `Spin` object.    spin = Spin(M0, T1, T2, f)    # Step 2: Create the RF pulse.    rf = InstantaneousRF()    # Step 3: Compute the steady-state signal.end

Step 3: Compute the Steady-State Signal

Now we are ready to tackle the core of the simulation. As mentioned earlier, there are two ways to go about computing the steady-state signal:

  1. we can manually simulate multiple TRs, or

  2. we can mathematically compute the steady-state magnetization.

We will discuss how to do both.

Version 1: Simulate Multiple TRs

The idea behind simulating multiple TRs is to simulate for long enough so that the magnetization reaches a steady-state. A good rule of thumb is to simulate for a few T1’s worth of time. Thus, the number of TRs we will simulate is 3T1 / TR, rounded up so we get an integer value.

nTR = ceil(Int, 3 * T1 / TR)

(Of course, the longer we simulate the closer to the true steady-state we will get.)

Then for each TR, we have excitation followed by free precession. But we also need to remember to increment the phase of the RF pulse. That means we will need to create a new InstantaneousRF object for each excitation. In this case, we will pass in the flip angle as we did previously, but we will also pass in the phase of the RF pulse (which will be the phase of the previous pulse plus the phase-cycling increment).

for i = 1:nTR    excite!(spin, rf)    freeprecess!(spin, TR)    rf = InstantaneousRF(, rf. + )end

Now we have the magnetization at the end of the TR, but we want the signal at the echo time. So we need to excite one more time, and then free precess until the echo time.

excite!(spin, rf)freeprecess!(spin, TE)

Then the steady-state signal can be obtained. One detail to remember, however, is that when doing phase-cycling, the phase of the receiver coil is matched to the phase of the RF pulse. To simulate this, we will multiply our signal by a phase factor.

sig = signal(spin) * cis(-rf.)

(cis() is an efficient way to compute \( e^{i \theta} \) .)

And there we have our steady-state signal!

Let’s update bssfp.

# bssfp.jlusing BlochSimfunction bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)    # Step 1: Create a `Spin` object.    spin = Spin(M0, T1, T2, f)    # Step 2: Create the RF pulse.    rf = InstantaneousRF()    # Step 3: Compute the steady-state signal.    if manual_ss        # Version 1: Simulate multiple TRs.        nTR = ceil(Int, 3 * T1 / TR)        for i = 1:nTR            excite!(spin, rf)            freeprecess!(spin, TR)            rf = InstantaneousRF(, rf. + )        end        excite!(spin, rf)        freeprecess!(spin, TE)        sig = signal(spin) * cis(-rf.)    else        # Version 2: Mathematically compute the steady-state.        # This version is used by default if `manual_ss` is not provided.    end    return sigend

Now let’s see how to compute the steady-state signal a different way.

Version 2: Mathematically Compute the Steady-State

In steady-state, a magnetization vector immediately after excitation (call it \( \mathbf{M}^{+} \) ) is the same as the one immediately after the previous excitation (call it \( \mathbf{M} \) ). Let \( \mathbf{A} \) and \( \mathbf{B} \) be such that \( \mathbf{A} \mathbf{M} + \mathbf{B} \) applies free precession to the magnetization vector, and let \( \mathbf{R} \) be the rotation matrix that represents the instantaneous excitation. Then in steady-state we have

\[\mathbf{M}^{+}=\mathbf{R}(\mathbf{A}\mathbf{M}+\mathbf{B})=\mathbf{M}.\]

Solving for \( \mathbf{M} \) gives us

\[\mathbf{M}=(\mathbf{I}\mathbf{R}\mathbf{A})^{1}\mathbf{R}\mathbf{B},\]

where \( \mathbf{I} \) is the identity matrix.

This derivation assumes a constant RF phase, so we need to account for phase-cycling somehow. Phase-cycling in bSSFP shifts the off-resonance profile, so it will suffice to use a Spin object with a suitably chosen off-resonance frequency to simulate phase-cycling.

Let’s use this derivation in our code.

spin_pc = Spin(M0, T1, T2, f - ( / 2 / (TR / 1000)))(R,) = excite(spin_pc, rf)(A, B) = freeprecess(spin_pc, TR)M = (I - R * A) \ (R * B)

Here, I comes from the LinearAlgebra standard library module, so we will need to make sure to load it.

Notice we didn’t modify spin, so we need to make sure to copy the computed steady-state magnetization to spin.

copyto!(spin.M, M)

The above gives us the steady-state magnetization immediately after excitation, so now we need to simulate free precession until the echo time.

freeprecess!(spin, TE)

Then we can get the steady-state signal.

sig = signal(spin)

Updating bssfp gives us the following.

# bssfp.jlusing BlochSim, LinearAlgebrafunction bssfp(M0, T1, T2, f, , TR, TE, , manual_ss = false)    # Step 1: Create a `Spin` object.    spin = Spin(M0, T1, T2, f)    # Step 2: Create the RF pulse.    rf = InstantaneousRF()    # Step 3: Compute the steady-state signal.    if manual_ss        # Version 1: Simulate multiple TRs.        nTR = ceil(Int, 3 * T1 / TR)        for i = 1:nTR            excite!(spin, rf)            freeprecess!(spin, TR)            rf = InstantaneousRF(, rf. + )        end        excite!(spin, rf)        freeprecess!(spin, TE)        sig = signal(spin) * cis(-rf.)    else        # Version 2: Mathematically compute the steady-state.        # This version is used by default if `manual_ss` is not provided.        spin_pc = Spin(M0, T1, T2, f - ( / 2 / (TR / 1000)))        (R,) = excite(spin_pc, rf)        (A, B) = freeprecess(spin_pc, TR)        M = (I - R * A) \ (R * B)        copyto!(spin.M, M)        freeprecess!(spin, TE)        sig = signal(spin)    end    return sigend

Now bssfp is complete and ready to use!

Using bssfp

To use bssfp, we need to run the bssfp.jl file that loads BlochSim.jl and defines bssfp.

julia> include("bssfp.jl")bssfp (generic function with 2 methods)

First, let’s make sure the two methods for computing the steady-state magnetization give approximately the same result.

julia> (M0, T1, T2, f, , TR, TE, ) = (1, 1000, 80, 0, /6, 5, 2.5, /2)(1, 1000, 80, 0, 0.5235987755982988, 5, 2.5, 1.5707963267948966)julia> sig_manual = bssfp(M0, T1, T2, f, , TR, TE, , true)0.09889476203144602 + 0.09290409266118088imjulia> sig_exact = bssfp(M0, T1, T2, f, , TR, TE, )0.09880282817239999 + 0.09281666742806782imjulia> isapprox(sig_exact, sig_manual; rtol = 0.001)true

The results are less than 0.1% different from each other, indicating the manually simulated steady-state does approach the mathematically computed steady-state, as expected.

Now let’s use bssfp to plot the characteristic bSSFP off-resonance profile. First we compute the steady-state signal for various off-resonance values.

julia> (M0, T1, T2, , TR, TE, , N) = (1, 1000, 80, /6, 5, 2.5, , 801)(1, 1000, 80, 0.5235987755982988, 5, 2.5, , 801)julia> f = range(-2 / (TR / 1000), 2 / (TR / 1000), N)-400.0:1.0:400.0julia> sig = bssfp.(M0, T1, T2, f, , TR, TE, );

Note the dot (.) when calling bssfp. This essentially is shorthand for

sig = zeros(ComplexF64, length(f))for i = 1:length(f)    sig[i] = bssfp(M0, T1, T2, f[i], , TR, TE, )end

(See our blog post on broadcasting for more information.)

And now we can plot the results.

julia> using Plotsjulia> pmag = plot(f, abs.(sig); label = "", title = "bSSFP Off-Resonance Profile", ylabel = "Magnitude", color = 1);julia> pphase = plot(f, angle.(sig); label = "", xlabel = "Off-Resonance (Hz)", ylabel = "Phase (rad)", yticks = ([-, 0, ], ["-", "0", ""]), color = 2);julia> plot(pmag, pphase; layout = (2, 1))

Summary

In this post, we have learned how to use BlochSim.jl to simulate the steady-state signal obtained from a bSSFP scan. We used the bssfp function we wrote to plot the characteristic bSSFP off-resonance profile to demonstrate that our code produces the expected results.

Additional Links

Mastering MRI Bloch Simulations with BlochSim.jl in Julia

By: Steven Whitaker

Re-posted from: https://blog.glcs.io/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_fp
10

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)
1000

julia> 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)
100

julia> 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 = /2
1.5707963267948966

julia> 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 Plots

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

julia> 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.M
Magnetization 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