2
$\begingroup$

I have a bloch simulation that I have implemented in Python. The implementation is very simple as:

import numpy as np

GYRO = 267538030.3797  # in radians / s / T
T1 = 0.3
inv_T1 = 1.0/0.3
T2 = 0.05
inv_T2 = 1.0/0.05
M0 = 1.0 

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B) - [M[0] * inv_T2,
                                    M[1] * inv_T2,
                                    inv_T1 * (M[2] - M0)]

Now, during the excitation, the effective field is given by: So, the pseudocode is:

\begin{align} B[0] &= \textrm{rf}(t) \cos(2 \pi \omega_{\rm eff} \ t) \\ B[1] &= \textrm{rf}(t) \sin(2 \pi \omega_{\rm eff} \ t) \\ B[2] &= 0 \end{align}

Assuming, that the RF excitation is on-resomance i.e. $\omega_{\rm eff} = 0$, I track the evolution of the magnetization of the spin across time and it works and after relaxing for T1 seconds, we have this curve:

enter image description here

So, the transverse magnetization starts to disappear and longitudinal magnetization recovers.

Now, I wanted to simulate a readout gradient which is applied just after the RF pulse. Now, the $X$ and $Y$ values for the effective magnetization is 0 as the RF pulse is turned off. The $Z$ component of the effective field at every time step is given by the product of gradient amplitude and the position. So,

\begin{align} B[0] &= B[1] = 0 \\ B[2] &= g * p \end{align}

Now, when I use the same equation to solve this, I get nonsensical results, which is shown by the curve below:

enter image description here

As you can see, the numbers are off the charts and I was expecting an echo at the end. The gradient is a simple trapezoidal gradient and I have checked the units and it all makes sense. I was wondering if my expression for the effective field in the presence of the gradient is correct? I have seen these being represented as rotation matrices as well and I wonder if there is something like that I might have missed.

Update

As per suggestion from M529, I did the following experiments:

# Change the method to ignore relaxation effects
import numpy as np

GYRO = 267538030.3797  # in radians / s / T

# This function is called at various time steps
def bloch(self, M, B):
    # Solve bloch equation in rotating frame.
    # M: Magnetization vector
    # B: Effective magnetic field
    return GYRO * np.cross(M, B) 

Now I call the method with the different flip angles (from 0-1200) and note the vector length and it does not maintain the length at larger flip angles. So, perhaps there is a numeric issue. I increased the RF bandwidth for computation reasons. Here is the graph of the magnetization at different flips. It is a weird non-linear effect. It first goes down and then climbs back up uncontrollably.

enter image description here

$\endgroup$
  • 1
    $\begingroup$ You can also use a rectangular RF pulse, if this eases the computational burden. Maybe this is even a simpler way to analyze the problem. $\endgroup$ – M529 Jul 12 '16 at 21:10
  • $\begingroup$ @M529 I am currently trying to see what happens with this "hard pulse" approximation so the little rotation matrices (which should commute as well). I am not sure how bad of an approximation this is especially when the FA is large like SE sequences. But this is whole thing is a learning exercise anyway.... I will try the rectangular pulse after that. Many thanks for your very useful insights. In the end, I guess I can use the spike RF approximation and just do an instantaneous excitation... $\endgroup$ – Luca Jul 12 '16 at 21:39
2
$\begingroup$

I am not sure if this is related, but this cross-product is part of the differential equation, something like $\frac{\textrm{d}M}{\textrm{d}t} = \gamma M \times B$... For a true infinitesimal $\textrm{d}t$, this makes sense. For nummerical considerations, this needn't be true. Think of the following:

  • neglect relaxation for a moment
  • $\vec{M}_0$ has a certain magnitude (i.e. length)
  • when you apply an rf pulse, you expect a rotation of $\vec{M}_0$, and a vector's length is preserved under rotations
  • when you calculate the rotation from the cross product, with several time-steps $\Delta t$, you probably do something like $\vec{M}_{0,n+1} = \vec{M}_{0,n} + \Delta t\;\gamma\vec{M}_{0,n} \times \vec{B}$ (i.e. you compute the change from the previous magnetization and the rf pulse and add it to the previous magnetization vector)
  • this operation does not preserve the length of the vector: $\vec{M}_{0,n} \times \vec{B}$ is perpendicular to $\vec{M}_{0,n}$. When you add the vectors, the resulting vector is longer than $\vec{M}_{0,n}$ (Pythagorean theorem)

Hence, if you do not take appropriate countermeasures, your magnetization vector will grow. Maybe this is what you see in your simulation: At the end you have an insanely huge signal that does not belong there. You can easily check this via a simple test: Neglect relaxation and apply a huge flip angle, lets say $90000^\circ$ and monitor your magnetization vector length. Certainly it should remain constant. If it doesn't, there might be a bug in your code.

Also try to reformulate your calculations with rotations matrices and check if the results get better.

$\endgroup$
  • 1
    $\begingroup$ Thanks for the response. I noticed that making the time steps smaller indeed makes things more reasonable but there are still errors. If I may ask, I read that with the rotations perhaps the RF rotations and the relaxation and free precession matrices do not commute. I am still looking into this but is this something to keep in mind? Thanks for suggesting these tests. I will update the thread as soon as I have some results. $\endgroup$ – Luca Jul 11 '16 at 18:50
  • 1
    $\begingroup$ You should certainly keep this non-commutativity in mind during the design of your simulations. It will follow (and guide?) you through the decisions whether you first do a relaxation and then a rotation due to an rf pulse, or the other way around. $\endgroup$ – M529 Jul 11 '16 at 19:00
  • 1
    $\begingroup$ It seems that there is a numerical issue as you suggested. I have updated the thread with magnetization evolution graph now. I am looking at the behavior for smaller time steps now and will try and formulate the whole thing with rotation matrices and update $\endgroup$ – Luca Jul 12 '16 at 14:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.