Bloch equation: adding gradients

Question

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.

Answer 1

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.