Magnetic resonance assessment of effective confinement anisotropy with orientationally-averaged single and double diffusion encoding

Porous or biological materials comprise a multitude of micro-domains containing water. Diffusion-weighted magnetic resonance measurements are sensitive to the anisotropy of the thermal motion of such water. This anisotropy can be due to the domain shape, as well as the (lack of) dispersion in their orientations. Averaging over measurements that span all orientations is a trick to suppress the latter, thereby untangling it from the influence of the domains' anisotropy on the signal. Here, we consider domains whose anisotropy is modeled as being the result of a Hookean (spring) force, which has the advantage of having a Gaussian diffusion propagator while still retaining the fact of finite spatial range for the diffusing particles. Analytical expressions for the powder-averaged signal under this assumption are given for so-called single and double diffusion encoding schemes, which sensitize the MR signal to the diffusive displacement of particles in, respectively, one or two consecutive time intervals.


I. INTRODUCTION
Magnetic resonance has proved to be an extremely effective tool to peer into materials and tissues noninvasively. It manipulates the magnetic orientation of molecules pervading the material into a natural precession which emits radio waves. By imposing a spatially-varying magnetic field (hence frequency of precession), the emitted radio frequency signal is made to encode in its spectrum the coordinates of the molecules that emit it. This way, the signal can be spectrally decomposed to trace back what proportion of it originates from where, that is, from which 'voxel.' Commonly achieved voxel sizes are in the neighborhood of a millimeter.
Another use of a spatially-varying precession rate involves sensitizing the signal to the motion of the molecules; diffusion in particular. When the molecules trace out random (Brownian) paths where different locations they visit impart different precessional angles on them, their precessions lose coherence, attenuating the sum of their emitted radio waves. This attenuation of MR signal, specifically its response to the direction of the gradient in precession rate, can reveal or quantify how mobile the molecules are along different directions. The difference of mobility can arise from pore boundaries, impurities, cell membranes, etc. Hence the diffusion-attenuated signal encodes the influence of structures that can be significantly smaller than voxel dimensions. Tracing out axon bundles in human brain white matter is for instance a widely employed application of this principle. This modality of magnetic resonance imaging, which we refer to as diffusion MR, is the subject of this contribution.
While the anisotropy (i.e., variance under a rotation transformation) of a single pore or a cell may be easily visualized, one would be mistaken to make a one-to-one connection with that and the anisotropy of the signal (i.e., the response of the signal to orientations of the specimen or the apparatus). For instance, the signal of a voxel consisting of an unaligned mixture of cylindrical aqueous compartments will be less sensitive to rotations than one consisting of an aligned bundle of cylinders. The anisotropy of the individual compartments is common in the two examples, but the aligned case has more ensemble anisotropy.
In some sense, then, ensemble anisotropy confounds compartment anisotropy at the signal level, and eliminating it pronounces features at the subvoxel level. One way to achieve this is to take an average of the signal over all orientations. If the material allows it, it can be ground into a powder to that effect; hence the term powder averaging. However this is generally impossible in bio-medical applications. Then, repeated applications of a measurement protocol in different orientations is the avenue to follow.
In this contribution, we are concerned with two particular diffusion MR schemes, single and double diffusion encoding (SDE and DDE), orientationally averaged in the aforementioned fashion to eliminate ensemble anisotropy. The single encoding scheme [1] is the bread and butter of most applications, employing a magnetic field gradient that remains on for a specified duration in one direction, and then the opposite direction after a specified delay (Fig. 1). The signal then encodes the probability of Brownian displacement between the application of the two pulses.
The double diffusion encoding (DDE) scheme employs two single-encoding blocks in succession [2], in different directions in general, before the signal is read (Fig. 2). This method has been studied extensively in recent years mostly because it allows the anisotropy of the microdomains to be quantified [3], which has important implications for medical imaging. The reader is referred to the reviews [4][5][6][7][8] for in-depth presentation of the method. In a nutshell, the DDE technique encodes into the signal the joint probability of two Brownian displacements taking place between the two pulses of each block. For freely diffusing molecules, the two displacements are uncorrelated [9]. However, when restrictions, and arguably inhomogeneities and forces, are present, this is no longer the case [10], hence imparting signatures of compartment size onto the signal. As mentioned above, DDE employing gradient blocks in different directions is sensitive to anisotropy of subdomains inside the voxel [3,11], but not independently of ensemble anisotropy [12].
While we do not consider restricted diffusion in the strict sense of restriction by hard walls in this article, we do respect the finite range of motion of the molecules by the aid of a harmonic attractive force [13][14][15][16][17]. This is an effective model which can mimic restrictive walls with reasonably tractable mathematics, and is actually approached when the gradient pulses are long [18].
This chapter is organized as follows. We first present the signal arising from an effectively confined domain under double diffusion encoding. Afterwards we derive analytical expressions for its powder average. We then treat the case of single diffusion encoding as an extreme case of double encoding and give its corresponding powder average expressions. The special cases of "stick" and "pancake" geometry are compared to their counterparts for free diffusion. The article is concluded after discussions based on the results of the previous sections.

II. DOUBLE DIFFUSION ENCODING AT THE COMPARTMENT LEVEL
Here we derive the double-encoded diffusion MR signal arising from a compartment which is characterized by an effective spring force attracting the molecules toward the center. Such a force mimics the confining effect of walls, membranes, etc., with minimal mathematical burden.
The spring (Hookean) force influences the motion of the molecules via a quadratic potential1 defining the confinement tensor C, which can be taken to have C = C without loss of generality. However, in order for the steady state molecule number density p st (x) ∝ e −V (x) to be normalizable, C must be positive (semi)definite.2 Under this (or any) potential, the magnetization density ρ(x, t) evolves according to Here, we have assumed that diffusion is governed by a spatially-uniform, possibly anisotropic diffusivity tensor D and that diffusion encoding is achieved by a magnetic field gradient waveform g(t). Note that we absorb the gyromagnetic ratio γ into arising from a single such confined compartment under a general encoding waveform g(t) can be found thanks to the Brownian paths having a Gaussian probability measure under the potential (2) [17]: with the generalized encoding wave vector Here, t e is the duration of the encoding protocol, and is a matrix of equilibration rates. As depicted in Fig. 2, double-diffusion encoding is achieved by a gradient waveform consisting of two pairs of bipolar rectangular pulses, each with a given duration δ and separation ∆, and magnitudes g 1 and −g 2 ; the minus sign is customary. The time between the leading edges of the second pulse of the first pair and the first pulse of the second pair, t m , is called the mixing time. The calculation of the signal (3) therefore entails very simple integrals, but in a cumbersome piecewise fashion. Upon significant simplification one finds One may refer to the tensors 1In units of the thermal energy scale k B T , where k B is the Boltzmann constant and T is the absolute temperature.
2Vanishing confinement (i.e., free diffusion) has an unnormalizable steady state, but it can be handled. 3Alternatively, g(t) can be called a precession rate gradient waveform.
respectively, as the self-coupling and cross-coupling tensors between encoding blocks.4 For the orientational average below, we consider |g 1 | = |g 2 | = g, with a fixed angle ψ between them. As the free diffusion limit is approached, one can see that the cross-coupling of encoding blocks vanishes, as T × ∼ (D/2)Ωδ 2 ∆ 2 (1 + Ωδ − Ωt m ), due to displacements in separate time intervals being uncorrelated in pure Brownian motion [9]. The dependence on the mixing time t m does not enter until second order in confinement, in line with approximate calculations done for a spherical wall [12].

III. DOUBLE DIFFUSION ENCODING: POWDER AVERAGE
The orientational (powder) average of the compartment signal (6) is performed as follows: The vectors g 1 and g 2 form a plane with unit normal vectorn One constructs the average by integrating over all orientations of (the plane normal to)n, and for each of these integrate over all orientations of the pair {g 1 , g 2 } within the plane, with their relative angle ψ fixed. The procedure hence described can be written as5S Here, gn(β) is a vector of magnitude g in the plane normal ton, whose orientation is parameterized by the (in-plane) azimuthal angle β, according to which we can identify g 1 → gn(β) and g 2 → gn(β + ψ).
We first take the in-plane integralS which serves as the definitions for the intermediate quantitiesS(n) and σn(β). According to Eq. (6), the latter is given by To anchor the angle coordinate β, we define the in-plane cartesian coordinates u, v such thatû gn(0) andv =n ×û, yielding gn(β) =ûg cos β +vg sin β .
After an exercise in trigonometric simplification, and denotingû T •v = T • uv etc., one finds σn(β) = ςn + ρn cos 2β − n sin 2β, with ςn ρn These depend onn throughû andv, of course. Hence via Eq. (10), where an integral representation of the modified Bessel function of order 0 was recognized [20]. The argument of the square root can be found after a semi-tedious calculation as and we rewrite the powder-averaged double-encoding signal (9) via Eq. (10) and (15) as For the actual evaluation of the integral, explicit expressions in terms of the polar and azimuthal angles (θ, ϕ) ofn need to be substituted,6 which are steps we omit here. Eyeballing how entries of T • and T × appear in Eq. (17) via Eq. (14a) and (16), one notes that it is useful to define the intermediate tensors Referring to their eigenvalues, in any preferred order (but the same for both), as m i andm i , and using the following shorthand the powder-averaged signal (17) attains the form upon substantial algebraic manipulation.
Here, i F j (. . .) are the confluent hypergeometric functions [22], tilde denoting regularization. Note that these expressions apply to the most general case, in which the confinement tensor (and thus the tensors that are functions of it) have three distinct eigenvalues. Special cases associated with coinciding eigenvalues are discussed in the appropriate occasion later. Which alternative among Eq. (22) yields better convergence is a matter of how the eigenvaluesm i are chosen to be ordered, by way of the sizes and signs ofm − i j . As a general guideline it would be wise to order the eigenvalues so as to avoid a sequence that alternates in sign, and with a large expansion parameter. Take Eq. (22b) for instance. Given that 1 F 1 (k+l+1; k+l+ 3 2 ;m − 13 ) > 0 and increasing for allm − 13 , it would be beneficial to makem − 13 negative (and large if possible), while keepingm − 21 positive (and small if possible). For a given set of eigenvaluesm i , ordering them in the fashionm 1 <m 2 <m 3 would be along this guideline, whereas the orderingm 2 <m 3 <m 1 would result in a sequence with larger terms and alternating sign.
Note however that while T • is a monotonic (decreasing) function of Ω, T × is not; see Eq. (7). Through Eq. (18) this means that their mixtures M andM are not necessarily monotonic in the confinement Ω. Hence what ordering of the confinement eigenvalues Ω i achieves what ordering in the eigenvalues m i andm i is a question which has an answer only on a case-by-case basis. Furthermore, the ordering of Ω i that yields a desirable ordering ofm i for a particular one of Eq. (22) may not produce an ordering of m i as desirable for the convergence of q mk in Eq. (21c).

A. Axisymmetric confinement
We refer to the condition when two of the eigenvalues Ω i of the confinement tensor coincide as axisymmetric confinement. Under this condition, the series expansions above undergo simplifications.
The most drastic simplification occurs when Ω 1 = Ω 2 . That is, given that two eigenvalues coincide, assigning first and second place to them is the wisest choice as far as the evaluation of the series expansions (21) is concerned.
First, the coefficient q mk simplifies as7 The coefficient Y mk , on the other hand, loses the interior summation in Eq. (21b) due to m − 12 vanishing. The remaining summation can be identified according to the definition of hypergeometric functions [22] as This expansion is not so sensitive to the ordering of eigenvalues number 1 and 3, as the expansion parameter (m − 13 ) is squared, and the hypergeometric function 1 F 1 (2n+1; 2n+3/2;m − 13 ) > 0 for all arguments. We depict in Fig. 3 the evaluation of the powder averaged signal (25) for axisymmetric confinement for representative values of encoding parameters. In Fig. 3(a), the confinement is anisotropic (prolate), whereas in Fig. 3(b), it is (nearly) isotropic. The bell-shaped dependence on the relative angle ψ between the gradient directions is seen in both cases, which is a sign that diffusion is not free [12]. This dependence is due mainly to the exponential prefactor in Eq. (25) that has nothing to do with the difference between the confinement eigenvalues (anisotropy). When the mixing time is increased, the bell-shaped modulation, indicating confinement regardless of anisotropy, stops overwhelming the relatively smaller influence of the rest of the expression (25): see Fig. 3(a) where the confinement is anisotropic. In an isotropic confinement, on the other hand, angular modulation simply disappears when the mixing time is increased (Fig. 3b), illustrating that the angular modulation that survives the increase in mixing time is due only to compartmental anisotropy, and not due to the fact of confinement (or to ensemble anisotropy, which was already eliminated by powder averaging).
7For this, it needs to be noted [21, supplementary information] that the way q mk arises in the calculation-before ever arriving at Eq. (21c)-is that it is the coefficient in the (double) series expansion of the Bessel function in Eq.

B. Insights from two dimensions
The signal expressions for double encoding are a bit unwieldy to get a conceptual handle on. However, some insight can be gleaned from considering the orientational averaging in two dimensions.
Obviously, in the spirit of Eq. (9), the 2D orientationally averaged signal can be written as That is, there is no normal vector to integrate over; everything takes place in an (x, y) plane. The same steps Eq. (10) through Eq. (15) apply, and one hasS The dependence on the relative angle ψ of the gradient directions occurs both in the exponential attenuation factor and in the Bessel function. The angular dependence in the Bessel function has to do with anisotropy (T • 1 − T • 2 ) which the exponential factor is insensitive to due to the trace. In the exponent, on the other hand, the angular dependence is controlled by Tr T × , whose presence is due to confinement (since T × → 0 in the free diffusion limit). For large mixing times (t m Ω i 1) the latter drops out and the angular modulation due to anisotropy is liberated. However, for smaller mixing time, it turns out that the angular dependence of the exponential factor dominates, which is due to confinement, suppressing the signature of anisotropy.
An interesting double-encoding scheme is its "symmetrized" version [23], which fixes ψ = π/2, but varies the magnitudes g 1 = g cos α and g 2 = g sin α as a function of some parameter α. In that scenario, one can easily calculate that showing that the modulation due to confinement in the exponential factor drops out, and one of pure anisotropy remains (with ). The result is a 'cleaner' version of the signal modulation wherein the confinement anisotropy is the only source of angular modulation characterized by the cos4α dependence. It should be remembered though that the confinement model is purely Gaussian, as such does not account for compartmental kurtosis. When truly restricted diffusion is considered, the compartment anisotropy and compartmental kurtosis both yield the same type of angular dependence compromising the interpretation of such angular dependence except when the compartments are isotropic [23].

IV. SINGLE DIFFUSION ENCODING
The powder-averaged single-encoding signal can be obtained by way of "hacking" its double-encoded counterpart. One notes that S c = e −g T • g is the single-encoding signal at the compartment level [17]. Then the features of the form (6) to get rid of are (i) the presence of a second vector g 2 unequal to the first, (ii) the presence of cross-coupling (T × ), and (iii) the double occurrence of self-coupling. The first is hacked away by setting ψ = 0. For the second, one simply sets T × = 0.8 The result is S c = e −2g T • g , suffering from the third problem above, which is fixed by replacing g → g/ √ 2. All this yields via Eq. (18) which is the substitution that converts the powder-average expressions for double-encoding into those of single-encoding.
With this condition applied, them − i j 's appearing in Eq. (21) turn into m − i j 's. However this alone does not produce drastic simplifications such as reducing summations. Rather than using Eq. (21), in fact, it is better to note that other results in Ref. [21] are quite suitable in the case of single-encoding. Since the compartment signal has the form S c = e − Tr(T • gg ) , with the matrix gg axisymmetric by dint of being rank-1, the results of Ref. [21] for axisymmetric diffusion or measurement tensor can be applied. Written in terms of the parameters of the present discussion, the relevant formulas of Ref. [21] indicate the following alternative expressionsS for fully anisotropic Ω (hence T • ); the axisymmetric case is taken up below. These summations are subject to the same guidelines that followed Eq. (22) except that the caveats about conflicting eigenvalue ordering do not apply here. The single-encoding expressions here only involve the self-coupling tensor T • which is a monotonically decreasing function of Ω; see Eq. (7). Therefore the ordering of the eigenvalues T • i and Ω i are certain to be in exactly the opposite sense of each other.
8Physically, this can be imagined as the limit t m → ∞. However, physics is not necessary. We simply have a set of expressions, e.g. Eq. (21), containing T • and T × via Eq. (18) via Eq. (19), and we want to remove instances of T × .

A. Axisymmetry and the power-laws for confined diffusion
The first alternative in Eq. (30) needs special care when T • 1 = T • 3 and the argument of the hypergeometric function 2 F 1 (. . .) diverges. Invoking a property of the hypergeometric function,9 and renaming T • 1 = T • ⊥ and T • 2 = T • , one finds10 Here, we have also explicitly denoted the dependence on the encoding protocol's timing parameters as per Eq. (7). Note that in the free diffusion limit (Ω → 0, see footnote 4), the known single-encoding powder-average signal [25][26][27] with q = gδ, is recovered. Two special limiting cases follow.

Stick:S ∝ g −1 scaling
When particles have negligible latitude to move in the transverse direction, the orientationally-averaged signal (31) assumes a special form. In terms of the confinement model, this corresponds to Ω ⊥ → ∞, which implies T • ⊥ → 0 via Eq. (7). The signal (31) then becomesS The large gradient regime has been important in identifying stick-like compartments via the q −1 scaling, which has been observed in white-matter areas of the brain and has been interpreted with the assumption of a free one-dimensional diffusion [28], which is adequate for channels of straight long channels of infinitesimal diameter. A notable exception is Ref. [29], which has incorporated the effects of finite size and curvature. Ref. [29] also pointed out that such a scaling is not the true asymptotic behavior of the signal; the latter is rather dictated by the Debye-Porod law for narrow pulses [30] and an even steeper attenuation is predicted when the pulses are wide. The above expression suggests that a similar decay is expected for the gradient magnitude rather than the q-value. The crucial difference is in the dependence on the timing parameters {δ, ∆} of the SDE sequence, see Figure 1. It would thus be interesting, e.g. in white-matter, to investigate whether the dependence on the timing parameters is more like 1/ ∆ − δ/3 (free difusion along the fiber) or 1/ T • (δ, ∆) (confined diffusion along the fiber), which would inform about the diffusion process along the axons. 9 lim x→0 x n 2 F 1 m+ 1 2 , −n; m+1; 10 Alternatively to using the results of Ref. [21] here, one may go back to the integral expression (20). First, we note that axisymmetry, with the choice Ω 1 = Ω 2 , makes the integrand independent of ϕ, yieldingS We are not aware of a closed-form evaluation of this integral. However, its special case (29) relevant here, making the arguments of the exponential and Bessel function match, has the result [24, 6.625

Pancake:S ∝ g −2 scaling
For completeness, we consider the opposite case where particles are able to spread in a plane, but not along the normal. Here, Ω → ∞, implying T • → 0 via Eq. (7). One then has11 Thus, the orientationally-averaged signal attenuates at a faster rate than in the case of sticks. Similarly though, the dependence on the timing parameters are different in free and confined diffusion scenarios.

V. DISCUSSION
We have provided, for the first time, explicit expressions for the orientationally-averaged SDE and DDE MR signal intensity for structures represented by confinement tensors [17]. The latter is the effective model of restricted diffusion when pulses are long enough for the diffusing particles to traverse distances larger than the pore size [18]. As such, our findings are relevant for a broad range of porous materials featuring isolated, small pores.
The counterpart of these results for compartments of free anisotropic diffusion were given in Ref. [21] for arbitrary encoding waveforms. For the time being, considering arbitrary waveforms for confined compartments seems extremely challenging. However, taking confinement into account is important, since the free diffusion model lacks features (such as the dependence on relative angle in double encoding) that a realistic signal will bear, see Figure 3. Such bell-shaped angular modulation was related to the radius of gyration of the pores [10] and has since been used to estimate the apparent size of pores [4,12,[31][32][33] via DDE measurements. The underlying reason was thought to be a restriction effect [5] that not only leads to an anisotropy of the diffusion process at a length scale smaller than the pore size [12,[34][35][36], but also makes the apparent diffusion coefficient depend on the diffusion time [37]. Thus, this effect is absent when free diffusion is thought to take place within individual pores, see Figure 3c. Our results indicate that the confinement tensor framework is capable of capturing such angular modulation, which makes it a suitable representation of diffusion within microdomains [17,18,38]. This is important in applications like q-space trajectory imaging [39] and diffusion tensor distribution imaging [40], which aim to characterize the structure of subdomains using general gradient waveforms.
Another angular modulation that is apparent in Figure 3 is w-shaped, which was pointed out by Mitra [10] for randomly distributed sticks. Later, it was proved that the dominant contribution to such angular modulation had the functional form cos 2ψ for fully restricted structures and cylinders of finite diameter [41]. Such modulation that manifests itself at twice the "angular frequency" [10,41,42] is present even when diffusion within the subdomains is envisioned to be free as long as it is anisotropic, see Figure 3c. Thus, such modulation is truly indicative of the anisotropy of the subdomains, be it free (Figure 3c), confined (Figure 3a&b), or truly restricted [35,41]. For more information on such anisotropy, the reader is referred to the review on this topic by Ianuş et al. [43] in this book series.
From a mathematical point of view, the expressions in Ref. [21] provided the Laplace transform of a tensor distribution, which includes rotated copies of a given diffusion tensor wherein all orientations are equally likely, thus extending the modeling approach that employs parametric diffusion tensor distributions [44][45][46][47] to a new type of tensor distribution. Evaluating the signal, in a similar fashion, for confinement rather than diffusion tensors can be regarded as the evaluation of a transform whose kernel is the compartmental signal given in (6) instead of the kernel of the matrix Laplace transform e −BD .

VI. CONCLUSION
We have given analytical expressions for the orientationally averaged diffusion MR signal originating from confined anisotropic compartments for two relatively simple encoding schemes. A number of observations related to signal modulation and power-law tails were made for such confined pores. These findings complement and extend the exact expressions for locally free diffusion provided in Ref. [21] to confined diffusion albeit for SDE and DDE measurements.