Gaussian Local Phase Approximation in a Cylindrical Tissue Model

In NMR or MRI, the measured signal is a function of the accumulated magnetization phase inside the measurement voxel, which itself depends on microstructural tissue parameters. Usually the phase distribution is assumed to be Gaussian and higher-order moments are neglected. Under this assumption, only the x-component of the total magnetization can be described correctly, and information about the local magnetization and the y-component of the total magnetization is lost. The Gaussian Local Phase (GLP) approximation overcomes these limitations by considering the distribution of the local phase in terms of a cumulant expansion. We derive the cumulants for a cylindrical muscle tissue model and show that an efficient numerical implementation of these terms is possible by writing their definitions as matrix differential equations. We demonstrate that the GLP approximation with two cumulants included has a better fit to the true magnetization than all the other options considered. It is able to capture both oscillatory and dampening behavior for different diffusion strengths. In addition, the introduced method can possibly be extended for models for which no explicit analytical solution for the magnetization behavior exists, such as spherical magnetic perturbers.


INTRODUCTION
The extraction of parameters describing the tissue architecture from MRI measurements is dependent on the correct description of susceptibility and diffusion effects during spin dephasing. While analytical solutions of the underlying Bloch-Torrey equation [1,2] are possible only for simple geometries, approximate methods are often sufficient in practice [3]. The important Gaussian phase approximation, introduced by Hahn [4] to describe the signal behavior inside a linear and uniform magnetic field gradient, assumes a Gaussian distribution of the accumulated spin phases. Recently, the Gaussian Local Phase (GLP) approximation was introduced [5]. Instead of the cumulative phase distribution inside a measurement voxel, the local phase distribution is expanded in a cumulant expansion and the integration over the voxel carried out afterward. Herein, we consider the case of a muscle tissue model which consists of a blood vessel surrounded by a cylindrical dephasing domain [6]. Due to the shape of this domain, the resulting diffusion propagator takes a rather complex form. The calculation of the magnetization time evolution is then dependent on zeros of cross products of Bessel functions, making it impractical to use [3,7]. To overcome this limitation, we present a convenient calculation method for the GLP moments based on a finite-difference approximation of the cylindrical Laplace operator. This permits formulating the first and second moments of the phase distribution as closed-form matrix exponentials, which are easy to implement numerically. We compare the new method to results obtained from integrating over the analytical diffusion propagator and to the numerically exact results from a finite-element simulation [8][9][10][11].
The manuscript is structured as follows. In section 2.1, we define the muscle tissue model and show the diffusion propagator. Integrating over products of the diffusion propagator and the local Larmor frequency allows the analytical calculation of the phase moments, as shown in sections 2.2 and 2.3. In section 3.1, the finite-difference scheme is introduced and used to write the governing differential equations as vector differential equations. The solutions are rather straightforward and result in the vector forms of the phase moments. These are then used to calculate the local magnetization in section 3.2 and the total magnetization in section 3.3. Finally, the introduced methods are compared to previously known solutions described in section 3.4. In Appendix A, the cumulant expansion and its statistical properties are reviewed. Alternative expressions for the cumulant functions are given in Appendix B, and a spectral expansion of the Laplace operator is shown in Appendix C. As the discretized Laplace operator is singular, its group inverse has to be used. This group inverse can be expressed in terms of simple functions, as detailed in Appendix D.

Model and Diffusion Process
We considered a muscle tissue model which consists of a radiallysymmetric dephasing domain parallel to the z-axis and bounded at R C and R D , as shown in Figure 1. As the model is invariant along the z-axis, it is sufficient to parametrize the dephasing domain in terms of polar coordinates (r, φ): r = (r, φ) ∈ [R C , R D ] × [0, 2π). The inner cylinder with radius R C > 0 represents a blood-filled vessel and therefore creates a spatially varying Larmor frequency [12] ω(r, φ) = δωR 2 C cos(2φ) with δω being the characteristic frequency shift at the vessel surface, which is proportional to the susceptibility difference between vessel and dephasing domain. The Larmor frequency is visualized in the upper part of Figure 1. In the following derivations, we assume without loss of generality that the external magnetic field is orthogonal to the cylinder axis. In the case of a different angle θ = π between the cylinder axis and the magnetic field, δω is modified by a factor of sin 2 (θ ) (see Equation (7) in [12]). The geometry can be described by the volume ratio Inside the dephasing domain, spin-1/2 particles can diffuse freely with an isotropic and spatially homogeneous diffusion constant D > 0. The characteristic diffusion time τ for this model is then given by As the muscle model is assumed to be surrounded by identical tissue cylinders, the crossing-over of particles into the domain of neighboring capillaries can be considered by introducing reflecting Neumann boundary conditions at R D . The boundary condition at the interior border R C is dependent on the assumed permeability of the vessel wall: negligible permeability corresponds to reflecting boundary conditions, while non-zero vessel permeability can be taken into account by choosing Robintype boundary conditions, as shown, for example, in Equation (7) in Ziener et al. [13]. The Gaussian diffusion process is described by the propagator p( r, r 0 , t), which is the solution of the diffusion equation with the initial condition Using the definition of the Laplace operator in cylinder coordinates the solution to the diffusion equation can be written as Due to the unboundedness of the radial Laplace operator r on the real line [0, ∞] under the L 2 -norm, the radial part of the operator exponential cannot be readily written as a Taylor series. It can, however, be expanded in terms of eigenfunctions of the Laplace operator [14][15][16][17][18]. The propagator can then be written in the classical form [19,20]: where the radial eigenfunctions n,ν (r) and the eigenvalues λ n,ν obey the cylindrical Bessel differential equation with the radial part of two-dimensional Laplace operator r given in Equation (6) inside an annular ring R C ≤ r ≤ R D .
Reflecting boundary conditions at r = R C and r = R D are assumed: The appropriate eigenfunctions that obey the reflecting boundary conditions (10) at r = R C are linear combinations of cylindrical Bessel functions J ν and Y ν of order ν and of the first and second kind, respectively: for all other combinations of indices. These eigenfunctions obey the orthogonality relation with the normalization constants The lowest eigenvalue takes the value zero, while all other eigenvalues λ n,ν can be obtained by introducing the eigenfunctions from Equation (12) into the reflecting boundary conditions (11) at r = R D as solutions of a transcendental equation The eigenvalues can be calculated in Mathematica (Wolfram Research, Inc., Champaign, IL, United States) using the package BesselZeros. For ν = 0 and n ≥ 2 the correct command is and for ν ≥ 1 and n ≥ 1, the command is Alternatively, the eigenvalues can be obtained as the eigenvalues of a discretized Laplace operator as demonstrated in Ziener et al. [21].

Moment Calculation
The signal generated during a NMR experiment can be calculated by approximating the local accumulated phase distribution by a cumulant expansion as detailed in Appendix A. Taking only the first cumulant into account, as seen in Equation (A7), the local magnetization can be written as which is further denoted as the first-order GLP approximation, with the first moment function denoted by α( r). For simplicity, the initial magnetization m 0 (r) is from now on set to the constant m 0 , corresponding to the voxel state immediately after application of the initial pulse. If both first and second cumulants are included, as seen in Equation (A8), this becomes further denoted as second-order GLP approximation, with the second moment function denoted by β( r, t). The real part of the local magnetization corresponds to the x-component of the transverse magnetization and the imaginary part to the ycomponent. Inserting m I ( r, t) and m II ( r, t) into the transverse Bloch-Torrey equation and collecting real and imaginary terms leads to two inhomogeneous partial differential equations for the moment functions: The initial conditions can be seen from the definitions in Equations (26) and (27): The boundary conditions on the moments similarly arise from the reflecting boundary conditions on m( r, t) for both real and complex parts: The solutions to the differential equations (21) and (22) are also given by integrals over products of the propagator and the Larmor frequency: where the spatial integrals go over the dephasing domain V. This relationship is visualized in the middle and lower subfigure of Figure 1. Here, the diffusion propagators take the form of ordinary functions, but later on they are used in a noncommutative operator form, in which case the specific ordering matters. The form presented here is valid in both cases. Inserting the propagator from Equation (B3) into the definitions of the moment functions in Equations (26) and (27) and performing the angular integration yields an explicit expression for the angular dependence of the moments: with The angular decomposition directly leads to a system of three inhomogeneous partial differential equations for α 2 (r, t), β 0 (r, t), and β 4 (r, t):

Analytical Moments
The moments (30)-(32) can be calculated analytically either by using an ansatz in terms of eigenfunctions of the Laplace operator or by integrating over the diffusion propagator as defined in Equation (8). The first moment is then and the coefficients of the second moment are where the function appearing in the expression for β 0 (r, t) in Equation (37) denotes the first derivative of the exceptional univariate Lommel function, which is given in Equations (2) or (7) in Glasser [22]. The derivation can be found in Equations (B4) and (B5) in the Appendix. The limiting behavior of the moments for long times is easily derived: It is obvious that only the moment β 0 (r, t) has a non-constant limit and instead linearly increases with time. This leads to exponential dampening of the local and total magnetization.  (36) in Bauer et al. [23]. The total magnetization can be calculated as the volume integral of the local magnetization inside the measurement voxel with the initial total magnetization In the limit of η → 1, i.e., dephasing on the surface of the cylinder, the spatial derivative in the governing differential equations is zero. Subsequently, a spatial discretization scheme is unnecessary and the coefficients can be found by use of a simple integrating factor. The resulting functions are Using the fact that all but the first eigenvalue λ n,ν for each ν diverge (see Figure 3 in and the special value of the eigenfunctions at the inner boundary (see Equation (63) in [25]) these limiting cases can also be obtained from Equations (36)-(38).
In order to extend the discretization to the differential equations, the effect of the radial Laplace operator r has to be considered. Applying a finite-differences scheme to the discretization grid from Equation (52) leads to the dimensionless matrix Laplacian The matrix elements r n+1/2 /r n and r n+1/2 /r n+1 can be directly calculated by evaluating the formula for the discretized radial grid points (52) at the respective locations n. In order to formulate the problem as a vector differential equation, it is prudent to consider the vectors being the constant vector of length N. The action of the Laplacian can then be described as product of the finite differencediscretized matrix r and the moments The discretized versions of Equations (33)-(35) are therefore This matrix differential equation system can be solved by first calculating α 2 (t) and subsequently β 0 (t) and β 4 (t) by the use of integrating factors. Multiplying Equation (66) from the left with exp(− t τ r − 4r −2 ) and collecting terms yield Integrating gives with the identity matrix 1 of size N × N. This can be simplified using the relation resulting in Similarly, the attenuation coefficients β 0 (t) and β 2 (t) can be derived by the use of integrating factors of exp(− t τ r ) and exp(− t τ r − 16r −2 ), respectively. For the non-angulardependent part β 0 (t), this leads to The second integral can be solved by using the identity Due to the singularity of r , the last integral cannot be solved by direct matrix inversion. As expanded upon in Appendix D, the necessary pseudoinverse can be expressed in terms of a simple function: Using the same identity from Equation (76), β 4 (t) becomes β 4 (t) [τ δω] 2 = 1 32 The time evolution of the vector moments at the inner and outer borders as well as one middle point of the dephasing domain is visualized in Figure 2. The numerical accuracy of the finite difference scheme and the spectral expansion method can be measured by the root mean squared error between the true values and the numerical results, as shown in Figure 3.

Local Magnetization
The discretized local magnetization is in the GLP first order and in the GLP second order These vector equations correspond to the scalar equations (18) and (19). As in the continuous case, the initial magnetization is set to the constant vector 1 of length N. The time course of the local magnetization reveals a complex relationship between the real and imaginary parts, as shown in Figures 4, 5.

Total Magnetization
For the total magnetization in the discretized case, the Jacobi determinant r in Equation (43) is replaced by the discretized differential area element (see Equation (64) in [26])

resulting in
In the first-order GLP, the angular interval can be solved and the total magnetization can be found by inserting Equation (80) into Equation (84) as with J 0 being the Bessel function of first kind and zeroth order.

Comparison With Other Solutions
In the Gaussian Phase approximation, the total magnetization is equal to where the definition of the eigenvalues λ n,2 is as usual, and the coefficients F n are defined by with the definitions of N n,2 given in Equation (14) and the definition of M − n,2 given in Equation (B8) in the Appendix. Using Equations (127)-(130) from Ziener et al. [25], this can be further simplified to The term linearly dependent on the time again has the same form as in Equation (42). Using the discretized version of the correlation function in Equation (A14), the total magnetization in the Gaussian phase approximation can also be written in terms of the Laplace operator: If the magnetization decay is modeled as monoexponential (ME), different approximations for different diffusion strengths exist. If the diffusion term dominates the susceptibility term, called the motional narrowing regime (ME MN), the formula is [27] M ME MN (t) and in the case of vanishing diffusion, called static dephasing regime (ME SD), it is [27] M ME SD (t) Both introduced solutions, i.e., the GLP first order and the GLP second order, are compared with the true total magnetization as obtained by a finite element simulation of the Bloch-Torrey equation and other approximations in Figure 6. In order to compare the numerical efficiency of the different approaches, 10 calculation runs of the total magnetization were timed for each method implemented in Mathematica using 24 computation threads in parallel. The results can be seen in Table 1. The matrix calculation methods is several times faster than the analytical calculation method. The finite element simulation is by far the slowest method, in part because of the high spatial and temporal resolution needed to avoid instability in cases of lower diffusion strength.

DISCUSSION AND CONCLUSION
Herein, we derive the GLP approximation for a cylindrical dephasing domain under the influence of the susceptibility effects of a blood-filled inner vessel [6]. We, further, present a solution method based on solving a matrix differential equation which arises from a finite-differences discretization scheme for the Laplace operator. As the phase distribution during magnetic resonance imaging is a complicated time-dependent function of the shape of the dephasing domain, of the diffusion strength and of susceptibility effects, it is often assumed to take a Gaussian form [4,28,29]. In the GLP, the ensemble distribution of the local phase is considered instead and the distribution approximated by its first two cumulants. This method of approximating a probability distribution in terms of cumulants rather than moments has been introduced by Kubo [30][31][32]. Derivations can be found in standard stochastic physics textbooks, e.g., in Chapter XVI section 4 in van Kampen [33] or in Cowan [28]. It has been shown that the GLP approximation leads to a better approximation of the total magnetization in the case of a constant gradient for cylindrical or spherical voxel shapes [5].
As can be seen in Figure 6, out of all options, the GLP approximation of second order has the least deviation from the true total magnetization time evolution. Interestingly, the GLP of first order shows greater deviations from the true total magnetization than the Gaussian Phase approximation. A possible explanation for this may be that in the GLP first order, the local magnetization does not converge to zero as expected, as Figures 4, 5 demonstrate. This leads to undampened oscillatory behavior that ultimately results in a large deviation both for the local and total magnetization. The Gaussian Phase approximation suffers from the opposite behavior: the total magnetization is overdampened and decays monotonically to zero over time. The GLP of second order does not suffer from either of these problems: it is able to follow oscillations to negative total magnetization while having proper dampening. This can also be seen from the definition of the local magnetization in Equation (19). If only the first cumulant is taken into account, the purely real coefficient α(r, t) is multiplied with the imaginary unit i in the exponent and its magnitude can only influence the angular frequency of the oscillations of m( r, t). The amplitude of the local magnetization cannot be modulated. This is only possible with the addition of a second cumulant: in this case, both α(r, t) and β(r, t) have non-zero real value in the exponent, which is then able to modulate its amplitude and phase. The local magnetization having constant absolute magnitude does not, however, mean that the total magnetization, corresponding to its domain integral, also has this property: due to destructive (or theoretically constructive) interference, its amplitude can, of course, vary.
From an algorithmic viewpoint, the calculation of the GLP cumulants either requires knowledge of the spectral expansion of the diffusion propagator or, alternatively, a proper finite difference scheme for the Laplacian on the dephasing domain. Both approaches are feasible for bounded domains with rotational symmetry. As Figure 3 demonstrates, the finitedifference approach has the advantage of a monotonically decreasing error in respect to the true solution. The increase in the error for the analytical solution for larger numbers of coefficients is most likely due to the idiosyncrasies of the numerical implementation and may very well not be present if one uses a different implementation of Bessel functions. Both methods, however, are not natively suited for extension to unbounded domains. For the spectral expansion, an increasing number of eigenvalues has to be considered as the dephasing domain becomes larger and larger, until ultimatively the eigenvalue spectrum becomes continuous for a truly unbounded domain. The finite difference scheme can be adapted by either increasing the number of grid points or by adapting the discretization scheme to allow a non-uniform grid with the last point in infinity [34][35][36][37][38]. These problems can most likely be overcome, and may allow extending the GLP approximation method for the cases of dephasing outside of cylinders or spheres [39][40][41][42]. Furthermore, the introduced method of comparing expressions for the cumulants for long time ranges allows writing the coefficient equations in a far simpler form. This may also be applicable for cumulant expansions in other models, such as Equations (70) and (82) in Ziener et al. [5].
In our model, we considered dephasing in the domain around a single cylindrical vessel. The influence of neighboring cylinders and possible crossing-over effects was included by imposing reflecting boundary conditions at the outer border. This assumption, however, is true only for an uniform arrangement of capillaries. The signal evolution from randomly arranged capillaries differs significantly. Yablonskiy and Haacke derived the relaxation rates of randomly arranged capillaries or spheres in the absence of diffusion and were able to find an expression for the relaxation rate that is based on a single capillary or sphere, similar to the uniform case [43]. It is unknown, however, if this holds for the non-negligible diffusion effects-numerical simulations of randomly arranged capillaries suggest that this is indeed the case [44], but no proof has been found.
Our model neglects the signal originating from inside the vessel. This is sensible if the volume ratio η is sufficiently low, as the total signal is, of course, the weighted sum of the intra-and extravascular signals. For very short experimental times, the flow inside the vessel can be neglected, and the protons can be assumed to be stationary. The intravascular magnetic field is then homogeneous [12], leading to a complete refocusing of the signal in spin echo experiments. For longer experimental times, the intravascular flow can no longer be neglected. Due to the low laminar flow speeds in the human capillaries, electrodynamic effects are most likely irrelevant, and the intravascular compartment can be modeled separately. To our knowledge, no systematic investigation of such effects has taken place.
The main advantage of GLP approximation and simultaneously its weak point arises from its handling of the interplay of diffusion and susceptibility effects: the diffusion propagator can be calculated separately, usually as the solution of the diffusion equation (or more generally the Fokker-Planck equation in case of a spatially variant diffusion constant [45]) and applied to the susceptibility term at each jump point, once for the first cumulant and twice for the second cumulant [46]. This places very few restrictions on the diffusion propagator, and its calculation does not have to be adapted to the form of the dipole field. Therefore, flow effects and anisotropic diffusion coefficients can be easily considered. Adapting the full Bloch-Torrey equation to such circumstances is rather cumbersome [47].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
LR and CZ performed the modeling and wrote the manuscript. CZ provided the initial concept of the article. EW contributed to the mathematical concepts described in the article. EW and H-PS provided the input on organization and presentation of models. All authors contributed to the manuscript revision, read, and approved the submitted version.