Conceptual parallels between stochastic geometry and diffusion-weighted MRI

Diffusion-weighted magnetic resonance imaging (MRI) is sensitive to ensemble-averaged molecular displacements, which provide valuable information on e.g. structural anisotropy in brain tissue. However, a concrete interpretation of diffusion-weighted MRI data in terms of physiological or structural parameters turns out to be extremely challenging. One of the main reasons for this is the multi-scale nature of the diffusionweighted signal, as it is sensitive to the microscopic motion of particles averaged over macroscopic volumes. In order to analyze the geometrical patterns that occur in (diffusion-weighted measurements of) biological tissue and many other structures, we may invoke tools from the field of stochastic geometry. Stochastic geometry describes statistical methods and models that apply to random geometrical patterns of which we may only know the distribution. Despite its many uses in geology, astronomy, telecommunications, etc., its application in diffusion-weighted MRI has so far remained limited. In this work we review some fundamental results in the field of diffusion-weighted MRI from a stochastic geometrical perspective, and discuss briefly for which other questions stochastic geometry may prove useful. The observations presented in this paper are partly inspired by the Workshop on Diffusion MRI and Stochastic Geometry held at Sandbjerg Estate (Denmark) in 2019, which aimed to foster communication and collaboration between the two fields of research.


Introduction
Diffusion-weighted magnetic resonance imaging (MRI) [23] is one of the few imaging modalities that is capable of mapping the immensely complex microstructural architecture of the human brain in a non-destructive manner. This is Supported by the Villum Foundation.
achieved by applying a specific sequence of diffusion-sensitizing magnetic field gradients during the MRI acquisition, producing a signal with a decay rate that is dependent on the relative mobility of water molecules in the tissue [30]. As the overall mobility of the molecules is decreased by the presence of any material barrier, the acquired signal effectively provides an indirect probe of the ambient structure. At the micrometer length scales accessible in current scanners, the dominant barriers to the diffusing molecules in the brain are the fiber-like neurites that transmit information between different regions of the brain [3], and whose properties are (naturally) of profound importance in neurology. Aside from the anticipated and present clinical value of this modality, diffusion-weighted MRI also stands to provide unique information about the evolution, morphogenesis, and function of the brain, already being pursued through the tracking of macroscopic fibers in tractography [2,4,20]. The challenge faced by the diffusion-weighted MRI community is to identify the relevant structural parameters determining the signal decay, and-to the extent that this is possible-to invert the relation between them. This is where we believe stochastic geometry could play a role.
Stochastic geometry [29] is an area of statistics that provides modeling and inference techniques for complex spatial objects whose structure can be described effectively as random patterns. One of the first cases studied in what is now called stochastic geometry-the problem of stereology [1]-included the inference of geometric properties of 3D objects from their intersections with a small number of 2D planes. This topic gained traction in the 1960s, as its solutions alleviated the challenging and computationally demanding task of actually computing 3D reconstructions. Stereology found a wealth of applications ranging from geology to for example microscopy of neuroanatomy. Since then stochastic geometry has evolved into a mature field of mathematics that offers a rich toolbox of rigorously developed techniques, including models with potential relevance for diffusionweighted MRI such as random tessellations and fiber processes. Skimming the stochastic geometry literature one quickly comes across a number of concepts that have obvious analogues independently developed in the diffusion-weighted MRI literature, and yet the two research domains have had very limited contact so far.
This chapter is written from the perspective of diffusion-weighted MRI as a field working on a set of challenging modeling problems, focusing on the potential utility of stochastic geometry in addressing them. Our aim is not to provide an exhaustive overview of relevant theory in either field, but simply to highlight some well-known results where the conceptual links between the two become apparent. A comprehensive introduction to stochastic geometry can be found in the book by Stoyan et al. [29], while the books by Jones [16] and by Johansen-Berg and Behrens [15] provide a good entryway to diffusion-weighted MRI.
To initiate an exchange of ideas between these two fields, the first Workshop on Diffusion MRI and Stochastic Geometry was co-organized by the authors and Eva B. Vedel Jensen at Sandbjerg Estate (Denmark, January 20-24, 2019). The observations presented in this chapter partly inspired this workshop, but the chapter likewise builds on the valuable discussions held at the lively and interactive workshop. We are not aware of any previous works pointing out the parallels described here, although some of them will undoubtedly have been noticed before.

Overview
The molecular dynamics relevant for diffusion-weighted MRI are conveniently described by a displacement probability density function-the diffusion propagator P (r, t). The propagator P (r, t) represents the probability of a displacement r at a diffusion time t, and the diffusion-weighted signal is related to the characteristic function S of the propagator given by where qr denotes the inner product between q and r. The Fourier parameter q can be considered an experimental parameter, determined in practice by the diffusion-sensitizing gradients of the acquisition. In the following sections we give examples of connections to stochastic geometry for three different 'limiting regimes' of the diffusion-weighted signal, where expressions for the relevant molecular dynamics can be simplified significantly. Secs. 2 and 3, which feature the short and long diffusion time limits respectively, can be well-understood in terms of the time-dependent diffusion coefficient D(t), which represents the mean squared molecular displacement at a time t. This diffusion coefficient appears as the first non-trivial coefficient in the Taylor series of log S, the cumulant expansion 3 From this expansion it can already be gathered that results based on the diffusion coefficient are mostly applicable when q is relatively small-i.e., when the gradients are weak enough for the D-dependent term in this expansion to dominateand in Sec. 4 we consider instead the limit where the gradient strength parameter ∼ q becomes large. In this regime the diffusion coefficient no longer provides an adequate vehicle for the description and analysis of the diffusion process, and we have to rely on other descriptors. Where necessary, additional details on these concepts will be given in the text, although technicalities will be skipped in favor of accessibility. To simplify the exposition further, we will make implicit use of the narrow pulse approximation [6,27] throughout. The paper concludes with an outlook on the possible future uses of stochastic geometry in Sec. 5.
2 Specific volumes and the short-time limit ". . . the time-dependent diffusion coefficient D(t) of mobile molecules confined in pores or cells carries information about the confining geometry. At early times, [a perturbative expansion of D(t)] gives, irrespective of details, the pore surface to volume ratio . . . , which is a measure of microscopic length." Sen (2004) The first regime we consider is the short-time limit, described in the seminal works by Mitra, Sen, and others [18,19,26], and with important insights dating back to Kac [17]. The situation analyzed in these works is the diffusion of molecules in the neighborhood of obstructive geometrical features, Fig. 1. At the boundaries between the diffusive medium and the geometry, the diffusion is prescribed by boundary conditions that can depend on e.g. the permeability, leaving the spatially averaged, time-dependent diffusion coefficient D(t) to be solved. While this is a very challenging problem for any non-zero, finite time t, the limiting behavior for t → 0 can be expressed in terms of a relatively small set of practically useful structural parameters. As D(t) can be measured in the scanner, we can use diffusion-weighted measurements to obtain estimates of these parameters. The emergence of these structural parameters in the diffusion-weighted signal can be understood as follows.
In the limit t → 0 the diffusing molecules do not have enough time to interact with the geometry, and the diffusion coefficient naturally approaches the medium's free diffusion coefficient D(0) = D 0 -the diffusion coefficient of the medium in the absence of any geometry. At times close to the limit, the particles move a typical distance in the order of √ D 0 t (by virtue of the definition of D in terms of the mean squared displacement), and so roughly speaking only the fraction of particles within some distance ∼ √ D 0 t can 'see' the geometry. As the geometry essentially acts as a barrier to the diffusion, the time-dependent diffusion coefficient decreases from its free diffusion limit, and this decrease is more significant if a larger fraction of the total number of particles can interact with the boundaries. For diffusion times approaching 0, this fraction becomes exactly proportional to the surface area of the geometry. Formalizing this notion, Mitra et al. [19] then showed that the first order correction describing the approach of the limit becomes proportional to the surface-to-volume ratio S/V of a smooth geometry, cf. Fig. 1, according to where d is the spatial dimension. From this relation, the surface-to-volume ratio can be estimated reliably from the diffusion-weighted measurements [13]. The order t terms in this expansion depend to varying degrees on the boundary condition parameters such as the permeability, as well as on the average curvature of the geometry, while higher order terms depend on even more intricate details of the environment. Although the complex interactions between particles undergoing random displacements and their surrounding structures have not been considered in stochastic geometry as such-this problem is closer to mathematical physics [5]-the surface-to-volume ratio uncovered in Eq. (3) is a commonly estimated quantity in for example stereology. The surface-to-volume ratio is also referred to as the surface density or as the specific surface area in stereology and stochastic geometry. An example of a question in this scenario for which stochastic geometry could be helpful is: "what is the smallest observation window (voxel) that has an acceptable error when estimating a given characteristic?" [29, Sec. 6.4.6].
3 Stationarity and the long-time limit "That is, we assume here that the voxel is statistically homogeneous. This macroscopic uniformity allows us to go from averaging over the contributions from all parts of the system . . . to ensemble averaging over all disorder realizations, leading to the description of the signal in terms of the statistical properties embodied in the correlators . . . "

Novikov and Kiselev (2010)
In the short-time limit of the previous section, particles could not explore their surroundings outside a vanishingly small window, leading to the simple relation between the surface-to-volume ratio and the time-dependent diffusion coefficient. In the long-time limit, on the other hand, we can assume that particles in fact see all the surrounding structures. Under fairly general conditions, the diffusion process in this limit can still be described in terms of a single diffusion coefficient-the limiting coefficient D ∞ = lim t→∞ D(t). While determining D ∞ for a given sample is very difficult, Novikov et al. [21,22,24] showed that once again there are basic geometrical properties of the structure that completely determine the limiting behavior of D(t). In this case, the diffusion coefficient approaches its limit as where c is a constant, and where the exponent ϑ = (p + d)/2 is determined by the spatial dimension d and a structural exponent p. The structural exponent p is the exponent of the limiting power law behavior of the structure correlation function, which roughly speaking quantifies the presence of long-range correlations in the diffusion barriers. While the constant c varies significantly between different configurations of barriers, we know of only a few possible values for the exponent p, with many vastly different configurations of barriers producing the same critical exponent in experiments. The structural exponent p effectively distinguishes between different structural universality classes [21]. A central concept in the cited works by Novikov et al. is that because the diffusing particles can explore their entire surroundings, the exact configuration of the barriers becomes irrelevant. Consequently, one can theoretically replace the complicated original structure with a simpler-to-analyze effective medium, permitting the derivation of e.g. Eq. (4). As explained in the quote above, the practical application of this simplifying methodology requires that we assume that the subsets of a sample explored by different particles can be viewed as random samples ('disorder realizations') generated by some basic statistical properties of the global structure. In the context of stochastic geometry, this same assumption is more commonly called (spatial) stationarity, and it is a key assumption in many classical proofs in the field. In particular, the assumption of stationarity enabled the first proofs of the fundamental formulas of stereology. In a number of stochastic geometry results it is now known that a weaker 'first order stationarity' assumption is sufficient [1], and it might be interesting to see if the same is possible for the long-time limit results discussed here.
4 Directional measures and the strong-gradient limit "The concept of disordered media and statistical averaging can be particularly valuable to deal with the geometric complexity of biological tissues. We believe that further progress in the field can be achieved by merging microscopic geometric models, statistical [descriptions] of disordered media, and high-gradient features of the signal formation." Grebenkov (2016) The gradient strength parameter set during a diffusion-weighted MRI experiment determines-roughly speaking-the scanner's sensitivity to diffusive motion. The mobility of the water molecules in the sample affects the decay rate of the diffusion-weighted signal, cf. Sec. 1, and stronger gradients make this effect correspondingly stronger. Although stronger gradients thus produce weaker, harder to detect signals, these signals contain information about interesting features of the diffusion that cannot be observed at lower gradient strengths. The results for the short-time and long-time limits described in the previous subsections are mainly used at low to moderate gradient strength acquisitions.
When we increase the gradient strength in a diffusion-weighted MRI experiment we first notice that the anisotropy, i.e., the orientation-dependence, becomes more significant. While we omitted this before, anisotropy is already a factor at the lower gradient strength experiments used for the concepts discussed in Secs. 2 and 3. In the general case, the long-time limit D ∞ of D(t) considered in Eq. (3) is, for example, in fact a tensorial quantity related to the mean squared molecular displacements along different orientations. This anisotropy in the diffusion reflects the anisotropy in the sample's micro-structure, so for example in the brain it is predominantly determined by the orientation distribution of neurites, cf. Sec. 1. At higher gradient strengths, the orientation-dependence of the diffusion can no longer be represented by a simple tensor, and a more complete set of orientational features becomes accessible.
The neurite orientations are generally characterized by a so-called fiber orientation distribution function (ODF), which specifies the likelihood that a neurite is locally tangent to a given orientation, and the estimation of this object from strong-gradient experiments is a common problem considered in diffusionweighted MRI. A large number of techniques have been proposed to deal with this question [14,31,32], but we will not discuss them in detail here. Instead we will point out that the ODF also has a stochastic geometrical analogue: the rose of directions [29]. The rose of directions is defined for fiber processes-a well-defined (locally finite) collection of randomly placed fibers. The rose of directions is the distribution of the direction tangent to a typical point on a fiber, where the meaning of a 'typical point' can be made explicit using Palm distribution theory. A stationary fiber process is an example of a stochastic geometrical model in the spirit of the work by Novikov et al. [22] that could be useful to model biological tissues. Furthermore, the rose of directions has a dual-the rose of intersections. These two objects are related through the Funk-Radon integral transform [9], which also appears naturally in for example the diffusion-weighted signal expression for narrow cylinders (the 'fiber ball' ODF) [14].
At even stronger gradients we move to another limiting regime-the stronggradient limit, or localization regime [11]. Here we can still make use of the intuition developed in Sec. 2: any structure in the sample will act as a barrier for particles in its vicinity, and thus slow down the average motion in that region. Slower diffusion in turn implies less signal decay, and so the main contributions to the signal now come from particles localized near the barriers [28]. The localization effect carries non-trivial consequences for the diffusion time-dependence of the signal, but for this chapter we are more interested in the orientation dependence in this regime.
A simple way to understand the impact of the geometry's anisotropy on the signal in the strong-gradient limit, is to look at a diffusion propagator with a compact support Ω. At strong gradients the integrand in Eq. (1) oscillates rapidly, and as a result the integral indeed tends to 0. However, the rate at which it decays is in fact governed by the support Ω, which is of course determined by the geometry. We can consider as an illustration the 'indicator' propagator which is 1 everywhere inside a sufficiently smooth and convex Ω, and 0 elsewhere. The signal decay is then exponential with an exponent sup r∈Ω qr [12], where the exponent, viewed as a function in q, determines Ω completely. We could for example use this observation to recover the shape of a convex pore enveloping a diffusing substance, provided that the pore is much smaller than the typical diffusion length √ D 0 t [6], but this approach is not limited to closed pores. The function q → sup r∈Ω qr is called the support function of Ω in stochastic geometry, where it can be associated with the rose of directions of a fiber process [29]. A similar interpretation can be given to the support function as it occurs here, leading to the definition of the barrier ODF [7,8]. It must be mentioned that as these developments are very recent, their practical utility is currently still being investigated.

Perspectives
We have given three examples of concepts in diffusion-weighted MRI that also occur in stochastic geometry: the surface-to-volume ratio, which is related to the specific surface area commonly estimated in stereology; statistical homogeneity, which is related to stationarity; and the fiber orientation distribution function, whose stochastic geometry analogue is the rose of directions. In the context of neuroimaging, the common thread between them appears to be the concept of stationary fiber and surface processes from stochastic geometry, which are completely characterized by the specific length/area and the rose of directions.
We believe further investigations in this direction could lead to new and useful methods for the analysis of diffusion-weighted MRI data, as well as novel research problems in stochastic geometry. In particular, we hope that the link to the mature statistical theory found in stochastic geometry may offer practical tools for the analysis of uncertainty and variation in applications of diffusion weighted MRI: How much faith can we put in estimated anatomical quantities? Which quantities can we expect to be able to derive from given data? What are the conditions we should put on our data acquisition in order to be able to draw sound conclusions for the hypotheses that we ultimately hope to investigate?