Magnetic Resonance Characterization of General Compartment Size Distributions

INTRODUCTION The influence of molecular diffusion on the MR signal can be exploited to estimate compartment size distributions in heterogeneous specimens [1]. There has been recent interest in adopting this idea to characterize axon diameter distributions in white matter [2]. To this end, previous studies have assumed a known statistical distribution of compartment sizes (such as log-normal or gamma distributions). However, neural tissue does not necessarily conform to such parametric distributions. In this study, we describe an alternative strategy to measure all moments of the compartment size distribution directly from a single diffusion-weighted MR signal decay curve, hence obviating the assumed existence of a known parametric size distribution.


Introduction
In heterogeneous media, one frequently encounters molecules of one kind trapped within isolated compartments of the host medium. For example in porous media, the vacant regions (pores) in a solid matrix may be filled by a fluid. In emulsions, globules of one fluid exist in another immiscible liquid, yielding a distribution of droplets. The physical properties of the resulting heterogeneous medium critically depend on the distribution of pore or globule sizes. Of the numerous methods for measuring compartment size distribution, nuclear magnetic resonance (NMR) offers significant advantages, probing pore geometry in large domains noninvasively. NMR's exquisite sensitivity to molecular diffusion makes it particularly useful in restricted domains. For example, one approach, originally proposed by Packer and Rees [1], assumes a known statistical distribution of compartment sizes, such as a log-normal distribution, and uses NMR-based molecular displacement measurements to infer the parameters of that distribution. Variants of this approach have been used in examining food products [2,3], oil suspensions [4,5] and most recently in characterizing axon diameter distributions in whitematter fibers of the brain [6].
Alternative experimental techniques, exploiting the internal magnetic field induced by susceptibility differences within the porous medium, have also been proposed [7,8]. However, our approach to determining compartment size distributions involves the sensitization of NMR acquisitions to diffusion via the existence or application of external magnetic field gradients; hence, the technique can be used in the absence of significant susceptibility differences. One such realization of the experiment is depicted in figure 1, in which a pair of pulsed-field gradients [9] of strength G, duration δ and separation are incorporated into a stimulated echo NMR experiment [10].  Figure 1. The pulsed-field gradient stimulated-echo NMR pulse sequence. The oscillatory-shaped radiofrequency pulses rotate the magnetization by 90 • . The rectangular pulsed-field gradients encode the displacement of spins.
In this paper, we describe a strategy to measure all moments of the compartment size distribution associated with a medium comprising a collection of non-connected pores directly from a single NMR signal attenuation profile, hence obviating the assumed existence [1]- [6] of a known parametric size distribution. Compared to other studies that have attempted to bypass the assumption of an underlying parametric distribution [11]- [14], our approach does not employ the Gaussian phase approximation; hence, it can be used when a significant attenuation of the NMR signal is present. Furthermore, our technique can be used for geometries comprising planar and cylindrical as well as spherical pore shapes, albeit a mixture of these pore geometries is not allowed. To this end, theoretical expressions linking the moments of the compartment size distribution to the NMR signal attenuation are derived. Next, we propose to use a complete basis of orthogonal functions to represent the general NMR signal attenuation profiles that makes it possible to estimate the moments of the compartment size distribution accurately and efficiently. We verify these derived expressions and the accuracy of the estimation method via simulations. Finally, we apply the method to real data obtained from controlled distributions of water-filled microcapillaries.

Theory
When an ensemble of pores with different sizes are examined via NMR, each pore's contribution to the overall signal is proportional to the number of spins that reside within that pore. In this paper, we consider ensembles of parallel plates, coherently oriented cylinders and spheres corresponding to the cases D = 1, 2 and 3, respectively, where D denotes the number of dimensions along which the diffusion process is restricted. The gradients are assumed to be applied perpendicular to the walls. If the size of each pore can be characterized by the length, a, then the total signal from the ensemble, as a function of q = (2π) −1 γ δG, where γ is the gyromagnetic ratio, is given by . (1) Here E(q, a) is the signal attenuation for a single pore of size a, and f (a) is the pore size distribution function whose features we are interested in extracting. The ensemble average of a quantity X (a) is defined via the relation 4 Clearly, the numerator and denominator in equation (1) can be written as a D E(q, a) e and a D e , respectively. Given the E(q) data for the ensemble, estimation of the nth order moment of the pore size distribution, a n e , is possible by applying an operator Q n that satisfies the conditions and Q n E(q, a) = ζ n a n−D , where ζ n is a constant. When these conditions are met, the nth moment of the pore size distribution is given by Therefore, the problem of calculating all moments of the pore size distribution is reduced to finding a class of operators Q n that satisfy equations (3) and (4) for n = 0, 1, 2, . . . and n = D.
In this paper, we shall consider the long diffusion time (D 0 a 2 max ) and the narrow pulse (D 0 δ a 2 min ) regimes, where D 0 is the bulk diffusivity of the fluid within the compartments, and a max and a min are, respectively, the sizes of the largest and smallest compartments in the ensemble. Under these experimental conditions, three different operators, denoted by I, D and K, provide all moments of a general pore size distribution. These operators are defined through their action on an arbitrary function F(q) to be where β 1 = 2, β 2 = 2π and β 3 = 4π , and To build a physical understanding of these operators, when applied on an arbitrary NMR signal attenuation profile E(q), we recall that E(q) can be envisioned to be the characteristic function of a displacement distribution because of the relationship [15]- [17] P(x) = dq e i2πq x E(q), where P(x) is an ensemble-average propagator indicating the likelihood that the spins will undergo a displacement, x, along the gradient direction. More generally, similar definitions 5 for higher dimensional spaces are possible, in which case q and x will be replaced by vectors. See [18] for discussions about the relations between propagators in different dimensional spaces. For example, using the three-dimensional version of equation (9), the action of the operator R 3 in isotropic environments can be shown to yield the return-to-origin probability [19]. Similarly, the result obtained by the application of R 2 on axially symmetric environments, when the gradient orientation is perpendicular to the symmetry axis, will yield a return-to-symmetry axis probability. Finally, the operator R 1 reveals the probability that the spins will return to their original plane whose normal is parallel to the gradient direction. The class of operators D l has a well-defined meaning as well. Taking the Fourier transform of both sides in equation (9) and expanding the exponential in a Taylor series, it is straightforward to show that the operator D l returns the 2lth-order moment of displacements, x 2l . It is not clear whether the values returned by the operators K m can be interpreted in a similar way. However, as we show below, they complement the information obtained from the preceding operators.
We now consider ensembles of parallel plates and cylindrical and spherical pores, and apply the above operators to the predicted NMR signal attenuation profiles. The strategy we employ involves the application of R D , which yields the denominator in equation (1). Subsequently, all of the remaining operators are applied, each of which returns one of the remaining moments.

Parallel plates (D = 1)
We shall start by considering an array of parallel slabs (D = 1) with variable spacings between consecutive plates whose normals are along the gradient direction. The NMR signal attenuation for a single pore with separation L is given as [20] where the total NMR signal can be obtained by inserting this expression into equation (1) with D = 1 and a = L. Applying the operator R 1 , followed by the other operators D and K, all moments of the distribution of separations are obtained as Next, we consider a pack of coherently oriented cylindrical tubes (D = 2) of varying radius. The NMR signal for a single tube of radius r 0 is given as [21] A general distribution of cylinder radii, f (r 0 ), can be incorporated into equation (1), with D = 2 and a = r 0 . Starting with the operator R 2 , all operators are applied on the resulting NMR signal attenuation, E 2 (q), revealing all moments of the distribution through the following relationships: The signal for a single spherical pore (D = 3) of radius R 0 is given as [20] For an ensemble of spheres whose radii are distributed according to f (R 0 ), all moments of this size distribution are obtained via the application of the operators starting with R 3 . The moments of f (R 0 ) are then given by The derivation of the above expressions for the moments is outlined in the appendix.

Estimation of the moments from sampled data
The above derivations establish the relationships between the moments of the compartment size distribution and a general NMR signal attenuation profile, E(q), which is inherently assumed to be a continuous function that extends to infinity. Since the real data are sparsely sampled and truncated along the q-axis, a method that analytically represents and approximates the discrete data would be useful. However, traditional methods, such as simple function fitting, fail to describe general E(q) profiles, while cumulant expansions suffer from a finite radius of convergence and are very inaccurate at large q-values. Instead, we employ the recently introduced one-dimensional simple harmonic oscillator (SHO)-based reconstruction and estimation (1D-SHORE) method [22], which overcomes these difficulties and is outlined here.

The one-dimensional simple harmonic oscillator-based reconstruction and estimation technique
In the 1D-SHORE approach, the diffusion-weighted NMR signal intensity is expressed as with where H n (x) is the nth-order Hermite polynomial and u is a characteristic length. The NMR signal attenuation, E(q) = S(q)/S(0), can also be expressed in this basis as where a n = a n S 0 (19) and S 0 = S(0) is the non-diffusion-weighted signal, which can be estimated from the coefficients a n : The functions φ n are based on the well-known eigenfunctions of the quantum-mechanical Hamiltonian for the SHO, which form a complete orthogonal basis for the space of squareintegrable functions [23]. However, our definitions of the eigenfunctions are slightly different from their forms as commonly used in quantum mechanics. Specifically, our basis is not normalized, but the scaling is such that when diffusion is Gaussian, and u 2 = 2D 0 ( − δ/3) = x 2 , the coefficients are given by a n = δ n0 , where δ i j is the Kronecker delta and D 0 is the diffusion coefficient. Moreover, the phase convention dictated by the factor i −n in equation (17) ensures that the real and imaginary parts of E(q) are symmetric and antisymmetric, respectively-a necessary condition for the propagator, P(x), to be real when a n are real. Despite these minor differences from the basis used in quantum mechanics, our basis functions still satisfy the relationship where A is the 'lowering operator' defined by We note that Assemlal et al [24] have represented the q-space signal in a basis of Gaussian-Laguerre functions and subsequently showed that these functions can be related to the basis functions we have used in this study.

Implementation
The problem of transforming a discrete set of measured signals into a set of coefficients can be cast as a set of linear equations. Given an M-dimensional vector of signal values, S, with components S m = S(q m ), an M × N dimensional matrix, Q, whose components are given by Q mn = φ n (u, q m ), can be computed. The problem is reduced to a matrix equation S = Qa , where a is the N -dimensional vector of a n coefficients. In our implementation, we solve this equation by computing the pseudoinverse of Q using singular-value decomposition [25]. Next, S 0 is computed using equation (20), which is inserted into equation (19) to determine the coefficients a n . Note that the above estimation scheme assumes a prior estimate of u. Although the completeness of the basis functions ensures the convergence of the estimates for any u, the rate of this convergence may depend on the particular choice of u. Therefore, in our implementation, we first estimate a maximum value for u obtained from the first few points of S(q) data under the assumption of Gaussianity of the signal. The basis functions are designed such that, when the signal is Gaussian, all coefficients except a 0 vanish yielding the expected behavior E(q) = exp(−2π 2 q 2 u 2 ). Starting from this estimate of the maximum u, we gradually reduce the u-value. At each value of u, we estimate the a n coefficients using the scheme described above and compute the signal attenuation values at the data points, which will be denoted by E est (u, q). Subsequently, the mean-squared error defined by where E data (q i ) are the original data points, is computed. The search for the optimal u continued until (u) falls below 1 × 10 −15 , or a local minimum is achieved, whichever comes first. The last set of u and a n values is retained for subsequent analysis.

Actions of the R, D and K operators in the simple harmonic oscillator basis
Once the a n coefficients and u are estimated as described above, the actions of the compartment size distribution operators are conveniently expressed in the SHORE basis via the relations n! a n n m=0,2,...
The derivations of the above expressions are straightforward but tedious and left out for brevity. All derivations involved writing the Hermite function in the definition of the basis functions as a power series and then carrying out the integrations term by term. Note that the SHO basis has the remarkable property that the lowering operator, A, naturally arises within the integrand of the definition for the K m operator. Therefore, the identity in equation (21) was used in the derivation of equation (24e). The action of K m for m > 1 was not evaluated because they concern highorder moments of the compartment size distribution, which are unlikely to be feasible or useful in practice.

Simulations
To validate the derived expressions for the moments and to assess the accuracy of the 1D-SHORE technique, we simulated distributions of planar, cylindrical and spherical pores. We assumed that the pore sizes were distributed according to a beta density function because, unlike in the case of more popular functions such as log-normal or gamma distributions, the beta density function is able to accommodate negatively skewed distributions. Here we include the results obtained for cylindrical pores for brevity. To ensure that the radii of cylinders are distributed according to a beta distribution, the interval between 0 and 20 µm was divided into The black curves and symbols correspond to the beta distribution with positive skewness (α = 3.5, β = 6.5), whereas the red color is used for the negatively skewed distribution (α = 6.5, β = 3.5).
10 000 segments. The E(q) profile corresponding to each of these radii was computed, then multiplied by the value of the beta distribution at that radii, where B(α, β) is the beta function and C is the maximum value of the radius, taken to be 20 µm in the simulations. Two such distributions of radii are considered: one with (α, β) = (3.5, 6.5) and the other with (α, β) = (6.5, 3.5). The inset in the left panel of figure 2 illustrates these two probability distribution functions, where the first distribution is drawn in black, while the second distribution is depicted in red. A total of 81 data points were generated for each ensemble of pores, and the resulting E(q) profiles are illustrated in the left panel of figure 2. The moments were computed via the 1D-SHORE methodology using equations (24a)-(24e) and (13a)-(13d). The estimated moments as well as the ground truth values obtained from equation (25) are plotted in the right panel of figure 2. It is clear that excellent agreement is achieved between the predicted and estimated values for the moments in both cases. Obviously, the first moment yields the mean value, i.e. µ = r 0 e , while the standard deviation of the ensemble can be obtained via the expression σ = ( r 2 0 e − µ 2 ) 1/2 . For the first distribution, the predicted values were 7.0 ± 2.9 µm, where the estimates yielded 7.3 ± 2.6 µm. For the second distribution, the predicted and estimated values were 13.0 ± 2.9 and 13.2 ± 2.3 µm, respectively.
It may be possible to estimate the higher-order descriptors of the underlying compartment size distribution as well. For example, the skewness of the distribution can be computed via Its expected values are 0.35 and −0.35 for the two distributions. The estimates from the moments yielded γ 1 values of 0.46 and −3.41. Clearly, although the moments are very accurate, measures that depend on the higher-order differences between these moments are not. Therefore, we expect that these higher-order descriptors cannot be reliably estimated from real data although the method performs reasonably well for the mean and standard deviation.

Experiments on size distribution phantoms
All experiments were performed on an 8.4 T Bruker NMR spectrometer equipped with a Micro5 probe capable of producing gradient strengths up to 1900 mT m −1 in each direction. The preparation of the size distribution phantoms is described in detail in [26]. Briefly, hollow microcapillaries with different inner diameters (ID) of 10 ± 1, 13 ± 1, 17 ± 1, 19 ± 1, 20 ± 1, 21 ± 1 or 29 ± 1 µm (Polymicro Technologies, Phoenix, AZ, USA) were immersed in water for several days for filling. Then, one phantom comprising only 19 ± 1 µm tubes (SD000) and three size distribution phantoms (SD001, SD002 and SD003) were prepared by mixing differentsized microcapillaries in one NMR tube. The size distribution phantoms were designed to have variations in both the mean diameter and the breadth of the distributions. The microcapillaries were carefully counted to ensure an accurate ratio of fibers needed in each specimen (the volumetric ratios in detail can be found in table 1 of [26]). Note that the microcapillaries were packed into a 4 mm glass sleeve, which was then inserted into a 5 mm NMR tube. The 5 mm NMR tube was aligned with the main axis parallel to the z-direction of the NMR magnet. Typical line widths of 4-20 Hz were obtained after shimming for all size distribution phantoms. Pulsed-field gradient experiments were performed using the stimulated echo sequence with the following parameters: 48 q-values were collected with a maximum gradient strength of 1600 mT m −1 and with /δ = 150/3 ms, resulting in a maximum q-value of 204.3 mm −1 and the number of scans of 32.
The top panel of figure 3 illustrates the inner diameter distributions weighted by the number of spins that reside within tubes of their respective IDs, and the predicted ID values. Note that these values are obtained by assuming that the nominal ID values provided by the manufacturer are correct. Previous experiments in similar microcapillaries [27]- [30] have demonstrated that there may be some global offset in the actual ID, while no evidence was found to suggest a significant spread within a specimen prepared from tubes of the same ID. Therefore, the standard deviation values of 1 µm reported by the manufacturer seem to represent the uncertainty rather than the variation within one specimen. On the other hand, the standard deviations that are tabulated in figure 3 quantify the spread of the distributions. Therefore, in the SD000 specimen, we report this value to be 0. We would like to stress that the predicted values may exhibit some inaccuracy due to any offset in the actual ID values and because some capillaries may be partially filled with water. Figure 3 also demonstrates the real data and the 1D-SHORE fit to them. All fits appear to be accurately approximating the signal attenuation curves, and the extrapolations seem to be acceptable. The mean and standard deviation estimates obtained using the theoretical framework introduced in this study are overlaid on the graphs. In the SD000 specimen, the µ 2 value was slightly larger than r 2 0 e . This observation is consistent with our predicted value of 0 for the standard deviation. The deviations between the predicted and estimated values may be due to the effects discussed above. However, comparing the results for different phantoms, it is clear that the trends in the mean and standard deviation values are consistent with the trends in their predicted values. Therefore, we conclude that it may be possible to create contrast based on the mean and standard deviation values of compartment size distributions using the methodology of this study. Estimated ID = 1 m µ 5.9 ± 2.9 Estimated ID = 1 m µ 5.0 ± 3.3 Figure 3. The top panel shows the histograms illustrating the composition of the four different specimens of water-filled microcapillaries. Each color represents one of these specimens. The SD000 specimen contains nothing but the 19 µm tubes. The corresponding real NMR data and the 1D-SHORE fit to them are depicted in the four panels below.

Discussion and conclusion
Accurate estimates of the moments are possible when the narrow diffusion gradient pulse duration (δ) and long pulse separation ( ) conditions are met. Further, a relatively dense sampling of q-values, large enough to create a significant attenuation of the NMR signal, is necessary. These conditions can be met, for example, in the fringe field of superconducting magnets or in currently available low-cost mobile NMR systems. When the signal is acquired in a more traditional setting that involves a homogeneous static magnetic field, pulsed-field gradients can be employed to sensitize the signal to diffusion. In such experiments, the longcondition is usually simple to fulfill for microscopic pores. However, it may be more challenging to meet the short-δ condition, the violation of which is expected to yield an underestimation of the moments. This deviation introduced by the diffusion pulse duration is expected to compete against the effect due to the finite range of q-values sampled. Therefore, the value of the mean compartment size may be overestimated or underestimated depending on the experimental design and the underlying pore size distribution.
A similar situation may arise when the surface relaxivity of the restricting walls, which is ignored by our method, is significant. In such environments, molecules may lose their magnetization when they get in contact with the walls. Since such collisions occur more frequently within smaller pores, their contribution to the aggregate signal may be reduced, hence introducing a bias towards larger estimates for the pore size. This effect is expected to be significant when the specimen comprises a broad size distribution and very small compartments. A competing effect becomes apparent when the signal attenuation for a single pore is considered. Although surface relaxation effects tend to attenuate the true signal, the signal attenuation value, E(q), may be increased, especially at smaller q-values [31]. In some cases, this effect may act against the effect of surface relaxation mentioned earlier.
As mentioned in the introduction, this formalism can be used to characterize compartment size distributions when the pores are isolated. As such, the method can be applied readily to study specimens such as emulsions and many food products. However, in more complex environments such as axonal fiber bundles in white matter in the nervous system, the accuracy of the results provided by the method needs to be investigated. One particular concern is the presence of extracellular medium in white matter. In some studies aiming to characterize axon diameter distributions, the diffusion of water molecules in the interstitial space has been assumed to be Gaussian [6]. If this assumption is to be employed, the Gaussian component can be removed from the signal profile [32] prior to the estimation of the moments using the introduced technique. For very densely packed cylindrical fibers, the space between the cells can be assumed to be isolated pores as well. Since the size of these 'pores' scales with the size of the fibers, the estimated moments can still be used to create contrast that depends on the features of the underlying fiber diameter distribution.
Note that we did not attempt to recover the distribution function from its moments in this study. However, we note that this problem, commonly referred to as the 'classical moment problem', is well studied [33]- [35] yet tricky [36], and care must be exercised in such an endeavor. Nonetheless, we recognize that much can be learned about the underlying pore size distribution from an examination of its first few moments.
In conclusion, we presented a method that yields all moments of a general compartment size distribution uniquely from a single NMR signal attenuation profile, E(q), without the need to solve an ill-conditioned problem, invoke an underlying parametric compartment size distribution or employ empirical statistical methods [37]. This formalism provides new insights into the relations between the features of a compartment size distribution and the NMR signal obtained from it. It is clear that different characteristics of diffusion, such as zero-displacement probabilities and the moments of molecular displacements, lead to different moments of the compartment size distribution. Further, these associations are different for pore spaces of different dimensionalities.
In the derivation of the expressions for the moments (equations (11a)-(11c), (13a)-(13d), (15a)-(15e)), the values returned by the operators R k are obtained by evaluating the corresponding integrals in equation (6). The action of the operators D l is given by calculating the moments of the average propagator evaluated using equation (7).
Unlike in the cases of R k and D l however, the action of the operators K m is non-trivial. Here, we provide a sketch of the derivations of equations (11c), (13d) and (15e), which are all obtained by applying the operator K m on the NMR signal attenuation functions.
For the purposes of the derivation, we redefine the pore size (a) as equal to L, 2r 0 and 2R 0 for the cases of parallel plates (D = 1), cylinders (D = 2) and spheres (D = 3), respectively. Then, we can define a useful dimensionless variable v = πqa. Note that both a and v are different for different values of D. However, we do not include D in a subscript for brevity. For all the three geometries considered, the NMR signal attenuation from a single pore of size a is given by the unified expression with C 1 = (π/2) 1/2 , C 2 = 2 and C 3 = (9π/2) 1/2 . Note that the following relationship holds [38]: 1 v