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:

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:

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.

• 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. – M529 Jul 12 '16 at 21:10
• @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... – Luca Jul 12 '16 at 21:39

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.

• 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. – Luca Jul 11 '16 at 18:50
• 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. – M529 Jul 11 '16 at 19:00
• 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 – Luca Jul 12 '16 at 14:52