Quantifying brain microstructure with diffusion MRI: Theory and parameter estimation
Dmitry S. Novikov, Els Fieremans, Sune N. Jespersen, Valerij G. Kiselev
QQuantifying brain microstructure with diffusion MRI: Theory and parameter estimation Dmitry S. Novikov, ∗ Els Fieremans, † Sune N. Jespersen, ‡ and Valerij G. Kiselev § Center for Biomedical Imaging, Department of Radiology, NYU School of Medicine, New York, NY, USA CFIN/MINDLab, Department of Clinical Medicine and Department of Physics and Astronomy, Aarhus University, Aarhus, Denmark Medical Physics, Deptartment of Radiology, Faculty of Medicine, University of Freiburg, Germany (Dated: April 12, 2018) We review, systematize and discuss models of diffusion in neuronal tissue, by putting them into an overarching physical context of coarse-graining over an increasing diffusion length scale. From this perspective, we view research on quantifying brain microstructure as occurring along the three major avenues. The first avenue focusses on the transient, or time-dependent, effects in diffusion. These effects signify the gradual coarse- graining of tissue structure, which occurs qualitatively differently in different brain tissue compartments. We show that studying the transient effects has the potential to quantify the relevant length scales for neuronal tissue, such as the packing correlation length for neuronal fibers, the degree of neuronal beading, and compartment sizes. The second avenue corresponds to the long-time limit, when the observed signal can be approximated as a sum of multiple non-exchanging anisotropic Gaussian components. Here the challenge lies in parameter estimation and in resolving its hidden degeneracies. The third avenue employs multiple diffusion encoding techniques, able to access information not contained in the conventional diffusion propagator. We conclude with our outlook on the future directions which can open exciting possibilities for designing quantitative markers of tissue physiology and pathology, based on methods of studying mesoscopic transport in disordered systems. CONTENTS
1. Diffusion MRI through a bird’s eye 2 theory 2 qt Imaging 5 coarse-graining: The three regimes 7
2. Time-dependent diffusion in neuronal tissue 12 effect? 13 (i) : Net surface area of restrictions 16 (ii) : Structural correlations via gradual coarse-graining 16 ∗ [email protected] † els.fi[email protected] ‡ sune@cfin.au.dk § [email protected]
3. The t → ∞ limit, regime (iii) : Multiple Gaussian compartments 21 Flat directions in the fitting landscape 26 Imaging (NODDI) 28 l = 0 invariant K ( b ) K l ( b ) for l = 2 , , . . .
4. Multiple diffusion encodings 31 O ( q ) O ( q ) and beyond. Microscopic anisotropy 33
5. Outlook and open questions 35 Acknowledgments 37 References 37 A. Causality and analytical properties in the frequency domain 47 B. OG with a finite number of pulses 48 C. Probing the
S/V limit with finite- N OG 50 a r X i v : . [ phy s i c s . b i o - ph ] A p r
1. DIFFUSION MRI THROUGH A BIRD’S EYE One of the most astonishing things about the world inwhich we live is that there seems to be interestingphysics at all scales. (cid:104) ... (cid:105)
To do physics amid thisremarkable richness, it is convenient to be able to isolatea set of phenomena from all the rest, so that we candescribe it without having to understand everything.Fortunately, this is often possible. We can divide theparameter space of the world into different regions, ineach of which there is a different appropriate descriptionof the important physics. Such an appropriatedescription of the important physics is an “effectivetheory” H. Georgi,
Effective Field Theory [1] Diffusion MRI (dMRI) is a macroscopic physical measure- ment of the voxel-averaged stochasticmotion of nuclear-spin- carrying molecules (typically water). This measurement oc- curs in a structurally complex tissue microenvironment such as the brain. Diffusion in complex media has been studied for about a century in a variety of fields, and is part of a broad class of transport phenomena in disordered systems. Our goal in this review article is to place biophysical dMRI modeling into a broader physical context. Our overarching theme will be that of coarse-graining and effective theory, which will allow us to present and discuss neuronal tissue models of diffusion from a unifying perspective. One of the key 20 th century advances in understanding the physics of complex systems was achieved by the develop- ment of effective theory , a paradigm to describe dynamics that involves only a handful of the so-called relevant degrees of freedom , or relevant parameters , thereby ignoring myriads of other, “irrelevant” ones [2–4]. This way of thinking was spurred by attempts to describe systems with an ever greater number of degrees of freedom, and a subsequent realization that it is plain impossible to keep track of all of them at once. The more complex the system, the more the challenge of building an adequate theory shifts towards identifying which (few) parameters to keep, and which ones (almost all!) to ignore. Over time, selecting relevant parameters and formu- lating an adequate effective theory has become synonymous with the notion of understanding the system’s behavior. Having NMR as an example, the quantum-mechanical cou- plings of a very complex multi-spin Hamiltonian, together with all molecular degrees of freedom describing rotations, vibrations and translations, relevant at the nm and ps level, average out to produce effective parameters such as the re- laxation rates R and R , and the diffusion coefficient D , at least for the most common NMR measurements. The param- eters R and R emerge in Bloembergen–Purcell–Pound the- ory and enter the Bloch equations describing the semiclassi- cal evolution of macroscopic magnetization [5, 6]. Reducing a myriad of variables describing molecular microenvironment to just a few relevant parameters has been a major scientific achievement of the 1940s–1960s NMR, and has formed the basis of effective theory of nuclear magnetization in liquids. The step from NMR in uniform liquids to biological tissues has brought along a new challenge, which our community is only beginning to fully embrace. This challenge is associated with the above effective parameters R ( r ) , R ( r ) and D ( r ) acquiring spatial dependence at the scale ∼ . − µ m, set by the cellular architecture, much coarser than molecular dimensions. These spatial variations become relevant at the corresponding ∼ − ms time scales of dMRI, — much slower than the ps time scales on which the local relaxation rates and the diffusion coefficient emerge. From the physics standpoint, the spatial variations of R ( r ) , R ( r ) and D ( r ) (with the latter including bound- ary conditions associated with cell membranes), occur at the mesoscopic scale , Fig. 1. The term “mesoscopic” originated in condensed matter physics some decades ago [7], signifying focussing on the intermediate scales (“meso”), in-between el- ementary (say, atomic or molecular), and macroscopic (asso- ciated with the sample size or the measurement resolution). By design, this term is relative, depending on which spatio- temporal scales are deemed small and large. For dMRI, the mesoscopic scale corresponds to tissue het- erogeneities at the scale defined by the MRI-controlled dif- fusion length, L ( t ) ∼ √ Dt ∼ − µ m, which is the root mean squared molecular displacement for times t ∼ − ms. In the dMRI literature, it is commonly referred to as the microstructure scale. This scale is commensurate with immense structural complexity of tissue architecture. At the mesoscale, quantum degrees of freedom become ir- relevant (at least for the dMRI purposes), and the dynamics of transverse magnetization m ( t, r ) (a two dimensional vec- tor, represented by a complex number) can be captured by the mesoscopic Bloch-Torrey equation ∂ t m ( t, r ) = ∂ r [ D ( r ) ∂ r m ( t, r )] − [ R ( r ) + i Ω( t, r )] m ( t, r ) . (1.1) Here Ω( t, r ) is the Larmor frequency offset that may include externally applied diffusion-sensitizing Larmor frequency gradients g ( t ) , Ω( t, r ) = Ω( r ) + g ( t ) r , and the static Ω( r ) arises from the intrinsic mesoscopic magnetic structure of tis- sues due to paramagnetic ions such as iron, myelin suscep- tibility in the white matter, or due to added contrast agent. While we focus on the transverse magnetization in what fol- lows, the full version of the above equation includes the lon- gitudinal magnetization components with m being a three- molecules, mesocells, macroresolution,elementary FIG. 1. The mesoscopic scale in brain dMRI, as an intermediatescale between the elementary (molecular) and the macroscopic (res-olution). dimensional vector. Further extension can incorporate multi- component m to describe the interplay between different pro- ton pools, e.g., to describe magnetization transfer [8, 9]. The mesoscopic Bloch-Torrey equation (1.1) is an ade- quate effective description at the µ m level, commensurate with typical diffusion length scales probed with dMRI. It is a mesoscopic equation in the sense that it involves scales in-between the quantum-mechanical molecular dynamics on the nm scale and the measurable signal in mm-sized MRI voxels. While the averaging up to the mesoscopic scale is already performed in its r -dependent parameters, it is our task to perform the remaining averaging over a macroscopic voxel V inherent to the observed (complex-valued) signal S [ t, g ( t )] ∝ (cid:82) V d r m ( t, r ) , for which the µ m-level spatially varying relaxation rates and diffusive properties produce the observable deviations from mono-exponential relaxation and Gaussian diffusion. It is because of this averaging that ad- dressing the mesoscopic tissue complexity requires bringing the tools and intuition from condensed matter and statistical physics, in contrast to the quantum-mechanical description at the molecular level [5, 6] and classical electrodynamics-based considerations used in designing MR hardware. The overarching goals of “microstructural”, or “meso- scopic” MRI modeling are (i). To identify the relevant tissue-specific parameters, which contribute to R ( r ) , R ( r ) , D ( r ) , Ω( r ) , and sur- vive in the voxel-averaged signal (i.e., to build an ap- propriate effective theory for the macroscopic signal); (ii). To suggest optimal ways to probe them (i.e., to solve the corresponding parameter estimation problem). Notice that to keep our terminology reasonably rigorous, we separated modeling into theory and parameter estimation (sometimes called “fitting”); hence our title. These two facets of modeling require very different tools and ways of thinking, as we will see below in Sections 2 and 3, respectively. Equation (1.1) is an example of an effective theory — i.e., an approximate description that emerges by averaging out the dynamics at the smaller spatial and temporal scales. It illus- trates a general principle: pretty much every dynamical equa- tion in physics is an effective theory (governed by an effective Hamiltonian or an effective action), i.e., it has emerged by identifying “collective” phenomena involving many-particle interactions at a more elementary level [1–4]. Over the past century, physicists have come to realize that, at each level of complexity, the effective theory and its rel- evant parameters can look very different [2], giving rise to the hierarchy of scales and of the corresponding emergent phenomena, from the most microscopic to the most macro- scopic. Interactions between quarks and gluons give rise to protons and neutrons, so that their charge and mass can be viewed as effective parameters emerging by averaging over the quark/gluon degrees of freedom. Interactions between protons and neutrons forms a nucleus; interactions between nuclei and electrons give rise to all of chemistry, whereby the details of interactions between protons and neutrons inside nuclei become irrelevant. Interactions between molecules, coarse-grained over nm scale, give rise to hydrodynamics, statistical mechanics and eventually, to biology, and so on. It is remarkable that, for instance, there is not a hint of clas- sical hydrodynamics at the level of the Schr¨odinger and Dirac equations describing the atomic structure; the large-scale hy- drodynamic description emerges after a highly nontrivial av- eraging over the corresponding quantum degrees of freedom of many molecules. Refined methods such as renormaliza- tion [3, 4] were crafted specifically to single out relevant pa- rameters from the rest upon iterative coarse-graining [10, 11], which is a procedure that averages the dynamics over finer- scale degrees of freedom to derive approximate effective dy- namics at the coarser scales involving a minimal number of parameters. This way of thinking reveals a fascinating hierar- chy of natural phenomena [2]. For quantifying tissue microstructure by measuring dif- fusion, transverse relaxation or magnetization transfer with MRI, the mesoscopic Bloch-Torrey equation (1.1) contains all relevant physical processes. This effective theory is as funda- mental for the mesoscopic MRI, as the Schr¨odinger equation is for the non-relativistic quantum mechanics, or the Navier- Stokes equation is for the classical hydrodynamics. It is al- ways the starting point for developing biophysical models re- lating the NMR signal to the mesoscopic tissue architecture. Diffusion in heterogeneous media is a beautiful and simple example of coarse-graining. It can be thought of as a gradual “forgetting”, or homogenizing over the increasing diffusion length. To illustrate this concept, consider a two-dimensional model example of a two-scale mesoscopic structure, repre- sented by randomly placed impermeable disks of two dif- ferent radii, embedded into an NMR-visible space with dif- fusion coefficient D , Fig. 2. To be specific, let us assign sizes, typical to cell dimensions: the small disks have ra- dius r small = 1 µ m and the large ones are 20 times larger. In a (hypothetical) tissue, this could describe diffusion in the extra-axonal space transverse to a fiber tract, hindered by two types of axons. Here we consider diffusion as a physical phe- nomenon; its relation to dMRI is discussed below in Sec. 1.4. At time t → , each water molecule only senses its own immediate environment; all molecules see the same “intrin- sic” diffusion coefficient D | t =0 = D , which is of the order ∼ µ m / ms . (For pure water at ◦ C, D = 3 µ m / ms .) As time increases (top row of Fig. 2), molecules get re- stricted by the walls of both small and large disks. As small disks have much higher net surface area than the large ones, the hindrance occurs mostly due to the small ones. Hence, the decrease of the resulting voxel-averaged diffusion coefficient would happen on the scale of a few ms, mostly dominated by the geometry of the small disks at the scale ∼ µ m. At t (cid:38) ms (bottom row), when the diffusion length FIG. 2.
Diffusion as coarse-graining.
An example of a medium where the mesoscopic structure is created by randomly placing black disksof two different radii, r small = 1 µ m and r large = 20 µ m, top left panel. To obtain snapshots of the medium as effectively seen by the diffusingmolecules at different time scales, we used a Gaussian filter with width L/ , where L ( t ) = √ Dt , and ignored the time dependence of D ( t ) in the definition of diffusion length, using a typical value D = 1 µ m /ms for the illustration purposes (cf. Sec. 2.4 below). L ( t ) strongly exceeds the small disk size, the effect of the small disks has become coarse-grained (while the effect of the large disks is not). Now, we can view the medium in-between the large disks as a homogeneous “effective medium”, with some effective diffusion coefficient D small < D given by the macroscopic (“tortuosity”) limit of a medium with the small disks only. It is important to note that if we did not have access to shorter times and could only resolve the dif- fusion times corresponding to the lower row of panels, there would be no way to identify the presence of the small disks — their effect has been homogenized, and their numerous parameters (e.g., size, coordinates) have become “irrelevant”, with their only role in renormalizing D down to D small . Hence, from time t (cid:38) ms on, we can adopt the coarse- grained description which only involves the large disks, im- mersed in a uniform medium with diffusion constant D small . The corresponding Eq. (1.1) would have the effective D ( r ) varying at the scale associated with large disks, with D ( r ) (cid:39) D small outside them, and the short-distance spatial harmonics of D ( r ) filtered out as it is obvious from Fig. 2; in Section 2, we will rigorously justify and use this intuitive picture. The measured diffusion coefficient would further decrease with t at the scale of a few hundred ms, corresponding to being hin- dered by the large disks — the remaining restrictions. Eventually, at even longer t (cid:38) ms, the effect of the large disks also becomes coarse-grained, and the whole sam- ple looks as if it were homogeneous with some macroscopic An experienced reader can recall the possibility to apply strong gradients q ∼ /r small to detect the small disks. This, however, practically requiressensitivity to short t ∼ / ( D small q ) , cf. Sec. 1.6 below. diffusion coefficient D ∞ , such that < D ∞ < D small < D . From this t onwards, one cannot distinguish this sample from a uniform medium with diffusion constant D ∞ . Our example shows that the hallmark of coarse-graining over larger- and larger-scale mesoscopic structure is the time- dependence of the overall diffusion coefficient. In the view of this time dependence, it will be convenient to work with the instantaneous diffusion coefficient, D inst ( t ) ≡ ∂∂t (cid:104) δx ( t ) (cid:105) , δx ( t ) = x ( t ) − x (0) , (1.2) defined as the rate of change of the mean squared molecular displacement (cid:104) δx ( t ) (cid:105) in a particular direction ˆx . (For sim- plicity, we assumed that our sample in Fig. 2 is statistically isotropic. For anisotropic samples, the diffusion tensor com- ponents will acquire the time dependence.) The average (cid:104) . . . (cid:105) in the definition (1.2) is actually a dou- ble average : (i) over the Brownian paths in the vicinity (of size ∼ L ( t ) ) of any given initial point r ≡ r ( t ) | t =0 , yielding the local coarse-grained diffusivity value D ( r ) ; and (ii) over the ensemble of random walkers (spins) originating from all possible initial points r . Because of the ensemble average, the measured diffusion characteristics, such as Eq. (1.2), de- scribe a macroscopic sample as a whole . They do not belong to any given Brownian path, but rather emerge as a result of averaging (i) over all possible Brownian paths that could be taken by a given molecule, and (ii) over the initial positions of all molecules in a sample. Upon taking into account increasing length scales, the ef- fective voxel-averaged D inst ( t ) flows towards the tortuosity limit D ∞ , starting, as in our example, from the “microscopic” D > D ∞ . We used the term “to flow” because the above picture mimics the renormalization group flow [3, 4] accord- ing to which the gradual evolution of a physically important parameter, such an elementary particle charge or an effective mass, occurs as a function of the coarse-graining length scale, from the high-energy = short-distance scale, down to the low- energy parameters relevant for the macroscopic description. Looking back, there was nothing special about requir- ing the disks to be impermeable (the black regions could have corresponded to some medium with diffusion coefficient D (cid:54) = D ); we could have used objects of a non-disk shape, and/or with non-sharp boundaries. Generally, as long as the random walkers can, in the limit t → ∞ , reach any point in a given “compartment”, the above coarse-graining picture applies to this compartment. If a voxel contains multiple non- exchanging compartments, it applies to them separately, with the net signal given by a sum of their contributions. A similar physical picture qualitatively applies to the ef- fects of spatially varying transverse relaxation rate R ( r ) — e.g., if the black and white regions in Fig. 2 instead repre- sented different local molecular-level R values, and spins were able to diffuse everywhere. The above argument would then lead to an effective R ( r ) entering Eq. (1.1) for times t exceeding the corresponding coarse-graining time scale. For instance, for t (cid:38) ms, the effect of small disks would homogenize to produce a uniform R , small rate in-between the large disks, and so on, leading to the time-dependent overall observed R ( t ) with an asymptotic macroscopic rate R | t →∞ observable at very long t . Likewise, if the meso- scopic structure in Fig. 2 represented spatially varying sus- ceptibility χ ( r ) , inducing the corresponding Ω( r ) , the result- ing R ∗ ( t ) rate would change — in this case, increase with t , approaching the R ∗ ( t ) | t →∞ macroscopic value as a result of gradual coarse-graining. We note that all the above mentioned quantities — D small and R , small ; D ∞ and R | t →∞ ; R ∗ ( t ) | t →∞ — are nonuni- versal , i.e., they depend on the numerous structural details, such as packing geometry (e.g., periodic versus random ar- rangement); they would change if the disks were instead squares, etc. Certainly, these quantities are not given by a simple averaging of the microscopic D ( r ) or R ( r ) over the sample. However, the initial values D inst (0) and R (0) are given by the sample-averaged (cid:104) D ( r ) (cid:105) and (cid:104) R ( r ) (cid:105) , corre- spondingly, since at t → (practically, at times just exceed- ing the ps time scale necessary for the local D ( r ) and R ( r ) to emerge), each spin senses only its immediate environment. The picture of gradual coarse-graining over an increasing diffusion length has a number of important consequences: We also note that the renormalization group flow can have fixed points ,i.e., microscopic parameter sets for which the effective parameters do notchange with the increasing coarse-graining scale. Approaching the asymp-totically normal diffusion with a finite D ∞ is a Gaussian fixed point,which is what happens for most structural arrangements; one can say thatdiffusion is most tissues is a continuous family of Gaussian fixed points(for each realization of microscopic tissue architecture). An example of anon-Gaussian fixed point is the so-called anomalous diffusion, cf. Sec. 1.9. The mesoscopic Bloch-Torrey equation (1.1) can be fully determined only after the relevant spatio-temporal scales are specified, since its parameters R ( r ) , Ω( r ) and D ( r ) are effective and, hence, scale-dependent. Generally, the observed voxel-averaged diffusion coef- ficient and the effective relaxation rate would depend on time t because of the presence of the mesoscopic structure (such as D inst ( t ) decaying from D | t =0 down to D ∞ in our example of Fig. 2). This time t can be set by the measurement sequence, and varying it provides a unique window into the tissue architecture at the scale of the corresponding diffusion length L ( t ) . This brings us to the fundamental challenge of inter- preting such time dependencies in terms of the meso- scopic structural complexity. Practically, we must fig- ure out which features in the effective R ( r ) , Ω( r ) and D ( r ) remain observable after the voxel-wise averaging as a result of a macroscopic acquisition (cf. Section 2). This is the overarching task — and justification — for the theoretical efforts in our community. If a measurement is too slow to track the transient processes, we are left (in each non-exchanging com- partment) with the t → ∞ macroscopic Bloch- Torrey equation, i.e., Eq. (1.1) with uniform effec- tive parameters D ( r ) → D ∞ , R ( r ) → R ( t ) | t →∞ , Ω( r ) → (cid:104) Ω( r ) (cid:105) . Its solution becomes trivial — mono- exponential relaxation and Gaussian diffusion (cf. Sec- tion 3); i.e., coarse-graining leads to the universal t → ∞ dynamics, albeit with nonuniversal macroscopic pa- rameters such as D ∞ . The mesoscopic information is now lost, as the signal is indistinguishable from that in a uniform medium. Effective macroscopic parameters are in general different from the intrinsic mesoscopic ones; for instance, D ∞ can be notably lower than the intrinsic water or axoplasmic diffusion coefficient. qt Imaging So far we managed to get away with looking at a single equa- tion (1.1) and wave hands based on drawing parallels with concepts developed in physics. It is now time to introduce basic notations; the content of this subsection should be fa- miliar to anyone actively working in dMRI. In what follows, for simplicity we will confine ourselves only to the mesoscopic structure as related to diffusion, and will assume the relaxation effects to be trivial (at least in each tissue compartment), setting R ( r ) → R , and a uniform voxel-wise Larmor frequency, Ω( r ) → (cid:104) Ω (cid:105) . (The nontrivial R ( r ) and Ω( r ) modify apparent diffusion metrics [12, 13]; this is beyond the scope of our review.) This allows us to factorize the magnetization m ( t, r ) ≡ e − R t − i (cid:104) Ω (cid:105) t ψ ( t, r ) , where ψ ( t, r ) is not subjected to the relaxation and frequency shift and obeys the following equation ∂ t ψ ( t, r ) = ∂ r [ D ( r ) ∂ r ψ ( t, r )] − i g ( t ) r ψ ( t, r ) . (1.3) FIG. 3. The parameter space of dMRI is at least two-dimensional: By increasing q one accesses the pro-gressively higher-order diffusion cumulants, Sec. 1.8,whereas the dependence along the t -axis reflects theirevolution over an increasing diffusion length scale L ∼√ Dt , Eq. (1.12). The b -value alone does not uniquelydescribe the measurement, unless diffusion in all tissuecompartments is Gaussian; contour lines of b = q t areschematically drawn in beige. Large- q limits: Top-leftis high-resolution limit L ( t ) (cid:28) l c , ql c (cid:29) , Sec. 1.5 (i) ;middle is the L ( t ) (cid:38) l c , ql c (cid:38) limit of probing thepore correlation function, Sec. 1.6. The hierarchy ofdMRI models (pictures), cf. Fig. 4, as well as the cumu-lant representation with different number of terms, cf.Fig. 5, are superimposed. The decrease of the signalfrom axonal bundles parallel to the increasing gradientis shown by their darkening (top right). In Section 2 wemove along the t -axis at low q , and in Section 3 we movealong the q -axis at asymptotically long t . Section 4 is de-voted to effects beyond this diagram, contained in voxel-averaged products of propagators at different t and q . We focus here on the most easily interpretable measure- ment with very narrow (i.e., short) gradient pulses. As we now discuss, serendipitously, this measurement accesses the propagator of the mesoscopic diffusion equation, which (cf. Sec. 1.3) describes evolution of particle density ρ ( t, r ) ∂ t ρ ( t, r ) = ∂ r [ D ( r ) ∂ r ρ ( t, r )] . (1.4) The fundamental solution of Eq. (1.4), or diffusion propa- gator G t ; r , r , satisfies this equation ∂ t G t ; r , r = ∂ r [ D ( r ) ∂ r G t ; r , r ] + δ ( t ) δ ( r − r ) (1.5) with the point-like and instant source at r = r . The source term corresponds to the solution with zero particle density for t < and with the initial condition δ ( r − r ) instantly appear- ing at t = 0 . The solution is thus proportional to the unit step function, θ | t> = 1 and θ | t< = 0 , such that ∂ t θ ( t ) = δ ( t ) . The propagator G t ; r , r is a fundamental quantity describing the diffusion process around the point r , with a meaning of the probability distribution function (PDF) of molecular dis- placements r − r over time t . (This PDF can be sampled using Monte Carlo simulations by releasing random walkers all at once from the point r .) Of course, since the local tissue structure is different around each initial point r , the propa- gator G t ; r , r depends on the points r and r separately. The fundamental connection between the diffusion process (1.4) and the NMR measurement stems from the gradient- dependent phase of ψ ( t, r ) as described by Eq. (1.3). In the The focus on narrow pulses helps one to gain physical intuition. Finitepulse-width δ for relatively weak gradients has an effect of a low-pass fil-ter, filtering out the frequencies (cid:38) /δ , acting on the narrow-pulse solution[14–16], cf. models of restricted [17–19] and hindered [20–23] diffusionrelevant for the brain. For strong gradients, long pulses lead to the local-ization regime, Sec. 1.10. Arbitrarily-shaped pulses in the Gaussian phaseapproximation will be considered in Sec. 2.2. limit of narrow pulses g ( τ ) = q · [ δ ( τ − t ) − δ ( τ )] and the initial condition as in Eq. (1.5), the magnetization ψ ( t, r ) dif- fers from G t ; r , r by the position-dependent phase e − i q ( r t − r ) acquired during the gradient application. The diffusion- weighted signal , which is a net magnetization (cid:82) d r m ( t, r ) in a voxel, S ( t, q ) S ( t, q ) | q =0 = (cid:90) d r d r t V e − i q ( r t − r ) G t ; r t , r ≡ G t, q (1.6) becomes equivalent to a spatial Fourier transform of the voxel-averaged propagator G t, r ≡ (cid:104)G t ; r + r , r (cid:105) r = (cid:90) d r V G t ; r + r , r . (1.7) In Eq. (1.6) we divided by the voxel volume V , such that the unweighted signal (the right-hand side) is normalized to unity. A thorough discussion can be found e.g., in ref. [24]. Note that exact “local” propagator G t ; r t , r is not transla- tion invariant, i.e., it depends on the absolute coordinates r t , r (and time t ). The voxel-averaging in Eq. (1.6) automat- ically restores translation invariance, which means that the measured propagator G is parameterized by the two variables: the spatial displacement r ≡ r t − r and the diffusion time interval t (equivalently, by q and t ). Hence, the parameter space of dMRI fundamentally con- sists of q and t , Fig. 3 (here we dropped the directionality in q to not overload the picture). Literally speaking, map- ping the diffusion propagator in the space of q and t can be referred to as qt Imaging. For multiple diffusion encoding, cu-tie imaging, or qtI (noun) : A noninvasive medical imagingtechnique for spatio-temporal mapping of the diffusion propagator in softtissues to quantify tissue structure below the nominal MRI resolution. Ofcourse, it is nothing but the familiar q -space imaging [14, 25, 26] sampledat various t , but don’t we all need a new acronym once in a while? which maps a more complex object than the diffusion prop- agator (Section 4), the parameter space in principle depends on the multiple q and t intervals. The so-called b -value [27] has historically become the of- ten single-quoted measurement parameter. However, it only defines the measurement if diffusion is Gaussian in every compartment , in which case the diffusion propagator G (0) t, q = θ ( t ) e − Dq t ≡ θ ( t ) e − bD , b ≡ q t (1.8) in each compartment is determined solely by the parameter combination q t . Schematically, the contour lines of constant b are outlined in Fig. 3. In general, for anisotropic tissues such as brain white matter, Gaussian diffusion in each com- partment is described by the diffusion tensor, bD → b ij D ij , where the b -matrix [28] b ij = q i q j t . The Gaussian limit (1.8), and its more general anisotropic Gaussian limit, are hallmarks of “full” coarse-graining, which occurs in the t → ∞ limit, cf. Fig. 2. In this case, no matter how structurally complex the tissue, it can be modeled as a sum of (anisotropic) Gaussian signals. Section 3 will be de- voted to the picture of multiple Gaussian compartments (the Standard Model), cf. the column of pictures at long t in Fig. 3. The three regimes From the unifying coarse-graining point of view, we can now categorize biophysical models of diffusion, Fig. 4, into the following three regimes. In either of the regimes, the the- oretical treatment simplifies. The regimes can be arranged according to the increasing diffusion length L ( t ) relative to characteristic mesoscopic tissue length scales: (i). No coarse-graining has yet occurred. If the local D ( r ) varies in space over the correlation length scale l c , then for L ( t ) (cid:28) l c and ql c (cid:29) , each molecule senses its own, locally homogeneous D ( r ) . In this high-resolution limit [24], the signal S ( b ) (cid:39) (cid:82) d D P ( D ) e − bD is a Laplace transform of the his- togram P ( D ) of all the local values D ( r ) . A more rele- vant to biology situation occurs when instead of smooth D ( r ) variations, there are sharp barriers. The relevant parameter is then the net surface-to-volume ratio S/V of all barriers (e.g., cell walls). For times such that L ( t ) S/V (cid:28) , one observes the S/V universal short- time limit of the diffusion coefficient [29]. (ii). Coarse-graining over the structural disorder [30] results in the power-law approach t − ϑ of the instantaneous dif- fusion coefficient D inst ( t ) towards the t → ∞ limit D ∞ . Here, the power-law exponent ϑ is connected to the large-scale behavior of the density correlation func- tion of the hindrances to diffusion, and to the spatial dimensionality, yielding qualitatively distinct behavior along [30, 31] and transverse [20, 31] to the neurites in the brain. In Section 2 we argue, following ref. [30], that the more heterogeneous, or “disordered”, the sam- ple is, the slower the approach (the smaller the expo- nent ϑ ). Conversely, in ordered media, such as in the model of perfectly ordered membranes [30, 32], the ap- proach of D inst ( t ) towards D ∞ is exponentially fast. (iii). Complete coarse-graining. Diffusion in each non- confining tissue compartment has approached its t → ∞ Gaussian (tortuosity) limit, as discussed above (cf. also a more detailed discussion in Sec. 1.9 below). If there is no exchange between compartments, we ob- tain the most common, “multi-exponential” model. For neuronal tissue, the compartments are anisotropic due to the presence of effectively one-dimensional neurites. In Section 3, we introduce the “Standard Model” of neuronal tissue that accounts for the neurites with as- sociated extra-neurite space, and with an orientation dispersion (Fig. 8). While known under a plethora of names and acronyms [33–47], from the physics stand- point, this is practically the same model, with differ- ences in the parameterization of the neurite orientation distribution function and variations in the descriptions of the extracellular space, as well as in the model pa- rameter estimation procedures and employed parame- ter constraints. The crossover between regimes (i) and (ii) occurs when the diffusion length, L ( t ) , is commensurate with the characteris- tic length scale of the structural disorder. The instantaneous diffusion coefficient D inst ( t ) decreases with time within this crossover; while no general results are available there, it can be studied using numerical simulations. Working in the t → ∞ limit (iii) can only give us com- partment volume fractions and their diffusion coefficients. Coarse-graining has already occurred and apparently washed out all traces of other microstructural parameters. Determining characteristic µ m-level length scale(s) l c , such as the correlation length of the arrangement of tissue building blocks (e.g., disk radii in Fig. 2), is in principle pos- sible using deviations from the Gaussian signal shape. In the spirit of Fig. 3, varying either t or q can yield the sensitivity of the diffusion signal (propagator) to the length scale, via the diffusion length (cid:112) D ( t ) t [cf. Eq. (1.12) below], and via /q , respectively. However, as we now discuss, these theoretically distinct ways are not that different in practice, because attain- ing q ∼ /l c at times t (cid:29) t c practically requires sensing the signal contributions that are small at least as some positive power of the small ratio t c /t (cid:28) , where t c ∼ l c /D ( t c ) . Stanisz 1997Assaf 2004 Assaf 2008Alexander 2010 Jespersen 2007, 2010 Fieremans 2010, 2011Sotiropoulos 2012Zhang 2012Mitra 1993 Burcaw 2015, Fieremans 2016Novikov 2014 Callaghan 1979Yablonskiy 2002Kroenke 2004Novikov 2011 Novikov 2014Sigmund 2014Fieremans 2016Tanner 1978, Powles 1992, Sukstanskii 2004Novikov 2014 Reisert 2014Jelescu 2016Reisert 2017Veraart 2017Novikov 2018
FIG. 4.
Models are pictures...
Here they are drawn with coarse-graining occurring, roughly, from left to right. References: Mitra 1993 [29],universal short- t limit; Novikov 2014 [30], universal long- t behavior; Burcaw 2015 [20] and Fieremans 2016 [31], long- t behavior transverseand along WM fibers; Tanner 1978 [48], Powles 1992 [49], Sukstanskii 2004 [32], periodic 1-dimensional lattice; Novikov 2011 [50], randompermeable barriers in any dimension, and its application to myofibers (Sigmund 2014 [51] and Fieremans 2016 [52]); Callaghan 1979 [53], firstmodel of diffusion inside random narrow cylinders; Yablonskiy 2002 [54], diffusion in finite-diameter cylinders modeling lung alveoli; Stanisz1997 [55], first model for WM fiber tracts made of ellipsoids; Assaf 2004 [56], CHARMED; Assaf 2008 [57], AxCaliber; Alexander 2010[58], ActiveAx; Kroenke 2004 [33], NAA diffusion inside neurites. The widely adopted t → ∞ picture of narrow “sticks” for the neurites,embeded in the extra-neurite space (the Standard Model): Jespersen 2007 [36], Jespersen 2010 [37], Fieremans 2010 [38], Fieremans 2011(WMTI) [39], Sotiropoulos 2012 (Ball and rackets) [40], Zhang 2012 (NODDI) [41], Reisert 2014 (MesoFT) [42], Jelescu 2016 (NODDIDA)[43], Reisert 2017 [44], Veraart 2017 (TEdDI) [47], Novikov 2018 (LEMONADE [45], RotInv [46]). Le
FIG. 5. ...while representations are formulas.
References: Le Bihan 1991 [59], first biexponential representation of dMRI signal frombrain; Basser 1994 [28], diffusion tensor imaging (DTI); Jensen 2005 [60], diffusion kurtosis imaging (DKI); Kiselev 2011 [15], cumulantexpansion; Novikov 2008 and 2010 [24, 61], effective medium theory (transverse relaxation and diffusion, correspondingly); ¨Ozarslan 2013[62], expansion in harmonic oscillator basis; Yablonskiy 2003 [63], inverse Laplace transform (multi-exponential representation). Anisotropicmulti-exponential representations: Wang 2011 [64], diffusion basis spectrum imaging (DBSI); White 2013 [65], restriction spectrum imaging(RSI); Scherrer 2016 [66], distribution of anisotropic microstructural environments in diffusion-compartment imaging (DIAMOND).
Varying t amounts to literally observing the diffusive dy- namics for short times, when the coarse-graining has not yet fully occurred, such as during the regimes (i) and (ii) above. In our example in Sec. 1.3, to identify the presence of the small disks, one could try, e.g., detecting time dependence in D inst ( t ) or in D ( t ) at t ∼ r /D . The random permeable barrier model [50, 52], a candidate for diffusion transverse to myofibers and for one-dimensional hindrances along the neu- rites [30], allows one to trace the effect of coarse-graining in D inst ( t ) or in D ( t ) across all the regimes (i) – (iii) . Varying q , by employing strong narrow gradients (with width δ (cid:28) t c so that q is well-defined [14–16, 18, 20–23]), can in principle allow one to unravel the coarse-graining, i.e., to observe features at q ∼ /l c even when t (cid:29) t c . How- ever, the price to pay for accessing such fine structures is the suppression of the signal. To give an example [67], consider diffusion in a porous medium with connected pores and iso- lated grains, cf. Fig. 3 (center-top drawing). For long dif- fusion times, the pore structure is effectively coarse-grained (similar to diffusion in-between the small grains in Fig. 2). While having the overall Gaussian-like envelope ∼ G (0) t, r with D ≈ D ∞ , the diffusion propagator ∼ Γ( r ) G (0) t, r (up to an overall normalization) replicates the pore shape on the fine scale of the order of l c , Fig. 3, with the density correlation function Γ( r ) of the pore space arising due to voxel averag- ing. Correspondingly, in q -space, G t, q ≈ G (0) t, q + 1 φ (cid:90) d d q (cid:48) (2 π ) d G (0) t, q − q (cid:48) Γ( q (cid:48) ) (1.9) where the first term originates from sample-averaged (cid:104) Γ( r ) (cid:105) ≡ φ , the pore water fraction, and represents the aver- age spin density spread, while the product [Γ( r ) − φ ] G (0) t, r be- comes a convolution (the second term) of this envelope with the correlation function in q -space (i.e., the pore space power spectrum) Γ( q ) . The longer the time, the sharper is the Gaus- sian propagator G (0) t, q in q -space, and the less is the blurring of the pore correlation function induced by the convolution. However, longer times result in a stronger power-law suppres- sion of this nontrivial convolution term, whose magnitude can be estimated as ∼ [ l c /L ( t )] d in d dimensions. This estimate comes from noting that we are essentially averaging the pore density fluctuations over the “diffusion volume” L d ( t ) ; for the short-range disorder in the grain placement, these fluctua- tions, ∼ at the scale ∼ l c , become uncorrelated at the scales much beyond l c . When a tissue consists of non-exchangeable compartments of different nature, one can tune the experiment to focus on one or the other. Consider an example of a tissue consist- ing of a non-confining compartment (e.g., the extra-cellular space), and a fully confining, i.e., restricted one (e.g., cells In the language of emergent phenomena, Sec. 1.2, this would be anal-ogous, e.g., to using neutron scattering (with large wave vectors q (cid:38)
10 nm − ) to resolve atomic structure of the fluid, for which the coarse-grained (large t and r (cid:29) /q ) continuous description is classical hydro-dynamics. of size a ). The diffusional coarse-graining in the latter stops at the time t a ∼ a /D . The whole medium possesses two relevant scales, L ( t ) and a , which are markedly different for long times when L ( t ) (cid:29) a . An experimentalist work- ing in this limit has a choice of selecting the wave number q of diffusion measurements (the strength of diffusion weight- ing). The choice q ∼ /a (diffusion diffraction [68]) en- ables measuring the size of the restricted compartment, but strongly suppresses the signal from the permeable one. The choice q ∼ /L ( t ) enables observation of diffusion dynam- ics in the permeable compartment. The signal from the re- stricted compartment remains unsuppressed, such that both the signal attenuation − ln S ∼ ( qa ) (cid:28) and the diffu- sivity D ( t ) ∼ a /t (cid:28) become negligible; one can then formally treat such a compartment as Gaussian with its dif- fusion coefficient D → . This leads to the picture of zero- radius “sticks” in the brain (neurites with zero diffusivity in the transverse direction), cf. Sec. 3.1 below. Models are pictures (Fig. 4), exemplifying a rough sketch of physical reality, specified by their assumptions meant to simplify nature’s complexity. This simplification relies on av- eraging over the irrelevant degrees of freedom, and keeping only a handful of relevant parameters describing the corre- sponding effective theory. Model assumptions are therefore a claim for the relevant parameters. They are more important than mathematical expressions, as they prescribe a parsimo- nious way to think about the complexity. Model validation is thereby validation of our frame of thinking. A representation could be defined as a model-independent mathematical expression used to store, to compress, or to compare measurements. It can be realized as a function with a few adjustable parameters or a set of coefficients for a decom- position in a basis (cf. Fig. 5 for a few most commonly used representations in dMRI). In contrast to models, representa- tions are as general as possible, and have very little assump- tions. As there are infinite ways to represent a continuous function, the choice of representation is often dictated by con- venience or tradition. Practically, not all representations are equivalent because one only uses a few basis functions rather than an infinite set; from this standpoint, sparser representa- tions are more favorable. By construction, representations do not carry any particular physical meaning and hence do not immediately invoke any picture of physical reality; one can say that representations are formulas. A detailed discussion on the choices between modeling and representing can be found in ref. [69]. In this Review, we mostly focus on models; however, there exists one fundamen- tally important representation that we will cover now. The ubiquitous nature of Gaussian diffusion, at least for suffi- ciently long times, has prompted a Taylor expansion [15, 70]: ln G ( t, q ) (cid:39) − D ij ( t ) q i q j + 16 ( ¯ Dt ) W ijkl ( t ) q i q j q k q l − . . . (1.10) in the powers of q , describing the deviation from the Gaus- sian form (1.8). The summation over the repeated coordinate indices i, j, ... = 1 ... is implied throughout; ¯ D = d D ii is the mean diffusivity used for normalization. The sym- metric tensor W ijkl is called the diffusional kurtosis ten- sor , while the kurtosis in a given direction, ˆn , is defined as K ( ˆn ) = ¯ D W ijkl n i n j n k n l / ( D ij n i n j ) [70]. The propagator expansion (1.10) stems from the corre- sponding cumulant expansion in probability theory noticed almost a century ago by Fisher and Wishart [71, 72]. For dif- fusion, only even orders in this series are nonzero due to the time-reversal symmetry in the absence of the bulk flow. Typi- cally, the Taylor series (1.10) converges within a finite radius in q which is model-dependent [73, 74]. A general diffusion propagator will have all even cumu- lant terms D ij ( t ) , W ijkl ( t ) , . . . nonzero and diffusion time- dependent [15, 24]. Experimentally we often access only a few first terms, especially when using low diffusion weight- ing on clinical systems. (We assume the narrow-pulse limit throughout. In Sec. 2.2 we discuss in detail how the lowest order of the ideal cumulant series (1.10) is modified by the arbitrary gradient shape). Upon coarse-graining, for a given tissue compartment the higher-order terms W ijkl , . . . flow to zero, such that the sig- nal approaches the Gaussian form (1.8) as t → ∞ . In this limit, the higher-order cumulant terms of the net diffusion propagator can originate only from the partial contributions from different tissue compartments (since a sum of Gaussians is non-Gaussian). For any t , the series (1.10) generates the cumulants (cid:104) x j x j . . . (cid:105) c (see e.g., refs. [15, 71, 72] for definition) of the PDF of molecular displacements (1.6), via taking derivatives at q = 0 , such as (cid:104) x i x j (cid:105) ≡ (cid:90) d r x i x j G t, r = − ∂ ∂q i ∂q j (cid:12)(cid:12)(cid:12)(cid:12) q =0 (cid:90) d r e − i qr G t, r . (1.11) Based on such averages, it is conventional to define the cu- mulative diffusion coefficient D ( t ) = (cid:104) x ( t ) (cid:105) t , (1.12) or, more generally, the cumulative diffusion tensor D ij ( t ) = (cid:104) x i ( t ) x j ( t ) (cid:105) t (1.13) (a symmetric × matrix with 6 independent parameters in 3 dimensions). These objects are defined in terms of the Since G t, r is written in terms of the relative displacements r = r t − r ,we re-denote δx i ( t ) → x i ( t ) in Eq. (1.12) to simplify the notation, anddrop the dependence on the initial position in the view of the translationalinvariance property (1.7). average rate of change of the mean-squared molecular dis- placement over the whole interval [0 , t ] (in contrast to the instantaneous rate of change (1.2) above). The linear estimation problem for D ij ( t ) , referred to as the diffusion tensor imaging (DTI) , has been solved by Basser et al. [28]. It requires a diffusion measurement along at least more, e.g., the b = 0 (unweighted) image. Likewise, the linear estimation problem for both the diffu- sion and kurtosis tensors, via the expansion up to ∼ q ∼ b , called diffusion kurtosis imaging (DKI) , has been introduced by Jensen et al. [60, 75]. It involves the 4 th order cumulant (cid:104) x i x j x k x l (cid:105) c related to W ijkl ( t ) . The number of parameters are now , hence one needs at least two b (cid:54) = 0 shells in the q -space, and at least 15 non-collinear directions. The weights for unbiased estimation of diffusion and kurto- sis tensors for non-Gaussian MRI noise were found recently [76]. A general method to calculate the number of parameters for a given order l c of the cumulant series (1.10) in 3 dimensions is based on the SO(3) representation theory (known in physics as theory of angular momentum in quantum mechanics). A term ∼ q l c of even rank l c is a fully symmetric tensor, which can be represented as a sum of the so-called symmetric trace- free (STF) tensors of ranks l c , l c − , . . . , , [77]. Each set of l + 1 STF tensors of rank l realizes an irreducible representation of the SO(3) group of rotations, equivalent to a set of l + 1 spherical harmonics Y lm [77]. Hence, the total number n c of nonequivalent components in the rank- l c cumulant tensor is n c ( l c ) = l c (cid:88) l =0 , ,... (2 l + 1) = 12 ( l c + 1)( l c + 2) , (1.14) so that n c = 6 for DTI ( l c = 2 ) and n c = 15 for DKI ( l c = 4 ). Suppose we truncate the cumulant series (1.10) at an (even) term of rank l c = l max . Hence we determine all the param- eters of cumulant tensors (diffusion, kurtosis, ...) of ranks , , . . . , l max . The total number of independent parameters in the truncated series N c ( l max ) = l max (cid:88) l =2 , ,... n c ( l ) = 112 l + 58 l + 1712 l max (1.15) corresponding to N c = 6 , , , . . . for l max = , , , . . . . Hence, DTI yields 6 parameters, DKI yields
21, etc. (Here we did not include the proton density S | b =0 in our counting.) The cumulants D ij , W ijkl , . . . of the signal obtained via Taylor-expanding its logarithm in the (even) powers of q i , or DTI, contrary to a widespread misconception, does not assume Gaussiandiffusion, as it merely provides the lowest-order cumulant term D ij , andtells nothing about the higher-order terms in the series (1.10). DTI applica-bility is thus dictated by the kurtosis term W to have negligible bias on theestimated D ij , and the employed b -range is practically set by balancingthe bias when b is too large and precision loss when b is too small. b , correspond to the cumulants of the genuine PDF of molecular displacements r = r t − r only in the narrow pulse limit, and in the absence of the mesoscopic magnetic structure (uniform R and Ω ). When the finite gradient pulse duration δ is comparable to the time scale of the transient processes, the measurement acts as a low-pass filter with a cutoff frequency ∼ /δ [14–23]. For finite t , the diffusion propagator in a heterogeneous medium is never Gaussian. The existence of domains with slightly different “local” D ( r ) at a given coarse-graining scale necessarily yields the time-dependent D inst ( t ) , as well as the higher-order terms in q , such as q , in the Taylor ex- pansion of ln G (0) ( t, q ) [24, 30]. Upon coarse-graining, these terms gradually flow to zero, and D inst ( t ) → D ∞ , such that diffusion becomes Gaussian asymptotically as t → ∞ in each separate non-confining tissue compartment. This was the picture of Sec. 1.3, cf. Fig. 2. In particular, we implied that the diffusion coefficient decreases, as a result of the coarse-graining, towards its finite tortuosity asymptote D inst ( t ) | t →∞ ≡ D ∞ > . How reliable is this picture? What does it take to destroy it? Existence of finite D ∞ is equivalent to mean squared dis- placement (cid:104) x ( t ) (cid:105) (cid:39) D ∞ t growing linearly with time for sufficiently long t , — this is a direct consequence of the defi- nition (1.2). One says that diffusion asymptotically becomes “normal”, i.e., the PDF of molecular displacements over a sufficiently large t approaches normal (Gaussian) distribu- tion, cf. Eq. (1.8) with D → D ∞ . Of course, if there are two or more non-exchanging tissue compartments, the total dis- tribution will be non-Gaussian (as a sum of Gaussians with different D ∞ ), but this non-Gaussianity is in a sense trivial; the total D ∞ would still exist (given by a weighted average for the corresponding compartment values) [38], and the scal- ing (cid:104) x ( t ) (cid:105) ∼ t at large t would hold. There exists a radical alternative, when (cid:104) x ( t ) (cid:105) ∼ t α for t → ∞ , with exponent α (cid:54) = 1 — the so-called anomalous diffusion [78]. According to the definition (1.2), D ∞ = 0 for α < (sub-diffusive behavior), and D ∞ = ∞ for α > (super-diffusive behavior). In other words, observation of anomalous diffusion is equivalent to stating that the macro- scopic diffusion coefficient D ∞ does not exist . (The trivial case D ( t ) ∼ a /t for a confining compartment of size a is not considered anomalous; (cid:104) x (cid:105) ∼ a , α = 0 .) The absence of D ∞ in a non-confining medium is always a drastic claim: it is potentially exciting yet should be thor- oughly validated, because the underlying physical assump- tions yielding α (cid:54) = 1 are generally quite peculiar and excep- tional, as we discuss below. In neuronal tissue, one always observes finite D ∞ in non-confining compartments (e.g., in the extra-cellular space), Section 2, hence diffusion is empir- ically never anomalous [20, 30, 31] for brain dMRI. We are not reviewing the MRI literature on anomalous diffusion, since our
From the point of coarse-graining, anomalous diffusion means that the sample never quite looks homogeneous — for example, a fractal has a self-similar structure, which im- plies similar statistics of static structural fluctuations at every length scale. In other words, when the coarse-graining over some scale has taken place, a larger scale looks statistically similar, so that the already averaged structural features are never forgotten, since they are reproduced again and again. In contrast, the structure in Fig. 2 implies that this memory is forgotten for each of the two length scales, correspondingly on the two well-defined time scales. Tissues empirically do not look self-similar; usually, when we look at a histological slide without a scale bar, we can still roughly say at which resolution the sample is imaged be- cause usually cell size is well defined (for a given tissue type) — otherwise, medical students would not pass their pathol- ogy exams. For instance, when we look at cross-sections of white matter tracts, the majority of axons are of the order of ∼ µ m in diameter [79–83], and the section does not look the same when magnified by factors of 10, 100, or 0.1, 0.01, etc. A more quantitative statement can be made by studying large-distance scaling behavior of the density-density correla- tion function of the tissue structure; recent investigation [20] confirms that the structural fluctuations in white matter tracts are short-range (and not diverging at large length scales). When can anomalous diffusion arise? In a broader context, this fundamental question has been extensively studied for the Fokker-Planck equation ∂ t ψ ( t, r ) = ∂ r [ D ( r ) ∂ r ψ ( t, r )] − ∂ r [ v ( r ) ψ ( t, r )] , (1.16) where in addition to the “diffusive” flow j ( t, r ) = − D ( r ) ∂ r ψ ( t, r ) (Fick’s law), one considers mesoscopic ran- dom flow because of some stationary local “velocity”, or “force” field v ( r ) (imagine active streams, such as vortices or currents in an ocean [84]). Equation (1.16) arises as a con- servation law ∂ t ψ = − ∂ r · j , where the total flow j = − D ( r ) ∂ r ψ ( t, r ) + v ( r ) ψ ( t, r ) . It turns out that the presence of the random flow term v ( r ) with short-range spatial correlations can drastically change the dynamics in dimensions d ≤ and drive the system away from the Gaussian diffusion. In dimension d = 1 , random force field causes sub-diffusive behavior (cid:104) δx (cid:105) ∼ ln t , a famous result by Sinai [85]. In d = 2 dimensions, super- diffusive behavior occurs when the flow v ( r ) is solenoidal, goal here is to discuss models which are relevant to observable diffusioneffects in neuronal tissue. A curious reader can find occasional claims ofanomalous diffusion, or dMRI signal as a stretched-exponential. We arenot aware of examples of a constructive derivation of the non-Gaussianfixed point [3, 4] starting from the stationary mesoscopic disorder withproperties relevant to the brain. Hence, these claims can merely be viewedas postulates “proven” by fitting in a finite range of t or q . If the model’sfunctional form contradicts the physics of the signal, the estimated pa-rameters will depend on the range of t and q , thereby characterizing theparticular measurement scheme, rather than the tissue [69]. div v ( r ) = 0 , and sub-diffusive if it is potential, curl v ( r ) = [84, 86–88]. In the absence of the random forces, v ≡ , small fluc- tuations in D ( r ) do not destroy the “trivial” Gaussian fixed point in dimensions d > [89, 90] (cf. footnote 2). In other words, for the spatial short-range disorder in D ( r ) to become relevant (i.e., to increase under the renormalization group flow), and for the anomalous diffusion to take over, the spatial dimension should formally be d = 0 . What this tells is that it is very difficult, without the flow term, to break the Gaussian fixed point of the finite D ∞ , at least starting from the weak disorder in D ( r ) . Extremely strong disorder, which is specially tuned, can induce the percolation transi- tion [91] D ∞ → ; another possibility for destroying finite D ∞ is to create the disorder in the mesoscopic D ( r ) with anomalously divergent spatial fluctuations [92, 93]. To the best of our knowledge, neuronal microstructure is compati- ble neither with a percolation transition, nor with diverging structural fluctuations [20]. Another class of phenomena where anomalous diffusion takes place corresponds to systems with slow dynamics, orig- inating from a broad distribution of time scales, such that the waiting time distribution p ( τ ) ∼ /τ µ between “hops” of random walkers has a power-law tail whose first moment di- verges, < µ < . Such broad distributions can emerge, e.g., in highly disordered amorphous solids, where escape times τ from various “traps” for electrons are distributed as a power law, first postulated by Scher and Montroll [94]. For the traps, the long tail in p ( τ ) can arise due to an ex- ponentially strong dependence of the activation rate τ − ∼ e − E/k B T on the energy barrier E at temperature T , such that a relatively flat distrubution ˜ p ( E ) can result in the L´evy-like p ( τ ) = ˜ p ( E ( τ )) /τ . Hopping with traps may lead to anoma- lous transport [95] and fluorescence [96]. Anomalously slow dynamics also occurs in viscoelastic systems where elemen- tary components are strongly coupled (Rouse polymer chain [97] of monomers tied to each other by elastic springs and undergoing Langevin dynamics). The simulated dynamics of single protein molecules [98] and of colloidal tracers re- stricted by crowded dynamical environments such as an F- actin network [99] can exhibit such a broad distribution of time scales [100]. While an active area of investigation, the anomalously slow dynamics is always characterized by strong disorder (e.g., broadly distributed traps) and/or interactions among the random walkers (e.g., parts of a polymer). To recap, coarse-graining over an increasing diffusion length L ( t ) provides a physical picture for time-dependent diffusion in mesoscopically disordered samples. This picture implies gradual “forgetting” of the memory about the struc- tural heterogeneities. In an overwhelming majority of sys- tems, the macroscopic dynamics is characterized by a Gaus- sian fixed point, the absence of long-term memory, and an asymptotically normal diffusion. In short, diffusion is almost always non-Gaussian, but almost never anomalous . In the brain, it is not anomalous specifically because the density fluctuations of brain structural units do not diverge at large scales, traps for water molecules with broad distribution of escape times do not seem to exist, and the “active” flow ef- fects (microstreaming, axonal transport) are negligible [101]. Before proceeding to review brain dMRI models, we would like to mention what we have left out, because of limited rel- evance to brain dMRI as of today, and/or due to exhaustive coverage elsewhere. On the methodological front, the leitmotif of the review is the language of coarse-graining, Fig. 2. It is most intuitive for modeling structurally disordered systems, typical for biol- ogy, cf. modeling the time-dependent diffusion in Section 2, based on including all the restrictions into the spatially vary- ing D ( r ) in Eq. (1.1). This took precedent to approaches to fully confining or periodic geometries, conventionally solved by formulating Eq. (1.1) as the Laplace equation with bound- ary conditions, thoroughly reviewed in ref. [102] in the con- text of diffusion in porous media. We also left out the physics of the localization regime , where diffusion in a strong constant gradient suppresses the signal everywhere except next to pore walls, within the gradient-dependent dephasing length L g = ( D /g ) / , which leads to signal decay [103–105] − ln S ∼ L ( δ ) /L g ∼ D / g / δ . This is an example where effects of non-narrow pulses lead to decoupling of the gradient magnitude g and the gradient pulse width δ in the narrow-pulse combination q = gδ . The “edge enhancement” also amplifies the role of the permeability of the walls [106]. Brain structures seem to be too small for the edge effects to be relevant, but such phenomena can become important in body dMRI. Playing with δ , e.g., using short-wide pulse combinations, we or one can map the Fourier transform of the shape of the confining pore [107, 108], which again requires prohibitively strong gradients for the narrow axons and dendrites in the brain, but is applicable in porous media NMR. The relevance of pulse width δ would add an extra dimension to the phase diagram in Fig. 3. Detailed review of practical aspects of dMRI measure- ments and biological applications are beyond our scope here. The reader is referred to the review [109] for implementation details of dMRI measurements, recent reviews of dMRI in white matter [110] and in cancer [111], as well as to other articles in this Special Issue.
2. TIME-DEPENDENT DIFFUSION IN NEURONAL
TISSUE
Everything should be made as simple as possible, butnot simpler Albert Einstein
The intuition of Sec. 1.3 suggests that the time-dependence of the diffusion coefficient defined as either Eq. (1.2) or
Eq. (1.12), is a hallmark of the mesoscopic structure, and the associated time scale can be translated into the corre- sponding mesoscopic length scale. Identifying µ m-level tis- sue length scales is the ultimate test of our ability to “quan- we are indeed sensitive to the micro -structure? The focus of this Section is on determining tissue properties on such length scales. Fundamentally, observation of the time dependent over- all D ( t ) is significant because it tells that diffusion is non- Gaussian in at least one of the tissue compartments . Indeed, at the lowest order O ( q ) of the cumulant expansion (1.10) of the signal S = (cid:80) f i S i , contributions from non-exchanging tissue compartments S i add up, such that the total diffusion coefficient is a weighted average: D ( t ) = (cid:88) f i D i ( t ) , (cid:88) f i = 1 . (2.1) An overall time-dependent D ( t ) necessarily means that at least one of D i depends on t . In turn, the time-dependent D i ( t ) must necessarily lead to a nonzero kurtosis and higher- order cumulants [24] in the i th compartment, arising from the same mesoscopic heterogeneity which has not yet been fully coarse-grained — and, hence, may still be possible to quan- tify. Conversely, if diffusion is Gaussian in all tissue compart- ments, all D i = const , and the overall diffusion coefficient is time-independent. The overall kurtosis is then a nonzero constant just because a sum of Gaussians is not a Gaussian. We begin this Section by reviewing experimental data on the time-dependent diffusion coefficient and kurtosis in brain, and then discuss the two physically distinct regimes of time- dependence, according to the hierarchy of Sec. 1.5: the short- time regime (i) , and the long-time regime (ii) approaching the asymptotically Gaussian diffusion in each non-exchanging compartment. Certainly, we are almost never in a pure limit experimen- tally — rather, we are typically in some crossover, e.g., in- between the regimes (i) and (ii) . However, it is still impor- tant to understand the behavior of the system in certain limits where it can be modeled with more confidence. Performing experiments in such limits provides a way to validate mod- els through observing definitive functional dependencies on the measurement parameters [69]; thus-identified relevant de- grees of freedom for tissues can then be incorporated into more complex theories of the crossover behavior relevant to a broader range of dMRI studies, and for clinical translation. Empirically, observing time-dependence of diffusion in brain tissue has been challenging because this effect occurs at time scales associated with diffusion across length scales featuring This argument can also be extended onto the regime of slow exchangebetween compartments, since Eq. (2.1) turns out to be valid in that regimein the long- t limit, following the coarse-graining argument [38] for gener-alizing the K¨arger model [112] to media with mesoscopic disorder. If theoverall D = const (i.e., the full coarse-graining is achieved in each ofthe compartments), and the overall kurtosis K ( t ) still depends on t , this t -dependence arises due to exchange [38, 60]. neurites (i.e., axons and dendrites). Typically, their diame- ters as well as the heterogeneities along them (e.g., spines, beads) are of ∼ µ m size, hence, the corresponding diffu- sion times are expected to be of the order of a few ms. Such short times are quite difficult to access, especially on human systems. Besides, the time-dependence is generally slow — which is theoretically expected due to its power-law character [30], as discussed below in Sec. 2.4 — therefore requiring a sufficiently broad range of times to detect. Time-dependence of the cumulative D ( t ) , Eq. (1.12), in brain tissue has been demonstrated using pulse gradient spin echo (PGSE) in several ex vivo studies for a range of diffu- sion times encompassing − ms [55, 113–115]. In vivo , time-dependent diffusion in both longitudinal and transverse directions was also observed in rat corpus callosum at t rang- ing from 9 to 24 ms [116], though another study yielded no change in the mean diffusivity of healthy and ischemic feline brain with respect to t between − ms [117]. In the human brain, it has been unclear for quite some time whether in vivo time-dependent diffusion properties can be observed. While several in vivo studies report no observ- able change over a broad range of diffusion times [118, 119], Horsfield et al. [120] reported time-dependent diffusion in several white matter regions at times ranging from 40 to 800 ms. Very recently, in vivo pronounced time-dependence in the longitudinal diffusivity and less pronounced time-dependence in the transverse diffusivity has been reported in several WM tracts of healthy human volunteers for relatively long diffu- sion times, t = 45 − ms, on a standard clinical scan- ner using stimulated echo acquisition mode (STEAM)-DTI [31]. Subsequently, a similar effect in the transverse direction to WM tracts was observed with STEAM-DTI in the range t = 48 − ms [121]. Oscillating gradient spin echo (OGSE) diffusion-weighted sequences are able to probe shorter diffusion time scales com- pared to conventional PGSE, and have clearly demonstrated time-dependent diffusion in the brain, including the obser- vation of time-dependent diffusivities in vivo in normal and ischemic rat brain cortex [122], as well as ex vivo in rat WM tracts [123]. By combining OGSE and PGSE, Pyatigorskaya et al. [124] observed time-dependent diffusion coefficient and a non-monotonic time-dependent kurtosis (with a maximum value K ≈ . at t ≈ ms) in healthy rat brain cortex at dependence in mouse cortex and hippocampus. In humans, Baron and Beaulieu [127] found eight major WM tracts and two deep gray matter areas to exhibit time-dependent diffu- sion using OGSE and PGSE, and Van et al. [128] reported a similar effect with OGSE in human corpus callosum. Fur- thermore, works using double diffusion encoding (cf. Section
4) indirectly point at the time-dependent nature of diffusion in brain tissue. Overall, while it is common to assume that diffusivities are approximately diffusion time-independent for t (cid:38) ms, the experimental data described above clearly demonstrates time- dependent diffusion both at short and long times. In what fol- lows, we describe the underlying theory for both limits and discuss the corresponding biophysical interpretation and po- In Sec. 1.4, we derived a general relation (1.6) between the dMRI signal and the diffusion propagator in the narrow-pulse limit. For gradient pulses g ( t ) of arbitrary shape, there is no such simple relation; the signal S [ g ( t )] is a functional of g ( t ) (i.e., a mapping of a function to a number). To obtain an explicit dependence of S on g ( t ) , one treats the gradient term in Eq. (1.3) perturbatively in g ( t ) , generalizing the cumulant series (1.10). Here, we will stay at the level of O ( g ) , the so-called Gaussian phase approximation (GPA), and describe the family of diffusion coefficients which define the second- order cumulant and carry the same information content, yet can be accessible using different techniques, Fig. 6. The GPA approximates the dMRI signal [14–16] S [ g ( t )] = (cid:104) e − iϕ (cid:105) ≈ e − (cid:104) ϕ (cid:105) (2.2) up to the second cumulant of the accumulated phase ϕ ( t ) = (cid:90) t d t (cid:48) g ( t (cid:48) ) r ( t (cid:48) ) = − (cid:90) t d t (cid:48) q ( t (cid:48) ) v ( t (cid:48) ) , (2.3) where we introduced the time-dependent wave vector q ( t ) = (cid:90) t d t (cid:48) g ( t (cid:48) ) (2.4) such that the Larmor frequency gradient g is given by its time derivative, g = ∂ t q . The first-order cumulant (cid:104) ϕ (cid:105) ≡ in the absence of the net flow. The balanced gradient condition sets q ( t ) | t>T ≡ at the end t = T of the gradient train interval. Eq. (2.4) generalizes the definition of q for narrow- pulse gradients, where q remained constant during an interval < t < T , cf. the propagator Eq. (1.6) with t = T . Writing (cid:104) ϕ (cid:105) as a double integral, and averaging over the Brownian paths, we obtain − ln S [ q ( t )] (cid:39) (cid:90) T d t d t q ( t ) (cid:104) v ( t ) v ( t ) (cid:105) q ( t ) , (2.5) where we from now on dropped the explicit vector notation of q and v (the corresponding tensor indices can be easily restored; one can think about isotropic media for simplicity). We can see that the diffusion process at the O ( q ) level is fully characterized by the autocorrelation function (cid:104) v ( t ) v ( t ) (cid:105) of molecular velocity, an even function of t − t in stationary media due to time translation invariance and time reversal symmetry of the Brownian motion. For uniform media, (cid:104) v ( t ) v ( t ) (cid:105) = 2 D δ ( t − t ) , which can be thought of as one of the equivalent definitions of the diffusion constant D . Technically, there is no such thing in nature as a zero-width δ ( t − t ) ; we can use this ap- proximation for simple liquids since the correlation time for forgetting the memory about the molecular collisions is of the order ∼ − ps, orders of magnitude below our ms- level time scales. We can say that diffusion in simple liq- uids is thereby Markovian (has no memory) on the relevant NMR time scales. This leads to the standard expression − ln S = bD with b = (cid:82) T q ( t ) d t , generalizing Eq. (1.8). For general mesoscopic media, microstructure introduces temporal correlations in positions and velocities of random walkers. For instance, if a walker just hit a wall, then its ve- locity will correlate negatively with the velocity just before the hit (since reflection and moving away from the wall is preferred), and this memory will last during the time depend- ing on the wall geometry and the presence of other restric- tions. To characterize such correlations, let us introduce the retarded velocity autocorrelation function [24] D ( t ) ≡ θ ( t ) (cid:104) v ( t ) v (0) (cid:105) , (2.6) where θ ( t ) is the unit step function, cf. Sec. 1.4. In terms of D ( t ) , Eq. (2.5) reads − ln S [ q ( t )] (cid:39) (cid:90) d t d t q ( t ) D ( t − t ) q ( t ) , (2.7) where the double integration can be extended over all real values of t , since q ( t ) is nonzero only on a finite interval. The time translation invariance of D allows us to rewrite the double intergral in the t -domain as a single integral in the ω -domain, by introducing the Fourier transform of D ( t ) , the dispersive diffusivity [24, 130] D ( ω ) = (cid:90) ∞ d t e iωt (cid:104) v ( t ) v (0) (cid:105) . (2.8) Eq. (2.7) can now be written in terms of the Fourier- transformed q ω = (cid:82) d t e iωt q ( t ) , as − ln S [ q ω ] (cid:39) (cid:90) d ω π q − ω D ( ω ) q ω . (2.9) Here, only Re D ( ω ) contributes, as Im D ( ω ) , odd in ω , yields zero after being integrated with an even function | q ω | . Equivalently, Im D ( ω ) does not contain extra information as it can be restored using the Kramers-Kronig relations [133]. The representation (2.9) underscores that, knowing the ve- locity autocorrelator D ( ω ) , one can evaluate the diffusion- weighted signal to O ( g ) for any gradient waveform g ( t ) . Conversely, by selecting a particular form of q ( t ) according to The real part, Re D ( ω ) , corresponds to the quantity called “ D ( ω ) ” in theNMR literature [129]. For anisotropic media, and for arbitrary q -space trajectories [131, 132], theintegrands in Eqs. (2.7)–(2.9) are q i ( t ) D ij ( t − t ) q j ( t ) , v i ( t ) v j (0) and q − ω,i D ij ( ω ) q ω,j respectively (with the sums over repeated indices). D ( ! )
S/V limit [67] was given in Sec. 1.5 (i) . Quantitatively, the short- t expansion of the diffu- sion coefficient (1.12) proceeds in powers of L ( t ) S/V , where L ( t ) = √ D t is the diffusion length: D ( t ) = D (cid:18) − √ π d SV (cid:112) D t + O ( t ) (cid:19) , (2.16) and S/V is the surface-to-volume ratio of the restrictions in d spatial dimensions. Identifying the √ t term practically in- volves very short diffusion times. Even for a red blood cell suspension, this limit was barely observable in the time do- main [139]; for the brain, with structural features even smaller than the red blood cell size, getting to this limit using PGSE is practically impossible due to very low b -values for short t . Hence, regime (i) is best accessed using OGSE. The corre- sponding functional form of D ( ω ) for Eq. (2.16) was recently derived in the N → ∞ limit [130] Re D ( ω ) (cid:39) D (cid:32) − d √ SV (cid:114) D ω (cid:33) . (2.17) For a finite total number N of oscillations, Eq. (2.17) is mod- ified, see Appendix C, by a correction factor c , Eq. (C5), in front of the / √ ω term. This factor approaches its N → ∞ limit c → rather fast, c − ∼ /N , such that c − < . as long as N ≥ for the cos , and N ≥ for the sin wave- forms. From directly comparing Eqs. (2.16) and (2.17), the rela- tion between OGSE frequency ν = ω/ π and diffusion time t = ∆ (in the narrow-pulse PGSE limit) is as follows [130]: S/V limit (i) : t = 964 · ν . (2.18) We note that Eq. (2.18) differs from the empirical relation (see, e.g., ref. [122]) wrong yet widely used: t = 14 ν , (2.19) which in fact is almost always incorrect [cf. Eq. (2.36) below]. Relation (2.19) originates from matching the b -value between one OGSE period and PGSE of the same duration. Since the whole notion of the b -value stems from Gaussian (i.e., time- independent) diffusion, it is not surprising that merely match- ing the diffusion attenuation between PGSE and OGSE for the constant D falls below the accuracy needed to define the diffusion time for the nontrivial, time-dependent case. Probing the short-time limit either in the time domain (Eq. (2.16)) or the frequency domain (Eq. (2.17)) potentially allows for decoupling the geometric effects of the surface-to- volume ratio S/V and the free diffusivity D . Recently this limit has been demonstrated in phantoms using both PGSE [141] as well as OGSE [142]. For in vivo brain measure- ment, OGSE provides the most practically feasible method, with the / √ ω dependency as the signature functional form (2.17) of this regime. In the healthy brain, this signature has so far never been observed, since presumably the achievable oscillation frequencies are still too low as compared to those needed to identify the effect of the restrictions from neurite walls with typical radius of curvature ∼ µ m, requiring dif- fusion times much below 1 ms (i.e., frequencies ν (cid:29) kHz). The search for the / √ ω regime has prompted using brain tumors with roughly spherical cells of larger size (about µ m), such that the required frequency range can be po- tentially accessible. Recently, the / √ ω functional form was observed by Reynaud et al. [143] in a mouse glioma model in the frequency range up to ω/ π = 225 Hz, which for the first time enabled the separation between the geometric ( S/V ) and “pure” diffusive ( D ) tumor features. Further combining the OGSE and PGSE methods has lead to the POMACE [144] and IMPULSED [145] methods for quantifying cell size and extra-cellular water fraction, cf. topical review [111]. Structural correlations via gradual coarse-graining Over time, random walkers probe the spatial organiza- tion of the sample’s microstructure, which makes the time- dependence of the diffusion metrics intricately tied to an in- creasingly large number of structural characteristics. Tech- nically, finding D ( t ) or D ( ω ) analytically in a realistic com- plex sample is nearly impossible as it amounts to including The often quoted relation t = ∆ − δ/ for the diffusion time of a finite-width PGSE is a myth for the same reasons. One can only say that themeasurement gives D ( t ) with t ≈ ∆ (with the accuracy of this approxi-mation controlled by δ ). More rigorously, the effect of finite pulse width δ creates a low-pass filter [14, 16] on D ( ω ) , whose effect is again model-dependent, see, e.g., Eq. (24) in ref. [20], as well as refs. [21–23, 140], forthe examples of this filter effect on the models of D ( ω ) relevant for brain. fusion coefficient D ( r ) and of the positions of all restrictions up to an infinitely high order. The intuition based on coarse-graining, Sec. 1.3, turns out to be helpful in solving this problem in the long time regime [30], when the diffusion coefficient (1.2) gradually ap- proaches its macroscopic (tortuosity) value (2.15). As men- tioned in Sec. 1.3, in the limit t → ∞ , any non-confining tissue compartment effectively looks completely uniform. Let us step back just a bit from t → ∞ and consider t long enough (yet finite) for the sample to look almost ho- mogeneous from the point of the diffusing molecules, Fig. 2, — no matter how heterogeneous it is in reality (e.g., at the cellular scale). In this limit, the problem of finding the diffu- sion propagator maps onto a much simpler problem of finding the diffusion propagator in a weakly heterogeneous medium (which is the corresponding effective theory), characterized by the diffusion equation (1.4) with | δD ( r ) | D ∞ (cid:28) , δD ( r ) ≡ D ( r ) − D ∞ . (2.20) This problem admits a perturbative solution [24, 30], with Eq. (2.20) defining a small parameter, as long as the macro- scopic (tortuosity) limit (2.15) exists, < D ∞ < ∞ (i.e., diffusion is not anomalous, which is practically always the case for dMRI in tissues, cf. Sec. 1.9). The lowest order in δD ( r ) vanishes, and the second order in the parameter (2.20) yields D inst ( t ) (cid:39) D ∞ + 1 d (cid:10) ( δD ( r )) (cid:11) | L ( t ) D ∞ (2.21) in d spatial dimensions. The last term in Eq. (2.21) involves the variance of the slowly-varying D ( r ) at a given coarse-graining length scale defined by the diffusion length L ( t ) . This variance decreases as a result of self-averaging , i.e., when different diffusing molecules on average begin to experience more and more similar mesoscopic structure with an increasing L ( t ) , such that any sample begins to approximate the ensemble of dif- ferent disorder realizations of D ( r ) more and more precisely. The always positive “fluctuation correction” to D ∞ (the last term) elucidates why the diffusion coefficient can only de- crease with t ; observation of its increase with diffusion time is a red flag for imaging artifacts. To be more rigorous, Eq. (2.21) can be expressed as [30] D inst ( t ) | t (cid:38) t (cid:39) D ∞ + 1 D ∞ d (cid:90) d d k (2 π ) d Γ D ( k ) e − D ∞ k t , (2.22) in terms of the power spectrum Γ D ( k ) = (cid:82) d r e − i kr Γ D ( r ) of the underlying effective D ( r ) | L ( t ) coarse-grained over the diffusion length L ( t ) corresponding to some sufficiently long time scale t , for which the relative deviation (2.20) from D ∞ is sufficiently small. The correlation function Γ D ( r ) = (cid:104) δD ( r + r ) δD ( r ) (cid:105) r (2.23) embodies the fluctuation correction in Eq. (2.21). We can see that diffusion indeed acts as a Gaussian filter (cf. Fig. 2 in Sec. 1.3), with a filter width ∼ L ( t ) ∼ √ D ∞ t , over the effective medium defined via the correlation function of the weakly heterogeneous D ( r ) . Hence, for sufficiently long t , Eqs. (2.21) and (2.22) be- come asymptotically exact with L ( t ) → ∞ , no matter how strongly heterogeneous the “true” (microscopic) D ( r ) is. From the renormalization group flow standpoint, we can say that the time-dependent corrections (last terms of Eqs. (2.21) and (2.22)) to the asymptotically Gaussian propagator be- come irrelevant as a result of integrating out the fluctuations of the locally varying D ( r ) over larger and larger scales. Likewise, the kurtosis and higher-order cumulants in this compartment will decay to zero, as governed by similar fluc- tuation terms. How to relate the time-dependence (2.21) and (2.22) to the mesoscopic structure? Here, one realizes [30] that the coarse- grained D ( r ) | L ( t ) depends on the similarly coarse-grained lo- cal density n ( r ) | L ( t ) of mesoscopic restrictions to diffusion (e.g., the disks in Fig. 2). Hence, the variance of D ( r ) | L ( t ) entering Eq. (2.21) is proportional to a typical density fluc- tuation (cid:10) ( δn ) (cid:11) | L ( t ) of the restrictions in a volume of size L d ( t ) in d dimensions (this becomes valid when the devia- tions δn ( r ) | L ( t ) = n ( r ) | L ( t ) − (cid:104) n (cid:105) from the mean sample density (cid:104) n (cid:105) become small). This proportionality, asymptot- ically exact at small k (i.e., after coarse-graining over large distances, cf. Eq. (2.22) for long t ), leads to the proportional- ity Γ D ( k ) ∝ Γ( k ) , k → (2.24) between the correlation functions (power spectra) of D ( r ) and of the underlying structure n ( r ) , Γ( r ) = (cid:104) n ( r + r ) n ( r ) (cid:105) r . (2.25) The structural correlation function can behave qualitatively distinctly at large distances, i.e., small k : Γ( k ) ∼ k p , k → . (2.26) The structural exponent p in Eq. (2.26) defines the structural universality class to which a sample belongs, according to its large-scale structural fluctuations embodied by its correlation function (2.25). The greater the exponent p , the more sup- pressed are the structural fluctuations at large distances (low k ); conversely, negative p signify strong disorder, where the fluctuations are stronger than Poissonian (for which p = 0 ). Hence, p characterizes global structural complexity, taking discrete values robust to local perturbations. This enables the classification of mesoscopic disorder [30], and its relation to the Brownian dynamics, as we now explain. From Eq. (2.22) it directly follows that the time-dependent instantaneous diffusion coefficient (1.2) approaches the finite bulk diffusion constant D ∞ as a power law D inst ( t ) (cid:39) D ∞ + const · t − ϑ , ϑ > , (2.27) with the dynamical exponent [30]: ϑ = ( p + d ) / (2.28) the structural exponent p , and to the spatial dimensionality d . To illustrate the above general relations, consider Pois- sonian disorder (uncorrelated restrictions, e.g., completely randomly placed disks in Fig. 2). Their density fluctuation within the “diffusion volume” L d scales as the inverse vol- ume, (cid:104) ( δn ( r )) (cid:105) ∼ /L d ( t ) ∼ t − d/ according to the central limit theorem. Equivalently, Γ( k ) → const ∼ k as k → , i.e., the exponent p = 0 . As a result, when restrictions are uncorrelated (or, more generally, short-range disordered , i.e., have finite correlation length in their placement), the instanta- neous diffusion coefficient approaches its macroscopic limit as D inst ( t ) (cid:39) D ∞ + const · t − d/ (2.29) in d dimensions, i.e., ϑ = (0 + d ) / . This is the intuitive picture behind the formal results [146, 147]. The approach described in ref. [30] generalizes this picture onto any universality class of structural disorder and enables identifying relevant structural fluctuations by measuring the dynamical exponent (2.28). This exponent manifests itself in the power law tail of the molecular velocity autocorrelation function (2.6) D ( t ) ∼ t − (1+ ϑ ) (2.30) and in the dispersive diffusivity, Eq. (2.8), D ( ω ) (cid:39) D ∞ + const · ( iω ) ϑ (2.31) whose real part is accessible with OGSE, Sec. 2.2.3. Relation (2.28) provides a way to determine the exponent p (or the effective dimensionality d ) and, thereby, the struc- tural universality class, using any type of macroscopic time- dependent diffusion measurement. Local properties, con- tributing to biological variability, affect the non-universal co- efficients, e.g., the values of D ∞ and the prefactor of t − ϑ in Eq. (2.27), but not the exponent ϑ . The latter exponent is universal , i.e., is independent of microscopic details, and is robust with respect to variations between samples of a similar nature. From the point of dMRI in biological tissues, the ex- ponent (2.28) is robust with respect to biological variability . We can also see that the stronger the fluctuations (the smaller the exponent p ), the smaller is the dynamical expo- nent ϑ , i.e., the slower is the approach to D ∞ . Physically, this happens because it takes longer for the coarse-graining to self-average the sample’s structural fluctuations. Conversely, if a sample is regular (a periodic lattice, formally equivalent to p → ∞ ), the approach of D ∞ will happen exponentially fast (i.e., faster than any finite inverse power law) [30]. The above approach exemplifies the power of an effective theory way of thinking, where, to make fairly general state- ments about the relation between the diffusive dynamics and The dispersive terms reads iω ln( − iω ) for the special case of ϑ = 1 ,hence Re D ( ω ) will depend on ω as | ω | , cf. ref. [20]. the structural disorder, we did not have to solve the full non- perturbative problem (starting from the microscopic restric- tions n ( r ) ), but instead ended up solving a relatively simple problem of finding lowest-order corrections [24, 30] to Gaus- sian diffusion in a weakly heterogeneous medium. Note that the undefined constants in Eqs. (2.27) and (2.31) are different. It is possible to find a more precise corre- spondence between the time-dependent terms in D inst ( t ) and D ( ω ) by using Eq. (2.12) followed by contour integration in the complex plane of ω , yielding t − ϑ ←→ π/ ϑ ) sin πϑ · ω ϑ , ϑ < (2.32) (here Γ( ϑ ) is Euler’s Γ -function; cf. also ref. [30], compare Supplementary Eqs. [S17] and [S18].) The above “conversion” between PGSE and OGSE works for D inst ( t ) . However, the cumulative D ( t ) follows the be- havior (2.27) only for ϑ < ; for ϑ ≥ , the PGSE D ( t ) expansion at long t will begin with the /t term, due to the integral in inverting the relation (2.10), D ( t ) = 1 t (cid:90) t d τ D inst ( τ ) (2.33) converging at short t for the t − ϑ term with ϑ > in D inst ( t ) . Therefore, the structure-specific dynamical exponent (2.28) is masked in PGSE if it exceeds unity; to reveal it, one has to use D inst ( t ) , which amounts to differentiating noisy experi- mental data [30, 141]. The borderline case ϑ = 1 has been considered in detail in ref. [20]; the PGSE diffusion coeffi- cient has a (ln t ) /t tail due to the logarithmically divergent integral in Eq. (2.33), D ( t ) (cid:39) D ∞ + A · ln( t/ ˜ t c ) t , t (cid:29) ˜ t c ∼ max { t c , δ } , (2.34) whereas the OGSE counterpart is given by Re D ( ω ) (cid:39) D ∞ + A · π | ω | , | ω | t c (cid:28) . (2.35) Here t c ∼ l c /D ∞ is the time to diffuse across the corre- lation length of the corresponding disordered environment (e.g., correlation length of the disordered axonal packing [20] in the case p = 0 and d = 2 considered below in Sec. 2.4.2). When the pulse duration δ exceeds t c , it starts to play the role of a cutoff time for the power-law tail [20–23]. Finally, we give one more illustration of the absence of any universal relation between PGSE diffusion time and OGSE frequency ω = 2 πν addressed in Sec. 2.2.3. The PGSE- OGSE correspondence, empirically, means that the constants in the tail t − ϑ in D ( t ) , and in the ω ϑ tail of Eq. (2.31) are equal. According to Eqs. (2.32) and (2.33), regime (ii) : t = (cid:34) ϑ ) sin πϑ π (1 − ϑ ) (cid:35) /ϑ · πν , ϑ < . (2.36) This relation is neither obvious, nor has it anything to do with the empirical Eq. (2.19). For example, for ϑ = 1 / (random ( ω / π ) / , kHz / D ( ω ) , µ m / m s rat, normalmouse, normalrat, injurymouse, injury FIG. 7. OGSE measurements in cortical GM: circles are data fromaverage of 5 rats [122] and squares from 6 neonatal mice at 24 hoursafter unilateral hypoxic ischemic injury [125]. Red: normal rat brainand contralateral side of mouse brain. Blue: globally ischemic ratand ipsilateral side of hypoxia-ischemia injured mouse brain. PGSEdata not shown. Dashed lines are fits from Fig. 4 of ref. [30], dot-ted lines are ω / fits (shown as guide to the eye; power-law ex-ponent fit for mouse data was not robust due to narrow frequencyrange). Note that while the absolute D ( ω ) values differ between ratand mouse, the general features are similar: data is well describedwith ω / behavior for normal and ischemic GM (except, possibly,the ischemic mouse, where major structural changes may have oc-curred in 24h); and the coefficient in front of √ ω (the slope) in-creases in ischemia, consistent with short-range structural disorderincrease along the neurites (e.g., due to beading). permeable barrier model [50] or short-range disorder in d = 1 relevant for the neurites, Sec. 2.4.2 below), we obtain t = (4 /π ) /ν . And besides, Eq. (2.36) only applies for ϑ < ; for media characterized by larger ϑ , the power-law PGSE and OGSE tails will not match. Probing the diffusivity time-dependence at long times poten- tially allows for identifying the disorder universality class of the mesoscopic structure, which then helps to parsimoniously model the relevant features of tissue architecture and extract the corresponding tissue length scale(s) and other parameters. Validation with Monte Carlo simulations:
Equation (2.28) was verified in d = 1 MC simulations using different placements of identical permeable barriers, according to peri- odic, short-range, hyperuniform, and strong disorder classes [30]. The same relation for the p = − universality class of random permeable membranes in d = 2 was verified in ref. [50]. The borderline “ log ” case of ϑ = 1 for p = 0 in d = 2 dimensions, Eq. (2.34), was verified in ref. [20] in the t -domain. In ref. [148], the same universality class was con- sidered in the ω -domain, verifying Eq. (2.35) and revealing the dependence of the prefactor A on the degree of disorder in fiber packing. Subsequently, the scaling (2.28) was verified with MC along synthetic model neurites, dimension d = 1 , featuring realistic spines, leaflets and beads placed randomly according to the Poissonian statistics, p = 0 [149]. The slight deviation from the ϑ = 1 / power law at the longest t is at- tributed to the periodic boundary conditions for a relatively short sample. In the same setting, Eq. (2.17) transverse to the neurites was verified for the sub-ms times, regime (i) . A more empirical approach [150] was employed for diffu- sion of cell-specific metabolites up to t = 2 s that was mea- sured by diffusion-weighted MR spectroscopy in vivo. Due to the broad time range, the D ( t ) -dependence was more pro- nounced, in comparison to the earlier measurements [151] for the narrower t -ranges. Distinct tissue morphologies were rec- ognized by comparing with large-scale MC simulations for particles diffusing in many synthetic cells generated as tree- like structures, by varying statistics of the number of pro- cesses, branches, and segment lengths. While the simula- tions matched the measurements, the relevant structural de- grees of freedom and the associated functional forms of the t -dependence were not unequivocally identified. Validation in phantoms:
The exponent (2.28) corre- sponding to the short-time disorder, p = 0 , has been demon- strated in d = 2 dimensions in an anisotropic fiber phantom mimicking the extra-axonal space [20] ( ϑ = 1 , leading to the ln t singularity in the PGSE D ( t ) , Eq. (2.34)), as well as more recently in a d = 1 -dimensional phantom [141] for which ϑ = 1 / for p = 0 and ϑ = 3 / for hyperuniform placement of permeable membranes ( p = 2 ). Cortical GM was probed with OGSE in rat [122] and mouse [125, 126], Fig. 7. It was suggested [30], that the ap- parent ω ϑ behavior with ϑ = 1 / can be explained by the dominance of the effectively d = 1 -dimensional diffusion along the narrow neurites with the short-range disorder (e.g., spines, beads, varicosities) along them [152–154]. The vari- cosities are known [155, 156] to become more pronounced in ischemia, which is consistent with the increase in the struc- tural disorder-induced coefficient in front of the ω -dependent term in Eq. (2.31). Conversely, assuming short-range disorder (e.g., in the varicosity placement), the power law exponent ϑ = 1 / then validates the effectively d = 1 -dimensional diffusion along the so-called “sticks” (narrow channels used to model the intra-neurite space), cf. Sec. 3.1 below. Human WM was probed in vivo in both the longitu- dinal [31] and transverse directions [31, 121]. The time- dependence of longitudinal diffusivity suggests restrictions are present along axons, which, similar to the GM case above, augments the commonly used “hollow tube” model for dif- fusion inside and outside neurites (cf. Sec. 3.1). The “hol- low tube” filled with some effective Gaussian medium with diffusion coefficient D ∞ becomes an effective theory tech- nically valid only in the t → ∞ regime (iii) ; for finite t , D ( t ) and higher-order cumulants) will be present. Recent quantitative analysis [31] based on Eq. (2.27) for d = 1 , revealed that this time-dependence is compatible with short-range disorder in the placement of restrictions along ax- ons. Intriguingly, the corresponding correlation lengths of about − µ m are similar to those reported in the litera- ture for varicosities along axons [152–154], suggesting them as potential sources for the reduction of the longitudinal dif- fusivity with time. Varicosities are often found to be rich in mitochondria and could therefore form obstacles for the diffu- sion along the fibers, or they could act as temporary traps for the longitudinal diffusion. Additional potential sources of the short range disorder could be axonal undulations [140, 157], or functional gap junctions unevenly spaced between 20 and µ m along the myelin sheath in sciatic nerve [158]. Note an interesting observation that the reduction in the diffusivity in acute stroke patients occurs predominantly along the axons when measured at the frequency ω/ π =
50 Hz , while the decrease in both the longitudinal and trans- verse directions directions is observable for the diffusion time
40 ms [159]. Originally explained by axonal beading, this ef- fect has to be taken into account in more general models of the diffusion response to tissue damage. Quantifying µ m-level structure of neuronal tracts in vivo has been brought to the forefront of neuroscience research pri- marily due to the axonal diameter mapping (ADM) concept, developed within the CHARMED and AxCaliber frameworks ([56, 57, 160]) and their extensions [58, 161]. Their common theme is the focus on the intra-axonal compartment (assum- ing no exchange) as the source of the diameter sensitivity, while typically approximating axons as impermeable cylin- ders, and building on exact solutions [14, 18, 102, 137]. Wa- ter diffusion in the extra-axonal space in all of the above ap- proaches is approximated as Gaussian (time-independent). Large overestimation of axonal diameters, by factors 3-15 in humans, cf., e.g., refs. [58, 162], provoked a debate [163– recognized that ADM may be confounded by two issues. The first ADM issue is the smallness of the signal attenu- ation for typical thin axons. The attenuation inside a cylin- der up to O ( g ) (GPA) is given by van Gelderen’s formula [117] that depends on the PGSE sequence timings ∆ and δ . However, in the practically relevant case δ (cid:29) r /D , the dependence on ∆ drops out, hence the diffusion time is not actually being used to “probe” the cylinder diameter. This is the Neuman’s limit [18], in which attenuation − ln SS (cid:39) g r · δD ≈ . · − (2.37) is proportional to the total time δ the gradients are on. The proportionality to the time δ can be understood in terms of the mapping onto the transverse relaxation in the diffusion-narrowing regime, where ∼ ( gr ) · r /D is the effective R ∗ rate [23]. The above attenuation was eval- uated for typical values of the Larmor frequency gradient g = 0 . µ m · ms) − corresponding to 40 mT/m; free axoplasmic D = 2 . µ m / ms ; pulse duration δ = 10 ms, and a typical inner diameter r = 1 µ m [79–83]. Likewise, the OGSE attenuation − ln SS (cid:39) b · Re D ( ω ) in the relevant ωr /D (cid:28) limit [cf. Eqs. (B5) and (B7) in Appendix B] − ln SS = 796 ( g / r · TD (2.38) becomes independent of the OGSE frequency ω and of the OGSE initial phase φ , since b ∼ T /ω and Re D ( ω ) ∼ ω . The analogy with Eq. (2.37) becomes obvious if we realize that, following the mapping onto the transverse relaxation in the presence of an (oscillating) gradient g = g cos( ωt − φ ) , what matters is the total time δ → T the gradients are on (here T = N · π/ω is the total OGSE train duration), and the time-averaged gradient power g → (cid:10) g ( t ) (cid:11) = g / . In other words, in the low-frequency limit, OGSE is just the Neuman’s limit (2.37) albeit with the reduced average gradi- ent amplitude, since the gradients are not at their peak value all the time. This yields that OGSE is not beneficial to map small compartment sizes, and does not provide any indepen- dent parameter combination, in the limit ωr /D (cid:28) (cf. Appendix B). Obviously, the most optimal setting is to keep the diffusion gradient at its maximum all the time δ (cid:46) T , cf. Eq. (2.37). Practical resolution limits for axonal radii de- pending on the SNR and fiber geometry were considered for both pulse-gradient and OGSE sequences [174–176]. The second ADM issue is potentially more significant. Had the above smallness been the only problem, we would just not see any dependencies of intra-axonal signal on experimen- tal parameters. Yet the fits of ADM model to data do show definitive trends in the estimated “apparent diameters” — for instance, with the gradient strength [162] — suggesting some unaccounted physical phenomenon. This has prompted tak- ing into account the coarse-graining outside [20] randomly packed axons. The (ln t ) /t term (2.34) from the extra-axonal space appears to completely overwhelm the weak attenuation This value is based on the observation [166] that axoplasmic diffusioncoefficient in squid giant axon is 20% below the water diffusion coeffi-cient at the same temperature, and is consistent with the recent estimate of D a ≈ . − . µ m / ms along axons in human WM at t = 50 ms ob-tained by suppressing extra-axonal compartment using either high b [167]or planar diffusion encoding [140], which sets a lower bound for D .Another large axon study in excised lamprey spinal cord [168] reporteda similar deviation of about 25% for the longitudinal diffusion coefficientfrom the free water diffusion coefficient. A somewhat larger value wasreported in excised pig spinal cord [169]. Alternatively, using NAA asan intracellular reporter molecule, the ratio for the in vivo measured par-allel diffusion coefficient in the corpus callosum relative to its diffusioncoefficient in dilute aqeous solutions, ranges from 0.4 up to 0.46 [33], cor-responding to estimates of water D a ∼ . − . µ m / ms . Here we consider brain; axons are about factor of 5 thicker in the spinalcord, and the ADM prospects are much better there [123, 170–173], dueto the r scaling in Eq. (2.37); see, however, the need for beyond-GPAcorrections discussed below. D ( t ) measure- ments transverse to human WM fiber tracts [31, 121]. This reveals an exciting unexpected mesoscopic effect: The struc- tural disorder in axonal packing within a WM fiber bundle completely changes the interpretation of ADM at low to mod- erate diffusion weightings. Furthermore, recently observed logarithmic dependence on the pulse duration δ in human WM [23], cf. Eq. (2.34), instead of the linear one, Eq. (2.37), and validated in a fiber phantom [22], confirms the meso- scopic extra-axonal origin of the “apparent” ADM effects. The decreasing apparent diameter trend with increasing gradient strength [162] is consistent with eventual suppres- sion, as ∼ e − b ( D ∞ + A ln t/t ) , of the extra-axonal contribu- tion; however, since the transverse D ∞ (cid:46) . µ m / ms (cf. Section 3), very large b -values are needed to fully suppress this effect [177]. However, when sufficiently strong gra- dients are used in animal settings, the GPA results of van Gelderen [117] and Neuman [18] should be corrected. Un- fortunately, no analytical solution exists beyond GPA for fi- nite δ . Lee et al. [23] estimated that the GPA will break down when g (cid:38) g ∗ = D /r for axons of radius r . This may become relevant for large axons: e.g., for r = 3 µ m and D = 2 µ m / ms , the critical gradient g ∗ corresponds to
277 mT / m . Furthermore, the next-order O ( g ) correction to the right-hand side of Eq. (2.37) will be of the same sign as the main effect, scaling as ∼ g r δ/D , implying that strong gradients cause extra attenuation relative to what GPA predicts. If the GPA is used instead of the exact solution, the GPA-derived radii will be overestimated, which may explain some residual overestimation of axonal radii in a recent ani- mal study with gradients as large as . / m [177]. We note that an even stronger dominance of the extra- axonal contribution occurs in the OGSE domain, since the fully confined water (within an impermeable cylinder trans- verse to its axis) yields a regular, ω contribution to Re D ( ω ) (Appendix B), whereas the extra-axonal water would con- tribute linearly , as | ω | , cf. Eq. (2.31) with ϑ = 1 . The linear term will dominate at low ω , in agreement with the linear dis- persion observed transverse to fibers with OGSE by Portnoy et al. [178] and analyzed in ref. [20]. Such linearity has been also observed in rat spinal cord by Xu et al. [123, 171]. The predominance of the | ω | scaling means that using OGSE is not optimal for probing inner axonal diameters — not just numerically as discussed after Eq. (2.38), but para- metrically! However, the | ω | scaling makes OGSE paramet- rically better for probing the extra-axonal space geometry. We warn that the gradient waveform optimization which does account for the mesoscopic effects in the extra-axonal space may sometimes give unfair preference to OGSE [174–176]. Overall, the above mesoscopic effects, measured in vivo on both animal and human scanners, may enable a novel kind of structural contrast at the micrometer scale (e.g., axonal loss and demyelination), and open up exciting possibilities of monitoring subtle changes of structural arrangements within GM and neuronal tracts in disease, aging, and development. When ADM is feasible (e.g., spinal cord, due to much larger r ), Neuman’s r scaling (2.37), together with volume- weighting ∼ r , gives a large weight to a small number of axons with largest diameters, effectively measuring [20] r Neuman = (cid:16) r /r (cid:17) / , (2.39) where the averages are taken over the voxel-wise axonal dis- tribution. Hence, the metric (2.39) may become susceptible to the mesoscopic fluctuations governed by the tail of the dis- tribution (practically, for sufficiently small voxels in which such fluctuations can be pronounced). Additionally, sam- pling fluctuations confound the comparison of dMRI mea- surements with histology, — where the metric derived based on, e.g., Eq. (2.39) can be strongly sample-dependent, espe- cially if small fields of view are utilized. This general phenomenon of rare structural configurations determining the measurement outcome has parallels with similar effects found in hopping conduction in disordered semiconductors, kinetics of reaction-diffusion systems, and other phenomena in disordered media [91, 179, 180]. In our case, an incidentally large number of thick axons may sig- nificantly skew the intra-axonal attenuation for a particular voxel. This could lead to strongly enhanced variations (rela- tive to those expected based on the measurement noise alone) in the corresponding parametric maps. The issue of the mesoscopic flucutations is fundamental, and the separation of the effects of biological variabiilty from the randomness in measurement outcomes due to the thermal noise requires model-independent ways of estimating local noise level [181, 182], as well as precisely quantifying the tails of the corresponding distributions of biophysical tissue parameters (e.g., of the axonal diameter distribution [79–83]).
3. THE t → ∞ LIMIT, REGIME (III): MULTIPLE GAUSSIAN COMPARTMENTS All science is either physics or stamp collectingErnest Rutherford The flamboyant century-old quote of a founder of the atomic age could be excusable, as scientific disciplines other than physics in his days were mostly collecting empirical informa- tion. Today, with so much more knowledge about the world and the associated abundance of data, Rutherford’s quote could as well sound “All science is either physics or fitting”. While the purpose of physics remains to identify relevant pa- rameters and to produce an explanation (an effective theory), and its instance for a particular measurement — a model [69], the complexity of models and the amount of data have turned parameter estimation into a field on its own, if not into a mul- titude of fields, employing a wealth of approaches, known under different names and incorporating advanced tools of statistics, machine learning and artificial intelligence. Within MRI, modern parameter estimation approaches are tied to the k -space data, which spurred the applications of compressed sensing [183, 184] and MR fingerprinting [185]. Tissue microstructure mapping presents its own set of pa- rameter estimation challenges. As we will illustrate in this Section, while from the physics standpoint, the dMRI mod- els in the t → ∞ regime become trivial (a sum of Gaus- sians = exponentials in b ), their number of parameters, and the inherent degeneracy of the fitting landscape in face of the typically low SNR of dMRI acquisitions, have turned param- eter estimation into an active area of investigation. In other words, the problem remains largely unsolved — even with a densely sampled q -space, and fully sampled k -space. So far, arguably, most intellectual efforts in the regime (iii) (as de- fined in Sec. 1.5) have been spent on the “fitting” rather than on the “physics”. This Section is hence primarily about the parameter estimation aspect of modeling (cf. Sec. 1.1). Below, after introducing the stick compartment in Sec. 3.1, we formulate the overarching Standard Model of diffusion in neuronal tissue as a sum of anisotropic Gaussian compart- ments (Sec. 3.2, Figs. 4 and 8), and then discuss challenges of its parameter estimation, Sec, 3.3, focussing on its degen- eracies. We subsequently review works involving constraints on the Standard Model parameters (Sec. 3.4), followed by the unconstrained, rotationally-invariant methods (Sec. 3.5), and conclude this Section with a summary of unresolved issues (Sec. 3.6). In this Section, we assume that the t → ∞ regime (iii) has been practically achieved, and neglect the time-dependent power-law “tails” describing the approach of the diffusion co- efficient to its tortuosity limit, discussed in Section 2. Full coarse-graining in the intra-neurite space then leads to the most anisotropic Gaussian compartment possible — the so-called “stick” compartment — first introduced by Kroenke et al. [33] and Jespersen et al. [36] in 2004–2007. Its main features are: A stick is a cylinder whose radius r is negligible com- pared with the “free” diffusion length ∼ √ D t at given t . Equivalently, the transverse diffusion coefficient inside neurites D ⊥ a (cid:39) r / t ∼ . µ m / ms for typical t ∼ ms is negligible compared to D ≈ . µ m / ms , and hence can be set to zero, D ⊥ a → . In other words, the measurement is insensitive to neu- rite radii (cf. discussion in Sec. 2.4.2). Mathematically speaking, a power-law approach, being scale-invariant,means that the t → ∞ regime is never fully achieved — there is no timescale that tells us where we can neglect the residual non-Gaussian effectsin each compartment. However, practically, their detection limit is set bya finite SNR. See footnote 15 in Section 2 The longitudinal diffusion inside a neurite becomes Gaussian, with the macroscopic (tortuosity) asymptote D a . Of course, D a , being the effective coarse-grained parameter, can be notably reduced relatively to the in- trinsic axoplasmic diffusion coefficient D , cf. Sec- tion 2. The parameter D a takes into account all re- strictions, such as varicosities (beads) and undulations, along the (average) neurite direction. Hence, it can have important biophysical and diagnostic value in the cases when the structure along neurites changes, e.g., in acute stroke [155, 156, 159] and in Alzheimer’s dis- ease [186]. Exchange between intra- and extra-neurite water can be neglected, at least at the time scales t used in clin- ical dMRI. Measuring exchange times in vivo is very difficult, making this assumption hard to validate. The consensus so far has been that this assumption holds for WM tracts, where sticks represent (myelinated) ax- ons (and possibly some glial processes). The filter- exchange study in [187] supports the presence of two compartments in healthy brain WM with the exchange time of about s, consistent with assuming negligible exchange on the clinical t ∼ ms time scale. At which t this assumption might break for glial cells, den- drites in GM, or for unmyelinated axons, is a subject of investigation. It was hypothesized [188] that tran- scytolemmal water exchange in astrocytes is fast, since inhibition of aquaporin-4 significantly reduced the dif- fusion coefficient already for t < ms, without mod- ifying tissue histology. However, recent work [189] of measuring T in the presence of fast extra-cellular flow for cultures of astrocytes and neurons grown on beads puts the intracellular residence time around and ms, correspondingly. Likewise, MR relaxation measurement [190] in the live rat brain organotypic cortical cultures yields the net cellular water efflux rate .
02 s − , with a significant fraction ( ∼ − %) of this exchange rate attributed to active transcytolemmal exchange related to the Na + -K + -ATPase activity. From the modeling standpoint, the stick compartment is the defining feature of dMRI inherent to the neuronal tissue, as compared to all other kinds of soft tissues. It is chiefly responsible for the anisotropy of the diffusion propagator in the brain (at least, in the white matter), and in spinal cord (where finite axonal radii can be detected, see footnote ?? ). The diffusion propagator for a stick pointing in the unit direction ˆn , measured in the unit gradient direction ˆg , G ˆn ( ˆg , b ) = e − bD a ( ˆg · ˆn ) (3.1) is determined by cos θ ≡ ˆg · ˆn , where θ is the angle be- tween ˆg and ˆn . The signal is not suppressed for ˆg ⊥ ˆn and decays fast with b when ˆg (cid:107) ˆn . Hence, when bD a (cid:29) , the stick dMRI response (3.1) becomes a thin “pancake”, non-negligible when | ˆn · ˆg | (cid:46) ( bD a ) − / nearly transverse to ˆg , whose angular thickness scales as δθ ∼ / √ bD a [167, 191, 192]. Both estimates follow from setting the ar- gument of the exponential to unity. n D ? e f ; D k e ≈ " voxel" = " n fiber""segment" P ( n ) n g θ f ; D a , D k e , D ? e FIG. 8.
The Standard Model of diffusion inneuronal tissue, Eq. (3.4).
In the t → ∞ regime (iii) , elementary fiber segments (fasci-cles), consisting of intra- and extra-neurite com-partments, are described by at least 4 indepen-dent parameters: f , D a , D (cid:107) e and D ⊥ e . CSF canbe further added as the third compartment, cf.Eq. (3.5). Within a macroscopic voxel, suchsegments contribute to the directional dMRI sig-nal according to their ODF P ( ˆn ) . An extensive review of dMRI validation studies is beyond the scope of this article. However, given the essential role sticks play in dMRI models, we mention the following two kinds of results. First, metabolites, such as N-Acetylaspartic acid (NAA), intrinsic to the intra-neurite space, can be used to identify the stick compartment, as reviewed by Ackerman and Neil [151]. A seminal NAA study was performed by Kroenke, Ackerman and Yablonskiy [33], who demonstrated a very good agree- ment between the dMRI signal from large voxels in rat brain averaged over three gradient orientations, at diffusion times t = 50 − ms, and the isotropically averaged stick signal, (cid:90) d ˆn G ˆn ( ˆg , b ) = (cid:114) π bD a erf (cid:16)(cid:112) bD a (cid:17) , (3.2) where erf is the error function. Taking a large voxel, which presumably has all neurite orientations, and subsequently av- eraging over 3 directions, maps the signal to that from a completely random stick arrangement, first considered by Callaghan [53] in 1979, and subsequently by Yablonskiy et al. [54] in 2002 for He diffusion in the lung, resulting in Eq. (3.2) (Fig. 4). The agreement with the theory (3.2) was very good in the whole range < b (cid:46)
20 ms /µ m [33]. Re- cent directional NAA imaging by Ronen et al. [193], quanti- fying the dMRI signal anisotropy, agrees well with the struc- ture tensor [194] derived from the axonal histology in the cor- pus callosum. Second, for the water dMRI, identifying a distinct func- tional form inherent to the stick compartment can also val- idate the pictures of sticks at sufficiently large b , when the extra-neurite signal becomes exponentially suppressed (be- cause the extra-neurite diffusion coefficient is nonzero in any direction), while the intra-neurite signal is only suppressed as a very slowly decaying power-law ∼ b − / , scaling as the width of the pancake-shaped stick response function (3.1), cf. Eq. (3.2) and refs. [167, 191, 192]. The recently observed [167, 192] power-law water signal attenuation in human WM in vivo , isotropically averaged over multiple gradient direc- tions ˆg , S | bD a (cid:29) (cid:39) β · b − α + γ , (3.3) with exponent α = 1 / , provides a unique signature of axons as sticks for water dMRI. The isotropic average (cf. Sec. 3.5.1 below) of the signal makes it equivalent to that from isotropic set of sticks, cf. Eq. (3.2) above, with erf approaching at large b . Either detectable axonal diameter values, or a notable ex- change rate between intra- and extra-axonal water, would de- stroy the very particular stick-related b − / scaling (3.3). In particular, the analysis in ref. [167] shows that, had the ax- onal diameters been notably higher than their histological es- timates of ∼ µ m [79, 80, 82], as the ADM results often show (cf. Sec. 2.4.2), the power law exponent α would differ from / for b ≤
10 ms /µ m . Hence, human dMRI mea- surements are practically insensitive to axonal diameters even with gradients of
80 mT / m employed, confirming the first stick assumption in Sec. 3.1. The human measurement [167] also revealed that the immobile water fraction, not decaying with b for any direction, is below detection limit. The same conclusion was made independently using isotropic diffusion weighting [195]. Small and slightly negative γ , of about ∼ , is a signature of the breakdown of zero-radius stick picture relevant at very high b [167]. A low, but measurable transverse intra-axonal diffusivity, even in vanishingly thin axons, can emerge from deviations of their form from perfectly straight sticks [140, 157]. Measure- ment with suppression of extra-axonal signal in the human brain suggests D ⊥ a = 0 . ± . µ m / ms for the diffu- sion time about
120 ms , which was obtained as the difference between the trace of intra-axonal diffusion tensor and the lon- gitudinal diffusivity, D a [140]. With the intra-neurite diffusion modeled as a collection of sticks, and the isotropically fully restricted water out of the picture, how should we model the remaining water? The an- swer depends on the coarse-graining length scale. If the diffusion time is as large as needed for water molecules to sample a statistically representative part of the extra-neurite space (ENS) within a voxel, then diffusion in this space should become Gaussian and be described by the overall ENS diffusion tensor, S ENS ∼ exp( − b ij D ENS ij ) , where the b -matrix b ij = q i q j t ≡ b g i g j depends on the com- ponents of the unit diffusion gradient ˆg , cf. Eqs. (1.8) and (1.10). The ENS tensor D ENS would then by definition de- scribe all ENS water, including cerebrospinal fluid (CSF) that could, e.g., contribute if a voxel contains part of a ventricle. a macroscopic voxel. For typical diffusion times t = 50 − ms, the corresponding diffusion length L ( t ) ∼ µ m de- fines the coarse-graining window, where the diffusion prop- erties locally become (almost) Gaussian. At this scale, the neuronal tissue, at least in the WM, looks as a highly aligned fiber segment (fascicle) , Fig. 8, the leftmost image. Because the fibers are locally coherent at the scale L ( t ) , the local ENS tensor D ENS ij ( r ) can differ at the different positions r within a voxel if they are separated by distance much exceeding L ( t ) (for instance, if a voxel contains fiber crossings). This local microscopic anisotropy of ENS suggests that, strictly speak- ing, the ENS diffusion is never Gaussian — but is rather a sum of local anisotropic Gaussian contributions, which are highly aligned with the corresponding local stick arrangements. These coarse-graining considerations lead to the general picture of anisotropic compartments, Fig. 8. The signal mea- sured in the unit direction ˆg , is a convolution S ˆg ( b ) = (cid:90) | ˆn | =1 d ˆn P ( ˆn ) K ( b, ˆg · ˆn ) (3.4) between the fiber orientation distribution function (ODF) P ( ˆn ) normalized to (cid:82) d ˆn P ( ˆn ) ≡ , and the response kernel K from a perfectly aligned fiber segment (fascicle) pointing in the direction ˆn . The kernel K ( b, ˆg · ˆn ) depends on the rel- ative angle θ , cos θ ≡ ˆg · ˆn (cf. Eq. (3.1). The general representation (3.4) gave rise to a number of methods for de- convolving the fiber ODF from the dMRI signal for a given | q | = q shell in q -space, using different empirical forms of the kernel [197–203]. Following the coarse-graining arguments above, the ker- nel’s functional form K ( b, ξ ) = S (cid:104) f e − bD a ξ + (1 − f − f CSF ) e − bD ⊥ e − b ( D (cid:107) e − D ⊥ e ) ξ + f CSF e − bD CSF (cid:105) , ξ = ˆg · ˆn , (3.5) is a sum of exponential (in b ) contributions from the aligned intra- neurite and extra-neurite spaces, respectively modeled by a stick compartment, by the axially symmetric Gaussian compartment with transverse and longitudinal diffusivities D ⊥ e and D (cid:107) e and the principal direction along the stick, and the CSF compartment, cf. refs. [44–46] and Fig. 8. The myelin water compartment is typically neglected due to its short T time [204] as compared to the echo times T E employed in clinical dMRI. We emphasize that the fractions f and − f are the relative signal fractions , and not the abso- lute water volume fractions, due to neglecting myelin water, cf. footnote 17; effectively, we neglect D ( t ) − D ∞ within each fascicle. This picture might be not as definite in GM, where the dendrites are moreintertwined, and the concept of the overall ENS tensor may be better justi-fied [36, 37]. But then, much less is known about diffusion in GM overall,and the precise knowledge of in vivo neurite residence times is still lack-ing. Therefore, here we mostly discuss SM in the context of WM, leavingGM modeling as one of the unresolved problems in Sec. 5.2. The convolution is on a unit sphere | ˆn | = 1 [196]. We normalize [46]d ˆn ≡ sin θ d θ d φ π such that (cid:82) d ˆn · ≡ , while K| b =0 = S ≡ S | b =0 . as well as due to generally different T values for the intra- and extra-neurite compartments [47, 205]. The isotropic CSF compartment has D CSF ≈ µ m / ms . Because of the ODF normalization (cid:82) d ˆn P ( ˆn ) ≡ , the CSF term can be included in the kernel or added separately to signal (3.4); we choose to include it in the kernel (3.5), as it makes the formulation (3.4) more elegant. In what follows, we will describe in depth the two-compartment kernel (3.5) with f CSF ≡ . The overarching model (3.4) and (3.5) includes nearly all previously used models [33, 36–47] (at least in the WM) as particular cases, also described in taxonomy studies in refs. [206, 207] where some of the “models” are actually rep- resentations, in the sense of Sec. 1.7. In other words, the numerous models (and acronyms), corresponding to the right part of Fig. 4 in Section 1, describe the same physics, and hence they are really the same model . Because of the overall popularity and inclusiveness of the above picture, here we suggest to call the model (3.4)–(3.5) the Standard Model (SM) . We note that a major limitation of the SM kernel (3.5) is using the same scalar parameter values for different fiber tracts passing through a voxel (noted, e.g., in refs. [208, 209]), which prompted assigning different (albeit constant) fiber re- sponses to different tracts to deconvolve the ODF [210–212], as an alternative. Over the past decade, it has become clear that the scalar parameters f , D a , D (cid:107) e and D ⊥ e , and the spherical tensor pa- rameters (the spherical harmonics coefficients of the ODF P ( ˆn ) ), carry distinct biophysical significance. Deconvolving the voxel-wise fiber ODF, instead of relying on the empirical directions obtained by Fourier-transforming the dMRI signal from the q to r space, in principle provides a more adequate starting point for fiber tractography, an essential tool for map- ping structural brain connectivity and for presurgical planning [213–217]. Furthermore, as illustrated further in Sec. 3.4.1, the ability to estimate scalar parameters of the kernel (3.5) would make dMRI measurements specific — rather than just sensitive — to µ m-level manifestations of disease processes, such as demyelination [218–220] ( D ⊥ e ), axonal/dendritic loss [220– f ), beading [224], inflammation and oedema ( f CSF , as well as, potentially, mostly D a for cytotoxic and mostly D (cid:107) e , D ⊥ e for vasogenic oedema [225]). Combining f with the extra-axonal volume fraction derived either from tortu- osity modeling based on D ⊥ e [218, 219] or from the myelin The name is suggested by the tongue-in-cheek analogies with the Stan-dard Model in particle physics. In both communities, SM represents theconsensus knowledge about the subject, satisfactorily describes (almost)everything, has been out there for a while, and yet one really hopes thatthere is exciting physics beyond it — which is far more difficult to access.Our community has a doubtless advantage in that investigations beyondSM are much cheaper (cf. Section 2) than building particle accelerators. to determine axonal g -ratio [226]. Since the precise nature and pathological changes in microarchitecture of restrictions leading to the scalar parameter values are unknown, ideally, to become specific to pathology, one needs to estimate f , f CSF , D a , D (cid:107) e and D ⊥ e separately. The first attempt to estimate the parameters and to vali- date the 2-compartment SM was performed by Jespersen et al. [37] using direct fitting to an extensive ex vivo data set covering both WM and GM, while parametrizing the ODF using spherical harmonics Y lm up to l ≤ . (The ENS diffu- sion in this work was described by the overall tensor, whose 6 components were estimated.) Data was acquired on a 16.4 T magnet using 16 shells with ≤ b ≤
15 ms /µ m and 144 di- rections, with acquisition taking over 15 hours. The fraction f correlated well with AMG staining for the neurites, and the ODF directionality agreed with the histology. Further quanti- tative comparison of the predictions to histology was carried out in [227], where the neurite ODF was determined from Golgi stained cortical neurons in immature ferret brains. While “brute-force” fitting of ∼ parameters could work for an extensive data set [37], clinical dMRI data is far noisier, with much less q -space coverage. Hence, because of the high dimensionality of parameter space and the unfavorable fitting landscape [43], SM parameter estimation for Eqs. (3.4)–(3.5) from realistic noisy clinical dMRI data has emerged as an overarching challenge (Sec. 3.3), which has until recently been addressed by introducing parameter constraints, as dis- cussed further in Sec. 3.4. The open challenge of parameter estimation also means that the literature [33, 36–47] differs largely in the ways SM parameters have been constrained. To quantify the problem’s complexity, we find here how many parameters N p the model (3.4) should have — and, hence, how many we have to estimate from (noisy) dMRI data. We follow here the treatment in ref. [46]. The answer depends on the maximal power l max of the dif- fusion weighting b l max / ∼ q l max to which an acquisition is sensitive, at a given SNR (by the time-reversal symmetry of diffusion, only even l max are considered). This can be seen either from the cumulant expansion (1.10) of ln S ˆg ( b ) , or, equivalently [15], from the Taylor expansion of the signal (3.4) S ˆg ( b ) S (0) = 1 − bM (2) i i g i g i + b M (4) i ...i g i . . . g i − . . . (3.6) in the fully symmetric moments M ( l ) i ...i l . These moments are proportional to angular averages (cid:104) n i . . . n i l (cid:105) over the ODF P ( ˆn ) , as it is evident from expanding the exponential terms containing ξ = n i g i in kernel (3.5), such that subsequent terms have the form b (cid:104) n i n j (cid:105) g i g j , b (cid:104) n i . . . n i (cid:105) g i . . . g i , etc. The maximal (even) order l of the product (cid:104) n i . . . n i l (cid:105) always appears with the corresponding power b l/ of the dif- fusion weighting. The symmetric tensors (cid:104) n i . . . n i l (cid:105) , after subtracting all possible traces, can be turned into the corresponding symmet- ric trace-free tensors (STF) of rank l [77], which are equiv- alent to the set Y lm of the spherical harmonics (SH) dis- cussed in Sec. 1.8 above. In other words, the ODF aver- ages (cid:104) n i . . . n i l (cid:105) correspond to the SH coefficients p lm of the ODF, P ( ˆn ) ≈ l max (cid:88) l =2 , ,... l (cid:88) m = − l p lm Y lm ( ˆn ) . (3.7) In particular, the highest-order cumulant C ( l max ) or the mo- ment M ( l max ) still practically resolvable from the signal, sets the maximal order l max for the even-order SH expansion (3.7). The correspondence between l max in Eq. (3.7) and the maximal order b l max / in expansion (3.6) embodies the per- turbative radial-angular connection in the q -space [46]. We note here an obvious corollary from the radial-angular connection in q -space: Oversampling the directions within the low- b shells does not improve angular resolution in es- timating P ( ˆn ) — in other words, optimal q -space coverage should match the sensitivity to the power b l max / ∼ q l max of the shell radius with the minimal number n c ( l max ) of di- rections per shell. Naive sampling, say, directions at b ≈ /µ m would not in principle yield better angular ODF resolution than, say, ∼ averages of directions. In- deed, the clinical dMRI signal at this b -value can be fully de- scribed using O ( b ) (DTI, l max = 2 ), or, at best, O ( b ) (DKI, l max = 4 ) cumulant terms, corresponding to being sensitive to the ODF expansion coefficients p lm up to l = 2 or l = 4 (containing 5 or 14 parameters). There is no way to deter- mine, say, p m and beyond, if the diffusion weighting is too weak for the b terms to be discernible at a given SNR. Coming back to counting the SM parameters, the (min- imum) N s = 4 scalar parameters from the kernel (3.5) in the absence of CSF (or N s = 5 if the CSF compart- ment is added), are complemented by the n c ( l max ) − l max ( l max + 3) / tensor parameters p lm , where n c ( l ) is the number of the even-order spherical harmonics coefficients up to the order l given by Eq. (1.14) in Sec. 1.8, and we sub- tracted one parameter because p ≡ √ π is set by the ODF normalization. This yields the total SM parameter count N p ( l max ) = N s + l max ( l max + 3) / (3.8) such that N p = 9 , , , , . . . for l max = 2 , , , , . . . already for the two-compartment kernel (3.5), without includ- ing S (0) and f CSF in the count. Equation (3.8) reveals that the model complexity grows fast, as l , if we are to account for the rich orientational content of realistic fiber ODFs in the brain. For an achievable l max ∼ − , the dMRI signal in principle “contains” a few dozen parameters, none of which are known a priori . We can now contrast the SM parameter count (3.8) with the number N c ( l max ) [Eq. (1.15) in Sec. 1.8] of independent pa- rameters “contained” in the cumulant or moment series trun- cated at l = l max . Since N c ( l ) ∼ l grows faster with l than N p ( l ) ∼ l , the moments/cumulants estimated from the signal should, beyond some l , contain more than enough in- formation to determine all the corresponding SM parameters. Direct comparison of Eqs. (1.15) and (3.8) yields that N c ( l max ) > N p ( l max ) for l max ≥ (3.9) for both 2- and 3-compartmental SM. This naive counting would let one believe that, having mastered sufficiently pre- cise DKI parameter estimation ( l max = 4 ), we would be able to find all the scalar SM parameters, as well as estimate arbi- trary fiber ODF up to p m , Eq. (3.7). This intuition, however, is deceptive, as the information content is not evenly distributed among all the N c ( l max ) com- ponents . It turns out that the minimal order l max for which the moments/cumulants contain enough information to deter- mine all SM parameters is l max = 6 , while at the l max = 4 level, there exists a one-dimensional manifold (for N s = 4 ), which can look as a single curve, or as two disjoint continu- ous “branches”, or families, of scalar parameters, which ex- actly match the signal’s moment tensors M (2) i i and M (4) i ...i (equivalently, the diffusion and kurtosis tensors) [45, 46]. The two families of solutions, or the two parts of the above one- dimensional curve (“bi-modality”) technically emerge as the two branches of a square root in a solution for a quadratic equation. In what follows, we illustrate this effect with a toy model of parallel fibers [38, 39] and then show that increasing the model complexity does not cure the problem. Let us now see how two branches appear as solutions of a quadratic equation involving directional diffusion and kur- tosis values for a very simple ODF of perfectly aligned fibers, for the 2-compartment SM case. Here we follow the cumulant-series DKI approach as in ref. [38]; an equivalent formulation in terms of the moments, cf. Eq. (3.6), can be found in ref. [46]. A similar approach was used to estimate white matter tract integrity (WMTI) metrics from DKI [39] and subsequently adapted [228] to enable their estimation with reduced data requirements using axially symmetric DKI [229]. Staying at the O ( b ) level (DKI), the overall radial and ax- ial components of the diffusion tensor, estimated from an ide- ally measured signal (the left-hand side), correspond to the following combinations of the scalar parameters (the right- hand side): D ⊥ = (1 − f ) D ⊥ e (3.10) D (cid:107) = f D a + (1 − f ) D (cid:107) e , (3.11) and for kurtosis components [38, 39] K ⊥ = 3 f − f , (3.12) K (cid:107) = 3 f (1 − f )( D a − D (cid:107) e ) D (cid:107) . (3.13) “Transverse” parameters f and D ⊥ e are uniquely determined from K ⊥ and D ⊥ : f = K ⊥ K ⊥ + 3 , D ⊥ e = (cid:18) K ⊥ (cid:19) D ⊥ . (3.14) However, there are two possible solutions in the parallel di- rection. The duality arises from choosing [39, 46] the ζ = ± branch of the square root in Eq. (3.13), D a − D (cid:107) e = ζ · (cid:115) K (cid:107) f (1 − f ) · D (cid:107) , ζ = ± . (3.15) Here √ K (cid:107) ∝ | D a − D (cid:107) e | ≡ ζ ( D a − D (cid:107) e ) , where ζ = sgn ( D a − D (cid:107) e ) . Note that, since the ground truth is unknown, our choice for the branch ζ may differ from the correct one. From Eqs. (3.11) and (3.15), we find that, not surprisingly, the correct choice of ζ yields the true values D a and D (cid:107) e . However, if the sign choice is wrong, then the “apparent” dif- fusivities do not agree with the true ones: D app a = (2 f − D a + 2(1 − f ) D (cid:107) e , (3.16) D (cid:107) e app = 2 f D a + (1 − f ) D (cid:107) e . (3.17) Note, that in this case, as expected, D (cid:107) e app − D a app = − ( D (cid:107) e − D a ) , i.e., the difference has the same absolute value and the wrong sign. In particular, for f = 1 / , the diffusivi- ties are swapped, — i.e we mistake D (cid:107) e for D a and vice-versa . Yet the above “apparent” values can seem completely bio- physically plausible, especially if f ≈ . . From the above derivation it is evident that the bi-modality of the parameter estimation originates from having two tissue compartments, and that the branch choice is certainly not obvious based on the parameter values estimated at low b . Flat directions in the fitting landscape The simplest nontrivial model revealing general degenera- cies in parameter estimation, NODDIDA (Neurite Orienta- tion Dispersion and Density Imaging with Diffusivities As- sessment) [43], is a two-compartment SM variant ( N s = 4 ) that assumes a 1-parameter Watson ODF shape [41] and sets f CSF ≡ . In this model, unconstrained nonlinear fitting has revealed two families (trenches) of biophysically plausi- ble solutions to fit optimization already in the relatively small, (4+1)-dimensional, parameter space, and flat directions along them [43]. Hence, NODDIDA exemplifies the two-fold nature of the parameter estimation challenge. Beyond the existence of the toy model above), each of them represents a shallow “trench” (a “continuous” degeneracy) in the parameter land- scape of nonlinear fitting. The flatness, or continuous degeneracy, can be formulated as having the number of estimated parameters exceeding the number of relations between the parameters obtainable from the data. In its simplest form, this problem exists already for the simplest single-directional fitting with a biexponential function [74]. A normalized biexponential S = f e − bD + (1 − f ) e − bD has 3 parameters; however, if b is low enough so that we are practically only sensitive to the terms ∼ b and ∼ b — i.e., when the DKI representation works well — we can only estimate 2 combinations of 3 model parameters, and will have a flat direction in the corresponding 3-dimensional fitting parameter landscape. (This problem persists in the presence of exchange, as discussed in ref. [230].) The expansion (3.6) of the 2-compartment SM with any ODF into moments has been analytically and numerically shown to possess a similar kind of degeneracy. This frame- work, called LEMONADE (Linearly Estimated Moments provide Orientations of Neurites And their Diffusivities Ex- actly) [45, 46], exactly relates the moment tensors M ( l ) to SM parameters. It turns out that, at the O ( b ) level, there are only 4 independent equations, which relate rotationally invariant combinations of moments M (2) and M (4) to 5 SM parameters — the 4 scalar ones: f , D a , D (cid:107) e , D ⊥ e , and the ODF invariant p (that characterizes the ODF anisotropy, de- fined in Sec. 3.5 below). Hence, the existence of the flat trenches in nonlinear fitting of NODDIDA is actually com- pletely general; both the discrete bi-modality and the continu- ous trenches follow from the exact relations [45, 46] between the moments and the SM parameters, and will be present for any fiber ODF. Hence, it is only capturing the moment M (6) that can lift both kinds of degeneracies — as we mentioned briefly after Eq. (3.9) above — which is practically quite dif- ficult to become sensitive to. As for the “discrete” degeneracy, the works [38, 39, 43, 45,
46] have collectively raised the fundamental question:
Which “branch” of parameters should be chosen, out of at least two biophysically plausible ones? In ref. [46], the discrete branch choice was formulated in terms of the the ratio β between ground truth compartment diffusivities falling or not falling within the interval − (cid:113) < β < (cid:113) , β = D a − D (cid:107) e D ⊥ e (3.18) determined by the discriminant of a quadratic equation. Note that this condition involves all three compartment diffusiv- ities, rather than just the two axial ones as in Sec. 3.3.3. Of course, only one branch corresponds to the truth; other(s) should be discarded. Obviously, selecting the wrong branch can radically change biophysical and diagnostic implications of the estimated parameters. Yet, branch choice is nontrivial, since often times, both parameter sets look equally biophysi- cally plausible [39, 46], and it is very difficult to have a pre- cise enough idea about the ground truth values, especially of D ⊥ e whose precision crucially affects β . In Sec. 3.6, we will comment on the branch selection. Due to the challenge explained above (Sec. 3.3) of fitting the SM to noisy dMRI data, especially those acquired in clini- cal setting, most attempts of SM parameter estimation so far are based on superimposing additional constraints on both the scalar SM parameters, as well as the fiber ODF, to improve robustness of the fitted SM parameters. An overview of employed models used so far is given in Figure 4. In what follows we consider two representative modeling approaches that have been popular because of their robustness, which potentially allows for clinical translation, and try to explain the quantitative differences in parameter estimates between them in the light of each model’s assump- tions and consequent biases. WMTI, as proposed by Fieremans et al. [38, 39], extracts the 2-compartment SM parameters by relating the DKI com- ponents to scalar parameters in the aligned-fiber framework [38], already explained above (Sec. 3.3.3). Subsequently, the perfectly aligned approximation was somewhat relaxed by al- lowing for some dispersion within the fiber bundle, as de- scribed by an intra-axonal diffusion tensor, while the diffu- sion in the extra-axonal space is still modeled as an overall Gaussian compartment [39]. While the advantage of using DKI eliminated the need for direct nonlinear fitting to the diffusion signal, two different biophysically plausible solu- tions still exist similar to the two branches as described above (Sec. 3.3.3). In WMTI, the branch was chosen as D a < D (cid:107) e based on the available data [39]. Parameter histograms corre- sponding to this choice, yielded f ≈ . , D a ≈ . µ m / ms and D (cid:107) e ≈ . µ m / ms in human corpus callosum. Since no specific model is assumed for the tortuosity D (cid:107) e /D ⊥ e , as D (cid:107) e and D ⊥ e are fitted separately, along with D a and f , it was suggested the WMTI parameters could be used to disentangle between acute damage such as neurite bead- ing, as reflected in D a [232], and chronic damage including different types of demyelination and axonal loss, reflected in changes in the tortuosity, D ⊥ e and f [218, 219]. As an in vivo validation, the age-related changes in the WMTI met- rics were studied during the first two years of healthy brain development [231] (Figure 9), showing significant nonlinear increases in f , and D ⊥ e , related to increased myelination and axonal density, while no changes in the longitudinal compart- ment diffusivities, D a and D (cid:107) e , as expected. Ex vivo animal validation studies provided reasonably accurate estimates of f in a mouse model of hypomyelination [233] and de-, and re- myelination [220, 234]. Furthermore, mouse validation stud- ies demonstrated that D a decreased during acute inflamma- tion, while the axonal water fraction f decreased during the chronic phase of cuprizone intoxication [235], whereby D ⊥ e age [years] age [years] FIG. 9. Comparison of NODDI (Sec. 3.4.2)and WMTI (Sec. 3.4.1) parameter evolutionwith age in human corpus callosum sple-nium [231]. A qualitatively (but not quanti-tatively) similar trend of continued increasein the intra-axonal water fraction f intra ≡ f was observed for both models, consis-tent with on-going myelination. WMTIdisplays a trend of increased fiber align-ment (expressed by the orientation disper-sion (cid:104) cos ψ (cid:105) , derived from the intra-axonaldiffusion tensor ), which could be a manifes-tation of continued pruning in the first yearof life, while NODDI does not. The CSFfraction is set, f iso ≡ f CSF = 0 , in WMTI. and f were found to be respectively more sensitive to global and patchy demyelination [220]. These validation studies suggest increased specificity of the WMTI parameters to mi- crostructural changes as compared to empirical diffusion met- rics. However, while the WMTI metrics correlate as expected with the concentration of (purely intra-axonal) NAA under the assumption D a < D (cid:107) e [236], it should be noted that the measured values of D a with low b dMRI protocols in the range of . − . µ m / ms are significantly lower compared to recently reported values for D a measured in the range of . − . µ m / ms using advanced diffusion protocols provid- ing additional information by varying TE [47], diffusion time [169], or using double diffusion encoding [237], cf. Sec. 4, planar [140], isotropic [238], or very strong multi-directional (linear) [167] diffusion encoding. Further research is war- ranted to understand whether the discrepancy is due to a wrong choice of branch ( D a > D (cid:107) e , cf. Sec. 3.6 below, which would affect primarily estimates of D a and D (cid:107) e , but not so much of f or D ⊥ e ), or, alternatively, due to a potential bias when estimating cumulants dMRI data over a to large b -range [239] (affecting estimates for all model parameters). We also note that branch selection (3.18) for the uncon- strained problem (3.4)–(3.5) is qualitatively similar but quan- titatively different from that in the WMTI highly-aligned tracts case [38, 39], cf. the toy model of Sec. 3.3.3. While qualitatively, the “wrong branch” in both the full model (3.4)– (3.5) and WMTI [39] corresponds, roughly, to swapping of intra- and extra-neurite parameters, there is no exact corre- spondence with the full model that includes dispersion; for instance, f and D ⊥ e are also different between the branches. The difference between WMTI and the full model comes from the fact that in the toy model (WMTI prototype), the perfectly-aligned fiber constraint p = p = 1 has been implemented, together with effectively mixing the LEMON- ADE equations with moments M (4) , m and M (4) , m . There- fore, the branch choice based on sgn ( D a − D (cid:107) e ) is sufficiently different from that of Eq. (3.18). An intermediate case be- tween the two including dispersion was obtained in [43, 169], by constraining the ODF to the Watson distribution, effec- tively mitigating the degeneracy by parameterizing all p l in terms of the Watson distribution concentration parameter κ . (NODDI) NODDI, proposed by Zhang et al. [41], is a 3-compartment SM that assumes a Gaussian-like (Watson) ODF shape char- acterized by one [41] or two [240] parameters. In addition, the three diffusivities are effectively fixed, in the following way: D (cid:107) e = D a D a = 1 . µ m / ms Mean-field tortuosity model [241], D ⊥ e /D (cid:107) e = 1 − f . The estimated parameters are f and f CSF , as well as the ODF parameters (one or two parameters, depending on the Watson [41] or Bingham [240] distribution used). Using high resolution ex vivo imaging, it was recently showed that Bingham-NODDI is able to capture the cortical fibers known to exhibit fanning/bending in human neocortex [242]. It was also shown that NODDI-derived dispersion agrees with histology measures in post-mortem normal and demyelinating lesions in spinal cord samples [243]. Further- more, recent work from Schilling et al. [244] shows a strong overall correlation between the fiber orientation dispersion in- dex (ODI) derived from NODDI versus derived from histol- ogy based on 3D confocal z-stacks in areas to the size of an MRI voxel in adult squirrel monkey brains. However, the same study [244] also showed a small, but systematic overestimation of the true histology-based ODI, as well as a correlation of NODDI-derived ODI with the es- timates f and f CSF . Furthermore, recent extensive human measurements up to b ≤
10 ms /µ m [46] also suggest that the above three parameter constraints generally do not hold, and therefore may bias the estimates of the fractions and fiber dispersion. While both NODDI and WMTI rely on the same over- arching SM, they have different constraints, particularly in ted in WMTI), the ODF (Watson in NODDI, single bun- dle in WMTI), and number of compartments (3 in NODDI, been evaluated by studying changes in the model parameters through normal human early development [231] (Figure 9). In this work, qualitatively similar trends were observed in f , in full agreement with expected on-going myelination, fiber classification and asynchrony of development. The quantita- tive estimates, however, are model-dependent, exhibiting bi- ases and limitations related to the models’ assumptions. Sim- ilarly, changes during the first two years in fiber dispersion in the splenium corpus callosum were qualitatively different be- tween NODDI and WMTl. This illustration clearly calls for extreme caution when interpreting modeling studies based on limited clinical dMRI data, where accuracy is typically sacri- ficed in favor of precision. Indeed, both WMTI and NODDI have made assumptions that allowed for a robust, rather than accurate estimation of the SM model parameters that are not fixed according to each model. This prompts both for im- proved SM parameter estimation methods (discussed next in Sec. 3.5), as well as for “orthogonal” and more comprehen- sive validation methods to gain better understanding of the relevant tissue features of modeling (discussed in Sec. 3.6), prior to applying them to clinical dMRI data. Let us now introduce the recently proposed family of ap- proaches to SM parameter estimation that do not rely on spec- ifying the ODF shape, by factoring it out in a rotationally invariant way. This will enable separation of estimating the scalar and the tensor (ODF) parameters. Of course, all the degeneracies of the parameter estimation will persist — and in fact, factorization has been used as a tool to prove that the above discussed degeneracies are completely general [46]. Much like convolutions become products in the Fourier do- main, the convolution (3.4) between the individual fiber re- sponse K and the ODF P becomes a product in the “spherical Fourier” domain (i.e., the SH basis) [196]: S lm ( b ) = p lm K l ( b ) (3.19) where K l ( b ) is the projection of the kernel K ( b, ξ ) onto the Legendre polynomial ( − l/ P l ( ξ ) [36, 44, 46]. Since any rotation transforms SH components S lm and p lm according to a unitary transformation belonging to the (2 l + − dimensional irreducible representation of SO(3) group labeled by “angular momentum” l , the 2-norms (cid:107) p lm (cid:107) ≡ (cid:113)(cid:80) lm = − l | p lm | and (cid:107) S lm (cid:107) (defined likewise) are conserved under rotations, i.e., are rotational invariants. It is thus conve- nient to introduce rotational invariants p l ≡ (cid:107) p lm (cid:107) / N l and The idea to operate with a single “energy” L norm per each “frequency”band l of SH has been previously applied, e.g., to the problem of shapematching in computer graphics [245] and recently for dMRI data harmo-nization [246]. S l ≡ (cid:107) S lm (cid:107) / N l , where normalization N l = (cid:112) π (2 l + 1) is chosen so that ≤ p l ≤ . Hence, equations (3.19) for the ( l, m ) SH components give rise to the corresponding equa- tions for the rotational invariants [44, 46], S l ( b, x ) = p l K l ( b, x ) , l = 0 , , . . . , (3.20) where we denoted by x the dependence on the kernel’s scalar parameters x = { f, D a , D (cid:107) e , . . . } to be estimated. The in- variant p ≡ is trivial (ODF normalization); the remaining ODF invariants p l , one for each l , characterize its anisotropy irrespective of the chosen basis. l = 0 invariant K ( b ) The l = 0 invariant for Eq. (3.20) has been independently in- troduced as “powder averaging” and “spherical mean” [247– swapping the order of integrations over ˆg and ˆn : S ∝ (cid:90) d ˆg (cid:90) d ˆn P ( ˆn ) K ( b, ˆg · ˆn ) = (cid:90) d ˆn P ( ˆn ) (cid:90) d ˆg K ( b, ˆg · ˆn ) ≡ (cid:90) d cos θ K ( b, cos θ ) , since (cid:82) d ˆg K ( b, ˆg · ˆn ) is independent of fiber direction ˆn due to the “translational invariance” on a unit sphere, and the ODF is normalized to (cid:82) d ˆn P ( ˆn ) ≡ . The last identity above gives the projection of kernel (3.5) onto the l = 0 Legendre poly- nomial P ( ξ ) ≡ , where ξ = cos θ in our case; for a stick compartment, this projection yields Eq. (3.2) above. K l ( b ) for l = 2 , , . . . Equation (3.20) formally yields an infinite family of rota- tional invariants K l ( b ) [44, 46], one for every l = 2 , , . . . . However, it turns out that by far the most useful is the next- order, l = 2 invariant, since the projections of e − bDξ onto the Legendre polynomials with l > , giving the compart- ment contributions to K l ( b, x ) , are too slowly varying [36] and thereby adversely affecting the sensitivity to the esti- mated parameters x . We also note that including the l > invariants in system (3.20) is only possible for anisotropic ODFs, with p l > . Physically, it is expected since the less symmetric the system, the more inequivalent ways it enables for probing it. In the brain, the ODF is at least somewhat anisotropic; its lowest- order invariant p is generally nonzero even in GM. Parameter estimation based on the ODF factorization via the rotational invariants amounts to inverting the nonlinear This intuition underlies theory of quantum-mechanical excitations of non-spherical nuclei [253], where analogs of our rotational invariants arethe corresponding irreducible tensor operators underpinning the Wigner-Eckart theorem. x and p l . Such inversion has so far been technically implemented in four distinct ways: a. Analytically inverting relations between their Taylor expansions — i.e., expressing model parameters in terms of the moments of the signal (LEMONADE) [45]. At typical bD ∼ , the biases in estimating the moments cause notable bias in the model parameters. b. Using the LEMONADE output as initialization for the RotInv solution of Eqs. (3.20) via nonlinear fitting using the gradient-descent optimization of the corre- sponding objective function [46]. This notably in- creases the accuracy of LEMONADE. c. The prevalence method [46]: To avoid the branch se- lection issue, initialize the fit objective function for Eq. (3.20) with a large number ( ∼ − ) of random starting points within the plausible parameter range (e.g., < f, p < , and < D < for all diffu- sivities), observe that the fit outcomes cluster around a few sets in the parameter space, and select the mean of the largest cluster (after excluding outcomes outside the bounds). The method works best for large b , say, b (cid:38) /µ m , since increasing b broadens the basin of attraction of the true minimum [43]. d. Machine learning framework: Generate distributions of the invariants based on the prior distributions of x , and numerically invert these relations based on the train- ing set [44]. The invertibility of these relations re- quires the resolution of the bi-modality problem (sec- tion 3.3). In particular, the constraint of close traces of intra- and extra-axonal tensors, | D (cid:107) e + 2 D ⊥ e − D a | < . µ m / ms , was applied according to results obtained using isotropic diffusion weighting [195]. While it is the fastest data processing method, it is sensitive to the way the training data are generated. Overall, our current experience tells that, no matter the implementation, the sensitivity to different scalar parame- ters varies dramatically (e.g., f is obtained reasonably well while the sensitivity to D (cid:107) e and D ⊥ e is much worse) [44, 46]; with decreasing SNR, methods (a)–(c) yield noisier param- eter maps while method (d) yields “too clean” maps com- pletely dominated by the mean values of parameter priors; branch multi-modality manifests itself in the need for the branch selection (3.18) in all these approaches. It is yet difficult to evaluate accuracy of these methods in vivo because of lack of understanding of what the ground truth is, and because all these methods are strongly depen- dent on the branch selection/initialization/priors. The lack of precision (due to the “continuous” degener- acy of shallow trenches) generally exists due to the multi- compartmental nature of the kernel (3.5) [46]. One can say that any standard (directional) dMRI measurement effectively under-samples the scalar part of the model (3.4), not pro- viding enough relations between the scalar parameters (cf. Sec. 3.3.4 above), and over-samples the tensor (ODF) part. In other words, the system’s true complexity lies within the kernel’s parameters hidden in functions K l ( b, x ) , Eq. (3.20) — while the ODF is in some sense “on the surface”. This prompts the need for “orthogonal” measurement schemes [47, 140, 195, 252, 254–257] which probe the scalar parameters in different combinations than entering the kernel projections K l ( b, x ) , as we are now going to discuss. Estimating precise maps of ground truth values, as well as the branch selection (3.18), remains an essential problem for quantifying neuronal microstructure, and is currently an ac- tive topic of research. Recent experiments using advanced dMRI protocols have been either employing very strong dif- fusion gradients (e.g., on unique Connectome scanners with gradients up to 300 mT/m) [167], or adding “orthogonal” ac- quisitions such as extra-neurite water suppression by strong unidirectional gradients [256] or planar diffusion weighting [140], isotropic diffusion weighting [238, 252, 254, 255, 257, [47] and the diffusion time [169, 259]. The choice of the branch, and an independent estimation of the compartment diffusivities D a and D (cid:107) e is of particular in- terest. Isotropic weighting (spherical tensor encoding) yields S ( b ) /S = f e − bD a + (1 − f ) e − b ( D (cid:107) e +2 D ⊥ e ) , (3.21) which seems to produce relations D a ≈ D (cid:107) e + 2 D ⊥ e due to an empirically small iso-weighted kurtosis of signal (3.21) [195, 257]. While this can be interpreted as favoring one of the branches, this relation cannot be used as a global con- straint: Szczepankiewicz et al. [254] show it failing in tha- lamus (note however that thalamus is a GM/WM mixture). Another possibility for using orthogonal measurements to re- solve the parameter estimation degeneracy is the application of double diffusion encoding (DDE), see Sec. 4, with promis- ing preliminary results [237]. In rat spinal cord, DDE seems to indicate the branch-merging case D a ≈ D (cid:107) e [256]. Note that such assumption is made in NODDI (section 3.4.2), al- beit this model fixes (rather than fits) the compartment dif- fusivities to equal values. This assumption does not seem to universally hold in the human brain [46]. As the ultimate goal of biophysical modeling is to study pathological and other changes (e.g., aging and development), it is imperative to es- timate the compartment diffusivities independently, because changes in one of them may indicate the earliest sign of a pathological or other process of interest. Overall, the Standard Model presents a microcosm of pa- rameter estimation challenges: a relatively low SNR in clin- ical dMRI coupled with both discrete and continuous degen- eracies, require careful validation and prompt employing the widest possible arsenal of measurements, to probe parameters from as many vantage points as possible. Achieving compart- mental specificity, crucial in studying pathological and other processes, remains a difficult but worthy goal.
4. MULTIPLE DIFFUSION ENCODINGS The whole is greater than the sum of its parts Aristotle Multiple diffusion encoding (MDE) generalizes the Stejskal- Tanner (Sec. 1.4) pulse sequence design by adding one or more extra diffusion weighting blocks, as illustrated in Fig. 10 for the case of double diffusion encoding (DDE) [131, 260–262]. Figure 10 also defines the main pulse se- quence parameters for DDE, which in addition to the familiar pulse-gradient parameters of each block, includes a mixing time τ . In the following, we will restrict our attention to the nar- row pulse limit. Generally, each diffusion weighting block is characterized by an independent diffusion wave vector q n and diffusion time t n , and the mixing times define delays between blocks. Thus, a rich set of experimentally control- lable parameters can in principle enable qualitatively differ- ent ways of probing the microstructure, as compared with the conventional, single diffusion encoded (SDE) sequences. The fundamental question of the information content of MDE signal S N relative to a set of independently acquired SDEs S can be formulated using an example of the DDE signal S (here we ignore the trivial e − R t factors, assuming the unweighted signal to be normalized to unity): S ( q , q , t , t , τ ) ≡ (cid:68) e i q · [ r (0) − r ( t )]+ i q · [ r ( t + τ ) − r ( t + t + τ )] (cid:69) = (cid:90) d r a d r b d r a d r b V e i q · ( r a − r b )+ i q · ( r a − r b ) G t ; r b , r a G τ ; r a , r b G t ; r b , r a ? = G t ; q G t ; q . (4.1) Technically, the above question is as follows: When does the convolution of the local propagators (1.5) defined in Sec. 1.4 contain more information than the product of the voxel- averaged translation-invariant SDE propagators (1.7)? Let us first outline the three cases when Eq. (4.1) holds, i.e., there is no extra information in MDE relative to SDE. a. Microscopic translation invariance: If G t ; r b , r a ≡ G t, r b − r a depends only on the relative displacement, the above equality holds for any G t, r . Of course, this is true for the Gaussian diffusion, when G t, q is de- scribed by Eq. (1.8), but the statement is much broader, since its proof (by change of integration variables r = r b − r a and r = r b − r a ) involves only the trans- lation invariance requirement. Practically, this means that the time scales involved in Eq. (4.1) exceed the time needed for the coarse-graining to restore sample’s translation invariance, whether this implies the Gaus- sian fixed point (1.8) or its anomalous counterpart, cf. Sec. 1.9. b. Long mixing time limit of a single pore:
If all spins are confined in the same pore of volume V , and τ exceeds the time to diffuse across the pore size, the “mixing” propagator G τ ; r a , r b → /V approaches a constant, and Eq. (4.1) again factorizes, irrespective of the (non- translation-invariant) functional form of G t ; r b , r a . c. Weak diffusion weighting:
Equality (4.1) holds for any G at the level of O ( q ) [263], cf. Sec. 4.2 below. For Set by the corresponding interval ∆ n between the fronts of the gradientpulses, Fig. 10. For finite pulse width δ n , see footnote 13 in Section 2. this statement to hold, it is only required that the cu- mulant expansion (1.10) has a nonzero convergence ra- dius. (This common property breaks down for a diffu- sion propagator of a stretched-exponential form, whose assumptions contradict experimental evidence [69].) Generally, the above requirements do not hold — tissues are microscopically not translation-invariant, a voxel can con- t δ g δ g Double Diffusion Encoding sequence
Resulting gradient waveform τ δ g δ g t Parallel wavevectors q q = -q Antiparallel wavevectors q q = q First diffusion encoding Second diffusion encoding
FIG. 10. (a) Example of a DDE sequence within the framework ofa double spin echo. (b) The resulting gradient waveform obtainedby multiplying each gradient by ( − n π , where n π is the numberof π pulses following the given gradient. In the text, we assumenarrow-pulse approximation, such that δ i → , with the Larmorfrequency gradients g i sufficiently large to yield finite q i = g i δ i (nosummation over i ). diffusion weighting can be strong, so that the O ( q ) terms are relevant. This justifies using MDE to obtain extra infor- mation. To get a feel for the difference between DDE and SDE, it is instructive to consider the long mixing time limit τ → ∞ in a system of disconnected pores. The average over Brownian paths splits into two parts: an average (cid:104) . . . (cid:105) over paths within a single pore, followed by an average (denoted by an overbar) over pores α = 1 , . . . , with volume fractions w α = V α /V adding to unity. For example, for the SDE signal (1.6): S ( q , t ) = (cid:10) e i q · ( r (0) − r ( t )) (cid:11) paths in pores ≡ (cid:88) pores: α w α (cid:104) e i q · ( r (0) − r ( t )) (cid:105) paths in α . For DDE, spin displacements in each of the diffusion-weighting blocks become independent of one another within each pore in the limit τ → ∞ , G τ ; r a , r b → /V α if r a and r b in Eq. (4.1) are from the same pore α with volume V α , while the probability to hop between pores is zero: G τ ; r a , r b ≡ if r a and r b belong to different pores α (cid:54) = β . This effective Kroeneker δ αβ eliminates the cross-terms between different pores that are present in the right-hand side of Eq. (4.1): S ( q , q , t , t , τ ) = (cid:10) e i q · ( r (0) − r ( t ))+ i q · [ r ( t + τ ) − r ( t + t + τ )] (cid:11) = (cid:88) α w α (cid:68) e i q · ( r (0) − r ( t )) (cid:69) α (cid:68) e i q · ( r ( t + τ ) − r ( t + t + τ )) (cid:69) α . (4.2) Here the subscript “paths in α ” was replaced with α for brevity. On the other hand, for the product of two SDE’s we have S ( q , t ) S ( q , t ) = (cid:10) e i q · ( r (0) − r ( t )) (cid:11) · (cid:10) e i q · ( r ( t + τ ) − r ( t + t + τ )) (cid:11) = (cid:88) α w α w β (cid:68) e i q · ( r (0) − r ( t )) (cid:69) α (cid:68) e i q · ( r ( t + τ ) − r ( t + t + τ )) (cid:69) β (cid:54) = S ( q , q , t , t , τ ) . (4.3) The physical meaning of the above equations is as follows: it is not possible in general to split the coherent averaging of the product (4.2) over pores into the product of the averages (4.3). The coherent disorder averaging of the propagators in equation (4.1) is also the reason that the effective medium theory [24, 30] for the disorder-averaged SDE propagator (1.6) has to be further augmented to incorporate the coarse- graining effects for MDE, relevant at finite τ . O ( q ) Historically, DDE was noted to provide a method for deter- mination of compartment dimensions [262] at low diffusion weighting in the limit of zero mixing time and long diffusion times. Taylor-expanding Eq. (4.2) in this limit, Mitra et al. [262] showed that in a system of identical pores S ( q ˆn , q ˆn , t, t, −→ t →∞ − q (cid:104) r (cid:105) (cid:20) θ (cid:21) , (4.4) where θ is the angle between the directions ˆn and ˆn of the diffusion wave vectors, and (cid:104) r (cid:105) = (cid:82) d r d r (cid:48) ( r − r (cid:48) ) / V ≡ (cid:82) d r ( r − r cm ) /V is the pore mean squared radius of gyration ( r cm is pore center-of-mass), a measure of pore size. Hence, a measure of the pore size can be determined from the signal dependence on diffusion wave vector angle in isotropic systems, or more simply from the signal difference between parallel and antiparallel diffusion wave vectors. Equation (4.4) has since been generalized to take into ac- count, e.g., partial volume, multiple concatenations, pulse se- quence timings (e.g., finite gradient width) for various ge- ometries [264–270]. This has later been demonstrated by sev- eral groups in model systems and biological samples ex vivo [271–277], and in vivo in humans [276, 278–280]. However, it was recently realized [263] that this property, i.e., the sensitivity (4.4) to pore gyration radius, is a general feature of any diffusion-weighted signal at the O ( q ) level, and hence it does not rely on information beyond that already contained in the SDE signal, which in the same regime be- haves as [281]: S ( q , t ) −→ t →∞ − q i q j (cid:104) ( x i ( t ) − x i (0))( x j ( t ) − x j (0)) (cid:105) = 1 − q (cid:104) r (cid:105) . (4.5) More generally, it was shown [263] that up to order O ( q ) , ln S ( q , q , t , t , τ ) = − q i q j D ij ( t ) t + q i q j D ij ( t ) t + q i q j [ D ij ( t + t + τ )( t + t + τ ) + D ij ( τ ) τ − D ij ( t + τ )( t + τ ) − D ij ( t + τ )( t + τ )] + O ( q ) , (4.6) D ij ( t ) is the cumulative diffusion tensor, Eq. (1.13). This explicitly demonstrates that the signal is fully character- ized by the time-dependent diffusion tensor, a quantity which is obtainable from the SDE acquired at a few diffusion times. This statement is valid for any diffusion sequence up to sec- ond order in the diffusion wave vector [263, 282], and is a consequence of the existence of the cumulant series, whose lowest order can be completely reproduced by knowing the diffusion tensor for all t or, equivalently, for all ω [15, 130], ln S = − (cid:90) (cid:104) v i ( t ) v j ( t ) (cid:105) q i ( t ) q j ( t ) d t d t , (4.7) where q ( t ) is the time integral of the arbitrary-shaped applied gradient, v is the molecular velocity, v = ∂ t r , and no bulk flow is assumed as usual. The (symmetric) autocorrelation function (cid:104) v i ( t ) v j ( t ) (cid:105) ≡ D ij ( | t − t | ) is constructed out of its retarded counterpart D ij ( t ) ≡ θ ( t ) (cid:104) v i ( t ) v j (0) (cid:105) = ∂ ∂t [ tD ij ( t )] (4.8) defined by generalizing Eqs. (2.6) and (2.14) to the anisotropic case. The function (4.8) is generally nonlocal in t [15, 16, 24,
30, 122, 130, 283, 284]. Figure 11 illustrates how this non- locality, integrated in Eq. (4.7), gives rise to the cross-term ∼ q i q j (second line of Eq. (4.6)); this term disappears in the Gaussian diffusion limit, when (cid:104) v i ( t ) v j ( t ) (cid:105) = 2 D ij δ ( t − t ) is infinitely narrow, and the function (4.8) is concentrated along the diagonal. Hence, the cross-term ∼ q i q j directly probes the time-dependence of the diffusion coefficient in Eq. (4.6), cf. Section 2. O ( q ) and beyond. Microscopic anisotropy At larger values of the diffusion weighting, double diffusion encoding was shown from the beginning to have the ability to characterize microscopic anisotropy ( µ A) in systems which are macroscopically isotropic, Fig. 12, see panel (c). Thus, in an early application of the sequence by Cory et al. [261], DDE was used to quantitatively measure the eccentricity of yeast cells, which was shown to be directly related to the difference in signals acquired with parallel and perpendicu- lar diffusion wave vectors. This has since been explored by many authors, e.g., in phantoms [285, 286], ex vivo tissues [287, 288], and in vivo [276, 289, 290]. The basic sensitivity to anisotropic pores can be understood already from Eq. (4.2) in the long diffusion time and long mixing time limit S ( q , q , t , t , τ ) = | χ α ( q ) | | χ α ( q ) | (4.9) where χ α ( q ) is the Fourier transform of the pore structure function, defined as χ α ( r ) ≡ /V α inside a pore and 0 other- wise [291]. For spherical pores, the structure function is isotropic, hence χ α ( q ) does not depend on the direction of q . For t t FIG. 11. Two-dimensional temporal integration involved in thesecond-order cumulant, Eq. (4.7) leading to Eq. (4.6) for the DDEmeasurement. Labels q and q indicate the time interval in which q ( t ) equals to q and q , respectively; for simplicity, the vector in-dices are not shown. The green shaded area along the diagonal sym-bolizes D ij ( | t − t | ) , Eq. (4.8), where it significantly deviates from , with the width of this region set by the correlation time t c . Thenontrivial cross-term q i q j in Eq. (4.6) arises from the off-diagonalquadrants. As this contribution is weighted with the velocity auto-correlation function, it tends to zero when the mixing time, τ (indi-cated by the thin lines along each dimension) becomes larger thanthe correlation time, τ (cid:29) t c . In particular, no non-trivial cross-termis present for Gaussian diffusion, for which t c → . anisotropic pores, say ellipsoids, the anisotropic structure functions χ α ( q ) in Eq. (4.9) ensure that the result depends on the directions of q and q . This is because diffusion in the two directions is generally correlated when pores are non- spherical. If the overall system is macroscopically isotropic, i.e., the orientations of the individual pores are randomly dis- tributed (Fig. 12c), the signal will be unaffected by rotations of the sample or of the laboratory system of reference, but the dependence on the relative angles between q and q will survive the pore averaging in Eq. (4.9). Mathemati- cally speaking, this is because the two terms in the product | χ α ( q ) | | χ α ( q ) | are not independent for a given pore, and hence the average of the product is different from the product of averages. A convenient measure of the eccentricity of the pore space (microscopic anisotropy) can therefore be found from the dif- ference of DDE signals acquired with parallel and perpen- dicular wave vectors q and q . In the presence of macro- scopic and microscopic anisotropy, the signal will depend on the orientations of both q and q , and microscopic diffu- sion anisotropy can no longer be extracted simply from the difference between parallel and perpendicular diffusion wave vectors. The rotationally invariant way to circumvent this (cf. Sec. 3.5 above) is to powder average the signal, analo- gously to how the K invariant was introduced in Sec. 3.5.1 (although, technically, the averaging here is over the SO(3) group instead of a 2-sphere), and practical recipes for doing this were proposed in refs. [287, 292, 293]. Microscopic dif- fusion anisotropy can then be defined as the difference be- FIG. 12. Examples of model systems considered in the text. In(a), a system of identical spherical pores is shown, whereas, in (b),the pores have a distribution of sizes. In (c), an approximatelyisotropic distribution of ellipsoidal pores is sketched and, in (d), thepores are coherently oriented. Systems (a)–(c) are macroscopicallyisotropic, system (d) is not. Systems (c) and (d) are microscopi-cally anisotropic . Ensemble heterogeneity is only seen in systems(b) (size) and (c) (orientation). Here, spins contributing to the signalare assumed to only reside within the pores. tween (log of) the powder averaged signals acquired with par- allel and perpendicular diffusion wave vectors [287]. Mitra’s original results were since generalized to arbitrary diffusion and pulse timings by several authors [287, 292, regime and for arbitrary diffusion times, the signal can be written as ln S ( q , q ) = − ( q i q j + q i q j ) D ij ( t ) t + D q i q j q k q l + q i q j q k q l ) W ijkl ( t ) + 14 q i q j q k q l Z ijkl ( t ) (4.10) where W is the kurtosis tensor as defined in [60] from the cumulant expansion (1.10) of the SDE propagator, whereas Z is a rank-4 tensor, unique to DDE, defined as Z ijkl = 4 t (cid:0) D αij D αkl − D αij D αkl (cid:1) . (4.11) The tensors D αij ( t ) refer to the microscopic t -dependent dif- fusion tensors characterizing diffusion within the pores, and the SDE-measured overall diffusion tensor D ij ( t ) entering the first line of Eq. (4.10) is an average over all pores, D ij ( t ) = D αij ( t ) ≡ (cid:80) α w α D αij ( t ) . The new tensor Z , Eq. (4.11), accessible with DDE (and inaccessible with SDE), is proportional to the covariance ten- sor of microscopic diffusion tensors. Microscopic diffusion anisotropy, defined as the difference between log of the pow- der averaged signals acquired with parallel and perpendicular diffusion wave vectors, can then be expressed as (see [287]) ε = 160 (cid:2) Z ijij − Z iijj + 2 t (3 D ij D ij − D ii D jj ) (cid:3) = t (cid:0) D αij D αij − D αii D αjj (cid:1) = 35 t var { σ α } . (4.12) In the last equality, the set σ α ≡ { σ α,i } i =1 denote the eigen- values of D α , and var { σ α } ≡ (cid:88) i =1 σ α,i − (cid:34) (cid:88) i =1 σ α,i (cid:35) . (4.13) With the above definition, D α ) − (tr D α ) = { σ α } in Eq. (4.12). The anisotropy metric ε has di- mensions of [length] . These somewhat awkward dimen- sions have a historical root in DDE eccentricity measure- ments [261]. While ε ≥ , in practice it is often estimated from the difference of signals, which can become negative due to noise. As an example, for randomly oriented (and identical) axially symmetric domains, such as fibers with (time- dependent) diffusivities D (cid:107) ( t ) and D ⊥ ( t ) , microscopic dif- fusion anisotropy becomes ε = 215 t (cid:2) D ⊥ ( t ) − D (cid:107) ( t ) (cid:3) . (4.14) If the domains are different, the corresponding Eq. (4.14) should be further averaged over them, cf. the var { σ α } term in Eq. (4.12). Microscopic diffusion anisotropy hence depends explicitly on diffusion time, but tends to the geometric mea- sures of pore shape anisotropy as the diffusion time increases, since for any confined region of size a , D ( t ) ∼ a /t , and t asymptotically drops out from Eq. (4.14). From the time dependence of the microscopic diffusion anisotropy, non- Gaussian effects of the individual compartments can be re- vealed by the time dependence of the compartmental (micro- scopic) diffusion tensors. Practically, the anisotropy metric (4.12) can be estimated from knowledge of the full Z tensor, or by the difference of the powder averaged log signals with parallel and perpendic- ular diffusion wave vectors. It has an advantage of being ad- ditive (cf. the pore average in Eq. (4.12)): if several distinct types of pore populations are present in the sample (e.g., a dis- tribution of D ⊥ and D (cid:107) in Eq. 4.14), ε simply becomes the volume-weighted mean over the corresponding ε from each of the populations. This is an advantage since it eases the interpretation; however, the disadvantage is the dependence on size of the pore in addition to its anisotropy. This has the additional consequence that ε is strongly biased by the larger pores: Since w α ∼ a and ε α ∼ a for a pore of size a , the population averaged eccentricity scales as ε ∼ a /a , heavily preferring the tail of the pore size distribution, — and hence susceptible to the mesoscopic fluctuations introduced in Sec. 2.5 above, cf. Eq. (2.39). To factor out the pore sizes, normalized dimensionless measures of microscopic diffusion anisotropy were intro- duced [287, 292, 295], such as the microscopic fractional µ FA: µ FA ≡ (cid:115)
32 ( σ − ¯ σ ) + ( σ − ¯ σ ) + ( σ − ¯ σ ) σ + σ + σ = (cid:115) εε + t (cid:0) tr D (cid:1) . (4.15) In the previous example with axially symmetric domains, mi- croscopic fractional anisotropy µ FA = (cid:114) (cid:12)(cid:12) D (cid:107) − D ⊥ (cid:12)(cid:12)(cid:113) D (cid:107) + 2 D ⊥ , (4.16) whereas fractional anisotropy FA is modulated also by the fiber orientation distribution function [287, 296], and only re- covers µ FA when the fibers are all coherently aligned. An- other metric which has been suggested to be of biological im- portance [252, 254], is the variance in isotropic diffusivity, V I ≡ ( D αii / − ( D αii / = 136 t Z iijj (4.17) which can also be inferred from the Z tensor [132, 297]. When diffusion within the individual pores is Gaussian, other methods such as the so-called magic angle spinning of the q-vector (q-MAS) [298] can also be used to estimate the dif- fusion tensor covariance [132, 299]. As we can see, MDE can potentially provide unique extra in- formation relative to SDE. However, this information content only starts at the level of O ( q ) , Eq. (4.10) (in addition to the standard SDE O ( q ) terms), and hence to claim the true novelty of the information, it has to be properly identified rel- ative to the SDE measurements with similar scan parameters (timings and gradients). Clarifying the advantages of MDE is practically essential in the view of much reduced SNR due to a notable increase of the echo time needed for the multiple gradients to play out. Overall, the main advantage of MDE so far seems to lie in its ability to detect and quantify microscopic diffusion anisotropy. In particular, the advantage of the DDE metrics of microscopic diffusion anisotropy (e.g., Eq. (4.12)) is that they do not rely on concrete assumptions regarding pore shapes. Of course, if a detailed model of microstructure is avail- able, e.g., of the Standard Model form (Section 3), micro- scopic diffusion anisotropy is directly accessible in terms of model parameters (which can in principle be determined us- ing SDE). However, even in this case, MDE adds value by providing an “orthogonal” way of being sensitive to model parameters — and, like in the SM case above, Sec. 3.6, it can provide rotationally invariant independent relations between parameters, which can help lift the parameter estimation de- generacies [237]. So far, the existing MDE models have been calculated in the limit of either Gaussian diffusion in all compartments, or in the t i → ∞ limit (e.g., closed pores). Importantly, the transient effects, cf. Section 2, have not yet been properly ac- counted for in the MDE framework. In particular, the struc- tural disorder-induced power-law tails in D inst ( t ) , or, equiva- lently, in D ( ω ) , such as the ones originating due to disordered axonal packing in the extra-axonal space [20, 31], will con- tribute to the “irreducible” MDE effects (that go beyond the product of a few SDE signals). Taking such transient pro- cesses into account seems a priori as crucial for the inter- pretation of MDE measurements, as it has been for the SDE — e.g., in the context of recent re-interpretation of the ax- onal diameter mapping results (cf. Sec. 2.4.2 above). The relevant coarse-graining formalism of the effective medium theory [20, 24, 30] seems perfectly suitable for the task — but it has not yet been developed.
5. OUTLOOK AND OPEN QUESTIONS There is nothing more practical than a good theoryL. Boltzmann We are writing this Review at a transformational moment, when our field of quantitative dMRI is experiencing a revolu- tion due to unprecedented quality of hardware and novel ac- quisition methods, enabling us to observe very subtle physical effects, even in human subjects and potentially in patients. Interpreting these effects in terms of the tissue microar- chitecture is highly nontrivial; it is safe to say that the the- oretical challenge has been so far greatly underappreciated. This, however, may swing the pendulum the other way, to- wards an “anti-modeling” point of view: Since, according to a widespread refrain, “biology is so much more complicated that anything physicists have ever studied”, there is little hope for the quantitative understanding of such effects, and the best we can do is to stay at the level of “representations” (cf. Sec. 1.7) and to draw empirical correlations between param- eters of such representations (e.g. mean diffusivity or frac- tional anisotropy) and the clinical disease scores. One of the messages of our Review is that the whole his- tory of Physics in the 20 th century offers the case for opti- mism. The quote from a nuclear physicist (before Section 1) has been a universal refrain for our sustained ability to un- derstand nature’s complexity, step-by-step, from the origins of elementary particles to the vast scopes of the Universe. The essence of the effective theory way of thinking is that one certainly does not need to understand everything about the world in order to understand some corner of its param- eter space really well. We certainly can quantify tissue mi- croarchitecture without uncovering the origins of the human conscience or mapping full details of the brain’s biochemical machinery. Too often, “biology is way more complex than all models. We reject this excuse [69]. We believe that having ap- propriate theoretical description of diffusion in tissues at the mesoscopic scale is not a luxury at this point — rather, this is an indispensable scientific method of investigation into pathological processes 2-3 orders of magnitude below nom- inally achievable resolution of MRI in any foreseeable fu- ture — or, in fact, ever, since the MRI resolution is strin- gently bounded by physical and physiological limitations that have been largely reached by now. The parallels with super- resolution microscopy [300] are quite obvious; that discipline took a century to develop, based on employing models and prior information. Our task is harder but, arguably, can lead to even more impactful advances. With that in mind, let us outline 10 exciting unresolved problems, focussing on which, to the best of our understand- ing, will propel our field forward.
1. Apparent vs. real diffusion metrics:
What are the confounding effects of mesoscopically varying R ( r ) , R ( r ) , and Ω( r ) in the mesoscopic Bloch-Torrey equa- tion (1.1) on the observed diffusion metrics, in the spirit of refs. [12, 13]? Can we develop a multi-modal meso- scopic imaging framework able to self-consistently quantify all these mesoscopic quantities and disentan- gle their effects in the apparent diffusion coefficients and higher-order metrics in each tissue compartment?
2. Relation between time-dependent D ⊥ ( t ) and its tor- tuosity limit D ∞ , and the geometric parameters of realistic axonal packings: As the time dependence of the diffusion transverse to fiber tracts is dominated by the extra-axonal water (Section 2), the natural question is what structural changes (e.g. demyelination, axonal loss) can affect this time dependence, as well as D ∞ . This is a difficult yet clinically impactful inverse prob- lem [218], whose approximate solution, relying on the ideas of coarse-graining and renormalization, has so far only been obtained in the t → ∞ limit [219].
3. Origin of structural disorder along the neurites: What causes the time dependence along the fibers or in the gray matter? Is it varicosities, beads, synap- tic boutons, undulations, or something else? Which of these structural units’ changes in pathology can be detectable?
4. Parameter estimation challenge for the Standard Model:
How many Gaussian compartments do we have to include? For increasing the precision, it looks like we need orthogonal measurements, such as MDE (e.g., isotropic diffusion weighting), and varying echo time. What is an optimal clinically feasible measure- ment protocol?
5. Time-dependent rotationally-invariant framework: Combining the ideas of Sections 2 and 3 can lead to de- scribing each fiber fascile in terms of the non-Gaussian propagators (inside and outside the neurites) with the corresponding time-dependent diffusion, kurtosis, etc, cumulant tensors; such fascicles then naturally com- bine into the SM-like signal based on the fiber ODF in a voxel. This difficult parameter estimation problem may offer the all-encompassing description of diffusion process measurable with dMRI in the brain.
6. Permeability/exchange time for the neurites:
How well we can approximate compartments as non- exchanging? At which time scales this assumption breaks? The answer most likely will be different for gray and white matter, and maybe even for different brain regions.
7. Standard Model for GM:
Can we apply SM as intro- duced in Section 3 to gray matter in vivo at clinical dif- fusion times, or should we modify the compartments? Do we have to include exchange, and if yes, at which level?
8. EMT for MDE:
Development of the effective medium theory framework [24] for the “disorder-averaging” in- volved in the multiple diffusion encoding signal (Sec- tion 4). Which physical effects, from the EMT stand- point, are best captured using MDE, or are completely absent in the SDE?
9. Signal vs. noise:
As the saying goes, “noise is sig- nal”. The fundamental question is to separate the ther- mal noise, imaging artifacts, as well as the genuine dif- ferences between parameters in voxels belonging to the same region of interest. Recently developed random matrix theory-based approaches [181, 182] offer an ex- citing prospect.
10. Mesoscopic fluctuations and biological variability: How different are the mesoscopic tissue parameters within a given region of interest? Their differences pro- vide the “natural” minimal width for the parameter dis- tributions within an ROI, in the limit of infinite SNR. Sometimes, relatively small differences in the meso- scopic parameters can translate into large differences of the dMRI metrics; the heavy sensitivity of the sig- nal from water inside axons to the tail of the axonal diameter distribution [20] (Secs. 2.5, 4.3) is an exam- ple of an effect of so-called “mesoscopic fluctuations” pioneered within condensed matter physics [7]. Study- ing these fluctuations can provide fundamental insights on the optimality and robustness of the organization of neuronal tissue microarchitecture, as well as offer prac- tical limits on our detection capabilities. ACKNOWLEDGMENTS It is a pleasure to thank our numerous colleagues and mem- bers of our research groups for stimulating discussions and collaborations which are reflected in our Review. In par- ticular, we thank Christian Beaulieu and Gene Kim for dis- cussions on experimental issues, and Dan Wu and Jiangyang Zhang for sharing their OGSE data displayed in Fig. 7. We also thank Ivana Drobnjak, Ileana Jelescu, Markus Nilsson and Marco Palombo, as well as the referees, for their thought- ful and constructive comments on the manuscript. Photo credit to Tom Deerinck and Mark Ellisman (National Center for Microscopy and Imaging Research) for the histol- ogy image illustrating a fiber fascicle in Figs. 1 and 8. E.F. and D.S.N. were supported by the National Institute of Neurological Disorders and Stroke of the NIH under award number R01NS088040. SNJ was supported by the Danish Ministry of Science, Technology and Innovations University Investment Grant (MINDLab, Grant no. 0601-01354B), and the Lundbeck Foundation R83-A7548. [1] Howard Georgi, “Effective field theory,” Annual review of nu- clear and particle science , 209–252 (1993). [2] Philip W Anderson, “More is different,” Science , 393–
396 (1972). [3] Kenneth G Wilson, “The renormalization group and critical phenomena,” Reviews of Modern Physics , 583 (1983). [4] John Cardy, Scaling and renormalization in statistical physics , Vol. 5 (Cambridge university press, 1996). [5] N Bloembergen, E M Purcell, and R V Pound, “Relaxation Effects in Nuclear Magnetic Resonance Absorption,” Physi- cal Review , 679–712 (1948). [6] A Abragam, Principles of Nuclear Magnetism (Oxford Uni- versity Press, New York, 1961). [7] Yoseph Imry, Introduction to mesoscopic physics (Oxford University Press, New York, 1997). [8] R. Mark Henkelman, Xuemei Huang, Qing-San Xiang, G. J. Stanisz, Scott D. Swanson, and Michael J. Bronskill, “Quan- titative interpretation of magnetization transfer,” Magnetic Resonance in Medicine , 759–766 (1993). [9] R. M. Henkelman, G. J. Stanisz, and S. J. Graham, “Magne- tization transfer in MRI: a review,” NMR in Biomedicine , [10] A A Migdal, “Phase transitions in gauge and spin-lattice sys- tems,” Soviet Physics JETP , 743–746 (1976). [11] Leo P Kadanoff, “Notes on Migdal’s recursion formulas,” An- nals of Physics , 359–394 (1976). [12] Jianhui Zhong, Richard P. Kennan, and John C. Gore, “Effects of susceptibility variations on NMR measurements of diffusion,” Journal of Magnetic Resonance , 267–280 (1991). [13] V. G. Kiselev, “Effect of magnetic field gradients induced by microvasculature on NMR measurements of molecular self- diffusion in biological tissues,” Journal of Magnetic Reso- nance , 228–235 (2004). [14] Paul T. Callaghan, Principles of nuclear magnetic resonance microscopy (Clarendon Press, 1993). [15] V. G. Kiselev, “The cumulant expansion: an overarch- ing mathematical framework for understanding diffusion { NMR } ,” in Diffusion { MRI } : Theory, methods and appli- cations , edited by Derek K. Jones (Oxford University Press, New York, 2010) Chap. 10. [16] Valerij G Kiselev, “Fundamentals of diffusion { MRI } physics,” { NMR } In Biomedicine
DOI:10.100 (2016). [17] J. S. Murday and R. M. Cotts, “Self-Diffusion Coefficient of Liquid Lithium,” The Journal of Chemical Physics , 4938 (1968). [18] C H Neuman, “Spin echo of spins diffusing in a bounded medium,” Journal of Chemical Physics , 4508–4511 (1974). [19] P van Gelderen, D DesPres, P C van Zijl, and C T Moonen, “Evaluation of restricted diffusion in cylinders. Phosphocrea- tine in rabbit leg muscle.” (1994). [20] Lauren M Burcaw, Els Fieremans, and Dmitry S. Novikov, “Mesoscopic structure of neuronal tracts from time- dependent diffusion,” NeuroImage , 18–37 (2015). [21] Hong-Hsi Lee, Lauren M Burcaw, Jelle Veraart, Els Fiere- mans, and Dmitry S. Novikov, “Low-Pass Filter Effect of Fi- nite Gradient Duration on Time-Dependent Diffusion in the Human Brain,” Proceedings of the International Society of Magnetic Resonance in Medicine 23, p. 2777 (2015). [22] Hong-Hsi Lee, Gregory Lemberskiy, Els Fieremans, and Dmitry S. Novikov, “Estimation of Fiber Packing Correlation Length by Varying Diffusion Gradient Pulse Duration,” Pro- ceedings of the International Society of Magnetic Resonance in Medicine 24, p. 2021 (2016). [23] Hong-Hsi Lee, Els Fieremans, and Dmitry S. Novikov, “What dominates the time dependence of diffusion transverse to axons: Intra- or extra-axonal water?” NeuroImage (2017), [24] Dmitry S. Novikov and Valerij G Kiselev, “Effective medium theory of a diffusion-weighted signal,” NMR Biomed , [25] P T Callaghan, C D Eccles, and Y Xia, “NMR microscopy of dynamic displacements: k-space and q-space imaging,” Jour- nal of Physics E: Scientific Instruments , 820–822 (1988). [26] D G Cory and A N Garroway, “Measurement of translational displacement probabilities by NMR: an indicator of compart- mentation,” Magn Reson Med , 435–444 (1990). [27] D Le Bihan, E Breton, D Lallemand, P Grenier, E Cabanis, and M Laval-Jeantet, “MR imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders.” Radiology , 401–7 (1986). [28] Pj Basser, J Mattiello, and D LeBihan, “Estimation Of The Effective Self-Diffusion Tensor From The { NMR } Spin-
Echo,” Journal of Magnetic Resonance Series B , 247–254 (1994). [29] Partha P. Mitra, Pabitra N. Sen, and Lawrence M. Schwartz, “Short-Time Behavior of the Diffusion Coefficient As a Geo- metrical Probe of Porous Media,” (1993). [30] Dmitry S. Novikov, Jens H Jensen, Joseph A Helpern, and
Els Fieremans, “Revealing mesoscopic structural universality with diffusion,” Proc Natl Acad Sci U S A , 5088–5093 (2014). [31] Els Fieremans, Lauren M Burcaw, Hong-Hsi Lee, Gregory
Lemberskiy, Jelle Veraart, and Dmitry S. Novikov, “In vivo observation and biophysical interpretation of time-dependent diffusion in human white matter,” NeuroImage , 414–427 (2016). [32] A L Sukstanskii, D A Yablonskiy, and J J H Ackerman, “Effects of permeable boundaries on the diffusion-attenuated { MR } signal: insights from a one-dimensional model,” Jour- nal of Magnetic Resonance , 56–66 (2004). [33] Christopher D Kroenke, Joseph J H Ackerman, and Dmitriy A Yablonskiy, “On the nature of the NAA diffusion attenuated MR signal in the central nervous system,” Mag- netic Resonance in Medicine , 1052–1059 (2004). [34] T E J Behrens, M W Woolrich, M Jenkinson, H Johansen- Berg, R G Nunes, S Clare, P M Matthews, J M Brady, and S M Smith, “Characterization and propagation of uncertainty in diffusion-weighted { MR } imaging,” Magn Reson Med , [35] T E J Behrens, H Johansen-Berg, M W Woolrich, S M Smith, C A M Wheeler-Kingshott, P A Boulby, G J Barker, E L Sillery, K Sheehan, O Ciccarelli, and Others, “Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging,” Nature Neuroscience , 750–757 (2003). [36] Sune N Jespersen, Christopher D Kroenke, Leif Ostergaard, Joseph J H Ackerman, and Dmitriy A Yablonskiy, “Modeling dendrite density from magnetic resonance diffusion measure- ments,” Neuroimage , 1473–1486 (2007). [37] Sune N Jespersen, Carsten R Bjarkam, Jens R Nyengaard, M Mallar Chakravarty, Brian Hansen, Thomas Vosegaard, Leif Ostergaard, Dmitriy Yablonskiy, Niels Chr. Nielsen, and Peter Vestergaard-Poulsen, “Neurite density from magnetic resonance diffusion measurements at ultrahigh field: Com- parison with light microscopy and electron microscopy,” Neu- roimage , 205–216 (2010). [38] Els Fieremans, Dmitry S. Novikov, Jens H Jensen, and Joseph A Helpern, “Monte { C } arlo study of a two- compartment exchange model of diffusion,” NMR in Biomedicine , 711–724 (2010). [39] Els Fieremans, Jens H Jensen, and Joseph A Helpern, “White matter characterization with diffusional kurtosis imaging,” Neuroimage , 177–188 (2011). [40] Stamatios N Sotiropoulos, Timothy E J Behrens, and Saad Jbabdi, “Ball and rackets: inferring fiber fanning from diffusion-weighted MRI,” NeuroImage , 1412–1425 (2012). [41] Hui Zhang, Torben Schneider, Claudia A Wheeler-Kingshott, and Daniel C Alexander, “ { NODDI } : practical in vivo neu- rite orientation dispersion and density imaging of the human brain,” Neuroimage , 1000–1016 (2012). [42] Marco Reisert, Valerij G Kiselev, Bibek Dihtal, Elias Kell- ner, and Dmitry S. Novikov, “MesoFT: unifying diffu- sion modelling and fiber tracking,” in Medical Image Com- puting and Computer-Assisted Intervention–MICCAI 2014 (Springer, 2014) pp. 201–208. [43] Ileana O Jelescu, Jelle Veraart, Els Fieremans, and Dmitry S. Novikov, “Degeneracy in model parameter estimation for multi-compartmental diffusion in neuronal tissue,” NMR in Biomedicine , 33–47 (2016). [44] Marco Reisert, Elias Kellner, Bibek Dhital, J¨urgen Hennig, and Valerij G. Kiselev, “Disentangling micro from mesostruc- ture by diffusion MRI: A Bayesian approach,” NeuroImage , 964–975 (2017). [45] Dmitry S. Novikov, Ileana O Jelescu, and Els Fieremans, “From diffusion signal moments to neurite diffusivities, vol- ume fraction and orientation distribution: An exact solution,” Proceedings of the International Society of Magnetic Reso- nance in Medicine 23, p. 469 (2015). [46] Dmitry S. Novikov, Jelle Veraart, Ileana O. Jelescu, and Els Fieremans, “Rotationally-invariant mapping of scalar and orientational metrics of neuronal mi- crostructure with diffusion MRI,” NeuroImage (2018), [47] Jelle Veraart, Dmitry S. Novikov, and Els Fieremans, “TE dependent Diffusion Imaging (TEdDI) distinguishes between compartmental T2 relaxation times,” NeuroImage (2017), [48] J E Tanner, “Transient diffusion in a system partitioned by permeable barriers. Application to NMR measurements with a pulsed field gradient,” The Journal of Chemical Physics , [49] J. G. Powles, M. J. D. Mallett, G. Rickayzen, and W. a. B. Evans, “Exact Analytic Solutions for Diffusion Impeded by an Infinite Array of Partially Permeable Barriers,” Proceed- ings of the Royal Society A: Mathematical, Physical and En- gineering Sciences , 391–403 (1992). [50] Dmitry S. Novikov, Els Fieremans, Jens H Jensen, and Joseph A Helpern, “Random walks with barriers,” Nature Physics , 508–514 (2011). [51] Eric E Sigmund, Dmitry S. Novikov, Dabang Sui, Obehi Ukpebor, Steven Baete, James S Babb, Kecheng Liu, Thorsten Feiweier, Jane Kwon, KellyAnne McGorty, and Others, “Time-dependent diffusion in skeletal muscle with the random permeable barrier model { (RPBM) } : application to normal controls and chronic exertional compartment syn- drome patients,” NMR in Biomedicine , 519–528 (2014). [52] Els Fieremans, Gregory Lemberskiy, Jelle Veraart, Eric E Sig- mund, Soterios Gyftopoulos, and Dmitry S. Novikov, “In vivo measurement of membrane permeability and myofiber size in human muscle using time-dependent diffusion ten- sor imaging and the random permeable barrier model,” NMR Biomed (2016), 10.1002/nbm.3612. [53] P T Callaghan, K W Jolley, and J Lelievre, “Diffusion of water in the endosperm tissue of wheat grains as studied by pulsed field gradient nuclear magnetic resonance.” Biophysi- cal journal , 133–41 (1979). [54] Dmitriy A Yablonskiy, Alexander L Sukstanskii, Jason C Leawoods, David S Gierada, G Larry Bretthorst, Stephen S
Lefrak, Joel D Cooper, and Mark S Conradi, “Quantitative in vivo assessment of lung microstructure at the alveolar level with hyperpolarized 3He diffusion MRI,” Proceedings of the
National Academy of Sciences , 3111–3116 (2002). [55] G J Stanisz, A Szafer, G A Wright, and R M Henkelman, “An Analyitical Model of Restricted Diffusion in Bovine Op- tic Nerve,” Magnetic Resonance in Medicine , 103–111 (1997). [56] Y Assaf, R Z Freidlin, G K Rohde, and P J Basser, “New modeling and experimental framework to characterize hin- dered and restricted water diffusion in brain white matter,” Magnetic Resonance In Medicine , 965–978 (2004). [57] Yaniv Assaf, Tamar Blumenfeld-Katzir, Yossi Yovel, and Pe- ter J. Basser, “Axcaliber: A method for measuring axon diam- eter distribution from diffusion MRI,” Magnetic Resonance in Medicine , 1347–1354 (2008). [58] D C Alexander, P L Hubbard, M G Hall, E A Moore, M Ptito, G J M Parker, and T B Dyrby, “Orientationally invariant in- dices of axon diameter and density from diffusion MRI,” Neu- roimage , 1374–1389 (2010). [59] D Le Bihan, C T Moonen, P C van Zijl, J Pekar, and D De- sPres, “Measuring random microscopic motion of water in tis- sues with MR imaging: a cat brain study,” J Comput Assist Tomogr , 19–25 (1991). [60] Jens H Jensen, Joseph A Helpern, Anita Ramani, Hanzhang Lu, and Kyle Kaczynski, “Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging,” Magn Reson Med , 1432– [61] Dmitry S. Novikov and V G Kiselev, “Transverse { NMR } relaxation in magnetically heterogeneous media,” J Magn Re- son , 33–39 (2008). [62] Evren ¨Ozarslan, Cheng Guan Koay, Timothy M Shepherd, Michal E Komlosh, M Okan \ .Irfano \ u glu, Carlo Pierpaoli, and Peter J Basser, “Mean apparent propagator (MAP) MRI: a novel diffusion imaging method for mapping tissue mi- crostructure,” Neuroimage , 16–32 (2013). [63] Dmitriy A Yablonskiy, G Larry Bretthorst, and Joseph J H Ackerman, “Statistical model for diffusion attenuated MR signal,” Magn Reson Med , 664–669 (2003). [64] Yong Wang, Qing Wang, Justin P Haldar, Fang-Cheng Yeh, Mingqiang Xie, Peng Sun, Tsang-Wei Tu, Kathryn Trinkaus, Robyn S Klein, Anne H Cross, and Others, “Quantification of increased cellularity during inflammatory demyelination,” Brain , 3590–3601 (2011). [65] Nathan S White, Trygve B Leergaard, Helen D’Arceuil, Jan G Bjaalie, and Anders M Dale, “Probing tissue microstructure with restriction spectrum imaging: histological and theoreti- cal validation,” Human brain mapping , 327–346 (2013). [66] Benoit Scherrer, Armin Schwartzman, Maxime Taquet, Mustafa Sahin, Sanjay P Prabhu, and Simon K Warfield, “Characterizing brain tissue by assessment of the distribu- tion of anisotropic microstructural environments in diffusion- compartment imaging (DIAMOND),” Magn Reson Med , [67] P P Mitra, P N Sen, L M Schwartz, and P Le Doussal, “Diffu- sion Propagator as a Probe of the Structure of Porous Media,” Physical Review Letters , 3555–3558 (1992). [68] Paul T. Callaghan, A. Coy, D. MacGowan, K. J. Packer, and F. O. Zelaya, “Diffraction-like effects in NMR diffusion stud- ies of fluids in porous solids,” (1991). [69] Dmitry S. Novikov, Valerij G. Kiselev, and Sune N. Jes- persen, “On modeling,” Magnetic Resonance in Medicine , [70] Jens H Jensen, Joseph A Helpern, Anita Ramani, Hanzhang Lu, and Kyle Kaczynski, “Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging,” Magn Reson Med , 1432– [71] R A Fisher and J Wishart, “The derivation of the pattern for- mulae of two-way partitions from those of simpler patterns,” Proc. London Math. Soc. (1) s2-33 , 195–208 (1932). [72] N G van Kampen, Stochastic Processes in Physics and Chem- istry , 1st ed. (Elsevier, Oxford, 1981). [73] A F Frøhlich, L Ostergaard, and V G Kiselev, “Effect of Im- permeable Boundaries on Diffusion-Attenuated { MR } Sig- nal,” J. Magn. Reson. , 223–233 (2006). [74] Valerij G Kiselev and Kamil A Il’yasov, “Is the “biexponen- tial diffusion” biexponential?” Magn Reson Med , 464–469 (2007). [75] Jens H Jensen and Joseph A Helpern, “MRI quantification of non-Gaussian water diffusion by kurtosis analysis,” NMR Biomed , 698–710 (2010). [76] Jelle Veraart, Jeny Rajan, Ronald R Peeters, Alexander Lee- mans, Stefan Sunaert, and Jan Sijbers, “Comprehensive framework for accurate diffusion MRI parameter estimation,” Magn Reson Med , 972–984 (2013). [77] Kip S Thorne, “Multipole expansions of gravitational radia- tion,” Reviews of Modern Physics , 299–339 (1980). [78] J-P Bouchaud and A Georges, “Anomalous diffusion in dis- ordered media - statistical mechanisms, models and physical applications,” Physics Reports - Review Section of Physics Letters , 127–293 (1990). [79] A S Lamantia and P Rakic, “Cytological And Quantitative Characteristics Of 4 Cerebral Commissures In The Rhesus- Monkey,” Journal of Comparative Neurology , 520–537 (1990). [80] F Aboitiz, A B Scheibel, R S Fisher, and E Zaidel, “Fiber composition of the human corpus callosum,” Brain Research , 143–153 (1992). [81] Y Tang, J R Nyengaard, B Pakkenberg, and H J G Gunder- sen, “Age-induced white matter changes in the human brain: a stereological investigation,” Neurobiol Aging , 609–615 (1997). [82] Roberto Caminiti, Hassan Ghaziri, Ralf Galuske, Patrick R Hof, and Giorgio M Innocenti, “Evolution amplified process- ing with temporally dispersed slow neuronal connectivity in primates,” Proceedings of the National Academy of Sciences , 19551–19556 (2009). [83] Daniel Liewald, Robert Miller, Nikos Logothetis, Hans Joachim Wagner, and Almut Sch¨uz, “Distribution of axon diameters in cortical white matter: an electron- microscopic study on three human brains and a macaque,” Biological Cybernetics , 541–557 (2014). [84] V E Kravtsov, I V Lerner, and V I Yudson, “Random walks in media with constrained disorder,” Journal of Physics A: Mathematical and General , L703 (1985). [85] Sinai Y G., “The limiting behavior of a one-dimensional ran- dom walk in a random medium,” Theory Probab. Appl. , [86] Daniel S Fisher, “Random walks in random environments,” Physical Review A , 960 (1984). [87] J. A Aronovitz and Nelson. David R, “Anomalous diffusion in steady fluid flow through a porous medium,” (1984). [88] Daniel S. Fisher, Daniel Friedan, Zongan Qiu, Scott J. Shenker, and Stephen H. Shenker, “Random walks in two- dimensional random environments with constrained drift forces,” (1985). [89] Yan-Chr Tsai and Yonathan Shapir, “On the scale-invariant distribution of the diffusion coefficient for classical particles diffusing in disordered media,” Journal of Physics A: Mathe- matical and General , 39 (1993). [90] Igor V. Lerner, “Distributions of the diffusion coefficient for the quantum and classical diffusion in disordered me- dia,” Nuclear Physics, Section A , 274–292 (1993), arXiv:9402028 [cond-mat]. [91] B I Shklovskii and A L Efros, Electronic properties of doped semiconductors (Springer, Heidelberg, 1984). [92] S Havlin, M Schwartz, R B Selinger, A Bunde, and H E Stan- ley, “Universality classes for diffusion in the presence of cor- related spatial disorder,” Physical Review A , 1717–1719 (1989). [93] S Havlin and D Ben-Avraham, “Diffusion in disordered me- dia,” Advances In Physics , 187–292 (2002). [94] Harvey Scher and Elliott W Montroll, “Anomalous transit- time dispersion in amorphous solids,” Physical Review B , [95] Dmitry S. Novikov, M. Drndic, L. S. Levitov, M. A. Kast- ner, M. V. Jarosz, and M. G. Bawendi, “Levy statistics and anomalous transport in quantum-dot arrays,” Physical Review B - Condensed Matter and Materials Physics , 1–7 (2005). [96] Siying Wang, Claudia Querner, Tali Dadosh, Catherine H Crouch, Dmitry S. Novikov, and Marija Drndic, “Collective fluorescence enhancement in nanoparticle clusters.” Nature Communications , 364 (2011). [97] Prince E. Rouse, “A Theory of the Linear Viscoelastic Prop- erties of Dilute Solutions of Coiling Polymers,” The Journal of Chemical Physics , 1272–1280 (1953). [98] Xiaohu Hu, Liang Hong, Micholas Dean Smith, Thomas Neu- sius, Xiaolin Cheng, and JeremyC. Smith, “The dynamics of single protein molecules is non-equilibrium and self-similar over thirteen decades in time,” Nature Physics , 171–174 (2016). [99] I. Y. Wong, M. L. Gardel, D. R. Reichman, Eric R. Weeks, M. T. Valentine, A. R. Bausch, and D. A. Weitz, “Anomalous Diffusion Probes Microstructure Dynamics of Entangled F- Actin Networks,” Physical Review Letters , 178101 (2004). [100] Ralf Metzler, Jae-Hyung Jeon, Andrey G Cherstvy, and Eli Barkai, “Anomalous diffusion models and their proper- ties: non-stationarity, non-ergodicity, and ageing at the cente- nary of single particle tracking.” Physical chemistry chemical physics : PCCP , 24128–64 (2014). [101] Matan Mussel, Lilah Inzelberg, and Uri Nevo, “Insignif- icance of active flow for neural diffusion weighted imag- ing: A negative result,” Magnetic Resonance in Medicine , DOI:10.1002/mrm.26375 (2016). [102] Denis S Grebenkov, “NMR survey of reflected Brownian mo- tion,” Rev. Mod. Phys. , 1077–1137 (2007). [103] S. D. Stoller, W. Happer, and F. J. Dyson, “Transverse spin relaxation in inhomogeneous magnetic fields,” (1991). [104] Thomas M. de Swiet and Pabitra N. Sen, “Decay of nuclear magnetization by bounded diffusion in a constant field gradi- ent,” The Journal of Chemical Physics , 5597 (1994). [105] M.D. Hurlimann, K.G. Helmer, T.M. Deswiet, and P.N. Sen, “Spin Echoes in a Constant Gradient and in the Presence of Simple Restriction,” Journal of Magnetic Resonance, Series A , 260–264 (1995). [106] Denis S. Grebenkov, “Exploring diffusion across permeable barriers at high gradients. II. Localization regime,” Journal of Magnetic Resonance , 164–176 (2014). [107] Frederik Bernd Laun, Tristan Anselm Kuder, Wolfhard Semmler, and Bram Stieltjes, “Determination of the defin- ing boundary in nuclear magnetic resonance diffusion experi- ments,” Physical Review Letters , 2–5 (2011). [108] Tristan Anselm Kuder and Frederik Bernd Laun, “NMR- based diffusion pore imaging by double wave vector mea- surements,” Magnetic Resonance in Medicine , 836–841 (2013). [109] Daniel C. Alexander, Tim B. Dyrby, Markus Nilsson, and Hui Zhang, “Imaging brain microstructure with diffusion MRI: practicality and applications,” NMR in Biomedicine (2018), [110] Ileana O. Jelescu and Matthew D. Budde, “Design and Vali- dation of Diffusion MRI Models of White Matter,” Frontiers in Physics (2017), 10.3389/fphy.2017.00061. [111] Olivier Reynaud, “Time-Dependent Diffusion MRI in Can- cer: Tissue Modeling and Applications,” Frontiers in Physics , 1–16 (2017). [112] J¨org K¨arger, “NMR self-diffusion studies in heterogeneous systems,” Advances in Colloid and Interface Science , 129–
148 (1985). [113] Christian Beaulieu and Peter S Allen, “An in vitro evaluation of the effects of local magnetic-susceptibility-induced gradi- ents on anisotropic water diffusion in nerve,” Magnetic reso- nance in medicine , 39–44 (1996). [114] Yaniv Assaf and Yoram Cohen, “Assignment of the water slow-diffusing component in the central nervous system using q-space diffusion MRS: Implications for fiber tract imaging,” Magnetic resonance in medicine , 191–199 (2000). [115] A Bar-Shir, L Avram, E ¨Ozarslan, P J Basser, and Y Co- hen, “The effect of the diffusion time and pulse gradient dura- tion ratio on the diffraction pattern and the structural informa- tion estimated from q-space diffusion MR: Experiments and simulations,” Journal of Magnetic Resonance , 230–236 (2008). [116] Nicolas Kunz, St´ephane V. Sizonenko, Petra S. H¨uppi, Rolf Gruetter, and Yohan Van de Looij, “Investigation of field and diffusion time dependence of the diffusion-weighted signal at ultrahigh magnetic fields,” NMR in Biomedicine , 1251– [117] Peter van Gelderen, Marloes H M de Vleeschouwer, Daryl DesPres, James Pekar, Peter van Zijl, and Chrit T W Moo- nen, “Water diffusion and acute stroke,” Magnetic resonance in medicine , 154–163 (1994). [118] Chris A Clark, Maj Hedehus, and Michael E Moseley, “Dif- fusion time dependence of the apparent diffusion tensor in healthy human brain and white matter disease,” Magnetic res- onance in medicine , 1126–1129 (2001). [119] Markus Nilsson, Jimmy L¨att, Emil Nordh, Ronnie Wirestam, Freddy St r ahlberg, and Sara Brockstedt, “On the effects of a varied diffusion time in vivo: is the diffusion in white matter restricted?” Magnetic resonance imaging , 176–187 (2009). [120] Mark A Horsfield, Gareth J Barker, and W Ian McDon- ald, “Self-diffusion in CNS tissue by volume-selective proton NMR,” Magnetic resonance in medicine , 637–644 (1994). [121] Silvia De Santis, Derek K. Jones, and Alard Roebroeck, “In- cluding diffusion time dependence in the extra-axonal space improves in vivo estimates of axonal diameter and density in human white matter,” NeuroImage , 91–103 (2016). [122] M D Does, E C Parsons, and J C Gore, “Oscillating Gradi- ent Measurements of Water Diffusion in Normal and Glob- ally Ischemic Rat Brain,” Magnetic Resonance In Medicine , 206–215 (2003). [123] Junzhong Xu, Hua Li, Kevin D. Harkins, Xiaoyu Jiang, Jing- ping Xie, Hakmook Kang, Mark D. Does, and John C. Gore, “Mapping mean axon diameter and axonal volume fraction by MRI using temporal diffusion spectroscopy,” NeuroImage , 10–19 (2014). [124] Nadya Pyatigorskaya, Denis Bihan, Olivier Reynaud, and Luisa Ciobanu, “Relationship between the diffusion time and the diffusion MRI signal observed at 17.2 tesla in the healthy rat brain cortex,” Magnetic resonance in medicine , 492–
500 (2014). [125] Dan Wu, Lee J. Martin, Frances J. Northington, and
Jiangyang Zhang, “Oscillating gradient diffusion MRI reveals unique microstructural information in normal and hypoxia- ischemia injured mouse brains,” Magnetic Resonance in
Medicine , 1366–1374 (2014), arXiv:NIHMS150003. [126] Dan Wu and Jiangyang Zhang, “The Effect of Microcircula- tory Flow on Oscillating Gradient Diffusion { MRI } and Dif- fusion Encoding with Dual-Frequency Orthogonal Gradients { (DEFOG) } ,” Magnetic Resonance in Medicine , n/a—-n/a (2016). [127] Corey A Baron and Christian Beaulieu, “Oscillating gradi- ent spin-echo (OGSE) diffusion tensor imaging of the human brain,” Magnetic resonance in medicine , 726–736 (2014). [128] Anh T. Van, Samantha J. Holdsworth, and Roland Bam- mer, “In vivo investigation of restricted diffusion in the hu- man brain with optimized oscillating diffusion gradient en- coding,” Magnetic Resonance in Medicine , 83–94 (2014), arXiv:NIHMS150003. [129] J. Stepi v snik, “Analysis of NMR self-diffusion measurements by a density matrix calculation,” Physica B+C , 350–364 (1981). [130] Dmitry S. Novikov and Valerij G Kiselev, “Surface-to-volume ratio with oscillating gradients,” Journal of Magnetic Reso- nance , 141–145 (2011). [131] N. Shemesh, S. N. Jespersen, D. C. Alexander, Y. Cohen, I. Drobnjak, T. B. Dyrby, J. Finsterbusch, M. A. Koch, T. Kuder, F. Laun, M. Lawrenz, H. Lundell, P. P. Mitra, M. Nilsson, E. Ozarslan, D. Topgaard, and C. F. Westin, “Conventions and nomenclature for double diffusion encod- ing nmr and mri,” Magn Reson Med , 82–7 (2016). [132] Daniel Topgaard, “Multidimensional diffusion mri,” J Magn Reson , 98–113 (2017). [133] L D Landau and E M Lifshitz, Statistical Physics, Part I (Course of Theoretical Physics, Vol. 5) (Pergamon Press, Ox- ford, 1969). [134] K D Merboldt, W Hanicke, and J Frahm, “Self-diffusion NMR Imaging Using Stimulated Echoes,” Journal Magnetic Resonance , 479–486 (1985). [135] B Gross and R Kosfeld, “Anwendung der spin-echo-methode der messung der selbstdiffusion,” Messtechnik , 171—-177 (1969). [136] J E Tanner, “Self diffusion of water in frog muscle.” Biophys- ical journal , 107 (1979). [137] J Stepisnik, “Time-dependent self-diffusion by NMR spin- echo,” Physica B , 343–350 (1993). [138] Paul T. Callaghan and Janez Stepisnik, “Frequency-Domain Analysis of Spin Motion Using Modulated-Gradient NMR,” Journal of Magnetic Resonance, Series A , 118–122 (1995). [139] L L Latour, K Svoboda, P P Mitra, and C H Sotak, “Time- Dependent Diffusion of Water in a Biological Model Sys- tem,” Proceedings of the National Academy of Sciences of the United States of America , 1229–1233 (1994). [140] Bibek Dhital, Marco Reisert, Elias Kellner, and Valerij G Kiselev, “Intra-axonal Diffusivity in Brain White Matter,” preprint , 1–9 (2017), arXiv:1712.04565v1. [141] Antonios Papaioannou, Dmitry S Novikov, Els Fieremans, and Gregory S Boutis, “Observation of structural universal- ity in disordered systems using bulk diffusion measurement,” Physical Review E , 61101 (2017). [142] Gregory Lemberskiy, Steven H. Baete, Martijn A. Cloos, Dmitry S. Novikov, and Els Fieremans, “Validation of surface-to-volume ratio measurements derived from oscillat- ing gradient spin echo on a clinical scanner using anisotropic fiber phantoms,” NMR in Biomedicine , e3708 (2017). [143] Olivier Reynaud, Kerryanne Veronica Winters, Dung Minh Hoang, Youssef Zaim Wadghiri, Dmitry S. Novikov, and Sungheon Gene Kim, “Surface-to-volume ratio mapping of tumor microstructure using oscillating gradient diffusion weighted imaging,” Magnetic Resonance in Medicine , [144] Olivier Reynaud, Kerryanne Veronica Winters, Dung Minh Hoang, Youssef Zaim Wadghiri, Dmitry S. Novikov, and Sungheon Gene Kim, “Pulsed and oscillating gradient { MRI } for assessment of cell size and extracellular space { (POMACE) } in mouse gliomas,” NMR in Biomedicine (2016), 10.1002/nbm.3577. [145] Xiaoyu Jiang, Hua Li, Jingping Xie, Eliot T Mckinley, Ping Zhao, John C Gore, and Junzhong Xu, “In Vivo Imaging of Cancer Cell Size and Cellularity Using Temporal Diffusion Spectroscopy,” Magnenetic Resonance in Medicine , 156–
164 (2017). [146] M H Ernst, J Machta, J R Dorfman, and H van Beijeren, “Long-time tails in stationary random media. 1. { Theory } ,” Journal of Statistical Physics , 477–495 (1984). [147] P B Visscher, “Universality in disordered diffusive systems - exact fixed-points in 1,2, and 3 dimensions,” Physical Review B , 5472–5485 (1984). [148] K´evin Ginsburger, Fabrice Poupon, Justine Beaujoin, Del- phine Estournet, Felix Matuschke, Jean-Franc¸ois Mangin, Markus Axer, and Cyril Poupon, “Improving the Realism of White Matter Numerical Phantoms: A Step toward a Bet- ter Understanding of the Influence of Structural Disorders in Diffusion MRI,” Frontiers in Physics , 12 (2018). [149] Marco Palombo, Clemence Ligneul, Edwin Hernandez- Garzon, and Julien Valette, “Can we detect the ef- fect of spines and leaflets on the diffusion of brain intracellular metabolites?” NeuroImage (2017), [150] Marco Palombo, Cl´emence Ligneul, Chlo´e Najac, Juliette Le Douce, Julien Flament, Carole Escartin, Philippe Hantraye, Emmanuel Brouillet, Gilles Bonvento, and Julien Valette, “New paradigm to assess brain cell morphology by diffusion- weighted MR spectroscopy in vivo,” Proceedings of the Na- tional Academy of Sciences , 201504327 (2016). [151] Joseph J H Ackerman and Jeffrey J Neil, “The use of MR- detectable reporter molecules and ions to evaluate diffusion in normal and ischemic brain,” NMR Biomed , 725–733 (2010). [152] Gordon M G Shepherd, Morten Raastad, and Per Andersen, “General and variable features of varicosity spacing along un- myelinated axons in the hippocampus and cerebellum,” Pro- ceedings of the National Academy of Sciences , 6340–6345 (2002). [153] Gordon M G Shepherd and Morten Raastad, “Axonal vari- cosity distributions along parallel fibers: a new angle on a cerebellar circuit,” The Cerebellum , 110–113 (2003). [154] Dominique Debanne, Emilie Campanac, Andrzej Bialowas, Edmond Carlier, and Gis`ele Alcaraz, “Axon Physiology,”
Physiological Reviews , 555–602 (2011). [155] Shengxiang Zhang, Jamie Boyd, Kerry Delaney, and Timo- thy H Murphy, “Rapid reversible changes in dendritic spine structure in vivo gated by the degree of ischemia,” J Neurosci , 5333–5338 (2005). [156] Ping Li and Timothy H Murphy, “Two-photon imaging during prolonged middle cerebral artery occlusion in mice reveals recovery of dendritic structure after reperfusion,” J Neurosci , 11970–11979 (2008). [157] Markus Nilsson, Jimmy L¨att, Freddy Stahlberg, Danielle van Westen, and Hakan Hagsl¨att, “The importance of axonal un- dulation in diffusion { MR } measurements: a Monte Carlo simulation study,” NMR in Biomedicine , 795–805 (2012). [158] Aaron J Schain, Robert A Hill, and Jaime Grutzendler, “Label-free in vivo imaging of myelinated axons in health and disease with spectral confocal reflectance microscopy,” Nat Med , 443–449 (2014). [159] Corey A llan Baron, Mahesh Kate, Laura Gioia, Kenneth Butcher, Derek Emery, Matthew Budde, and Christian
Beaulieu, “Reduction of Diffusion-Weighted Imaging Con- trast of Acute Ischemic Stroke at Short Diffusion Times,”
Stroke; a journal of cerebral circulation , 2136–2141 (2015). [160] D Barazany, P J Basser, and Y Assaf, “In vivo measurement of axon diameter distribution in the corpus callosum of rat brain,” Brain , 1210–1220 (2009). [161] H Zhang, P L Hubbard, G J M Parker, and D C Alexan- der, “Axon diameter mapping in the presence of orientation dispersion with diffusion MRI,” Neuroimage , 1301–1315 (2011). [162] Susie Y Huang, Aapo Nummenmaa, Thomas Witzel, Tanguy Duval, Julien Cohen-Adad, Lawrence L Wald, and Jennifer A McNab, “The impact of gradient strength on in vivo diffusion MRI estimates of axon diameter,” Neuroimage , 464–472 (2015). [163] Assaf Horowitz, Daniel Barazany, Ido Tavor, Moran Bern- stein, Galit Yovel, and Yaniv Assaf, “In vivo correlation be- tween axon diameter and conduction velocity in the human brain,” Brain Structure and Function , 1777–1788 (2015). [164] G M Innocenti, R Caminiti, and F Aboitiz, “Comments on the paper by Horowitz et al.(2014),” Brain structure and function , 1789–1790 (2015). [165] Assaf Horowitz, Daniel Barazany, Ido Tavor, Galit Yovel, and Yaniv Assaf, “Response to the comments on the paper by Horowitz et al. (2014),” Brain Structure and Function , [166] C Beaulieu, “The basis of anisotropic water diffusion in the nervous system - a technical review,” NMR in Biomedicine , 435–455 (2002). [167] Jelle Veraart, Els Fieremans, and Dmitry S. Novikov, “Uni- versal power-law scaling of water diffusion in human brain defines what we see with { MRI } ,” preprint arXiv:1609.09145 https://ar (2016). [168] Masaya Takahashi, David B Hackney, Guixin Zhang, Suzanne L Wehrli, Alex C Wright, William T O’Brien, Hide- masa Uematsu, Felix W Wehrli, and Michael E Selzer, “Mag- netic Resonance Microimaging of Intraaxonal Water Diffu- sion in Live Excised Lamprey Spinal Cord,” Proceedings of the National Academy of Sciences of the United States of America , 16192–16196 (2002). [169] Sune Nørhøj Jespersen, Jonas Lynge Olesen, Brian Hansen, and Noam Shemesh, “Diffusion time dependence of mi- crostructural parameters in fixed spinal cord,” NeuroImage , [170] M.E. Komlosh, E. ¨Ozarslan, M.J. Lizak, I. Horkayne-Szakaly, R.Z. Freidlin, F. Horkay, and P.J. Basser, “Mapping average axon diameters in porcine spinal cord white matter and rat corpus callosum using d-PFG MRI,” NeuroImage , 210–
216 (2013). [171] Junzhong Xu, Hua Li, Ke Li, Kevin D. Harkins, Xiaoyu Jiang, Jingping Xie, Hakmook Kang, Richard D. Dortch, Adam W. Anderson, Mark D. Does, and John C. Gore, “Fast and simplified mapping of mean axon diameter using temporal diffusion spectroscopy,” NMR in Biomedicine , 400–410 (2016). [172] Tanguy Duval, Jennifer A. McNab, Kawin Setsompop, Thomas Witzel, Torben Schneider, Susie Yi Huang, Boris Keil, Eric C. Klawiter, Lawrence L. Wald, and Julien Cohen- Adad, “In vivo mapping of human spinal cord microstructure at 300mT/m,” NeuroImage , 494–507 (2015). [173] Dan Benjamini, Michal E. Komlosh, Lynne A. Holtzclaw, Uri Nevo, and Peter J. Basser, “White matter microstructure from nonparametric axon diameter distribution mapping,” Neu- roImage , 333–344 (2016). [174] Ivana Drobnjak, Bernard Siow, and Daniel C Alexander, “Optimizing gradient waveforms for microstructure sensitiv- ity in diffusion-weighted MR,” Journal of Magnetic Reso- nance , 41–51 (2010). [175] Ivana Drobnjak, Hui Zhang, Andrada Ianu??, Enrico Kaden, and Daniel C. Alexander, “PGSE, OGSE, and sensitivity to axon diameter in diffusion MRI: Insight from a simulation study,” Magnetic Resonance in Medicine , 688–700 (2016). [176] Markus Nilsson, Carl-fredrik Westin, Ivana Drobnjak, and T Daniel, “Resolution limit of cylinder diameter estimation by diffusion MRI : The impact of gradient waveform and ori- entation dispersion,” NMR Biomed , 1–13 (2017). [177] Farshid Sepehrband, Daniel C Alexander, Nyoman D Kur- niawan, David C Reutens, and Zhengyi Yang, “Towards higher sensitivity and stability of axon diameter estimation with diffusion-weighted MRI,” NMR Biomed , 293–308 (2016). [178] S Portnoy, J J Flint, S J Blackband, and G J Stanisz, “Os- cillating and Pulsed Gradient Diffusion Magnetic resonance Microscopy Over an Extended b-Value Range: Implications for the Characterization of Tissue Microstructure,” Magnetic Resonance in Medicine , 1131–1145 (2013). [179] I M Lifshitz, “Energy Spectrum of Disordered Systems,” Ad- vances In Physics , 483–536 (1964). [180] Nevill F Mott and Edward A Davies, Electronic processes in non-crystalline materials (Oxford University Press, New York, 1971). [181] Jelle Veraart, Els Fieremans, and Dmitry S. Novikov, “Dif- fusion { MRI } noise mapping using random matrix theory,” Magnetic resonance in medicine , 1582–1593 (2016). [182] Jelle Veraart, Els Fieremans, and Dmitry S. Novikov, “Uni- versal power-law scaling of water diffusion in human brain defines what we see with MRI,” preprint arXiv:1609.09145 (2016), arXiv:1609.09145. [183] Emmanuel J Candes, Justin K Romberg, and Terence Tao, “Stable signal recovery from incomplete and inaccurate mea- surements,” Communications on Pure and Applied Mathe- matics , 1207–1223 (2006). [184] Michael Lustig, David Donoho, and John M Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine , 1182–1195 (2007). [185] Dan Ma, Vikas Gulani, Nicole Seiberlich, Kecheng Liu, Jef- frey L Sunshine, Jeffrey L Duerk, and Mark A Griswold, “Magnetic resonance fingerprinting,” Nature , 187–192 (2013). [186] Milos D Ikonomovic, Eric E Abrahamson, Barbara A Isan- ski, Joanne Wuu, Elliott J Mufson, and Steven T DeKosky, “Superior frontal cortex cholinergic axon density in mild cog- nitive impairment and early Alzheimer disease,” Arch Neurol , 1312–1317 (2007). [187] Bj¨orn Lampinen, Filip Szczepankiewicz, Danielle van Westen, Elisabet Englund, Pia C Sundgren, Jimmy L¨att,
Freddy St r ahlberg, and Markus Nilsson, “Optimal experi- mental design for filter exchange imaging: { A } pparent ex- change rate measurements in the healthy brain and in intracra- nial tumors,” Magn Reson Med (2016), 10.1002/mrm.26195. [188] J´erˆome Badaut, Stephen Ashwal, Arash Adami, Beatriz Tone, Rebecca Recker, David Spagnoli, B´eatrice Ternon, and An- dre Obenaus, “Brain Water Mobility Decreases after Astro- cytic Aquaporin-4 Inhibition Using RNA Interference,” Jour- nal of Cerebral Blood Flow and Metabolism , 819–831 (2011). [189] Donghan M Yang, James E Huettner, G Larry Bretthorst, Jeffrey J Neil, Joel R Garbow, and Joseph J.H. Acker- man, “Intracellular water preexchange lifetime in neurons and astrocytes,” Magnetic Resonance in Medicine (2017), [190] Ruiliang Bai, Charles S. Springer, Dietmar Plenz, and Peter J. Basser, “Fast, Na < sup > + < /sup > /K < sup > + < /sup > pump driven, steady-state transcytolemmal water exchange in neu- ronal tissue: A study of rat brain cortical cultures,” Magnetic Resonance in Medicine , 3207–3217 (2018). [191] Jens H. Jensen, G. Russell Glenn, and Joseph A. Helpern, “Fiber ball imaging,” NeuroImage , 824–833 (2016). [192] Emilie T McKinnon, Jens H Jensen, G Russell Glenn, and Joseph A Helpern, “Dependence on b-value of the direction- averaged diffusion-weighted imaging signal in brain,” Mag- netic resonance imaging , 121–127 (2017). [193] Itamar Ronen, Matthew Budde, Ece Ercan, Jacopo Annese, Aranee Techawiboonwong, and Andrew Webb, “Microstruc- tural organization of axons in the human corpus callosum quantified by diffusion-weighted magnetic resonance spec- troscopy of N-acetylaspartate and post-mortem histology,” Brain Struct Funct , 1773–1785 (2014). [194] Matthew D Budde and Joseph A Frank, “Examining brain microstructure using structure tensor analysis of histological sections,” NeuroImage , 1–10 (2012). [195] Bibek Dhital, Elias Kellner, Valerij G. Kiselev, and Marco Reisert, “The absence of restricted water pool in brain white matter,” NeuroImage (2017), [196] Dennis M. Healy, Harrie Hendriks, and Peter T. Pt Kim, “Spherical deconvolution,” Journal of Multivariate Analysis , 1–22 (1998). [197] J. Donald Tournier, Fernando Calamante, David G. Gadian, and Alan Connelly, “Direct estimation of the fiber orien- tation density function from diffusion-weighted MRI data using spherical deconvolution,” NeuroImage , 1176–1185 (2004). [198] Adam W. Anderson, “Measurement of fiber orientation dis- tributions using high angular resolution diffusion imaging,” Magnetic Resonance in Medicine , 1194–1206 (2005). [199] J-Donald Tournier, Fernando Calamante, and Alan Connelly, “Robust determination of the fibre orientation distribution in diffusion { MRI } : non-negativity constrained super-resolved spherical deconvolution,” Neuroimage , 1459–1472 (2007). [200] F Dell’Acqua, G Rizzo, P Scifo, R A Clarke, G Scotti, and F Fazio, “A Model-Based Deconvolution Approach to Solve Fiber Crossing in Diffusion-Weighted MR Imaging,” (2007). [201] B Jian and B C Vemuri, “A Unified Computational Frame- work for Deconvolution to Reconstruct Multiple Fibers From Diffusion Weighted MRI,” (2007). [202] Enrico Kaden, Thomas R. Kn¨osche, and Alfred Anwan- der, “Parametric spherical deconvolution: Inferring anatom- ical connectivity using diffusion MR imaging,” NeuroImage , 474–488 (2007). [203] Nathan S. White and Anders M. Dale, “Optimal diffusion MRI acquisition for fiber orientation density estimation: An analytic approach,” Human Brain Mapping , 3696–3703 (2009). [204] Alex Mackay, Kenneth Whittall, Julian Adler, David Li, Don- ald Paty, and Douglas Graeb, “In vivo visualization of myelin water in brain by magnetic resonance,” Magnetic Resonance in Medicine , 673–677 (1994). [205] Richard D. Dortch, Kevin D. Harkins, Meher R. Juttukonda, John C. Gore, and Mark D. Does, “Characterizing inter- compartmental water exchange in myelinated tissue using relaxation exchange spectroscopy,” Magnetic Resonance in Medicine , 1450–1459 (2013). [206] Eleftheria Panagiotaki, Torben Schneider, Bernard Siow, Matt G Hall, Mark F Lythgoe, and Daniel C Alexander, “Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison,” Neuroimage , [207] Uran Ferizi, Torben Schneider, Thomas Witzel, Lawrence L Wald, Hui Zhang, Claudia A M Wheeler-Kingshott, and Daniel C Alexander, “White matter compartment models for in vivo diffusion MRI at 300mT/m,” Neuroimage , 468–
483 (2015). [208] Yaniv Assaf and Peter J. Basser, “Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain,” NeuroImage , 48–58 (2005). [209] Silvia De Santis, Daniel Barazany, Derek K. Jones, and Yaniv Assaf, “Resolving relaxometry and diffusion proper- ties within the same voxel in the presence of crossing fibres by combining inversion recovery and diffusion-weighted ac- quisitions,” Magnetic Resonance in Medicine , 372–380 (2016). [210] Anthony J. Sherbondy, Matthew C. Rowe, and Daniel C. Alexander, “MicroTrack: An algorithm for concurrent projec- tome and microstructure estimation,” Lecture Notes in Com- puter Science (including subseries Lecture Notes in Artifi- cial Intelligence and Lecture Notes in Bioinformatics) LNCS , 183–190 (2010). [211] Jacques-Donald JD Tournier, Susumu Mori, and Alexander Leemans, “Diffusion tensor imaging and beyond.” Magnetic Resonance in Medicine , 1532–56 (2011). [212] Gabriel Girard, Alessandro Daducci, Laurent Petit, Jean- Philippe Thiran, Kevin Whittingstall, Rachid Deriche, Demian Wassermann, and Maxime Descoteaux, “AxTract: Toward microstructure informed tractography,” Human Brain Mapping , 5485–5500 (2017). [213] T. E J Behrens, H. Johansen Berg, S. Jbabdi, M. F S Rush- worth, and M. W. Woolrich, “Probabilistic diffusion tractog- raphy with multiple fibre orientations: What can we gain?” NeuroImage , 144–155 (2007). [214] Maxime Descoteaux, Rachid Deriche, Thomas R. Kn¨osche, and Alfred Anwander, “Deterministic and probabilistic trac- tography based on complex fibre orientation distributions,” IEEE Transactions on Medical Imaging , 269–286 (2009). [215] Shawna Farquharson, J-Donald Tournier, Fernando Cala- mante, Gavin Fabinyi, Michal Schneider-Kolsky, Graeme D Jackson, and Alan Connelly, “White matter fiber tractogra- phy: why we need to move beyond DTI.” Journal of neuro- surgery , 1367–77 (2013). [216] Stamatios N. Sotiropoulos, Saad Jbabdi, Junqian Xu, Jes- per L. Andersson, Steen Moeller, Edward J. Auerbach,
Matthew F. Glasser, Moises Hernandez, Guillermo Sapiro,
Mark Jenkinson, David A. Feinberg, Essa Yacoub, Christophe
Lenglet, David C. Van Essen, Kamil Ugurbil, and Timo- thy E.J. Behrens, “Advances in diffusion MRI acquisition and processing in the Human Connectome Project,” NeuroImage , 125–143 (2013), arXiv:NIHMS150003. [217] Bryce Wilkins, Namgyun Lee, Niharika Gajawelli, Meng Law, and Natasha Lepore, “Fiber estimation and tractogra- phy in diffusion MRI: Development of simulated brain im- ages and comparison of multi-fiber analysis methods at clini- cal b-values,” NeuroImage , 341–356 (2015). [218] Els Fieremans, Jens H Jensen, Joseph A Helpern, Sungheon
Kim, Robert I Grossman, Matilde Inglese, and Dmitry S.
Novikov, “Diffusion Distinguishes Between Axonal Loss and
Demyelination in Brain White Matter,” Proceedings of the In- ternational Society of Magnetic Resonance in Medicine 20, p.714 (2012). [219] Dmitry S. Novikov and Els Fieremans, “Relating Extracel- lular Diffusivity to Cell Size Distribution and Packing Den- sity as Applied to White Matter,” Proceedings of the Interna- tional Society of Magnetic Resonance in Medicine 20, p.1829 (2012). [220] Ileana O Jelescu, Magdalena Zurek, Kerryanne V Winters, Jelle Veraart, Anjali Rajaratnam, Nathanael S Kim, James S Babb, Timothy M Shepherd, Dmitry S. Novikov, Sungheon G Kim, and Els Fieremans, “In vivo quantification of de- myelination and recovery using compartment-specific diffu- sion MRI metrics validated by electron microscopy,” Neu- roimage , 104–114 (2016). [221] Ahmad Raza Khan, Andrey Chuhutin, Ove Wiborg, Christo- pher D Kroenke, Jens R Nyengaard, Brian Hansen, and Sune Nørhøj Jespersen, “Biophysical modeling of high field diffusion MRI demonstrates micro-structural aberra- tion in chronic mild stress rat brain,” NeuroImage (2016), [222] Ahmad Raza Khan, Andrey Chuhutin, Ove Wiborg, Christo- pher D Kroenke, Jens R Nyengaard, Brian Hansen, and Sune Nørhøj Jespersen, “Summary of high field diffusion MRI and microscopy data demonstrate microstructural aber- ration in chronic mild stress rat brain,” Data in brief , 934 (2016). [223] Peter Vestergaard-Poulsen, Gregers Wegener, Brian Hansen, Carsten R Bjarkam, Stephen J Blackband, Niels Christian Nielsen, and Sune N Jespersen, “Diffusion-weighted MRI and quantitative biophysical modeling of hippocampal neu- rite loss in chronic stress,” P L o S One , e20653 (2011). [224] Matthew D Budde and Joseph A Frank, “Neurite beading is sufficient to decrease the apparent diffusion coef fi cient af- ter ischemic stroke,” Proc Natl Acad Sci U S A , 14472– [225] A. W. Unterberg, J. Stover, B. Kress, and K. L. Kiening, “Edema and brain trauma,” Neuroscience , 1021–1029 (2004). [226] Nikola Stikov, Jennifer S W Campbell, Thomas Stroh, Mariette Lavel´ee, Stephen Frey, Jennifer Novek, Stephen Nuara, Ming-Kai Ho, Barry J Bedell, Robert F Dougherty, Ilana R Leppert, Mathieu Boudreau, Sridar Narayanan, Tan- guy Duval, Julien Cohen-Adad, Paul-Alexandre Picard, Al- icja Gasecka, Daniel Cˆot´e, and G Bruce Pike, “In vivo histol- ogy of the myelin g-ratio with magnetic resonance imaging,” Neuroimage , 397–405 (2015). [227] Sune Nørhøj Jespersen, “Equivalence of double and single wave vector diffusion contrast at low diffusion weighting,” NMR Biomed , 813–818 (2012). [228] B Hansen, A.˜R. Khan, N Shemesh, T.˜E. Lund, R Sangill, L Ostergaard, and S.˜N. Jespersen, “White matter biomarkers from fast protocols using axially symmetric diffusion kurtosis imaging,” preprint arXiv:1610.02783 https://ar (2016). [229] Brian Hansen, Noam Shemesh, and Sune Nørhøj Jespersen, “Fast imaging of mean, axial and radial diffusion kurtosis,” NeuroImage http://dx. (2016). [230] Hang Tuan Nguyen, Denis Grebenkov, Dang Van Nguyen, Cyril Poupon, Denis Le Bihan, and Jing-Rebecca Li, “Param- eter estimation using macroscopic diffusion MRI signal mod- els,” Physics in Medicine and Biology , 3389–3413 (2015). [231] Ileana O Jelescu, Jelle Veraart, Vitria Adisetiyo, Sarah S Milla, Dmitry S. Novikov, and Els Fieremans, “One diffu- sion acquisition and different white matter models: How does microstructure change in human early development based on { WMTI } and NODDI?” NeuroImage , 242–256 (2015). [232] Edward S Hui, Els Fieremans, Jens H Jensen, Ali Tabesh, Wuwei Feng, Leonardo Bonilha, Maria V Spampinato, Robert Adams, and Joseph A Helpern, “Stroke assessment with diffusional kurtosis imaging,” Stroke , 2968–2973 (2012). [233] Nathaniel D Kelm, Kathryn L West, Robert P Carson, Daniel F Gochberg, Kevin C Ess, and Mark D Does, “Evalu- ation of diffusion kurtosis imaging in ex vivo hypomyelinated mouse brains,” NeuroImage , 612–626 (2016). [234] Maria F Falangola, David N Guilfoyle, Ali Tabesh, Edward S Hui, Xingju Nie, Jens H Jensen, Scott V Gerum, Caixia Hu, John LaFrancois, Heather R Collins, and Joseph A Helpern, “Histological correlation of diffusional kurtosis and white matter modeling metrics in cuprizone-induced corpus callosum demyelination,” NMR in Biomedicine , 948–957 (2014). [235] C Guglielmetti, J Veraart, E Roelant, Z Mai, J Daans, J Van Audekerke, M Naeyaert, G Vanhoutte, R Delgado y Palacios, J Praet, E Fieremans, P Ponsaerts, J Sijbers, A Van der Lin- den, and M Verhoye, “Diffusion kurtosis imaging probes cor- tical alterations and white matter pathology following cupri- zone induced demyelination and spontaneous remyelination,” NeuroImage , 363–377 (2016). [236] Elan J Grossman, Ivan I Kirov, Oded Gonen, Dmitry S. Novikov, Matthew S Davitz, Yvonne W Lui, Robert I Gross- man, Matilde Inglese, and Els Fieremans, “N-acetyl-aspartate levels correlate with intra-axonal compartment parameters from diffusion { MRI } ,” NeuroImage , 334–343 (2015). [237] Santiago Coelho, Leandro Beltrachini, Jose M. Pozo1, and Alejandro F. Frangi, “Double diffusion encoding vs single dif- fusion encoding parameter estimation of biophysical models in diffusion weighted mri,” in Proc. Int. Soc. Magn. Reson. Med. , Vol. 25 (2017) p. 3383. [238] Els Fieremans, Jelle Veraart, Benjamin Ades-Aron, Filip Szczepankiewicz, Markus Nilsson, and Dmitry S Novikov, “Effect of combining linear with spherical tensor encoding on estimating brain microstructural parameters,” Proc. 26th Annual Meeting ISMRM, Paris, France (2018). [239] Andrey Chuhutin, Brian Hansen, and Sune Nørhøj Jespersen, “Precision and accuracy of diffusion kurtosis estimation and the influence of b-value selection,” NMR in Biomedicine , [240] Maira Tariq, Torben Schneider, Daniel C Alexander, Claudia A Gandini Wheeler-Kingshott, and Hui Zhang, “Bingham–
NODDI: Mapping anisotropic orientation dispersion of neu- rites using diffusion MRI,” NeuroImage (2016). [241] A Szafer, J Zhong, and J C Gore, “Theoretical model for water diffusion in tissues,” Magn Reson Med , 697–712 (1995). [242] Maria Tariq, Michiel Kleinnijenhuis, Anne-Marie van Cap- pellen van Walsum, and Hui Zhang, “Validation of NODDI estimation of dispersion anisotropy in V1 of the human neo- cortex,” Proceedings of the International Society of Magnetic Resonance in Medicine 23, p. 0477 (2015). [243] Francesco Grussu, Torben Schneider, Carmen Tur, Richard L.
Yates, Mohamed Tachrount, Andrada Ianu, Marios C. Yian- nakas, Jia Newcombe, Hui Zhang, Daniel C. Alexander,
Gabriele C. DeLuca, and Claudia A. M. Gandini Wheeler-
Kingshott, “Neurite dispersion: a new marker of multiple sclerosis spinal cord pathology?” Annals of Clinical and
Translational Neurology , 663–679 (2017). [244] Kurt G. Schilling, Vaibhav Janve, Yurui Gao, Iwona Step- niewska, Bennett A. Landman, and Adam W. Ander- son, “Histological validation of diffusion MRI fiber orienta- tion distributions and dispersion,” NeuroImage , 200–221 (2018). [245] Michael Kazhdan, Thomas Funkhouser, and Szymon Rusinkiewicz, “Rotation invariant spherical harmonic repre- sentation of 3 d shape descriptors,” in Symposium on geome- try processing , Vol. 6 (2003) pp. 156–164. [246] H Mirzaalian, L Ning, P Savadjiev, O Pasternak, S Bouix, O Michailovich, G Grant, C E Marx, R A Morey, L A Flash- man, M S George, T W Mcallister, N Andaluz, L Shutter, R Coimbra, R D Zafonte, M J Coleman, M Kubicki, C F Westin, M B Stein, M E Shenton, and Y Rathi, “Inter-site and inter-scanner diffusion MRI data harmonization,” Neu- roImage , 311–323 (2016). [247] Sune Nørhøj Jespersen, Henrik Lundell, Casper Kaae Sønderby, and Tim B. Dyrby, “Orientationally invariant met- rics of apparent compartment eccentricity from double pulsed field gradient diffusion experiments,” NMR in Biomedicine , 1647–1662 (2013). [248] Samo Lasic, Filip Szczepankiewicz, Stefanie Eriksson, Markus Nilsson, and Daniel Topgaard, “Microanisotropy imaging: quantification of microscopic diffusion anisotropy and orientational order parameter by diffusion MRI with magic-angle spinning of the q-vector,” Frontiers in Physics , 1–14 (2014). [249] Enrico Kaden, Frithjof Kruggel, and Daniel C Alexander, “Quantitative mapping of the per-axon diffusion coefficients in brain white matter,” Magnetic Resonance in Medicine , [250] Brian Hansen, Torben Ellegaard Lund, Ryan Sangill, and Sune Nørhøj Jespersen, “Neurite density from an isotropic diffusion model,” Proceedings of the International Society of Magnetic Resonance in Medicine 23, p. 3043 (2015). [251] Enrico Kaden, Nathaniel D Kelm, Robert P Carson, Mark D Does, and Daniel C Alexander, “Multi-compartment micro- scopic diffusion imaging,” NeuroImage , 346–359 (2016). [252] Filip Szczepankiewicz, Danielle van Westen, Elisabet En- glund, Carl-Fredrik Westin, Freddy St r ahlberg, Jimmy L¨att, Pia C Sundgren, and Markus Nilsson, “The link between diffusion { MRI } and tumor heterogeneity: Mapping cell ec- centricity and density by diffusional variance decomposition { (DIVIDE) } ,” NeuroImage , 522–532 (2016). [253] A Bohr and Ben R Mottelson, Nuclear Structure, Vol. II: Nu- clear Deformations (World Scientific, 1998). [254] Filip Szczepankiewicz, Samo Lasi v c, Danielle van Westen, Pia C. Sundgren, Elisabet Englund, Carl Fredrik Westin, Freddy St r ahlberg, Jimmy L¨att, Daniel Topgaard, and Markus Nilsson, “Quantification of microscopic diffusion anisotropy disentangles effects of orientation dispersion from microstructure: Applications in healthy volunteers and in brain tumors,” NeuroImage , 241–252 (2015). [255] Joao P de Almeida Martins and Daniel Topgaard, “Two- Dimensional Correlation of Isotropic and Directional Diffu- sion Using { NMR } ,” Phys. Rev. Lett. , 87601 (2016). [256] Nathan P Skinner, Shekar N Kurpad, Brian D Schmit, L Tu- gan Muftuler, and Matthew D Budde, “Rapid In Vivo De- tection of Rat Spinal Cord Injury With Double-Diffusion- Encoded Magnetic Resonance Spectroscopy,” Magnenetic Resonance in Medicine , 1639–1649 (2017). [257] Bjorn Lampinen, Filip Szczepankiewicz, Johan Martensson, Danielle van Westen, Pia C. Sundgren, and Markus Nils- son, “Neurite density imaging versus imaging of micro- scopic anisotropy in diffusion MRI: A model comparison us- ing spherical tensor encoding,” NeuroImage , 517–531 (2017). [258] Bibek Dhital, Elias Kellner, Marco Reisert, and Valerij G Kiselev, “Isotropic diffusion weighting provides insight on diffusion compartments in human brain white matter in vivo,” Proceedings of the International Society of Magnetic Reso- nance in Medicine 23, p. 2788. (2015). [259] Hong-hsi Lee, Els Fieremans, and Dmitry S Novikov, “LEMONADE(t): Exact relation of time-dependent diffusion signal moments to neuronal microstructure,” Proc Intl Soc Mag Reson Med (2018). [260] Y. Cheng and D. G. Cory, “Multiple scattering by nmr,” J. Am. Chem. Soc. , 7935–7936 (1999). [261] D. G. Cory, A. N. Garroway, and J. B. Miller, “Applications of spin transport as a probe of local geometry,” Polym. Prep. Am. Chem. Soc. Div. Polym. Chem. , 149 (1990). [262] P. P. Mitra, L. L. Latour, R. L. Kleinberg, and C. H. Sotak, “Pulsed-field-gradient nmr measurements of restricted diffu- sion and the return-to-the-origin probability,” J. Magn. Reson. Series A , 47–58 (1995). [263] Sune Nørhøj Jespersen, “Equivalence of double and single wave vector diffusion contrast at low diffusion weighting,” NMR Biomed , 813–818 (2012). [264] J. Finsterbusch, “Extension of the double-wave-vector diffusion-weighting experiment to multiple concatenations,” J. Magn. Reson. , 174–182 (2009). [265] J. Finsterbusch, “Numerical simulations of short-mixing-time double-wave-vector diffusion-weighting experiments with multiple concatenations on whole-body mr systems,” Journal of Magnetic Resonance , 274–282 (2010). [266] Sune Nørhøj Jespersen and Niels Buhl, “The displacement correlation tensor: microstructure, ensemble anisotropy and curving fibers.” Journal of magnetic resonance (San Diego, Calif. : 1997) , 34–43 (2011). [267] Martin A Koch and J¨urgen Finsterbusch, “Numerical simu- lation of double-wave vector experiments investigating diffu- sion in randomly oriented ellipsoidal pores.” Magnetic reso- nance in medicine , 247–254 (2009). [268] E. Ozarslan and P. J. Basser, “Microscopic anisotropy re- vealed by nmr double pulsed field gradient experiments with arbitrary timing parameters,” J Chem Phys , 154511 (2008). [269] E. Ozarslan, N. Shemesh, and P. J. Basser, “A general frame- work to quantify the effect of restricted diffusion on the nmr signal with applications to double pulsed field gradient nmr experiments,” J Chem Phys , 104702 (2009). [270] Noam Shemesh, Evren Ozarslan, Peter J Basser, and Yoram Cohen, “Measuring small compartmental dimensions with low-q angular double-pgse nmr: The effect of experimental parameters on signal decay.” Journal of magnetic resonance (San Diego, Calif. : 1997) , 15–23 (2009). [271] M. A. Koch and J. Finsterbusch, “Compartment size estima- tion with double wave vector diffusion-weighted imaging,”
Magn. Reson. Med. , 90–101 (2008). [272] M. E. Komlosh, E. ¨Ozarslan, M. J. Lizak, F. Horkay, V. Schram, N. Shemesh, Y. Cohen, and P. J. Basser, “Pore diameter mapping using double pulsed-field gradient mri and its validation using a novel glass capillary array phantom,” J
Magn Reson , 128–35 (2011). [273] E. Ozarslan, M. E. Komlosh, M. J. Lizak, F. Horkay, and
P. J. Basser, “Double pulsed field gradient (double-pfg) mr imaging (mri) as a means to measure the size of plant cells,”
Magn Reson Chem
49 Suppl 1 , S79–84 (2011). [274] N. Shemesh, E. Ozarslan, P. J. Basser, and Y. Cohen, “Accu- rate noninvasive measurement of cell size and compartment shape anisotropy in yeast cells using double-pulsed field gra- dient mr,” NMR Biomed , 236–46 (2012). [275] N. Shemesh, E. ¨Ozarslan, A. Bar-Shir, P. J. Basser, and Y. Co- hen, “Observation of restricted diffusion in the presence of a free diffusion compartment: single- and double-pfg experi- ments,” J. Magn. Reson. , 214–225 (2009). [276] Marco Lawrenz and J¨urgen Finsterbusch, “Double-wave- vector diffusion-weighted imaging reveals microscopic diffu- sion anisotropy in the living human brain,” Magn Reson Med , 1072–82 (2013). [277] T. Weber, C. H. Ziener, T. Kampf, V. Herold, W. R. Bauer, and P. M. Jakob, “Measurement of apparent cell radii using a multiple wave vector diffusion experiment,” Magn. Reson. Med. , 1001–1006 (2009). [278] A. V. Avram, E. Ozarslan, J. E. Sarlls, and P. J. Basser, “In vivo detection of microscopic anisotropy using quadruple pulsed-field gradient (qpfg) diffusion mri on a clinical scan- ner,” Neuroimage , 229–39 (2013). [279] M. A. Koch and J. Finsterbusch, “Towards compartment size estimation in vivo based on double wave vector diffusion weighting,” NMR Biomed , 1422–32 (2011). [280] M. Lawrenz and J. Finsterbusch, “Mapping measures of mi- croscopic diffusion anisotropy in human brain white matter in vivo with double-wave-vector diffusion-weighted imaging,” Magn Reson Med , 773–83 (2015). [281] S. N. Jespersen, “Comment on ”measuring small compart- ments with relatively weak gradients by angular double- pulsed-field-gradient nmr” by morozov bar, sochen, and co- hen,” Magn Reson Imaging , 1643–4 (2013). [282] J. H. Jensen, “Sufficiency of diffusion tensor in characterizing the diffusion mri signal to leading order in diffusion weight- ing,” NMR Biomed , 1005–7 (2014). [283] F Frøhlich A., L Ostergaard, and G Kiselev V., “Effect of Impermeable Interfaces on Apparent Diffusion Coefficient in Heterogeneous Media,” App. Magn. Reson. , 123–137 (2005). [284] Astrid F Frøhlich, Sune N Jespersen, Leif Ostergaard, and Valerij G Kiselev, “The effect of impermeable boundaries of arbitrary geometry on the apparent diffusion coefficient,” Journal of Magnetic Resonance , 128–135 (2008). [285] N. Shemesh, E. Ozarslan, T. Adiri, P. J. Basser, and Y. Co- hen, “Noninvasive bipolar double-pulsed-field-gradient nmr reveals signatures for pore size and shape in polydisperse, ran- domly oriented, inhomogeneous porous media,” J Chem Phys , 044705 (2010). [286] Noam Shemesh, Tal Adiri, and Yoram Cohen, “Probing mi- croscopic architecture of opaque heterogeneous systems us- ing double-pulsed-field-gradient nmr,” Journal of the Ameri- can Chemical Society , 6028–6035 (2011). [287] S. N. Jespersen, H. Lundell, C. K. Sonderby, and T. B. Dyrby, “Orientationally invariant metrics of apparent compartment eccentricity from double pulsed field gradient diffusion ex- periments,” NMR Biomed , 1647–62 (2013). [288] Noam Shemesh and Yoram Cohen, “Microscopic and com- partment shape anisotropies in gray and white matter revealed by angular bipolar double-PFG MR,” Magnetic resonance in medicine , 1216–1227 (2011). [289] N. Shemesh, D. Barazany, O. Sadan, L. Bar, Y. Zur, Y. Barhum, N. Sochen, D. Offen, Y. Assaf, and Y. Co- hen, “Mapping apparent eccentricity and residual ensemble anisotropy in the gray matter using angular double-pulsed- field-gradient mri,” Magn Reson Med , 794–806 (2012). [290] Grant Yang, Qiyuan Tian, Christoph Leuze, Max Winter- mark, and Jennifer A McNab, “Double diffusion encoding mri for the clinic.” Magnetic resonance in medicine (2017), [291] P. P. Mitra, “Multiple wave-vector extensions of the nmr pulsed-field-gradient spin-echo diffusion measurement,” Phys. Rev. B , 15074–15078 (1995). [292] M. Lawrenz, M. A. Koch, and J. Finsterbusch, “A tensor model and measures of microscopic anisotropy for double- wave-vector diffusion-weighting experiments with long mix- ing times,” J Magn Reson , 43–56 (2010). [293] H. Lundell, T. B. Dyrby, P. L. Hubbard, Feng-Lei Zhou, G. J. Parker, and S.N. Jespersen, “Validation of double diffusion schemes of microscopic fractional anisotropy,” in Proc. Int. Soc. Magn. Reson. Med. , Vol. 23 (2015) p. 155. [294] E. ¨Ozarslan, “Compartment shape anisotropy (csa) revealed by double pulsed field gradient mr,” J. Magn. Reson. , 56–
67 (2009). [295] Sune Nørhøj Jespersen, Henrik Lundell, Casper Kaae Sønderby, and Tim B. Dyrby, “Commentary on ”mi- croanisotropy imaging: quantification of microscopic diffu- sion anisotropy and orientation of order parameter by diffu- sion mri with magic-angle spinning of the q-vector”,” Fron- tiers in Physics (2014), 10.3389/fphy.2014.00028. [296] S. N. Jespersen, L. A. Leigland, A. Cornea, and C. D. Kroenke, “Determination of axonal and dendritic orientation distributions within the developing cerebral cortex by diffu- sion tensor imaging,” Ieee Transactions on Medical Imaging , 16–32 (2012). [297] Sune Nørhøj Jespersen, “White matter biomarkers from dif- fusion mri,” Journal of magnetic resonance (San Diego, Calif. : 1997) (2018). [298] Stefanie Eriksson, Samo Lasic, and Daniel Topgaard, “Isotropic diffusion weighting in pgse nmr by magic-angle spinning of the q-vector.” Journal of magnetic resonance (San Diego, Calif. : 1997) , 13–18 (2013). [299] Carl-Fredrik Westin, Hans Knutsson, Ofer Pasternak, Filip Szczepankiewicz, Evren ¨Ozarslan, Danielle van Westen, Ce- cilia Mattisson, Mats Bogren, Lauren J O’Donnell, Marek Kubicki, Daniel Topgaard, and Markus Nilsson, “Q-space trajectory imaging for multidimensional diffusion mri of the human brain.” NeuroImage , 345–362 (2016). [300] Eric Betzig, “Nobel Lecture: Single molecules, cells, and super-resolution optics,” Reviews of Modern Physics , [301] R Kubo, “The fluctuation-dissipation theorem,” Rep. Prog. Phys. , 255–284 (1966). [302] Alexander L Sukstanskii, “Exact analytical results for { ADC } with oscillating diffusion sensitizing gradients,” Journal of Magnetic Resonance , 135–140 (2013). [303] Keh-Jim Dunn and David J Bergman, “Self diffusion of nu- clear spins in a porous medium with a periodic microstruc- ture,” The Journal of Chemical Physics , 3041–3054 (1995). [304] Junzhong Xu, Mark D. Does, and John C. Gore, “Quanti- tative characterization of tissue microstructure with temporal diffusion spectroscopy,” Journal of Magnetic Resonance , [305] Janez Stepisnik, Samo Lasic, Ales Mohoric, Igor Sersa, and Ana Sepe, “Velocity autocorrelation spectra of fluid in porous media measured by the CPMG sequence and constant mag- netic field gradient,” Magnetic Resonance Imaging , 517–
520 (2007). Appendix A: Causality and analytical properties in the frequency domain As stated after Eq. (1.4), the diffusion propagator G t ; r , r is the density of particles released at a point r at the moment t = 0 . While the diffusion equation is often formulated for positive times only, t > , with an explicit initial condition, writing Eq. (1.4) for all times t is more convenient due to reasons that will become clear below. Being the response to an instant point source of particles, the diffusion propa- gator helps finding the particle density for an arbitrary source f ( t, r ) (“particle injection”), using the linearity of diffusion equation, [ ∂ t − ∂ r D ( r ) ∂ r ] ρ ( t, r ) = f ( t, r ) . (A1) The solution takes the form of a t - and r -convolution, ρ ( t, r ) = (cid:90) d t d d r G t − t ; r , r f ( t , r ) , (A2) which is straightforward to prove by acting with the bracketed operator from Eq. (A1) on G t ; r , r under the integral. This so- lution explains the notion of causality that implies that the response ρ ( t, r ) follows the source, f ( t , r ) , and it cannot precede it . This means that G t − t ; r , r ≡ for t < t . This is guaranteed by the proportionality of G t ; r , r to the step func- tion θ ( t ) , as stated after Eq. (1.5). Synonymous to causality is the notion that G t ; r , r is a retarded propagator, which implies that any perturbation, f , of the system propagates into the future, which is opposed to the formally possible advanced propagator for which perturbations propagate into the past. Likewise, quantities D ( t ) , D inst ( t ) and D ( t ) , entering Eqs. (2.10)–(2.14), are retarded, since they are identically zero for t < . In particular, the retarded velocity autocorre- lator (2.6) has a physical meaning of a response of the current (sometimes called flux) J ( t, r ) of diffusing particles to that of a lump of particle density ρ ( t, r ) (the generalized Fick’s law), cf. ref. [24], J ( t, r ) = − (cid:90) d t D ( t − t ) ∂ r ρ ( t , r ) + O ( ∂ r ρ ) . (A3) Equivalently, in the Fourier domain, the convolution becomes a multiplication, cf. Eq. (2.8): J ω, r = −D ( ω ) ∂ r ρ ω, r + O ( ∂ r ρ ) . (A4) In physics, the above retarded response functions are known as particular cases of the general linear response theory and the fluctuation-dissipation theorem [133, 301]. Our main goal in this Appendix is to investigate how causality, i.e., the retarded character of any response func- tion that identically vanishes for t < , manifests itself in the frequency domain. We will now show that causality imposes a strict constraint on the analytical properties of its Fourier transform in the complex plane of ω . Namely, a retarded response must be an analytic function, i.e., it must have no singularities (e.g., poles or branch points), in the upper half- plane Im ω ≥ , Fig. 13 [133]. x xx FIG. 13. Analytical structure of a causal (retarded) response func-tion on the complex plane of ω . When calculating the inverse Fouriertransform such as Eq. (A5), the original integration contour over thereal axis can be closed in the infinite semicircle with Im ω > (lightblue dashed line) when t < , according to the Jordan’s lemma.Causality then requires that no singularities are present in the upperhalf of the complex plane, in which case the integration contour canbe shrunk to a point. For t > , the contour can be closed whereIm ω < (light red dashed line). This contour can be shrunk to en-circle the singularities of the transformed function (red solid lines).The shown examples include a simple pole as in Eq. (A7), two otherpoles and a branch cut along a part of the imaginary axis, to illustratea few typical options. To show that, consider the inverse Fourier transform back to the time domain (having D ( ω ) as an example): D ( t ) = (cid:90) d ω π e − iωt D ( ω ) , (A5) and demand that the resulting D ( t ) ≡ for t < . The integration in Eq. (A5) is performed along the Re ω axis of the complex plane of ω (i.e., over all the frequen- cies). Because of the Fourier exponential, one has to close the integration contour along an infinite semi-circle on which e − iωt → (Jordan’s lemma), and then shrink this contour to single out contributions of all singularities, according to Cauchy’s theorem. Equivalently, we can view the Fourier integration as proceeding along the equator of the Riemann sphere (topologically equivalent to the complex plane with an added point at infinity); in this case, the fact that the Re ω axis corresponds to a closed contour is more obvious. Cauchy the- orem again applies, and the Jordan lemma dictates in which hemisphere — top or bottom — of the Riemann sphere the contour should be shrunk from the equator to encompass the singularities. For negative times, t < , e − iωt diverges when Im ω → −∞ , and vanishes when Im ω → + ∞ , which dictates clos- ing the integration contour in the upper half of the ω plane, In Eq. (A5), the sign in the exponential e − iωt is chosen in the traditionof physics; the opposite sign would invoke the interchanging of the upperand lower halves of the ω plane. in the upper half-plane, in which case the integration contour is constricted to a point, yielding D ( t ) | t< ≡ . All singular- ities of D ( ω ) must then be present in the lower half of the ω plane, where the contour is closed for t > . The presence of singularities is necessary, since a function without any singu- larities on a Riemann sphere is a constant. Let us illustrate the above general considerations by an- alyzing the analytical structure of the retarded propagator of a uniform diffusion equation. Eq. (1.5) with a constant D ( r ) = D takes the form (cid:2) ∂ t − D ∂ r (cid:3) G (0) t ; r = δ ( t ) δ ( r ) , (A6) where we selected r = 0 due to translation invariance. The above equation in the Fourier domain G (0) ( t, r ) = (cid:82) d ω π d d q (2 π ) d e i qr − iωt G (0) ω, q becomes algebraic, as the differen- tial operators ∂ t → − iω and ∂ r → i q become diagonal: (cid:2) − iω + D q (cid:3) G (0) ω ; q = 1 , with the solution of a Lorentzian form G (0) ω ; q = 1 − iω + D q . (A7) This solution preserves causality, since its only singularity, at ω = − iD q , resides in the lower half-plane of ω . For t > , closing the integration contour in the lower half of the complex plane and using the residue theorem gives the Gaussian propagator in the qt representation, Eq. (1.8). The above consideration shows that causality is tightly re- lated to the integration in the time domain. We now inspect this relation closely, first without explicit reference to diffu- sion propagator, and then applying it to Eqs. (2.12) and (2.13). While the differentiation in the time domain of any function f ( t ) corresponds to the multiplication with − iω in the Fourier domain, ∂ t f ( t ) ⇔ − iωf ( ω ) , the inverse of the differential operator — i.e., the factor / ( − iω ) — corresponds to an in- definite integration (the antiderivative) in the time domain. However, infinitesimal shifts of the pole at ω = 0 result in different integration limits of definite time integrals. Con- sidering both possible shift directions, the same integration technique as above gives the Fourier transformation: (cid:90) d ω π e − iωt − i ( ω ± iε ) = ± θ ( ± t ) , (A8) where θ ( t ) is the unit step function, and ε → +0 . For an arbi- trary function f ( ω ) , which is integrable on the real axis of ω , the product f ( ω ) / ( − iω ± ε ) is Fourier-transformed according to the convolution theorem, f ( ω ) − iω + ε ⇔ (cid:90) d t (cid:48) f ( t (cid:48) ) θ ( t − t (cid:48) ) = (cid:90) t −∞ d t (cid:48) f ( t (cid:48) ) , (A9) f ( ω ) − iω − ε ⇔ − (cid:90) d t (cid:48) f ( t (cid:48) ) θ ( t (cid:48) − t ) = − (cid:90) ∞ t d t (cid:48) f ( t (cid:48) ) . (A10) Note that differentiating both Eqs. (A9) and (A10) with re- spect to t yields back f ( t ) . Eqs. (A9) and (A10) show that while the addition of ε is unimportant for the Fourier trans- form of derivatives, it is crucial for inverting differential op- erators. Shifting the pole of f ( ω ) / ( − iω ) downwards from the real axis results in the causal integration (A9), for which the resulted integral up to any time moment t depends on the integrand in the past, t (cid:48) < t . The opposite shift results in the dependence on the future, t (cid:48) > t . (The results differ by f | ω =0 = (cid:82) d t f ( t ) , and coincide if f | ω =0 = 0 , when f ( ω ) /ω is not singular.) Obviously, the first choice is adequate for the majority of solutions to equations describing the time evolu- tion of physical quantities such as Eq. (1.5), or Eqs. (2.12) and (2.13). The notation ε with ε → +0 is often abbreviated to simply +0 as it is done in Eqs. (2.12) and (2.13). Appendix B: OG with a finite number of pulses Consider the OG gradient wave form g ( t ) = g cos( ω t − φ ) with arbitrary initial phase φ and N oscillations, such that the total gradient duration T = N · π/ω . The corresponding g ( ω ) = g (cid:20) e − iφ · e i ( ω + ω ) T − i ( ω + ω ) + e iφ · e i ( ω − ω ) T − i ( ω − ω ) (cid:21) , e ± iω T = 1 , (B1) results in q ω = g ( ω ) / ( − iω ) , such that the wave form acts as the following “filter” for D ( ω ) in Eq. (2.9): q − ω q ω = g (1 − cos ωT )2 ω (cid:20) ω − ω ) + 1( ω + ω ) + 2 cos 2 φ ( ω − ω )( ω + ω ) (cid:21) . (B2) As q ω and Eq. (B2) are not singular when ω → and ω → ± ω , we do not need to specify how the zeroes of denominators in Eq. (B2) are shifted. Hence, one can directly substitute Eq. (B2) into Eq. (2.9) and integrate with any D ( ω ) along the real axis. However, to reveal the analytical structure, we find it useful to shift the frequency poles inside the square brackets by an infinitesimal positive imaginary part below the real axis, ω → ω + i in all the denominators, cf. Appendix A. As D ( ω ) is also analytic in the upper half-plane of the complex ω , Appendix A, the whole integrand in Eq. (2.9) remains analytic there. Hence, in the prefactor − cos ωT , the terms and − e iωT can be dropped for T > , as they identically vanish when closing the T → − T symmetry of Eq. (2.9) with q − ω q ω from Eq. (B2), which is forgone by this procedure. Hence, − cos ωT → − e − iωT , yielding the causal expression ln S = g (cid:90) d ω π D ( ω ) e − iωT ω (cid:20) ω + − ω ) + 1( ω + + ω ) + 2 − φ ( ω + − ω )( ω + + ω ) (cid:21) , ω + = ω + i , (B3) which starts to mimic the functional form of Eq. (2.13) — and that’s the goal! Dropping and − e iωT made the integrand in Eq. (B3) singular; however, the prescription how to go around its poles regularizes the result, which we obtain by a transformation into a sum of simple fractions. This yields the general relation between the intrinsic D ( ω ) and OG with N = ω T / π pulses: ln S = g ω (cid:90) d ω π D ( ω ) e − iωT (cid:20) (cid:18) ω + − ω ) + 1( ω + + ω ) (cid:19) + 2 sin φω − φ ω (cid:18) ω + − ω − ω + + ω (cid:19)(cid:21) . (B4) We can now see that the singularities in Eq. (B4) occur separately at ω = ± ω , and at ω = 0 for finite φ (when the gradient is not a pure cos wave form). Hence, we expect the response of D ( ω ) on these two frequencies. This response has the contributions of three distinct physical origins. The first two terms (in the braces) in Eq. (B4) yield the “pure OG” effect, that of the cos wave form with φ = 0 in the limit N → ∞ . The second term describes the ω = 0 singularity due to the finite time-average component in q ( t ) , ¯ q ≡ q ω → /T = ( g /ω ) sin φ present for φ (cid:54) = 0 , which is similar to the PG measurement. Not surprisingly, it yields the b PG D ( T ) ∝ D ( T ) sin φ contribution via the exact relation (2.13), i.e., the narrow-pulse D ( t ) with the diffusion time t = T weighted by the PG b -value contribution b PG = ¯ q T = ( g /ω ) T sin φ , which vanishes for the cos wave form (cf. Eq. (7) in ref. [130]). Finally, the last term, representing a finite OG-linewidth effect, is small as ∼ /N , as we will now show. While the remaining calculations can be done using the correspondence of e − i ( ω − ω ) T / ( ω − ω ) in the limit T → ∞ with the delta-function, we proceed via a more transparent transformation to the time domain. We use the relations (cid:90) d ω π e − i ( ω − ω ) T ω − ω + i D ( ω ) = − i (cid:90) T d t e iω t D ( t ) , (cid:90) d ω π e − i ( ω − ω ) T ( ω − ω + i D ( ω ) = (cid:90) T d t e iω t ( t − T ) D ( t ) valid for any retarded response functions D ( t ) and D ( ω ) related by the Fourier transformation (derived based on the property that a convolution in ω is a product in t , similar to how Eq. (A9) was derived). We can see that causality means that only t < T are integrated over, as expected. As a result, we transform Eq. (B4) into − ln S = b · (cid:34) (cid:82) T d t (cid:0) − tT (cid:1) D ( t ) cos ω t + 2 D ( T ) sin φ φ + 12 πN (cid:90) T d t D ( t ) sin ω t (cid:35) , b = πg Nω (1 + 2 sin φ ) , (B5) where the b -value is calculated [130, 302] assuming constant D ( ω ) ≡ D , such that D ( t ) = D δ ( t − (cid:15) ) | (cid:15) → +0 [24]. We can see that the total b = b | φ =0 + b PG is a sum of the pure OG and the pure PG contributions, since sin ω t is orthogonal to a constant. In the N (cid:29) limit of a clear separation of time scales π/ω (cid:28) T , we can extend the upper integration limit T → ∞ ; using t cos ω t = ∂ ω sin ω t in the first term of Eq. (B5), we can separate the main contribution and the ∼ /N correction: − ln S N (cid:29) (cid:39) b · (cid:34) Re D ( ω ) + 2 D ( T ) sin φ φ + 12 πN (cid:0) φ − ω ∂ ω (cid:1) Im D ( ω )1 + 2 sin φ (cid:35) , D ( ω ) = (cid:90) ∞ d t D ( t ) e iω t . (B6) Here, the term ∝ D ( T ) is exact due to Eq. (2.13), while setting ω = ω in D ( ω ) in the other two terms is precise up to ∼ ω /N . We can see that the first term in the square brackets of Eq. (B6) is a leading effect (it is not small when N → ∞ ); it gives the balance of the ω = ω and ω = 0 contributions with the “filter” weights defined by the phase φ in agreement with the general property of the second-order cumulant term, discussed after Eq. (2.9). The second term, ∼ Im D ( ω ) /N , is suppressed as ∼ /N . In the limit N → ∞ , the imaginary part of D ( ω ) does not contribute to the OG measurement [24]. There exists a case in which the /N term is comparable with the first one, namely, of a closed pore with a characteristic size ∼ a in the limit of low frequencies (long-times), ω a /D (cid:28) . In this limit, D ( T ) ∼ a /T happens to be parametrically as small as the second term. More precisely, the second term ∼ / πN of Eq. (B6) exactly cancels the D ( T ) sin φ contribution to the first term, and the net result is determined by Re D ( ω ) ∼ ω . This was first noticed by Sukstanskii [302] via the time- domain calculation for a one-dimensional impermeable box. Such behavior, however, is completely general. Indeed, for a pore of arbitrary shape, the dispersive diffusivity is a sum over the eigenmodes of the Laplace operator, D ( ω ) = D (cid:88) k C k − iωβ k D /a − iω ⇒ Re D ( ω ) | ω → (cid:39) (cid:88) k C k β k · ω a D (B7) with C k and β k determined by the pore geometry [14, 102]; to the order O ( ω ) , D ( ω ) (cid:39) − iωµa , D ( T ) = µa /T , where the dimensionless parameter µ = (cid:80) k ( C k /β k ) , and so both the D ( T ) and the Im D ( ω ) contributions to Eq. (B6) cancel each D ( ω ) ∼ ω a /D , the leading effect, vanishes very fast (quadratically) in the low- frequency limit [14, 20, 102, 137, 176, 303–305] (this scaling was first observed by Stepisnik et al. [305] in porous media). For the relevant case of a cylinder of radius a , C k = 2 / ( β k − and the sum (cid:80) k β k ( β k − = [18]. The estimate (B7) tells that OG (with any phase φ ) is less efficient than PGSE in creating adequate diffusion weighting if one wants to measure sizes of small fully confining compartments, cf. text after Eq. (2.38). Appendix C: Probing the
S/V limit with finite- N OG Let us now use the general finite- N relation (B4) for the S/V model (2.16); this setup is practically relevant to study cell density in brain tumors [143, 144]. The intrinsic D ( ω ) in this limit can be obtained performing Fourier transform of the outcome of Eq. (2.14) [130]: D ( ω ) (cid:39) D (cid:18) − S √ D V d · e iπ/ √ ω (cid:19) ⇔ D ( t ) (cid:39) D θ ( t ) (cid:18) δ ( t ) − S √ D V d · √ πt (cid:19) . (C1) Substituting Eq. (C1) into Eq. (B5), using (cid:90) T d t √ πt e iω t = (cid:114) ω (cid:104) C (2 √ N ) + i S (2 √ N ) (cid:105) , (cid:90) T d t √ πt tT cos ω t = − πN (cid:114) ω S (2 √ N ) , (C2) where the Fresnel integrals are defined in a standard way, S ( x ) = (cid:90) x d u sin πu , C ( x ) = (cid:90) x d u cos πu , (C3) we obtain − ln S = b · D (cid:104) − c ( φ, N ) · SV d √ (cid:113) D ω (cid:105) + 2 D ( T ) sin φ φ , (C4) where b is given in Eq. (B6), and the finite- N correction factor c ( φ, N ) = 2 C (2 √ N ) + 3 + 4 sin φ πN S (2 √ N ) . (C5) Here, the /ω term in Eq. (B4) yields the exact D ( T ) according to Eq. (2.13) as discussed in Appendix B. In the “ideal OG” limit N → ∞ , using C ( ∞ ) = , we obtain c ( φ, N ) → for any φ , such that Eq. (C4) yields Eq. (7) of ref. [130]. When, additionally, φ = 0 , we obtain − ln S = b · Re D ( ω ) , cf. Eq. (2.17) in the main text. We emphasize that performing the calculation in the frequency domain allows us to separate the time scales and identify the contributions to Eq. (C4) of two distinct physical origins. The D ( T ) term is completely general — i.e., the applicability of Eqs. (C4)–(C5) only requires the period of the oscillation to be short enough so that the model (C1) applies; the whole gradient train T can be long and (practically, always) falls out of the short-time S/V limit. Often times, at that point one can set D ( T ) ≈ D ∞ . As discussed in Appendix B , the D ( T ) term appears because at finite φ , the OG wave form can be thought of as a pure cos wave form and a PG with diffusion time t = T [130]. Hence, Eq. (C4) allows one to probe the S/V ratio by keeping the period short, yet the total gradient train as long as needed, to accumulate the diffusion weighting b ∝ N , Eq. (B6). If, additionally, the whole OG train T falls into the short-time limit, under a more stringent condition ( S/V ) √ D T (cid:28) , then D ( T ) is given by Eq. (2.16). In this limit, Eq. (C4) yields − ln S = bD (cid:34) − ˜ c ( φ, N ) · SV d √ (cid:114) D ω (cid:35) , ˜ c ( φ, N ) = c ( φ, N ) + √ N sin φ φ , (C6) where the re-defined correction factor ˜ c ( φ, N ) corresponds to Eq. (14) of ref. [302], where the problem of finite- N correction in the S/V limit was first considered. It was noted there, that ˜ c ( φ, N ) nominally diverges for N → ∞ as √ N . It is clear that this divergence occurs due to the √ T scaling from D ( T ) , and it eventually gets cut off when N becomes so large than D ( T ) falls out of the validity regime of Eq. (2.16). Hence, this spurious divergence is a result of defining ˜ c ( φ, N ) in ref. [302] by forcing Eq. (C6) to mimic the form of Eq. (2.16), instead of separating the physics at the two time scales, π/ω and T . The separation of scales identified in Eq. (C4) based on the general expression (B4) extends the validity of OG in the S/V limit far beyond the claim of ref. [302], “the high-frequency regime can be achieved only when the total diffusion time is smaller than the characteristic diffusion time” (implying ( S/V ) √ D T (cid:28) ), onto the practically relevant domain ( S/V ) (cid:112) D /ω (cid:28) , for any φ . Note that for pure cos gradient, ˜ c (0 , N ) = c (0 , N ) , due to the absence of the ω = 0 singularity in Eq. (B4), and Eq. (C6) agrees with Eqs. (C4) and (C5).agrees with Eqs. (C4) and (C5).