Signal Generation
The basic principle of signal generation in MPI relies on the non-linear magnetization response M(H) of ferromagnetic particles to an applied magnetic field H, cf. figure 1. An oscillating drive field H
D(t) of sufficient amplitude leads to a magnetization response M(t) of the particles, which has a different spectrum of higher harmonics than the drive field. If, for instance, a harmonic drive field is used, the drive field spectrum only contains the base frequency, whereas the particle response also contains multiples thereof. The information contained in these higher harmonics is used for MPI. Experimentally, the time-dependent change in particle magnetization is measured via the induced voltage in receive coils. Assuming a single receive coil with sensitivity σ
r
(r) at spatial position r, the changing magnetization induces a voltage
according to Faraday's law. μ
0 is the magnetic permeability of vacuum. The receive coil sensitivity σ
r
(r) = H
r
(r)/I
0 derives from the field H
r
(r) the coil would produce if driven with a unit current I
0 [5].
In the following, the sensitivity of the receive coil is approximated to be homogeneous over the region of interest, i.e., σ
r
(r) is constant. If M
x
(r, t) is the magnetization component picked up by a receive coil in x direction, i.e., having the sensitivity σ
r
= (σ
x
, 0, 0)T, the detected signal can be written as
Neglecting constant factors, we introduce the notation s(r, t) for signal generated by a point-like distribution of particles at position r. If the particle distribution is approximated by a δ distribution, the volume integral vanishes and the particle magnetization M
x
(r, t) is determined by the local field H(r, t). For the moment, the field is assumed to have only one spatial component H
x
(r, t), which is pointing in receive-coil direction x. The signal can then be written as
Since this equation holds for a general 1D setup where the field is aligned with the direction of the acquired magnetization component, the subscript x has been omitted. Equation 3 shows that high signal results from the combination of a steep magnetization curve with rapid field variations. Fourier expansion of the periodic signal s(t) generated by applying a homogeneous drive field H(r, t) = H
D(t) yields the signal spectrum S
n
, as shown in figure 1. Intensity and weight of higher harmonics in the spectrum are related to the shape of the magnetization curve M(H), and to the waveform and amplitude of the drive field H
D(t). To illustrate their influence on the spectrum, a number of representative cases are shown in figure 2.
The step function relates to an immediate particle response and creates a spectrum that is rich in high harmonics. The spectral components have constant magnitude at odd multiples of the drive frequency. The even harmonics are missing due to the sine-type pattern of the time signal s(t). The step function corresponds to an ideal particle response and represents the limiting case for the achievable weight of higher harmonics. For this magnetization curve, triangle and sine drive fields yield the same result.
If the particle response to the drive field is slowed down by introducing a linear range in the magnetization curve, the relative weight of higher harmonics is reduced. Thereby, the harmonic drive field performs better than the triangular excitation, since it sweeps faster over the linear range.
The last row in figure 2 shows a more realistic particle magnetization as given by the Langevin function [6]
M(ξ) = M
0 (coth ξ - 1/ξ), ((4)
where M
0 is the saturation magnetization and ξ is the ratio between magnetic energy of a particle with magnetic moment m in an external field H, and thermal energy given by the Boltzmann constant k
B and temperature T:
A higher magnetic moment results in a steeper magnetization curve and creates more higher harmonics for a given drive field amplitude. Alternatively, high harmonics can be generated from a shallow curve using faster field variations, e.g., induced by a higher drive field amplitude. It should be noted that MPI uses ferromagnetic particles to obtain a sufficiently steep magnetization curve [1]. For low concentrations, however, their mutual interactions can be neglected and they can be treated like a gas of paramagnetic particles with extremely large magnetic moment, a phenomenon also known as "super-paramagnetism" [7].
1D Spatial Encoding
To encode spatial information in the signal, a static magnetic gradient field H
S(r), also called selection field, is introduced. For 1D encoding, the selection field has a non-zero gradient only in x direction, G
x
= dH
S/dx. If the gradient is non-zero over the complete field of view (FOV), the selection field is monotonically rising or falling, thus it can cross zero only in a single point, the above-mentioned FFP. In regions far away from the FFP, the particle magnetization is driven into saturation by the selection field.
Application of a spatially homogeneous and temporally periodic drive field H
D(t) in addition to the selection field H
S(r) corresponds to a periodic displacement of the FFP along the gradient direction. The particles experience a local field
H(x, t) = H
S(x) - H
D(t). (6)
The minus sign has been chosen to make later calculations more convenient. Since the FFP sweeps over each spatial position x at a different point in time, each position can be identified by its characteristic spectral response.
Harmonic Drive Field – Ideal Particles
Figure 3 shows spectra at three different spatial positions generated by ideal particles exposed to a selection field H
S of constant gradient strength G and a harmonic drive field H
D of frequency ω
0 = 2π/T and amplitude A. A derivation of their functional form is given in appendix A.1 and A.2. For the nth harmonic, corresponding to the nth component of a Fourier series expansion, one finds the following dependence on particle position x:
where the U
n
(x) represent Chebyshev polynomials of the second kind (42). The functions are defined in the range -1 < Gx/A < 1. A cosine drive field has been used instead of the sine drive field to arrive at a simpler expression. The spatial dependence for the first harmonics is plotted in the left part of figure 4. One finds an increasing number of oscillations with increasing frequency components n. This relates to the fact that Chebyshev polynomials form a complete orthogonal basis set (cf. appendix A.2), so that any particle distribution C(x) can be expanded into these functions. Successive frequency components have alternating spatial parity with respect to the center of the FOV (even/odd).
The S
n
(x) can be seen as a sensitivity map, describing the spatial sensitivity profile of each frequency component n. In an MPI experiment, a 1D particle distribution C(x) would generate spectral signal components
Thus, the S
n
(x) represent the system function introduced by [3]. The system function not only describes the spatial signal dependence, but also contains information about the particles' magnetization curve and system parameters (e.g. drive field amplitude A and frequency ω
0 = 2π/T, selection field gradient G).
Using (7), the spectral signal components (8) for ideal particles can be written as
In this notation, the V
n
correspond to coefficients of a Chebyshev series, cf. equation (58). It follows that the particle concentration can be reconstructed by doing a Chebyshev transform of the measured V
n
, i.e., by evaluating the Chebyshev series
Hence, for ideal particles under the influence of a harmonic drive field and a constant selection field gradient, reconstruction of the spatial particle distribution simply corresponds to calculating the sum over the measured harmonics V
n
weighted with Chebyshev polynomials of the second kind. In terms of the system function
, the particle concentration can be written as
Harmonic Drive Field – Langevin Particles
For more realistic particles as described by (4), the system function S
n
(x) is given by a spatial convolution between the derivative of the magnetization curve, M'(H
S), and the Chebyshev components, as derived in A.2. Using equation (7) one can write:
Depending on the steepness of M(H), the S
n
(x) will be a blurred version of the
, extending slightly beyond the interval which is covered by the FFP motion and to which the
are confined. Thus, particles that are located slightly outside the range accessed by the FFP also generate signal. The right part of figure 4 displays components of the system function for particles following the Langevin magnetization curve in a constant selection field gradient.
In the measurement process according to (8), the FOV now refers to the range where the S
n
(x) are non-zero. A sufficiently steep magnetization curve can provide confinement to a region not much larger than the range covered by the FFP, i.e., -A/G < x <A/G.
Since the system function components cannot be sharper than the convolution kernel, an MPI experiment with Langevin particles will run into a resolution limit correlating with the width of M'(x). Since the derivative of the magnetization curve is a symmetric function M'(x) = M'(-x), one can use (12) to show that
where
corresponds to the expression in square brackets in (13). Since (14) corresponds to (9), reconstruction for the harmonic drive field is given by (11), i.e.,
This means, that in the interval where the ideal particle system function
is defined, i.e., -A/G < x <A/G,
can be directly reconstructed. If the particle concentration C(x) is confined to the FOV,
is just the convolution of C(x) with M'(x):
From this equation, one can infer that the resolution of the reconstructed image is limited by the width of M'(x). If the particle magnetization curve is known and the measurement provides sufficient SNR, deconvolution can be used to overcome this limitation. However, in practical applications with distributions of different particles sizes and magnetization curves, deconvolution may be difficult. Therefore, the following paragraph gives an estimation of the achievable spatial resolution without deconvolution. The derivative of the Langevin curve in equation (4) is
which corresponds to the blue curve plotted in figure 1. The full width at half maximum (FWHM) of this curve can be determined numerically as Δξ
FWHM ≈ 4.16. If the particle magnetization m and the selection field gradient strength G are known, this can be translated to a spatial resolution using equation (5):
The particle magnetization depends on particle diameter d according to the following relation [8]:
Assuming magnetite particles (Fe3O4) with a saturation magnetization of M
S = 0.6 T/μ
0, the resolution limit Δx imposed by the magnetization curve can be calculated. Figure 5 shows Δx as a function of gradient strength for different particle diameters. The cross-hairs indicate the gradient strength of 5.5 T/m/μ
0 and the dominating particle diameter of 30 nm used in real-time in vivo MPI [4]. This gradient/particle combination theoretically allows a resolution better than 0.5 mm. However, due to the wide distribution of particle sizes and the regularization applied in reconstruction [3] to mitigate limited SNR, the observed resolution was not better than Δx ≈ 1.5 mm.
Triangular Drive Field
An illustrative case is to use a triangular drive field instead of the harmonic field, cf. appendix A.3. The system function for ideal particles then has the form
for an FFP motion covering the range 0 < x < 2A/G. Now, instead of the Chebyshev series, a Fourier series can be used to reconstruct a particle concentration
The measured frequency components V
n
are proportional to components
in k space, which are related to the spatial distribution C(x) by Fourier transformation. In terms of the system function, (21) becomes
For realistic particles, the system function has to be convolved with M'(H
S). Since equations (12–14) derived for harmonic drive field excitation hold for the triangular system function as well, a modified/convolved concentration
can be reconstructed in the range 0 < x < 2A/G.
Matrix Formulation
For MPI image reconstruction, the continuous spatial distribution of particles will be mapped to a grid, where each grid location represents a small spatial region. Furthermore, only a limited number n of frequency components is recorded. If the spatial positions are indexed with m, (8) becomes
or, in vector/matrix formulation,
v= Sc. (24)
Calculation of the concentration vector then basically corresponds to an inversion of matrix S:
c= S
-1
v. (25)
This notation will be used for 2D or 3D imaging as well, which requires collapsing spatial indices into the single index m. Thus, concentration is always a vector, independent of spatial dimension.
Going back to the 1D case for a harmonic drive field, introduction of a scalar
and a diagonal matrix
allows derivation of the following identity by comparing (25) with (11):
S
-1 = α βS
T. (28)
Thus, in the case of 1D imaging of ideal particles, the inverse matrix can simply be obtained by multiplication of the transpose with a scalar and a diagonal matrix.
Using only a limited number of frequency components corresponds to working with a truncated Chebyshev series. The Chebyshev truncation theorem then states that the error in approximating the real concentration distribution is bounded by the sum of the absolute values of the neglected coefficients. More importantly, for reasonably smooth distributions, the error is on the order of the last retained Chebyshev coefficient [9].
2D and 3D Spatial Encoding
1D Drive Field
A first step towards describing 2D and 3D imaging is to look at the 3D system function of particles in a 3D selection field H
S(r) combined with a 1D drive field H
D(t). Using a harmonic drive field and choosing a Maxwell coil setup to create a selection field as described in equation (68), the total field can be approximated by
The system function can be written as a convolution over the z component of the selection field (cf. (65))
In this vector, each component refers to the signal induced by the respective x/y/z magnetization component. Detection of these components requires three orthogonal (sets of) receive coils. For ideal particles (cf. (66)), the explicit spatial dependence becomes (cf. (70))
where the asterisk denotes convolution over the z component, i.e., the direction of the FFP motion resulting from the drive field. Thus, an expression describing the 3D spatial dependence of the respective magnetization component is convolved in drive field direction with the set of 1D Chebyshev functions.
The shape of the convolution kernel is determined by ∂ M/∂ H
z
, which describes how the magnetization responds to the drive field change. For ideal particles, it is singular at the origin. Figure 6 shows the xz plane of the 3D kernel for the signal components detected in z and x direction, S
n,z
(r) and S
n,x
(r), respectively. Along the center line in drive field direction, the kernel for the M
z
magnetization corresponds to the delta distribution, just as in the 1D situation (48). With increasing distance from the center line, the kernel broadens and its amplitude decreases rapidly. For M
x
, and for symmetry reasons also for M
y
, the kernel is zero on the symmetry axes. It has high amplitude close to the singularity at the origin.
To form the 3D ideal particle system function, the 3D kernel is convolved along the drive field direction with the 1D Chebyshev polynomials (31). Figure 7 shows central 2D slices extracted from selected harmonics for the above case of 1D drive field motion in z direction. Directly on the line covered by the FFP trajectory, the system function is given by Chebyshev polynomials and therefore can encode an arbitrary particle distribution, as discussed for the 1D situation. With increasing distance to the center line, the convolution kernel has an increasing blurring effect, so that finer structures of the higher Chebyshev polynomials are averaged to zero. Therefore, the signal in higher system function components condenses to the line of the FFP trajectory (cf. figure 7, harmonic 12 and 25), where the blurring effect is low. This can be explained by the fact that only in the close vicinity of the FFP, the field change is rapid enough to stimulate a particle response that generates high frequency components.
In an MPI experiment, it can be useful to exploit symmetries in the system function to partially synthesize the system function and thus speed up its acquisition and reduce memory requirements. From the 3D response to the 1D FFP motion, two basic rules can be derived for the parity of the system function in the spatial direction indexed with i ∈ {x, y, z}.
1. The "base" parity is given by the parity of the convolution kernel shown in figure 6. It is even, if the receive direction j ∈ {x, y, z} is aligned with the drive direction k ∈ {x, y, z}. This corresponds to the magnetization derivative component ∂M
j
∂H
k
for j = k. Otherwise kernel parity is odd:
2. If the spatial direction of interest is a drive field direction, i.e., i = k, then parity alternates between successive harmonics h of that drive field component:
The reason is the alternating parity of the Chebyshev polynomials in the 1D system function.
Parity observed for harmonic h in a spatial direction i then is p
i,j,k,h
= p
Cheb·p
kernel.
2D and 3D Drive Field
As displayed in figure 7, the spatial pattern of the particle response for higher harmonics is confined to a narrow region close to the line of the FFP motion. Therefore, to induce a particle response with high harmonics extending over the whole 2D plane, a lateral displacement of the FFP is necessary. This is achieved by adding a second drive field which displaces the FFP in x direction. For 3D encoding, a third drive field for FFP displacement in y direction is necessary. Figure 7 furthermore shows that high resolution is obtained only in the close vicinity of the FFP trajectory line. Thus, the 2D or 3D trajectory should be sufficiently dense to achieve homogeneous resolution over the imaging plane or volume. For a simple implementation, one can choose harmonic drive fields with slight frequency differences in the orthogonal spatial directions, causing the FFP motion to follow a 2D or 3D Lissajous pattern. In the following, a 2D system function requiring two drive frequencies is investigated. The treatment of a 3D system function using three drive frequencies would be analogous. Figure 8 displays a 2D Lissajous pattern generated by the superposition of two orthogonal harmonic drive fields with frequency ratio ω
x
/ω
z
= 24/25:
Using a 3D selection field according to (29), the doubled selection field gradient in z direction requires A
z
= 2A
x
to cover a quadratic FOV with the FFP motion. Figure 8 displays the first components of a simulated 2D system function for ideal particles exposed to the superposition of the 2D Lissajous drive field and the 3D selection field. Each receive direction has its own set of system functions, denoted by 'receive x' and 'receive z'. Components corresponding to higher harmonics of the respective drive frequency are indicated by red frames. On the x channel, they have a spacing of 24 components. In the spatial x direction, they closely resemble the 1D Chebyshev series, while in z direction, they show no spatial variation. On the z channel, harmonics of the drive frequency exhibit a spacing of 25 components with a spatial pattern that is basically rotated by 90 degrees with the respect to the x components.
While components corresponding to harmonics of the drive field frequencies only allow 1D encoding in the respective drive field direction, components arising from a mixture of both drive frequencies provide spatial variation in both directions at the same time. For instance, moving to the left from the first x drive-field harmonic on the x channel (component 24) corresponds to mix frequencies mω
x
+ n(ω
x
- ω
z
) with increasing integer n and m = 1. For larger m, one starts at a higher harmonic m. Moving to the right corresponds to negative n. Thus, pure drive field harmonics and their vicinity relate to low mixing orders, while increasing distance goes along with larger n and higher mixing orders.
It should be noted that the system function component observed for mω
x
+ n(ω
x
- ω
z
) appears a second time at frequency mω
x
+ n(ω
x
+ ω
z
). Thus every component corresponding to frequency mixes appears twice. Examples are components 23 and 73 (m = 1, n = 1, green frames) or 47 and 97 (m = 2, n = 1, blue frames), but also 26 and 74 (m = 1, n = -2, orange frames) on the x channel.
Figure 8 also plots the maximum intensities (weights) of the generated system function components. Highest intensities are found in pure multiples of the drive frequencies, however with a decrease towards higher frequencies. Components corresponding to mix frequencies have much lower intensity than pure harmonics. If the system function is acquired experimentally, these components will be the first to fall below the noise level. Thus, low SNR in the system function acquisition will reduce the achievable resolution.
The higher the order of a system function component, the finer its spatial structure. This behavior and the general spatial patterns closely resemble 2D Chebyshev polynomials, which can be written as a tensor product of the 1D polynomials for each direction: U
n
(x) ⊗ U
m
(z). Figure 9 plots the first components of these functions. The 2D Chebyshev functions satisfy an orthogonality relation similar to (50). A graphical representation of this relation for the first 256 components is shown in the left part of figure 10. The inner product between orthogonal functions vanishes, so that only the product of a function with itself is non-zero, leading to the diagonal line in figure 10. In the right part, the corresponding plot is shown for the normalized rows of an ideal particle 2D Lissajous system function. Bright spots and lines off the diagonal indicate that some system function components are not orthogonal with respect to each other. However, black regions prevail and one can infer that most components are orthogonal. Therefore, there is only little redundancy in the system function.
To demonstrate this, a phantom image (figure 11, left) is expanded into an equal number of 2D Chebyshev and Lissajous system function components, respectively. The number of components has been chosen to equal the the number of pixels in the image (64 × 64). The image obtained from the Chebyshev transformation exhibits reduced resolution compared to the original image. The reason is that the Chebyshev functions provide higher resolution at the edges of the FOV but reduced resolution at the center. To keep the high resolution at the image center, higher Chebyshev components would have to be included in the expansion. The image obtained from the system function components has been reconstructed by inverting the system function matrix using minimal regularization to suppress noise [10]. Half the system function components were taken from the receive x system function, the other half from the z function (as displayed in figure 8). Resolution of the image is better than observed for the Chebyshev expansion, but the image has small artifacts that make it appear less homogeneous. Considering the fact that some system function components are not orthogonal and therefore redundant, the image quality is quite good. The reconstructions from only the x or z components show significantly worse image quality, indicating that these subsets are not sufficient to homogeneously represent the image information.