Original Research ARTICLE
Generalized Moment Analysis of Magnetic Field Correlations for Accumulations of Spherical and Cylindrical Magnetic Perturbers
- 1Department of Neuroradiology, Heidelberg University, Heidelberg, Germany
- 2Department of Radiology, German Cancer Research Center, Heidelberg, Germany
- 3Department of Experimental Physics 5, University of Würzburg, Würzburg, Germany
In biological tissue, an accumulation of similarly shaped objects with a susceptibility difference to the surrounding tissue generates a local distortion of the external magnetic field in magnetic resonance imaging. It induces stochastic field fluctuations that characteristically influence proton spin dephasing in the vicinity of these magnetic perturbers. The magnetic field correlation that is associated with such local magnetic field inhomogeneities can be expressed in the form of a dynamic frequency autocorrelation function that is related to the time evolution of the measured magnetization. Here, an eigenfunction expansion for two simple magnetic perturber shapes, that of spheres and cylinders, is considered for restricted spin diffusion in a simple model geometry. Then, the concept of generalized moment analysis, an approximation technique that is applied in the study of (non-)reactive processes that involve Brownian motion, allows deriving analytical expressions of the correlation function for different exponential decay forms. Results for the biexponential decay for both spherical and cylindrical magnetized objects are derived and compared with the frequently used (less accurate) monoexponential decay forms. They are in asymptotic agreement with the numerically exact value of the correlation function for long and short times.
When exposed to an external magnetic field as in magnetic resonance imaging (MRI), spatial variations of magnetic susceptibility in heterogeneous systems such as biological tissue induce local inhomogeneities of the magnetic field that are usually visible for macroscopic structures. Susceptibility-weighted imaging and quantitative susceptibility mapping are MRI sequences that make use of this effect, e.g., to better visualize a cerebral thromboembolism or to quantify cerebral blood oxygen saturation [1, 2], and are therefore important in neuroradiological imaging. Susceptibility differences in biological tissue occur at the interfaces between different tissue structures and/or different types of tissue, thus generating local magnetic field gradients that lead to a dephasing of the magnetization. This effect can also be used to obtain quantifiable information on microscopic structures that are not resolvable in routine MRI due to technical limitations or optimization of measurement time. Such structures might be near-spherical as, for instance, iron accumulations in brain cells in neurodegenerative disease  or ultrasmall paramagnetic iron oxide particles that are ingested from macrophages in cell labeling techniques to allow assessment of the extent of cardiovascular disease , nerve injury , or to simply track immune cells . Also, due to the BOLD effect , capillaries containing erythrocytes with deoxygenated hemoglobine possess an inherent paramagnetic susceptibility that is different to that of the surrounding tissue. For such magnetized spherical and cylindrical subvoxel structures, diffusion effects due to spin movements are not negligible. In fact, an increase of spin dephasing through diffusion effects may lead to a characteristic loss of inter-spin coherence: a spin that diffuses around these magnetic field perturbers experiences the local magnetic field inhomogeneities even though the applied external field is homogeneous. Diffusion effects for an ensemble of diffusing spins therefore lead to a more rapidly evolving phase decoherence than that expected from diffusion-free spin relaxation . In addition, combined effects on dephasing of magnetic susceptibility and diffusion gradients are actively investigated, e.g., for diffusional kurtosis imaging . In fact, some authors suggest that anomalous or non-Gaussian diffusion in biological tissue is due to dephasing effects that are promoted by both diffusion multi-compartmentalization and magnetic susceptibility [10, 11].
Previously, local susceptibility gradients and their effect on MR-signal behavior have been examined through spin-echo relaxation rates R2 = 1/T2 and gradient-echo relaxation rates [12–14], but, since magnetization does not necessarily decay monoexponentially ~ exp(−t/T2) or in time , a rightful interpretation of the signal evolution becomes cumbersome. Another approach is given by the frequency autocorrelation function K(t) that describes the transition probability of a spin from one local Larmor frequency state to the next in a specific time . It is both sensitive to proton spin diffusion and the shape of the local magnetic field inhomogeneities. When spin dephasing around the local magnetic field inhomogeneity is Gaussian, a relation between magnetization time evolution and the frequency correlation function can be established  and it could be shown that information about correlation functions provide a more direct and measurable link to the local magnetic field inhomogeneity than relaxation rates . Magnetic field correlation imaging techniques have already been utilized for MR measurements of iron accumulation in brain parenchyma  and for the MR-analysis of porous media . In addition, since Carr-Purcell-Meiboom-Gill (CPMG) sequence relaxation rates are connected to the correlation function , a better knowledge of the correlation function may aid in determining microstrucutral parameters that quantify local capillary size, density, and oxygen extraction fraction in muscle tissue [21, 22].
Spin diffusion leads to a decrease of correlation over time that is mostly assumed to occur monoexponentially with a characteristic time constant or correlation time [23, 24], a simplification that is hardly valid for large parameter regimes. With this in mind, a detailed account of the intricate relation between magnetic field inhomogeneities and the frequency correlation function that characterizes the MRI relaxation process has been provided recently for restricted diffusion in Ziener et al. , and for unrestricted diffusion in Ayant et al.  and Sukstanskii et al. .
To account for the time dependence of microstructural quantities of the diffusion process, generalized moment expansion analysis can be employed to approximate their long time behavior. Generalized moment analysis is an extension of the first passage time approximation  and has been used to study, for instance, (non-)reactive processes that involve Brownian motion [29, 30], finite Ising systems , single-mode lasers , or to generally approximate autocorrelation functions . In fact, the mean relaxation time approximation represents the lowest order approximation of the generalized moment expansion method and could be shown to better approximate the exact correlation function than the high-frequency i.e., short-time approximation used in conventional Padé approximations [30, 33].
In this work, the generalized moment analysis will be introduced and utilized to determine Padé approximants in the biexponential approximation for correlation functions of spherical and cylindrical magnetic objects. This analysis extends and furthers a preceding analysis that only considered the monoexponential decay forms for cylindrical magnetic perturbers .
2. Materials and Methods
The following analysis is based on the consideration of compact and impermeable objects G inside a voxel of size VVoxel that are exposed to an external static magnetic field with strength B0 (see Figure 1A). Assuming homogeneous diffusion properties in the dephasing volume V = VVoxel − G, i.e., no barriers such as membranes or other materials, diffusing proton spins around object G are subject to stochastic field fluctuations induced by the distortion of the magnetic field in the vicinity of G.
Figure 1. Proton spin trajectory around a magnetized object and measurement of the frequency correlation function through asymmetric spin echos. (A) Magnetized object G with susceptibility χi embedded in a surrounding medium V with susceptibility χe; both volumes constitute the voxel volume VVoxel and the volume fraction is η = G/VVoxel. The spin is subjected to stochastic field fluctuations in V induced by a distortion of the magnetic field from object G. (B) The corresponding frequency correlation function K(t) = 〈ω(t)ω(0)〉, where the bra-ket-formalism indicates averaging over all protons inside the voxel, determines the MR-signal S: through measurement of asymmetric spin echos, i.e., spin echos S(t)|ts+t/2 that are shifted with time ts away from the standard spin echo value of t/2 (see main text) we have .
2.1. Frequency Autocorrelation Function
Generally, the distortion of the magnetic field is connected to the susceptibility difference Δχ = χi − χe between the magnetized object and its surrounding medium and leads to a variation in the local Larmor frequency
with frequency shift δω = γB0Δχ and gyromagnetic ratio γ = 2.675 × 108s−1T−1 . Let p(r, r0, t) be the transition probability of a proton spin at point r and time t that started at position r0 and time t = 0. Since, presumably, spin diffusion inside the dephasing volume encounters no barriers, transition probability p(r, r0, t) is nothing less than the Green's function of the diffusion equation
where Δ denotes the Laplace operator and D is the diffusion coefficient. Thus,
and it is evident that p(r, r0, t) only depends on the spatial difference r − r0.
Then, to investigate stochastic field fluctuations around object G, a two-point correlation function K(t) can be introduced
where we made use of Equation (3) in the last step. The probability density of the equilibrium p(r0) = 1/V due to the assumption of homogeneous spin density. The frequency autocorrelation function K(t) is linked to the measured MR-signal S(t) and can be determined through measurement of asymmetric spin echos (see also Figure 1B and ). An asymmetric spin echo S(t)|tp is a spin echo that is measured after a 180° refocusing pulse given at a time tp ≠ t/2 where a refocusing pulse given at time t/2 corresponds to a symmetric spin echo S(t)|t/2. When |tp − t/2| = ts, the following equality in our notation holds for low orders of ts according to Jensen et al. :
Usually, Equation (6) is used in practice to be fitted to MRI data to obtain the correlation function.
2.2. Spherical and Cylindrical Magnetic Objects
In this work, two kinds of magnetic perturbers are considered: spheres and cylinders (see Figures 2A,B, respectively). For instance, an ensemble of microscopic spherical magnetic objects can represent the intracellular accumulation of iron particles in neurodegenerative disease  or the accumulation of ultrasmall paramagnetic iron oxide as a contrast agent to track macrophages in nerve injury or atherosclerotic disease [4, 5]. Collectives of cylindrical magnetic perturbers can be represented through a microvascular network in skeletal or cardiac muscle tissue [36, 37] where capillaries contain paramagnetic deoxygenated hemoglobine .
Figure 2. Spherical and cylindrical object in an external static magnetic field B0. (A) Spherical magnetic object with radius Ri (red), that is surrounded by a dephasing volume sphere with radius R. The corresponding local Larmor frequency can be found in Equation (7). (B) Cylindrical object (red) with radius Ri of the cross section circle perpendicular to the cylinder axis. The cylinder is surrounded by a cylindrical dephasing volume with radius R. The angle θ between external magnetic field and cylinder axis is important in the determination of the local Lamor frequency in Equation (8).
The magnetic field around a spherical magnetic perturber is that of a magnetic dipole, i.e., with the external magnetic field B0 being parallel to the z-axis of the coordinate system and utilization of spherical coordinates (r, θ, φ) (see also Figure 2A), the Larmor frequency reads
Here, represents the frequency shift on the equatorial surface of the sphere with radius Ri. Its dephasing volume V ranges from the sphere surfaces to the surface of the mean relaxation volume. The latter depends on the form of the magnetized object and, therefore, the outer surface of V can be considered as a sphere with radius R. Consequently, the volume fraction of the magnetic material .
For a cylinder, the angle θ between its axis and the external magnetic field B0 has to be acknowledged in the frequency shift to obtain the local frequency (see Figure 2B)
where is the external frequency shift. In analogy to spheres, the local volume fraction for a outer cylindrical surface of the dephasing volume with radius R. In both cases, the rationale of the approach has already been discussed in detail elsewhere [13, 38, 39].
For ensembles of several (microscopic) magnetic objects, it is helpful to avail oneself of the principle of supply areas as in the Krogh capillary model . If the volume hosting the magnetic perturbers can be split into smaller volumes of similar geometry and if the local resonance frequency around the magnetic perturber is point symmetric, i.e., ω(r) = ω(−r) which is valid for both spheres and cylinders, it is possible to think of spin trajectories that leave their initial supply area or dephasing volume as being reflected at its boundary. Also, the same applies for spins whose trajectories hit the surface of the magnetic perturber that is assumed impermeable. Then, the Krogh model allows reducing the problem to one single object, i.e., diffusion is not considered in the whole tissue but, instead, for one single magnetic perturber only and reflective boundary conditions are assumed
The case of unrestricted diffusion for a dephasing volume V extending from the surface of the magnetic perturber to infinity has already been analyzed in detail for spheres [26, 27, 41] and cylinders  and shall not be the subject of this investigation. It corresponds to very small volume fractions η ≈ G/VVoxel where spin diffusion around the magnetic object does not reach the outer surface of V during the measurement. A detailed account in a similar notation can be found in Ziener et al. .
2.3. Eigenfunction Expansion
The diffusion equation (Equation 2) can be solved with the help of an eigenfunction expansion of the transition probability p(r, r0, t):
The eigenfunctions ϕn are solutions of the eigenvalue equation
with eigenvalues κn, and, they fulfill the orthogonality condition
where Ri is the radius of the inner sphere or cylinder, respectively. Thereby, the lowest eigenvalue κ0 = 0 corresponds to the lowest eigenfunction . The influence of diffusion is taken into account by the correlation time τ
By substituting the spectral expansion of the transition probability p(r, r0, t) into the definition of the autocorrelation function (Equation 4), we obtain
with expansion coefficients
Since , the expansion coefficient F0 = 0 does not contribute to the correlation function.
Explicit expressions can be given for both spherical and cylindrical geometries; in the case of spheres, reads
with κn obeying the eigenvalue equation
Functions and are the first derivatives of the spherical Bessel functions of the first and second kind, respectively . An approximation for the first zero κ1 is provided in Appendix A.
For cylinders, we have with cylindrical Bessel functions J2 and Y2 of the first and second kind, respectively,
and eigenvalues κn are the solution of
The long time behavior of K(t) can be determined by considering the first eigenvalue κ1 only:
Evidently, KL(t) converges asymptotically against zero for large t or K(∞) = 0. On the other hand, the value of the correlation function at t = 0 can be obtained by making use of Equation (3) with the initial condition p(r, r0, 0) = δ(r − r0):
The last relation is useful for numerical computation and corresponds to Parseval's theorem in Fourier analysis.
2.4. Generalized Moment Analysis
The generalized moment expansion method or generalized moment approximation (GMA) has been used recently as an effective algorithm to approximate the time dependence of observables that are connected to reactive and non-reactive processes involving Brownian motion [29, 30, 43, 44]. It can be applied for monitoring diffusive redistributions and barrier crossing problems as well as to approximate autocorrelation functions . In our case, the generalized moment analysis will be used to derive analytical expressions for the correlation function for biexponentially approximated decay forms that will prove to be more accurate than monoexponential approximations. The following subsections will introduce generalized moments and the Padé approximation and results from monoexponential approximations, as detailed in Ziener et al. , will be briefly recapitulated to give a buildup to biexponential approximations.
2.4.1. Generalized Moments
In using the Laplace transform of the correlation function K(t) in Equation (5),
it is convenient to expand in a power series of 1/s and s for the short and long time behavior, respectively, as
The high-frequency moments μn hereby correspond to the short time behavior and the low-frequency moments μ−(n+1) to the long time behavior. They are defined as
respectively. Equation (5) allows combining the last two equalities into one expression for positive and negative n as
Therefore, the moments μn only depend on diffusion coefficient D and the form of the distortion of the magnetic field from magnetic perturber G. To avoid confusion, only one expression μn will be used below that refers to high-frequency moments for n ≥ 0 and to low-frequency moments for n < 0. For n = 0, the moment is identical to K(0) in Equation (21) and, therefore, corresponds to . The evaluation of Equation (28) for positive indices is straight-forward, whereas moments with a negative index n can be determined through quadratures, i.e.,
where the moment function μ−n(r) for n ≥ 0 was introduced that satisfies
Thus, moments with negative index can be determined by solving the differential equation Equation (30) where μ−n(r) obeys the same reflecting boundary conditions in Equation (9).
2.4.2. Padé Approximation
The frequency autocorrelation function K(t) can be approximated by the sum of N exponential decays as
with prefactors f(Nh, Nl), n and decay constants Γ(Nh, Nl), n that are related to the moments μn. The sum of the number of high frequency moments Nh and that of low frequency moments Nl equals the 2N parameters to be determined: Nh+Nl = 2N. Then, the Laplace transform of K(t) in Equation (23) can be written as a sum of Lorentzians as
Padé approximants are required to accurately describe the high- and low-frequency behavior of in Equations (24) and (25) to a desired degree and to correctly reproduce the Nh high- and Nl low-frequency moments. The resulting description of the (Nh, Nl)-generalized-moment approximation, denoted as (Nh, Nl)-GMA, is a two-sided Padé approximation around s = 0 and s = ∞ . When introducing the Padé approximation of the frequency autocorrelation function according to Equation (31) into the definition of the high and low frequency moments according to Equations (26) and (27), the parameters f(Nh, Nl), n and Γ(Nh, Nl), n can be obtained through
The generalized moments μn are closely related to the shape of the local field inhomogeneity as shown above and parameters f(Nh, Nl), n and Γ(Nh, Nl), n can be determined by the multi-exponential approximation of the correlation function in Equation (31). Therefore, an approximation of the correlation function is given through the solution of Equation (33), provided that the generalized moments μn are known.
In the following, the correlation function within the GMA is derived for the biexponential approximation (N = 2) in the spherical and cylindrical model geometry, since it is then possible to provide an algebraic solution of Equation (33). However, the numerical solution for N > 2 is still possible by means of an equivalent eigenvalue problem [44, 46]. The case N = 1 was treated in Ziener et al.  and is presented in Appendix B for the sake of completeness.
For N = 2, the approximation of the correlation function can be extended to the sum of two exponential decays and Equation (33) leads to the sum
In this work, two cases are being considered. The first case is the (2, 2)-GMA reproducing the moments μ−2, μ−1, μ0, and μ1 of the exact correlation function. To find the best approximation for the long time behavior, we choose the (1, 3)-GMA that reproduces the moments μ−3, μ−2, μ−1, and μ0, i.e., more low-frequency moments are being taken into account. Introducing the following abbreviations
we find for the (2, 2)-GMA the parameters
Consequently, the correlation function for short times can be approximated as
Similar results have been obtained by Bauer et al.  and Barsky et al. : Equations (3.19–3.22) in Bauer et al.  as well as Equations (B15–B19) in Barsky et al.  are in agreement with the coefficients in Equation (35).
The second case of the (1, 3)-GMA yields parameters
and, consequently, the correlation function for long times reads
Similar expression could be derived by Bauer and Schulten .
Introducing the above results for the biexponential approximation in the (2, 2)- and (1, 3)-GMA into Equation (31) determines the corresponding correlation function at t = 0 :
which is in agreement with the generally valid Equation (21). In the following sections, the time evolution of the two biexponential approximations in the GMA for the correlation function are exemplified for both spherical and cylindrical objects, respectively (see below).
In order to solve differential equation (Equation 30) for spherical magnetic objects, it is helpful to start with the ansatz
that separates the angle dependent part of the local Larmor frequency given in Equation (7). Also, the radial part of the Laplace operator can be written as
with d = 3 for spheres and d = 2 for infinitely long cylinders. This leads to an iterative set of equations for the functions μ−n(r) of the form
The differential equation (Equation 43) can be transformed to an Euler type differential equation. Using the general expression (Equation 41) for the solution of the differential equation (Equation 30), where the local Larmor frequency corresponds to Equation (7), the low-frequency moments given in Equation (28) can be found as
The evaluation of the integral yields the general form
For spherical geometries, the functions μn(r) can be found as the solution of Equation (43). The recursion starts with μ0(r) which results from Equation (44). The functions are given by
We have used the abbreviation
and the constants c1, …, c6 that are given by
where . The functions ln(η) are necessary to determine the moments μn in Equation (46) for negative indices from Equation (45). They are given by
It is convenient to define a ratio τn of successive moments as
As in Ziener et al. , we make use of Koenig's theorem [48, 49] on the zeros of Padé approximants [48, 49]: the theorem states that the series 1/τn converges to the lowest eigenvalue of Equation (11) and , since then
where τ is the correlation time from Equation (13). The lowest eigenvalue κ1 is the solution to the geometry-dependent eigenvalue equation, i.e., Equation (17) for spheres and Equation (19) for cylinders. In Figure 3, a sequence of successive moment ratios τn ranging from n = −3, …, 0 is compared with the exact numerical solution for κ1 for spherical magnetic objects and Koenig's theorem is nicely exemplified for η → 1. The asymptotic value of the first eigenvalue for large volume fractions could be determined as for spheres and for cylinders (from Equations [A14] and [A5] in Ziener et al. , respectively), simply stating that, in this limit, the correlation function decays monoexponentially with time constant . An approximation for κ1 is provided in Appendix A. Figure 3 also validifies the concept of generalized moments to define a search interval for numerical computation of the exact solution from Equations (17) and (19). However, the successive moments curves diverge from the exact solution for small volume fractions.
Figure 3. Dependence of successive moment ratios on the volume fraction for spherical magnetic objects. The ratio of successive moments τμ−n(η)/μ−(n+1)(η) = l−n(η)/l−(n+1)(η) in dependence on volume fraction η is obtained from Equation (46) [functions l−n(η) are given in Equation (50)]. The numerically exact eigenvalue (solid line) can be obtained by solving Equation (17). An approximate solution for (dashed line), based on a Taylor expansion, is given in Equation (A3) in Appendix A. For η → 1, all graphs asymptotically converge to the value . A visualization of the dependence of the ratio of successive moments on the volume fraction for cylinders with a slightly different notation is given in Figure 9 in Ziener et al. . It can be obtained from Equation (46) with quantities δω, τ, and ln(η) determined from the cylindrical geometry.
Different cases of the monoexponential functions Ki(t) (see Appendix B and Equation B3) are being shown in Figure 4A to demonstrate the dependence of the monoexponential approximation on the chosen moments μn. Evidently, the initial value of Ki(t), i ≥ 1, decreases with increasing i while the exponential cutoff begins at longer times for increasing i. Also, the approximants Ki(t) converge to the long time behavior given in Equation (20), as expected from Koenig's theorem. To evaluate the accuracy of the biexponential approximation, a double logarithmic plot has been chosen in Figure 4B. The (2, 2)-GMA with the two lowest high and low frequency moments approximates the short as well as the long time limit sufficiently well. In contrast, the (1, 3)-GMA with the three lowest low and the first high frequency moment approximates the long time limit significantly better than the short time limit. An application to biological tissue is provided in Figure 4C where mono- and bi-exponential approximations in the GMA of the correlation function are visualized for values that are typical for susceptibility gradients based on nonheme iron-storage in brain tissue (through iron-rich oligodendrocytes in the Globus pallidus that approximately have radii Ri of 3–4 μm, a frequency shift of δω = 184 s−1, a tissue density of η = 0.04 and D = 0.76·10−9m2s−1, leading to τ = 16.12 ms; see [50, 51] for more details). It can be seen that the biexponential approximations K(1, 3)(t) and K(2, 2)(t) are in excellent agreement with the actual correlation function K(t), while the monoexponential functions Ki(t) only approximate K(t) sufficiently for t ≳ 30 ms. However, changes in the size Ri of the spherical perturbers affect the correlation time τ as well according to Equation (13). The correlation function in dependence on correlation time τ is therefore shown in Figure 4D at the characteristic time t = 10 ms for the same set of mono- and bi-exponential approximations as in Figure 4C. In this case, the biexponential approximations show an equally well agreement with the numerically obtained correlation function K(τ) for all correlation times, whereas the monoexponential approximations exhibit significant deviations from K(τ) for correlation times that range from 40 ms to 5 s. For very short and very long correlation times, however, the monoexponential approximations K0(τ) and K1(τ) coincide with K(τ). Mono- and bi-exponential approximations as compared to K(t) for nanoparticle experiments based on parameters in Muller et al. , and for lung tissue based on parameters given in Stone et al.  and Baete et al. , are shown in Appendix C.
Figure 4. Mono- and bi-exponential approximation in the GMA of spherical objects. (A) Monoexponential approximations Ki(t) of the correlation function K(t) (solid line) as obtained from Equation (B3) with moments given in Schulten et al. . The long time approximation KL(t) (dashed line) is given by Equation (20) where the first eigenvalue κ1 is determined by Equation (17) (η = 0.001 in A,B). (B) Biexponential approximations for long times K(1, 3)(t) [dashed line obtained from Equation (39)] and for short times K(2, 2)(t) [dotted line obtained from Equation (37)] of the correlation function K(t) [solid line obtained from Equation (14)]. Evidently, for spherical objects, the biexponential approximations are more accurate than monoexponential approximations and show a better agreement for large and short times with the actual correlation function K(t). (C) Mono- and bi-exponential approximations of K(t) for typical parameters of iron-storage in oligodendrocytes of the Globus pallidus, see main text for details. K(1, 3)(t) and K(2, 2)(t) both show an excellent agreement with K(t) whereas the monoexponential approximations exhibit significant deviations from K(t) for t ≲ 30 ms. (D) Correlation function approximations in dependence on correlation time τ at time-point t = 10 ms for otherwise equal parameters as in (C). The agreement of biexponential approximations with K(t) is excellent for all correlation times, whereas the monoexponential approximations mainly agree with K(t) for very short or very long correlation times.
In analogy to spherical magnetic objects, the correlation functions for spins that diffuse between two concentric cylinders can be obtained within the GMA by employing a similar ansatz as the one for a spherical geometry and by using Equation (42) for d = 2 as shown in Ziener et al. . Then, generalized moments are determined by the differential equation
Again, for n = 0, 1, 2, 3, the functions μ−n(r) for cylindrical geometries that obey Equation (53) can be determined as
and the constants c1, …, c6 are given by
The functions ln(η) are obtained in analogy to the spherical case and are given as
As in the case for spheres, the low-frequency moments μ−n for cylinders from Equation (28) can be expressed as
that correspond to Equation (77) in Ziener et al. . Evaluation of the integral leads to Equation (46), but with the respective values of δω, τ, and ln(η) for cylinders (see Equation 56).
In analogy to Figure 4A, the time evolution of the functions Ki(t) and their proximity to the long time approximation KL(t) from Equation (20) is shown in Figure 5A. As in the spherical case, the initial value of Ki(t), i ≥ 1, decreases with increasing i and successive moments converge to the long time approximation KL(t) for d = 2.
Figure 5. Mono- and bi-exponential approximation in the GMA of cylindrical objects. (A) Monoexponential approximations Ki(t) of the correlation function K(t) (solid line) as obtained from Equation (B3) with moments given in Schulten et al. . The long time approximation KL(t) (dashed line) is given by Equation (20) where the first eigenvalue κ1 is determined by Equation (19) (η = 0.01 in A,B). (B) Biexponential approximations for long times K(1, 3)(t) [dashed line obtained from Equation (39)] and for short times K(2, 2)(t) [dotted line obtained from Equation (37)] of the correlation function K(t) [solid line obtained from Equation (14)]. As in the case for spherical objects, the biexponential approximations of the correlation function for cylindrical objects show a better asymptotic agreement with the actual correlation function K(t) than monoexponential approximations. (C) Mono- and bi-exponential approximations of the correlation function for typical parameters of capillaries in brain tissue, see main text for details. The deviations of the biexponential approximations from K(t) are visibly smaller than that of the monoexponential approximations. (D) Mono- and bi-exponential approximations of the correlation function in dependence on correlation time τ at time-point t = 10 ms, for otherwise equal parameters as in (C). The approximations K0(τ), K1(τ), K(1, 3)(τ), and K(2, 2)(τ) coincide with K(t) for very short and very long correlation times, whereas correlation function K(τ) for intermediate correlation times between 1 ms and 1 s is best approximated by the biexponential approximations.
In comparison, the biexponential approximation is shown in Figure 5B: for the (2, 2)- and the (1, 3)-GMA it can be seen that the first approximation describes short and long time limit sufficiently well whereas the second approximation is more accurate in the long time limit. An application to biological tissue is provided in Figure 5C where mono- and bi-exponential approximations in the GMA of the correlation function are visualized for values that are typical for susceptibility gradients between capillaries, that contain deoxygenated hemoglobine, and the surrounding tissue [with approximate values Ri = 4 μm , δω = 200 s−1 at 3 Tesla [12, 56], a tissue density of η = 0.03 , and D = 1·10−9m2s−1 , leading to τ = 16 ms]. It can be seen that the biexponential approximations K(1, 3)(t) and K(2, 2)(t) are in excellent agreement with the actual correlation function K(t) for t ≳ 20 ms, while only the monoexponential functions K2(t) and K3(t) sufficiently approximate K(t) for t≳ 80 ms. Still, the deviation of the biexponential approximations to K(t) for t < 30 ms are visibly smaller than for the monoexponential approximations. However, changes in the size Ri of the cylindrical perturbers affect the correlation time τ as well according to Equation (13). The correlation function in dependence on correlation time τ is therefore shown in Figure 5D at the characteristic time t = 10 ms for the same set of mono- and bi-exponential approximations as in Figure 5C. In this case, the biexponential approximations show an equally well agreement with the numerically obtained correlation function K(τ), especially for correlation times τ ≲ 50 ms, whereas the monoexponential approximations exhibit significant deviations from K(τ) for correlation times that range from 1 ms to 5 s. For very short and very long correlation times, however, the monoexponential approximations K0(τ) and K1(τ) coincide with K(τ), as in the case of spherical perturbers.
4. Discussion and Conclusions
The incentive of this work was to gain a better understanding of the dynamic frequency correlation function of spins that diffuse around microscopic magnetized objects. The analysis is based on an eigenfunction expansion of the correlation function where eigenfunctions and eigenvalues can be determined in the case of restricted diffusion through reflective boundary conditions. A two-sided Padé approximation is then utilized to derive analytical expressions in the short and long time limit. The different regimes of the functional form of the correlation function as well as the case of unrestricted diffusion have been studied in a preceding analysis that considered a monoexponentially decaying correlation function in the generalized moment approximation (GMA) for cylindrical magnetic objects . The here-presented analysis extends and furthers this prior work to include biexponentially decaying correlation functions in the GMA for both high and low frequency moments and both spherical and cylindrical magnetic field perturbers, respectively. Expressions for four generalized moments in both model geometries are provided. They allow to determine the corresponding parameters of Padé aproximants for small and large times. Also, successive generalized moment ratios are shown to asymptotically converge to the exact first eigenvalue of the underlying spherical and cylindrical model geometry eigenvalue (Equations 17, 19), respectively, in agreement with Koenig's theorem.
The model with its boundary conditions is based on a presumed regular arrangement of accumulations of similarly shaped objects that can then be considered in analogy to Krogh's model for supply volumes around capillaries . Thus, large volume fractions or the long time limit of spin diffusion processes can be examined for a point-symmetric local Larmor frequency where diffusion trajectories into a neighboring supply volume are thought of as being reflected at the boundary to said volume. More detailed accounts about the case of unrestricted diffusion can be found in Yablonskiy et al. , Ziener et al. , Ayant et al. , and Ziener et al. . With Figures 4, 5 it can be seen that biexponential approximations in the GMA for K(1, 3)(t) and K(2, 2)(t) and for both spherical and cylindrical model geometries are in much better agreement with the actual correlation function K(t) than monoexponential approximations such as the frequently used “mean correlation time approximation” . Comparisons of the approximations are also provided for iron-accumulations in brain tissue and for cerebral cortex capillaries in Figures 4C,D, 5C,D, respectively. Both cases demonstrate that there are pronounced differences between mono- and bi-exponential approximations for t ~ 1 − 100 ms, with the biexponential approximations providing the better approximations to the correlation function for all time-points. These results are therefore relevant for diffusion MRI experiments applied to these tissue types since the typical diffusion time during which the diffusion process can be detected in such an experiment ranges between 1 ms and 1 s. However, for very short correlation times (e.g., iron-oxide nanoparticle suspensions, see Figure 6A) or very long correlation times (e.g., diffusion in peripheral lung tissue, see Figure 6B), mono- and bi-exponential approximations practically coincide.
Figure 6. Mono- and bi-exponential approximation in the GMA for nanoparticles and lung tissue. (A) Mono- and bi-exponential approximations are compared to the correlation function K(t) [solid line, obtained from Equation (14)]; parameter set obtained from Muller et al. , see also main text. For this extremely short correlation time τ, the monoexponential approximations K0(t) and K1(t) coincide with the biexponential approximations K(2, 2)(t) and K(1, 3)(t), respectively. (B) Mono- and bi-exponential approximations for peripheral lung tissue, based on parameters given in Stone et al.  and Baete et al. . In this case of a very long correlation time τ, most approximations coincide with K(t).
Recently, it could be shown that Carr-Purcell-Meiboom-Gill (CPMG) relaxation rates are connected to the correlation function [16, 59]. Therefore, knowledge of the correlation function allows obtaining microstructural parameters through CPMG measurements for accumulations of similarly shaped subvoxel structures, for example in the case of lung alveoli , microscopic blood products around or in malign cerebral tumors or for microvascular changes in skeletal muscle tissue . Also, based on a recently developed measurement method, magnetic field correlation imaging has been shown promising in its potential to evaluate iron-associated neuropathologies [18, 19]. Knowledge of the exact time course and limit case behavior of the correlation function will aid in further application of magnetic field correlation imaging to determine microstructural parameters of susceptibility gradient based local magnetic field inhomogeneities.
This work was carried out by the seven authors, in collaboration. CZ designed research; CZ, FK, LB, and TK performed research; CZ and FK contributed numerical tools; and FK, TK, LB, HS, MB, SH, and CZ wrote the paper. Moreover, all authors have read and approved the final manuscript.
This work was supported by grants from the Deutsche Forschungsgemeinschaft (DFG ZI 1295/2-1 and DFG KU 3555/1-1). FK was also supported by a postdoctoral fellowship from the medical faculty of Heidelberg University and the Hoffmann-Klose Foundation (Heidelberg University).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
3. Bartzokis G, Tishler TA, Lu PH, Villablanca P, Altshuler LL, Carter M, et al. Brain ferritin iron may influence age- and gender-related risks of neurodegeneration. Neurobiol Aging (2007) 28:414–23. doi: 10.1016/j.neurobiolaging.2006.02.005
4. Richards JM, Shaw CA, Lang NN, Williams MC, Semple SI, MacGillivray TJ, et al. In vivo mononuclear cell tracking using superparamagnetic particles of iron oxide: feasibility and safety in humans. Circ Cardiovasc Imaging (2012) 5:509–17. doi: 10.1161/CIRCIMAGING.112.972596
9. Palombo M, Gentili S, Bozzali M, MacAluso E, Capuani S. New insight into the contrast in diffusional kurtosis images: does it depend on magnetic susceptibility? Magn Reson Med. (2015) 73:2015–24. doi: 10.1002/mrm.25308
10. Palombo M, Gabrielli A, De Santis S, Capuani S. The γ parameter of the stretched-exponential model is influenced by internal gradients: validation in phantoms. J Magn Reson. (2012) 216:28–36. doi: 10.1016/j.jmr.2011.12.023
11. GadElkarim JJ, Magin RL, Meerschaert MM, Capuani S, Palombo M, Kumar A, et al. Fractional order generalization of anomalous diffusion as a multidimensional extension of the transmission line equation. IEEE J Emerg Sel Topics Circuits Syst. (2013) 3:432–41. doi: 10.1109/JETCAS.2013.2265795
14. Kurz FT, Kampf T, Heiland S, Bendszus M, Schlemmer HP, Ziener CH. Theoretical model of the single spin-echo relaxation time for spherical magnetic perturbers. Magn Reson Med. (2014) 71:1888–95. doi: 10.1002/mrm.25196
19. Jensen JH, Szulc K, Hu C, Ramani A, Lu H, Xuan L, et al. Magnetic field correlation as a measure of iron-generated magnetic field inhomogeneities in the brain. Magn Reson Med. (2009) 61:481–5. doi: 10.1002/mrm.21823
21. Kurz FT, Kampf T, Buschle LR, Heiland S, Schlemmer HP, Bendszus M, et al. CPMG relaxation rate dispersion in dipole fields around capillaries. Magn Reson Imaging (2016) 34:875–88. doi: 10.1016/j.mri.2016.03.016
22. Kurz FT, Buschle LR, Kampf T, Zhang K, Schlemmer HP, Heiland S, et al. Spin dephasing in a magnetic dipole field around large capillaries: approximative and exact results. J Magn Reson. (2016) 273:83–97. doi: 10.1016/j.jmr.2016.10.012
23. Ziener CH, Bauer WR, Melkus G, Weber T, Herold V, Jakob PM. Structure-specific magnetic field inhomogeneities and its effect on the correlation time. Magn Reson Imaging. (2006) 24:1341–7. doi: 10.1016/j.mri.2006.08.005
25. Ziener CH, Kampf T, Herold V, Jakob PM, Bauer WR, Nadler W. Frequency autocorrelation function of stochastically fluctuating fields caused by specific magnetic field inhomogeneities. J Chem Phys. (2008) 129:014507. doi: 10.1063/1.2949097
26. Ayant Y, Belorizky E, Aluzon J, Gallice J. Calcul des densités spectrales résultant d'un mouvement aléatoire de translation en relaxation par interaction dipolaire magnétique dans les liquides. J Phys France (1975) 36:991–1004.
27. Sukstanskii AL, Yablonskiy DA. Gaussian approximation in the theory of MR signal formation in the presence of structure-specific magnetic field inhomogeneities. Effects of impermeable susceptibility inclusions. J Magn Reson. (2004) 167:56–67. doi: 10.1016/j.jmr.2003.11.006
34. Salomir R, de Senneville BD, Moonen CTW. A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility. Concepts Magn Reson B (2003) 19B:26. doi: 10.1002/cmr.b.10083
37. Bauer WR, Nadler W, Bock M, Schad LR, Wacker C, Hartlep A, et al. The relationship between the BOLD-induced T2 and T2*: a theoretical approach for the vasculature of myocardium. Magn Reson Med. (1999) 42:1004–10.
38. Bauer WR, Schulten K. Nuclear spin dynamics (l=1/2) under the influence of random perturbation fields in the “strong collision approximation”. Berichte der Bunsengesellschaft - Phys Chem Chem Phys. (1992) 96:721–7.
55. Lauwers F, Cassot F, Lauwers-Cances V, Puwanarajah P, Duvernoy H. Morphometry of the human cerebral cortex microcirculation: general characteristics and space-related profiles. Neuroimage (2008) 39:936–48. doi: 10.1016/j.neuroimage.2007.09.024
56. Dickson JD, Ash TW, Williams GB, Sukstanskii AL, Ansorge RE, Yablonskiy DA. Quantitative phenomenological model of the BOLD contrast mechanism. J Magn Reson. (2011) 212:17–25. doi: 10.1016/j.jmr.2011.06.003
60. Kurz FT, Kampf T, Buschle LR, Schlemmer HP, Heiland S, Bendszus M, et al. Microstructural analysis of peripheral lung tissue through CPMG inter-echo time R2 dispersion. PLoS ONE (2015) 10:e0141894. doi: 10.1371/journal.pone.0141894
61. Ziener CH, Kampf T, Melkus G, Herold V, Weber T, Reents G, et al. Local frequency density of states around field inhomogeneities in magnetic resonance imaging: effects of diffusion. Phys Rev E (2007) 76:031915. doi: 10.1103/PhysRevE.76.031915
62. Kampf T, Helluy X, Gutjahr FT, Winter P, Meyer CB, Jakob PM, et al. Myocardial perfusion quantification using the T1-based FAIR-ASL method: the influence of heart anatomy, cardiopulmonary blood flow and look-locker readout. Magn Reson Med. (2014) 71:1784–97. doi: 10.1002/mrm.24843
63. Buschle LR, Kurz FT, Kampf T, Triphan SM, Schlemmer HP, Ziener CH. Diffusion-mediated dephasing in the dipole field around a single spherical magnetic object. Magn Reson Imaging (2015) 33:1126–45. doi: 10.1016/j.mri.2015.06.001
To determine the eigenvalues in Equation (17), we must find the roots of the function
The derivatives of the spherical Bessel functions j2 and y2 can be expanded in a Taylor series for small arguments of κn and, thus, an expression for small arguments of f(κn) is given as
The first zero of Equation (A1) then approximates the first zero of Equation (17) as
In the monoexponential approximation (N = 1) or (1, 1)-GMA, only two generalized moments have to be taken into account and Equation (33) reduces to
The first ratio τμ0/μ−1 = τ/τ−1 is determined in the (1, 1)-GMA that can be regarded as the “mean correlation time approximation,” since, for m = −1, 0, f(1, 1) = μ0, Γ(1, 1) = 1/τ−1 and with Equation (27)
where Pc(t) = −∂tK(t)/K(0) and τ−1 represents a “mean correlation time.” This is in analogy to the “mean relaxation time approximation” in nuclear spin dephasing [38, 62]. While considering the two high-frequency moments μ0 and μ1, the corresponding (2, 0)-GMA with f(2, 0) = μ0 and Γ(2, 0) = μ1/μ0 = 1/τ0 describes the short time behavior of K(t) whereas the (0, 2)-GMA with the two low-frequency moments μ−2 and μ−1 and parameters and Γ(0, 2) = μ−1/μ−2 = 1/τ−2 is a measure for the long time behavior. Since the MCTA is an approximation of the dynamic monoexponentially decaying correlation function, the deviation in Figure 3 of the first ratio τ/τ−1 reflects a strongly decaying non-exponential correlation function.
It is useful to define monoexponential functions Ki(t) for arbitrary successive moments as
where the index i = 1 uses only the lowest high and low frequency moment and indices i ≤ 0 only high and indices i ≥ 2 only low frequency moments, respectively. In Figures 4A, 5A, the time evolution of the functions Ki(t) is exemplified for both spherical and cylindrical objects, respectively (see below). With Koenig's theorem we find that where KL is given in Equation (20).
The relaxivity of particulate MR contrast media was examined in Muller et al.  to obtain the set of approximate parameters Ri = 50nm, δω = 34·106s−1, η = 2·10−6, D = 2.3·10−9m2s−1, and τ = 1.09·10−6s for iron-oxide nanoparticle suspensions. This extreme example of a very short correlation time was chosen to demonstrate the coincidence of mono- and bi-exponential approximations of the correlation function, as shown in Figure 6A. The biexponential approximations K(1, 3)(t) and K(2, 2)(t) then practically coincide with the respective monoexponential approximations K1(t) and K0(t), as anticipated.
The influence of magnetic susceptibility properties of peripheral lung tissue with near-spherical air-containing alveoli on MR signal decay was examined recently, see e.g. [60, 63]. For a typical set of parameters with Ri = 150 μm, δω = 1500 s−1, η = 0.85, D = 1·10−9m2s−1 (τ = 22.5 s) [53, 54], the correlation function K(t) is shown with mono- and bi-exponential approximations in Figure 6B. In this case of a very long correlation time, mono- and bi-exponential approximations coincide completely with the numerically obtained correlation function K(t), with the exception of K0(t).
Keywords: magnetized sphere/cylinder, magnetic susceptibility, correlation function, diffusion, magnetic resonance imaging
Citation: Kurz FT, Kampf T, Buschle LR, Schlemmer H-P, Bendszus M, Heiland S and Ziener CH (2016) Generalized Moment Analysis of Magnetic Field Correlations for Accumulations of Spherical and Cylindrical Magnetic Perturbers. Front. Phys. 4:46. doi: 10.3389/fphy.2016.00046
Received: 12 September 2016; Accepted: 14 November 2016;
Published: 02 December 2016.
Edited by:Federico Giove, Centro Fermi, Italy
Reviewed by:Silvia Capuani, National Research Council, Italy
Jonathan Doucette, University of British Columbia MRI Science Lab, Canada
Copyright © 2016 Kurz, Kampf, Buschle, Schlemmer, Bendszus, Heiland and Ziener. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Felix T. Kurz, email@example.com