# ASSESSING CELLULAR MICROSTRUCTURE IN BIOLOGICAL TISSUES USING *IN VIVO* DIFFUSION-WEIGHTED MAGNETIC RESONANCE

EDITED BY : Julien Valette, Itamar Ronen and Sune Nørhøj Jespersen PUBLISHED IN : Frontiers in Physics and Frontiers in Neuroscience

#### Frontiers Copyright Statement

© Copyright 2007-2019 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.

The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.

Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.

Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.

As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.

All copyright, and all rights therein, are protected by national and international copyright laws.

The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use. ISSN 1664-8714 ISBN 978-2-88945-841-7 DOI 10.3389/978-2-88945-841-7

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

# Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# ASSESSING CELLULAR MICROSTRUCTURE IN BIOLOGICAL TISSUES USING *IN VIVO* DIFFUSION-WEIGHTED MAGNETIC RESONANCE

Topic Editors: Julien Valette, Commissariat à l'Energie Atomique et aux Energies Alternatives, France Itamar Ronen, Leiden University, Netherlands Sune Nørhøj Jespersen, Aarhus University, Denmark

Magnetic resonance imaging (MRI) and spectroscopy (MRS) techniques have opened new doors for examining biological tissues in vivo. By combining sensitization to diffusion using magnetic field gradients with a variety of imaging and localization schemes, diffusion-weighted MRI and diffusion-weighted MRS allow investigating translational diffusion of endogenous molecules, such as water or metabolites, in biological tissues, most commonly the brain but also other organs such as the prostate.

The typical voxel resolution of MRI or MRS is in the millimeter to centimeter range, much lower than the cellular scale. However, as molecules are typically diffusing over just a few µm during the duration of the measurement (the "diffusion time") and encounter numerous biological membranes at these scales, the average cellular microstructure has a critical influence on the measured diffusion signal. Hence, diffusion-weighted MRI and diffusion-weighted MRS are sensitive to tissue microstructure at a scale well below the nominal imaging resolution. However, the connection between diffusion properties and tissue microstructure remains indirect, so any attempt to quantify microstructure will rely on modeling.

The goal of this Research Topic was to gather experts in various acquisition and modeling strategies and show how these approaches, despite their own strengths and weaknesses, can yield unique information about cellular microstructure, and sometimes complement each other.

Citation: Valette, J., Ronen, I., Jespersen, S. N., eds. (2019). Assessing Cellular Microstructure in Biological Tissues using *In Vivo* Diffusion-Weighted Magnetic Resonance. Lausanne: Frontiers Media. doi: 10.3389/978-2-88945-841-7

# Table of Contents

Frederik B. Laun, Kerstin Demberg, Armin M. Nagel, Micheal Uder and Tristan A. Kuder

*25 Effective Potential for Magnetic Resonance Measurements of Restricted Diffusion*

Evren Özarslan, Cem Yolcu, Magnus Herberthson, Carl-Fredrik Westin and Hans Knutsson

*32 Improving the Realism of White Matter Numerical Phantoms: A Step Toward a Better Understanding of the Influence of Structural Disorders in Diffusion MRI*

Kévin Ginsburger, Fabrice Poupon, Justine Beaujoin, Delphine Estournet, Felix Matuschke, Jean-François Mangin, Markus Axer and Cyril Poupon

Gregory Lemberskiy, Els Fieremans, Jelle Veraart, Fang-Ming Deng, Andrew B. Rosenkrantz and Dmitry S. Novikov

Qiang Zeng, Feina Shi, Jianmin Zhang, Chenhan Ling, Fei Dong and Biao Jiang

*136 Assessing Diffusion in the Extra-Cellular Space of Brain Tissue by Dynamic MRI Mapping of Contrast Agent Concentrations*

Sébastien Mériaux, Allegra Conti and Benoît Larrat

*144 Time-Dependent Diffusion MRI in Cancer: Tissue Modeling and Applications*

Olivier Reynaud

*160 Optimizing Filter-Probe Diffusion Weighting in the Rat Spinal Cord for Human Translation*

Matthew D. Budde, Nathan P. Skinner, L. Tugan Muftuler, Brian D. Schmit and Shekar N. Kurpad

*173 On the Viability of Diffusion MRI-Based Microstructural Biomarkers in Ischemic Stroke*

Ilaria Boscolo Galazzo, Lorenza Brusini, Silvia Obertino, Mauro Zucchelli, Cristina Granziera and Gloria Menegaz

# Editorial: Assessing Cellular Microstructure in Biological Tissues Using in vivo Diffusion-Weighted Magnetic Resonance

Julien Valette1,2 \*, Itamar Ronen3 and Sune Nørhøj Jespersen4,5

1 MIRCen, Commissariat à l'Energie Atomique et aux Energies Alternatives, Fontenay-aux-Roses, France, 2 Neurodegenerative Diseases Laboratory, Centre National de la Recherche Scientifique, Université Paris-Sud, Université Paris-Saclay, UMR 9199, Fontenay-aux-Roses, France, 3 C. J. Gorter Center for High Field MRI, Department of Radiology, Leiden University Medical Center, Leiden, Netherlands, 4 Center of Functionally Integrative Neuroscience and MINDLab, Department of Clinical Medicine, Aarhus University, Aarhus, Denmark, 5 Department of Physics and Astronomy, Aarhus University, Aarhus, Denmark

Keywords: medical imaging, cell morphology, diffusion-weighted (DW), MRI, NMR, diffusion, microstructure

#### **Editorial on the Research Topic**

### **Assessing Cellular Microstructure in Biological Tissues Using in vivo Diffusion-Weighted Magnetic Resonance**

Since the seminal works of Stejskal and Tanner in the 1960s [1], the study of molecular self-diffusion by diffusion-weighted (DW) nuclear magnetic resonance (NMR) has emerged as a dynamic and rich research field. Initially introduced as a tool to measure diffusion in materials and to estimate microstructural properties of materials by modeling deviations from free diffusion caused by microstructural barriers [2], DW-NMR was soon proposed as a tool to characterize cellular microstructure in biological tissues [3]. It was then quite a natural—although challenging—step to endow magnetic resonance imaging (MRI) with diffusion-weighting capabilities to map water diffusion in vivo [4], which almost immediately revolutionized the diagnosis of brain ischemic stroke [5–7]. This partial historical account of the early days highlights the different "flavors" still characterizing our research field at the interface of physics, chemistry, biology, computer science, and medical imaging: a unique blend of powerful concepts and theories, cutting-edge methodological, and technical innovations, ideally leading to highly-relevant biological insights or clinical applications.

Various approaches, each with their own strengths and weaknesses, have been developed over the years to probe microstructure with ever higher degrees of sensitivity and specificity. The goal of this Research Topic was to gather experts in various acquisition and modeling strategies and show how these approaches can yield unique and complementary information about cellular microstructure. In parallel to the topics discussed here, other developments are also driving in vivo DW-NMR forward. Important examples are the ever increasing field strength and gradient performance of MRI scanners, boosting signal-to-noise, and enabling larger and better temporally resolved diffusion-weighting, respectively [8, 9]. However, these technical innovations also come with new difficulties or limitations, e.g., increased tissue heating, decreased transverse relaxation times, increased sensitivity to motion and field fluctuations, decreased field homogeneity and increased geometrical distortions at high fields [10, 11], and increased peripheral nerve stimulation by powerful gradients [9]. The drawbacks resulting from these technical innovations are partially offset by parallel developments in coil and pulse sequence technology. Such technical innovations and solutions are not discussed in this Research Topic.

Edited and reviewed by: Ewald Moser, Medical University of Vienna, Austria

> \*Correspondence: Julien Valette<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 12 February 2019 Accepted: 27 February 2019 Published: 21 March 2019

#### Citation:

Valette J, Ronen I and Jespersen SN (2019) Editorial: Assessing Cellular Microstructure in Biological Tissues Using in vivo Diffusion-Weighted Magnetic Resonance. Front. Phys. 7:37. doi: 10.3389/fphy.2019.00037

Many fundamental questions remain open about the relationships existing between the diffusion-weighted signal and its derived diffusion metrics on one hand, and the underlying microstructural features on the other hand. Approximately half of the papers of this research topic aim at addressing such fundamental questions, using different approaches.

Some rely primarily on theoretical considerations. Özarslan et al. revisit the q −1 decay in narrow fibers, taking curvature into account. The problem is examined in three distinct temporal regimes of the Stejskal-Tanner pulse-field gradient experiment. For smaller cells, the q −1 decay of the orientationally-averaged signal is only predicted for straight fibers. This decay is more general for cells with longer projections, but fades away for curvy structures as gradient pulse duration increases. The authors stress that the q −1 decay could represent an intermediate range, as the true asymptotic behavior is governed by a steeper attenuation. Laun et al. tackle another kind of diffusion experiments, where the time-dependence of the apparent diffusion coefficient (ADC) is measured. In the short diffusion time (td) regime, the ADC is known to decrease as the square root of td multiplied by the surface to volume ratio (S/V). The next term in the expansion is proportional to td multiplied by a factor depending on permeability and curvature. Laun et al. show that this term is identically zero for oscillating gradients, meaning that permeability and curvature cannot be easily measured with these sequences, while S/V term may be easier to quantify than with pulsed-field gradients. Ozarslan et al. present a new perspective for the representation of restricted diffusion that is highly suitable for diffusion experiments performed with commonly available gradient hardware, for example on clinical scanners. Diffusion experiments under such conditions are typically performed with long gradient pulses, thus blurring the information necessary for the characterization of the restricting pore. The approach presented here assumes that diffusion occurs in an effective quadratic potential landscape, thus resulting in a linear (Hookean) dependence of the stochastic effective force on the center-of-mass position. The implications that this approach may have on the interpretation of diffusion experiments performed with long diffusion gradients are discussed, together with possible extensions of this theoretical framework to potential landscapes other than quadratic.

Theoretical analysis typically relies on idealized assumptions for mathematical tractability, and therefore struggle to evaluate the quantitative impact of complex structural features in realistic acquisition conditions. In that case, numerical simulations may help. Ginsburger et al. describe simulations in synthetic white matter mimicking axonal bundles (Ginsburger et al.), with a small set of parameters controlling the level of tissue complexity and disorder. Their results show that angular dispersion and most importantly axonal beading have non-negligible effects on extracellular diffusion which may be measured using DW-MRI clinical protocols with oscillating gradients.

A major limitation of both theoretical analysis and numerical simulation may be the insufficient knowledge about the precise nature or magnitude of some features potentially affecting diffusion (such as the properties and detailed distribution of membranes or the existence of molecular pools with different relaxation times), thus precluding solid theoretical or numerical approaches. In such situations, we must rely on experiments for gaining further fundamental insights about what diffusion measurements reflect. Shemesh investigates the microstructural correlates of fractional anisotropy (FA) and microscopic FA (µFA) using high resolution double-diffusion encoding (DDE) and double oscillating gradient diffusion encoding (DODE) experiments in fixed spinal cords. The results indicate that µFA(DODE) mostly depend on axon sizes and myelin water fraction, whereas µFA(DDE), FA(DODE), and FA(DDE) correlate to a much lesser extent with those cellular features. This is promising for future interpretations of these metrics in health and disease. Valette et al. recapitulate how intracellular metabolite diffusion has been measured in the brain, using DW-NMR spectroscopy in vivo, from td<1 ms to td>1 sec, and how the picture emerges that metabolites are not confined into subcellular compartments, but are instead diffusing in long fibers with some long-range complexity. Some possible future directions to extend the range of diffusion times, and what information about cell structure can be gleaned, are also discussed. As a final illustration of how well-thought experiments may yield fundamental insight, Lemberskiy et al. exploit data obtained in the prostate at multiple td and echo times to select and combine two different models, which allows disentangling the cellular compartment (exhibiting shorter relaxation time T2 and associated with diffusion through extended disorder) from the lumen (exhibiting longer T2 and associated with the short-time S/V limit). Beyond fundamental understanding, data modeling allows extracting lumen diameter, fiber diameter and membrane permeability, with possible application to noninvasive diagnosis and staging of prostate cancer with MRI.

The potential relevance for applications of DW-NMR in vivo certainly underlies all works in our field, even those aiming at gaining fundamental knowledge, but some works more directly aim toward applications in animals or humans. Such more "application-driven" papers, either describing optimized methods to measure biologically/clinically relevant parameters, or evaluating the value of these methods under various conditions (e.g., in disease), represent approximately another half of this research topic.

Hansen and Jespersen review fast kurtosis imaging protocols, which are capable of estimating the mean, axial and radial kurtosis (i.e., deviations from Gaussian diffusion when q is varied at constant td) with only a few diffusion directions. The review describes the theoretical foundations for these methods, technical and practical considerations, and potential applications mostly related to peripheral nerve or white matter microstructure. In a similar vein, Jelescu and Budde review popular models of diffusion in white matter to extract parameters such as axonal orientation dispersion, axon diameter, and extracellular volume fraction, from the q-dependence at constant td. Assumptions and required data are detailed for each model, before examining how these models can be validated using numerical simulations, phantoms, or ex vivo microscopy. Zeng et al. consider the fitting quality of multi-exponential signal expressions to diffusion data from white matter in healthy volunteers. Specifically, they propose to describe water diffusion in white matter in terms of a modified tri-exponential expression; a biexponential plus a "dot" compartment, i.e., a water pool having zero diffusivity. The latter component is interpreted as restricted diffusion, and the modified tri-exponential is found to have lower corrected Akaike criterion than both free tri-exponential and biexponential signal representations in several white matter ROIs.

Mériaux et al. present an original method which, unlike other papers of this research topic, does not assess diffusion by gradient-induced diffusion-weighting, but by dynamically acquiring concentration maps of contrast agents as they diffuse through the extracellular space (as more traditionally performed with optical imaging) to evaluate the ADC and tortuosity in the extracellular space. One original point of their approach is that the agent is delivered to the brain after i.v. injection and ultrasound-induced blood-brain barrier opening, avoiding invasive intracerebral injection. This technique may be useful to characterize the alterations of extracellular space diffusion properties in pathologies such as Alzheimer's disease, and also to predict drug biodistribution.

Some methods may prove particularly relevant for specific pathologies, such as tumors. Reynaud starts by reviewing the fundamentals of time-dependent diffusion in the short-time regime and near the tortuosity limit. The main three geometrical models (IMPULSED, POMACE, VERDICT) implemented so far to describe tumor cells as spheres are then reviewed, and their applications for estimating tumor intra- and extracellular volume fraction and cell size is then presented and discussed. Budde et al. present an optimization of pulse sequences for imaging of spinal cord injury. The authors consider orientationally invariant DDE, filter probe DDE and filter probe single-diffusion encoding (SDE), in which the idea of a filter probe is basically to suppress extra-axonal water with a diffusion gradient perpendicular to

# REFERENCES

the cord. They find that filter probe SDE coupled with reduced field-of-view excitation is the best protocol with highest signalto-noise ratio and with high sensitivity to injury, without contamination from edema and cerebrospinal fluid. This is demonstrated with simulations, and in vivo rat and human spinal cord imaging. Finally, one important clinical use of MRI is in assessing tissue damage and prospective recovery in ischemic stroke. Question remains as to what are the most meaningful diffusion metrics in ischemic stroke, that is, those that most sensitively reflect microstructural characteristics of tissue affected by ischemia. Galazzo et al. consider a set of indices obtained from the Simple Harmonic Oscillator Based Reconstruction and Estimation (SHORE). Assessment of quantities associated with the SHORE formalism are investigated in a tract-based fashion in white matter and region-of-interest analysis in gray matter in patients and controls, providing insight into the potential benefits of adopting the SHORE model and some of its metrics in the clinic.

In conclusion, this research topic illustrates the multifaceted nature of cellular microstructure assessment using DW-NMR. While no particular method can grasp the full complexity of cellular microstructure, each of them can shed light on different aspects and find specific applications. In the long term, some approaches will certainly prove more insightful or useful than others. For now, the present collection of works pays tribute to the vitality and diversity of our research field.

# AUTHOR CONTRIBUTIONS

JV drafted the first version of the editorial. IR and SJ revised the first draft and made contributions about papers they edited.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Valette, Ronen and Jespersen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Influence of the Size and Curvedness of Neural Projections on the Orientationally Averaged Diffusion MR Signal

Evren Özarslan1,2 \*, Cem Yolcu1 , Magnus Herberthson3 , Hans Knutsson1 and Carl-Fredrik Westin1,4

1 Department of Biomedical Engineering, Linköping University, Linköping, Sweden, 2 Center for Medical Image Science and Visualization, Linköping University, Linköping, Sweden, 3 Division of Mathematics and Applied Mathematics, Department of Mathematics, Linköping University, Linköping, Sweden, 4 Laboratory for Mathematics in Imaging, Department of Radiology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA, United States

#### Edited by:

Julien Valette, Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA), France

#### Reviewed by:

Silvia Capuani, Consiglio Nazionale Delle Ricerche (CNR), Italy Gernot Reishofer, Medical University of Graz, Austria Sune Nørhøj Jespersen, Aarhus University, Denmark

> \*Correspondence: Evren Özarslan<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 29 September 2017 Accepted: 09 February 2018 Published: 02 March 2018

#### Citation:

Özarslan E, Yolcu C, Herberthson M, Knutsson H and Westin C-F (2018) Influence of the Size and Curvedness of Neural Projections on the Orientationally Averaged Diffusion MR Signal. Front. Phys. 6:17. doi: 10.3389/fphy.2018.00017 Neuronal and glial projections can be envisioned to be tubes of infinitesimal diameter as far as diffusion magnetic resonance (MR) measurements via clinical scanners are concerned. Recent experimental studies indicate that the decay of the orientationally-averaged signal in white-matter may be characterized by the power-law, E¯(q) ∝ q −1 , where q is the wavenumber determined by the parameters of the pulsed field gradient measurements. One particular study by McKinnon et al. [1] reports a distinctively faster decay in gray-matter. Here, we assess the role of the size and curvature of the neurites and glial arborizations in these experimental findings. To this end, we studied the signal decay for diffusion along general curves at all three temporal regimes of the traditional pulsed field gradient measurements. We show that for curvy projections, employment of longer pulse durations leads to a disappearance of the q −1 decay, while such decay is robust when narrow gradient pulses are used. Thus, in clinical acquisitions, the lack of such a decay for a fibrous specimen can be seen as indicative of fibers that are curved. We note that the above discussion is valid for an intermediate range of q-values as the true asymptotic behavior of the signal decay is E¯(q) ∝ q −4 for narrow pulses (through Debye-Porod law) or steeper for longer pulses. This study is expected to provide insights for interpreting the diffusion-weighted images of the central nervous system and aid in the design of acquisition strategies.

Keywords: diffusion, magnetic resonance, anisotropy, Stejskal-Tanner, curvature, curvilinear, power-law, powder

# 1. INTRODUCTION

Diffusion-sensitized magnetic resonance acquisitions have been employed to recover the microscopic building blocks of complex nervous tissue. Simplified models exploiting the compartmentalized structure of the tissue are instrumental in this endeavor. Water molecules in the intra- and extra-cellular spaces have been envisioned to form separate compartments with different signal characteristics [2]. The intracellular signal is also thought to represent the superposition of contributions from cells of different types, shapes, and orientations [3]. The same argument has been employed even for a single neuron wherein each neurite has been considered to comprise

**8**

a collection of straight compartments [4]. Such representation of neurites as slender cylinders distributed in random orientations within the voxel is perhaps the model most relevant to the current study.

The diameter of neurites, and in fact all neural projections, is so small that diffusion in the transverse plane may be negligible. More explicitly, a cylinder with the same diameter as a neurite would not suffer any signal loss in a typical clinical diffusion MRI measurement when the diffusion gradients are applied in the direction perpendicular to the cylinder's axis. This justifies assigning zero value to transverse diffusivity for molecules confined in the cylinder. Such behavior [5] was indeed observed for N-acetyl-L-aspartate (NAA) diffusion in the brain [6] and have been employed for water in recent models [7].

In this work, we consider the pulsed field gradient measurement introduced by Stejskal and Tanner [8] featuring diffusion encoding gradients **G** of duration δ, whose leading edges are separated from each other by duration 1 (see **Figure 1A** for the effective gradient waveform). We define **q** = γ δ**G**, where γ is the gyromagnetic ratio and note that for sufficiently small values of q = |**q**|, the signal for each compartment can be approximated with a Gaussian, i.e.,

$$E(q) \approx e^{-q \dagger V q} \,. \tag{1}$$

Considering the form (Equation 1) of the signal, **V** can be referred to as the signal decay tensor1 . The geometric parameters of the compartment have typically a complicated relation to **V**; the exact form of such relation is dictated by the temporal parameters (δ and 1) of the diffusion encoding pulse sequence. For axially symmetric **V**, we shall denote by vk and v⊥ the eigenvalues of **V** associated with directions parallel and perpendicular to the symmetry axis, respectively.

Our focus in this work is the orientationally-averaged signal, which can be obtained by computing the "isotropic component" of the signal, as was referred to in Özarslan and Basser [10] and actually estimated as a byproduct of the q-space signal representation in Özarslan et al. [11]. Alternatively, the signal values measured over all gradient directions at a particular q-value can be averaged [12, 13] so that any dependence on the direction of the gradient vector is lost. Repeating this procedure for all q-values reduces the data collected over the three-dimensional q-space into a one-dimensional profile, which does not contain any information on ensemble (macroscopic) anisotropy2 . The estimated signal profile represents the decay for the so-called "powdered" specimen, which contains an isotropic distribution of each and every compartment in the original specimen [19, 20]. The orientationally-averaged signal for axisymmetric compartments, each of which contributes

FIGURE 1 | (A) Stejskal and Tanner's pulsed field gradient experiment features a pair of rectangular gradient pulses of duration δ whose leading edges are separated by 1. (B) The timing parameters of this experiment lie on or above the dashed line since the separation of the two pulses (1) has to be at least as long as the duration of each of them (δ). Thus, the experiment has three distinct regimes (labeled A, B, and C) based on whether these timing parameters are short or long. These are the three regimes that exhibit all interesting features of the signal. The same features are expected to be observed to various degrees in the intermediate region between these three regimes.

according to Equation (1), is thus the same as the signal for an isotropic ensemble of such compartments, and is given by [6, 21, 22]

$$\bar{E}(q) = \frac{\sqrt{\pi}e^{-q^2\nu\_{\perp}} \operatorname{erf}(q\sqrt{\nu\_{\parallel} - \nu\_{\perp}})}{2q\sqrt{\nu\_{\parallel} - \nu\_{\perp}}} \,. \tag{2}$$

This expression predicts a squared exponential decay in general. However, if signal loss is limited to only the fiber direction (v⊥ = 0), a much slower decay emerges. A power law of the form E¯ ∝ q −1 , to be specific. Therefore, the appearance of this particular power-law can be used as an indicator of vanishing transverse diffusivity, in other words, the signal decay tensor **V** being of rank 1. We note that the problem of characterizing the orientationally averaged signal is considerably more complicated when Equation (1) cannot be used to represent the compartmental signal; addressing this issue is one of our goals in this study.

The E¯ ∝ q −1 decay alluded to above was recently reported in white matter [1, 23]. These studies have found that in

1We note that the decay tensor is closely related to an apparent diffusion tensor (ADT) whose time-dependence has been shown to be sufficient for describing (approximately-)Gaussian diffusion via general gradient waveforms [9]. In this study, it proves convenient to employ the **V**-tensor, which encapsulates all dependencies other than that on the **q**-vector.

2Yet another approach would involve employing alternative gradient waveforms for isotropic diffusion weighting [14–18]. However, we do not discuss such sequences here because of the complicated dependence of the signal intensity on the gradient waveform.

white-matter-dominated regions of the brain, the orientationally averaged signal exhibits a decay ∼ q −c with an exponent c close to 1, in support of cylindrical neural projections as remarked above. In gray-matter-dominated regions, McKinnon et al. [1] have observed a larger exponent c ≈ 1.8 ± 0.2. They proposed that the apparent breakdown of the cylinder model may indicate a significantly larger permeability of the cellular membranes in gray matter vs. white matter. Here, we investigate an alternative hypothesis. Namely, the departure of the exponent c from 1 can well be due to impermeable but curved projections. To assess this point, we studied the influence of neural projections' size and shape on the orientationally-averaged diffusion MR signal. Due to the large variability in the geometric features of the neural cells, all temporal regimes of the Stejskal-Tanner sequence were considered, and the problem was studied both in the small-q regime as well as at larger q-values for which Equation (1) and thus Equation (2) are inaccurate.

Investigation of power-like tails in the diffusion MR signal goes back to Köpf et al. [24], where large values of the wave vector were achieved using a fringe field method. A range of exponents (roughly between −1.8 and −4.6) were observed across various nonneural tissue types, as well as stretched exponential behavior, which was ascribed to fractional Brownian motion. In a subsequent study, Yablonskiy et al. [25] predicted an exponent of −2 for specimens featuring compartments with a distribution of diffusivities. Jian et al. [26] considered a parametric tensor distribution, which suggested a signal decay with a general power-law tail. These studies, however, observe or predict the signal decay for measurements along a single direction. As for the orientationally (powder) averaged signal, a quite general statement regarding an asymptotic power-law decay in the diffusion MR signal is the Debye-Porod law [27]. Here, the orientationally averaged signal measured using narrow pulses is predicted to follow a q −4 tail under quite general considerations. In our discussion, we take up apparent violations of this.

Although the influence of fiber curvature on the MR signal has been considered [28–33] in various contexts, to our knowledge, this is the first study to provide explicit expressions for the signal decay for diffusion along a general parametrized curve and study its effect on the orientationally averaged signal.

In the next section, we provide explicit expressions for the MR signal for diffusion on curves in three distinct temporal regimes of the Stejskal-Tanner measurement. In the subsequent section, we discuss the implications of our theoretical findings as they relate to the morphology of neural cells and recent experimental observations. The article is concluded following a brief discussion of what observed power-law tails in the powder averaged signal may represent in that context, as well as from the perspective of the Debye-Porod law.

# 2. COMPARTMENTAL AND ORIENTATIONALLY-AVERAGED SIGNAL FOR DIFFUSION ALONG CURVES

The effective gradient waveform of a traditional Stejskal-Tanner measurement is shown in **Figure 1A**. In this work, we consider

(ignoring branchings and other features such as beading patterns [34]) are: its contour length denoted by ℓ, its radius of gyration Rg, the characteristic radius of curvature Rc, and the radius of the projection's cross-section R0.

three distinct regimes of this pulse sequence based on whether δ and 1 are short or long. These regimes are indicated by the letters A, B, and C on the δ-1 plane in **Figure 1B**. We note that the essential features of the signal at these three extreme situations are exhibited to some extent for more general timing values, i.e., within the interior of the triangle whose vertices are at A, B, and C.

The relevant size parameters of a simplified neural projection (ignoring branchings and other features such as beading patterns [34]) are: its contour length denoted by ℓ, its radius of gyration Rg, its characteristic curvature radius Rc, and the radius of the projection's cross-section R0. These parameters are illustrated for a representative projection in **Figure 2**. As mentioned in Introduction, R0 is typically so small that qR0 ≪ 1 for clinical MRI; this justifies representing the neurites and glial projections via one-dimensional curves. Thus, we shall consider diffusion taking place on a curve **r**(s) = (r1(s),r2(s),r3(s))⊺ parameterized by its arclength s, where 0 6 s 6 ℓ.

# 2.1. Regime A: Short Diffusion-Time

In the first case, diffusion is observed for such a short time that the hindrances have not been encountered. This condition implies that the pulse durations are short as well. In fact, when D1 ≪ R 2 c , spins spread so little that even the most curved point along the projection seems like a straight segment3 . Thus, the

3The "most curved point" should be taken with a grain of salt. Even though, mathematically, the approximations in this regime would require a diffusion length √ D1 much smaller than the smallest radius of curvature along the curve, Rc is not

compartmental signal occurs as the average, over the curve, of the signals originating from the tangents of the curve:

$$E(\mathbf{q}) = \frac{1}{\ell} \int\_0^\ell \mathrm{d}s \, e^{-D\Delta \mathbf{q} \cdot \mathbf{\hat{r}} \hat{\mathbf{r}}(s) \hat{\mathbf{r}}^\mathsf{T}(s) \mathbf{q}} \,, \tag{3}$$

where

$$
\hat{t}(\mathbf{s}) = \frac{\mathbf{d}r(\mathbf{s})}{\mathbf{d}s} \tag{4}
$$

is the unit tangent to the curve.

As shown in Appendix A (see the Supplementary Material), the very small-q behavior of the orientationally-averaged signal (q 2D1 ≪ 1) is given by E¯(q) ≈ e −q 2D1/3 . Thus, the decay rate is determined solely by the (intracellular) diffusivity, and bears no geometric feature of the neural arborization in this regime.

For a more general analysis involving larger q-values, one can still employ Equations (1, 2) as follows. We are interested in the orientational average of the compartmental signal in Equation (3), i.e.,

$$\bar{E}(q) = \frac{1}{\ell} \int\_0^\ell \mathrm{d}s \left\langle e^{-D\Delta q \mathsf{T} \hat{I}(s) \hat{t}^\mathsf{T}(s) q} \right\rangle , \tag{5}$$

where the order of integration and orientational averaging, indicated by the angular brackets, was changed. However, the expression within the angular brackets is of the Gaussian form (Equation 1) with a rank-1 decay tensor whose non-vanishing eigenvalue is D1. Consequently, the signal for the powdered specimen is given through Equation (2) by setting v⊥ = 0, and vk = D1 to be

$$\bar{E}(q) = \frac{\sqrt{\pi}\,\mathrm{erf}(q\sqrt{D\Delta})}{2q\sqrt{D\Delta}}\,. \tag{6}$$

Clearly, for larger q-values, the orientationally-averaged signal decay in regime A is proportional to q −1 irrespective of the shape of the curve as long as R 2 c ≫ D1. It can be observed that Equation (5) has the form of a signal arising from a uniform orientational distribution of "sticks". Hence the emergence of q −1 at large qvalues can be justified alternatively by Veraart et al.'s arguments [23].

#### 2.1.1. Incorporating Curvature Effects

The above expression holds when the diffusion distance is much smaller than Rc as pointed out earlier. Here, we would like to generalize this expression to larger timing parameters, 1 and δ, to allow for the possibility that the diffusion distance during the course of the experiment is long enough for the molecules to traverse an approximately circular arc along the curve. Moreover, we assume that there is a single characteristic radius of curvature that represents the effective curvedness of the entire projection. This characteristic curvature is denoted by Rc .

Let Earc(**n**ˆ, ϕ, **q**) denote the signal for a single such arc, where **n**ˆ is the unit vector normal to its plane and ϕ is the polar coordinate of the center of the arc in a cylindrical reference frame oriented along **n**ˆ. The orientationally averaged signal can then be written as the average of Earc(**n**ˆ, ϕ, **q**), over all possible realizations of a single arc, i.e.,

$$\bar{E}(q) = \frac{1}{4\pi} \int\_{S\_2} \mathrm{d}\hat{\mathfrak{n}} \, \frac{1}{2\pi} \int\_0^{2\pi} \mathrm{d}\varphi \, E\_{\mathrm{arc}}(\hat{\mathfrak{n}}, \varphi, q) \,, \tag{7}$$

where S2 denotes the unit sphere. The second average simply defines the signal for a full circle of radius Rc . If we denote by Ecirc(**n**ˆ, **q**) the signal for such a circle whose plane has the normal vector **n**ˆ, the orientationally averaged signal can be expressed as

$$\bar{E}(q) = \frac{1}{4\pi} \int\_{S\_2} \mathrm{d}\hat{n} \, E\_{\mathrm{circ}}(\hat{n}, q) \,. \tag{8}$$

Due to its axial symmetry, the signal for the circle has the functional dependence Ecirc(**n**ˆ · ˆ**q**, q), where **q**ˆ = **q**/q. Since **n**ˆ · ˆ**q** is invariant under exchange of the unit vectors **n**ˆ and **q**ˆ, one is free to replace the integration variable above with **q**ˆ and fix **n**ˆ instead [19] . With the variable θ defined through **n**ˆ· ˆ**q** = cos θ, we moreover note that since there is no motion and hence no signal attenuation in the direction along **n**ˆ, the integrand's functional dependence may further be reduced to Ecirc(q sin θ), namely, the signal obtained when a **q**-vector of magnitude q sin θ is applied in the plane of the circle. Upon taking these observations into account, we obtain

$$\bar{E}(q) = \int\_0^{\pi/2} \mathbf{d}\theta \,\,\sin\theta \,\, E\_{\text{circ}}(q\,\,\sin\theta)\,\,. \tag{9}$$

#### **<IP_ADDRESS>. Narrow Pulses**

First, we shall consider the scenario involving narrow pulses, i.e., Dδ ≪ R 2 c , but allow for the possibility that D1 ≈ R 2 c . For arbitrary 1, the signal for a full circle is given by Özarslan et al. [28]

$$E\_{\rm circ}(q) = J\_0(q\mathcal{R}\_\epsilon)^2 + 2\sum\_{n=1}^\infty e^{-n^2D\Delta/\mathcal{R}\_\epsilon^2}J\_n(q\mathcal{R}\_\epsilon)^2,\tag{10}$$

where Jn denotes the nth order Bessel function. This expression yields, via Equation (9), the orientationally-averaged signal to be

$$\begin{split} \bar{E}(q) &= f\_0(2q\mathcal{R}\_c) + \frac{\pi}{2} \left[ f\_1(2q\mathcal{R}\_c)\mathbf{H}\_0(2q\mathcal{R}\_c) - f\_0(2q\mathcal{R}\_c)\mathbf{H}\_1(2q\mathcal{R}\_c) \right] \\ &+ 2\sum\_{n=1}^{\infty} e^{-n^2D\Delta/R\_c^2} \frac{(q\mathcal{R}\_c)^{2n}}{(2n+1)!} \end{split}$$

$${}\_1F\_2\left(n+\frac{1}{2};n+\frac{3}{2},2n+1;-q^2R\_c^2\right),\tag{11}$$

where **H**α is the Struve function of order α and 1F2 represents the generalized hypergeometric function. We note that this orientationally averaged signal in the narrow pulse regime has the asymptotic behavior

$$\bar{E}(q) \sim \frac{\sqrt{\pi} \,\mathrm{erf}(q\sqrt{D\Delta})}{2q\sqrt{D\Delta}} + \frac{(D\Delta)^2 q^2 e^{-D\Delta q^2}}{6R\_c^2} - \mathcal{O}(\mathcal{R}\_c^{-4}) . \tag{12}$$

such a strict measure. Rather, it is the minimal radius of curvature that the curve exhibits along a portion of it significant enough to influence the signal.

This expression contains the curvature-related correction to Equation (6), which is valid for straight fibers, and suggests that curvature induced effects are rather limited for the case of narrow pulses, as the E¯(q) ∝ q −1 dependence is unaffected at larger q-values.

#### **<IP_ADDRESS>. Longer pulse durations**

To investigate the influence of longer pulse durations, we numerically evaluated the integral in Equation (9) using Simpson's rule [35]. We adapted the multiple correlation (MCF) framework to the problem of diffusion on a circle. The details of this procedure are provided in Appendix B (see the Supplementary Material).

In **Figure 3**, we illustrate the signal decay curves at different pulse durations. In these simulations, the D1 value was set equal to R 2 c ; during a time interval of duration 1, the spread of molecules on the circle is about 80◦ . When the pulse duration is short, we verified that the decay obeys the expressions in Equations (11) and (12) (results not shown). The power-law, E(q) ∝ q −1 , which is valid in the narrow pulse regime does not prevail when the pulses are prolonged.

# 2.2. Regime B: Long Diffusion-Time and Narrow Pulses

When diffusion is probed via a pair of impulses (or, when Dδ ≪ R 2 c ) separated from each other by a duration long enough for the molecules to reach all points within the curve, i.e., D1 ≫ ℓ 2 , the signal is given simply by the expression E(**q**) = | ˜ρ(**q**)| 2 , where ρ˜(**q**) is the Fourier transform of the density

$$\rho(r) = \frac{1}{\ell} \int\_0^\ell \mathrm{d}s \,\delta(r - r(s))\,,\tag{13}$$

which is uniform along the curve. The orientationally averaged signal decay is thus obtained by averaging | ˜ρ(**q**)| 2 , and is given by

$$\bar{E}(q) = \frac{1}{\ell^2} \int\_0^\ell \mathrm{d}s \int\_0^\ell \mathrm{d}s' \frac{\sin\left[q \left|r(s) - r(s')\right|\right]}{q \left|r(s) - r(s')\right|} \,. \tag{14}$$

This expression is referred to as the Debye scattering equation [36], which is widely utilized in studies employing small angle scattering experiments for characterizing the structure of polymers. The essential features of the resulting E¯(q) curve are well-understood [37, 38]. As shown in Appendix A (see the Supplementary Material), the small-q regime (qRg ≪ 1, the "Guinier regime" of scattering experiments) is described by the relationship E¯(q) ≈ e −(qRg) 2 /3 . As for larger q-values, such studies indicate that depending on the structure of the polymer, different sections of the curve could be characterized by different power-laws. For example, Gaussian chains undergo q −2 decay, while fractional exponents are obtained for curves exhibiting fractality. Perhaps the most relevant finding, however, is that wormlike structures (i.e., those characterized by a so-called persistence length over which the polymer is likely to retain its direction) are characterized by a decay ∝ q −1 at q-values about the reciprocal of the chain's persistence length [39]. Thus, even a class of non-straight structures exhibit q −1 decay in Regime B.

FIGURE 3 | Signal attenuation curves for different values of the pulse duration. The separation of pulses was taken so that D1 = R 2 c . The q −1 decay, visible at short pulse durations, disappears as the pulses are prolonged.

## 2.3. Regime C: Long Pulse-Duration

For the traditional Stejskal-Tanner measurement utilizing a pair of identical pulses in opposite directions, the compartmental signal has the form

$$E = \left\langle e^{-iq \cdot (\mathfrak{k}\_2 - \mathfrak{k}\_1)} \right\rangle\_{\text{paths}},\tag{15}$$

where the averaging is performed over all trajectories, with

$$\mathfrak{k}\_n = \frac{1}{\delta} \int\_{t\_n}^{t\_n + \delta} r(t) \, dt \tag{16}$$

being the center of mass coordinate [40] of the fragment of trajectory coinciding with each pulse (t1 = 0 and t2 = 1). Therefore, the MR signal (Equation 15) elicited by flat gradient pulses is not sensitive to the Brownian trajectories instant by instant but only in a time averaged sense.

In the limit Dδ/ℓ2 → ∞, the two random variables lose correlation and the explicit dependence on 1 disappears, leading to the signal intensity [40]

$$E = \left| \tilde{p}\_{\rm cm}(\mathbf{q}, \delta) \right|^2,\tag{17}$$

where

$$\tilde{p}\_{\rm cm}(q,\delta) = \int p\_{\rm cm}(\xi,\delta) e^{-iq\cdot\xi} \,\mathrm{d}\xi \tag{18}$$

is the Fourier transform of the center of mass distribution. Moreover, in this limit, the distribution for the random variable ξ approaches a Gaussian due to the central limit theorem [41] 4 . Consequently, the signal (Equation 17) also approaches a quadratic exponential form, encoding no more structural information than the variance of pcm(ξ , δ). Whether the domain

4This non-obvious statement, whose rigorous proof can be found in mathematics literature [42], was instrumental in our identification of an effective potential for restricted diffusion [43].

is an irregularly curved fiber or a much more regular shape like a sphere, its fine structural features will find no representation in the signal acquired this way, since the length of the time averaging (pulse) has suppressed all short-scale (high q) variation encoded in the cumulants higher than the second.

Hence, the compartmental signal has the form Equation (1). Here, though, **V** is the variance tensor for the center of mass coordinate (Equation 16) for a trajectory of (long) duration δ. This variance can be calculated for diffusion along a general continuous curve following Mitra and Halperin's [40] derivation in the case of slab geometry, with slight modifications. One finds

$$\begin{split} V\_{ij} &= \frac{2}{D\delta} \int\_0^\ell \mathrm{d}s \, r\_i(s) \int\_0^\ell \mathrm{d}s' \, r\_j(s') \left\{ B\_2 \left( \frac{|s - s'|}{2\ell} \right) + B\_2 \left( \frac{s + s'}{2\ell} \right) \right\} \\ &+ \frac{\ell^2}{3D\delta} \left[ B\_4 \left( \frac{|s - s'|}{2\ell} \right) + B\_4 \left( \frac{s + s'}{2\ell} \right) \right] \right), \end{split} \tag{19}$$

where Bn(·) denotes the nth order Bernoulli polynomial, and Dδ/ℓ2 ≫ 1. Here, the exponentially decaying terms are ignored as their contribution is negligible at even moderate durations. The details of the derivation of the above expression is provided in Appendix A (see the Supplementary Material). We note that alternative representations of the final result (Equation 19) can be given in terms of polylogarithmic or Hurwitz zeta functions. We verified that Equation (19) correctly reproduces Mitra and Halperin's expression [40] in the same regime for the slab geometry, i.e., for a straight line.

As were in the previous cases, the signal decay tensor for long pulse duration is not rank-1 for a general curve. However, unlike in previous regimes, the compartmental signal decay is truly Gaussian. Consequently, the orientationally averaged signal in this regime suffers q −1 decay if and only if the fibers are straight.

# 3. RESULTS AND DISCUSSION

## 3.1. Summary

The above findings can be summarized as follows: For narrow pulses (Dδ ≪ R 2 c , ℓ 2 ) in regimes A and B, the orientationally averaged signal exhibits a slow q −1 tail for large q for diffusion along curved as well as straight 1 dimensional structures. As the pulse duration is prolonged, as we have studied for regime A, the slow q −1 tail gives way to a steeper drop for substantially curved fibers, by which we mean Rc . √ Dδ. Indeed, the signal may be expected to bear less features of diffusion along a 1D structure and more of diffusion in a 3D domain (i.e., a steep decay), since a long pulse serves to average the motion of the spin carriers over a length ∼ √ Dδ along the path curved in 3D space. Fibers much shorter than the averaging length (ℓ ≪ √ Dδ) fall into regime C and may contribute to an exponent of −1 in the orientationally averaged signal only if they are straight (Rc ≫ ℓ).

# 3.2. Clinical Relevance

An important question to ask is: What regime is the most relevant for clinical MR examinations of the brain? There is no clear answer to this question, essentially because of the extremely wide variability in the size and shapes of cells within the brain [44] 5 . Consequently, it is impossible to suggest that diffusion within all neural cells takes place in a single experimental regime. We can, nonetheless, argue that regime B is the least relevant one as it is impossible to meet the narrow pulse condition along with the long diffusion time condition when δ ≈ 1 as in clinical acquisitions performed at larger (in a practical sense) b-values, where b = q 2 (1−δ/3). The diffusion distance √ 2D1 is expected to be about 10–20 µm in such acquisitions, which will place longer (ℓ ≫ √ D1) arborizations toward regime A while the shorter (ℓ ≪ √ D1) ones will tend to exhibit features of regime C. If the longer structures furthermore exhibit curvature radii safely above the "averaging length" (Rc≫ √ Dδ), the q −1 signature will be possible to observe over a range of large q-values without requiring strict straightness, as demonstrated in **Figure 3**. For structures whose radii of curvature are small (Rc . √ Dδ) the appearance of the slow decay q −1 becomes sensitive to curvature, and is not retained for as wide a range of q-values (see **Figure 4**). At the extreme end of the spectrum, for structures so short that ℓ ≪ √ Dδ, entering into regime C, no curvature is tolerated if the tail q −1 is to appear. These considerations will have to be revisited if one measures the diffusion of molecules other than water, due to differences in their diffusion characteristics [48–50].

# 3.3. On Experimental Findings

In light of the above deliberations, we can revisit the experimental observations of McKinnon et al. [1] who have reported the powder averaged signals stemming from various brain regions dominated by white- or gray-matter. In white-matter regions, they have measured a decay q −c for large q with c values in the neighborhood c ≈ 1.1 ± 0.1, whereas in gray-matter regions, the decay was significantly faster, with c ≈ 1.8 ± 0.2. Concerning the faster decay observed in gray-matter, they propose an explanation based on permeability differences in white and gray-matter regions. If one assumes impermeable membranes, as we did in this study, the following alternative interpretation seems adequate: In white-matter, the observation of the tail q −1 is compatible with regimes A and C, suggesting a substantial presence of fibers that fall into these regimes. The former implies long fibers (ℓ ≫ √ D1) that could in fact be modestly curved (as long as Rc ≫ √ Dδ remains valid). The latter implies short fibers (ℓ ≪ √ Dδ) that are straight. In graymatter regions, the loss of the slow q −1 decay suggests that the signal originates predominantly from fibers that fall outside these descriptions, i.e., exhibiting strong curvature (Rc . √ Dδ). Indeed, gray-matter is rich in dendrites and unmyelinated axons, which will exhibit a fair amount of bending distributed across any voxel. The gray box in **Figure 4** depicts roughly the q range used in McKinnon et al.'s study, where it is seen that while modest curvatures do exhibit the exponent −1 for a while, strong curvature yields a steeper decay for most of the range.

5A collection of neuron images for various species and anatomical regions can be found in the Neuromorpho database, which can be accessed through its web site, http://www.neuromorpho.org. For a recent review on the findings based on this database, see Parekh and Ascoli [45]. We also note [46] wherein the authors employ the approach taken in Jespersen et al. [47] on this database to relate the neuronal morphology in gray-matter to the MR signal at very low diffusion sensitivity.

Another potential explanation involves glial cells, which constitute a substantial portion of all neural cells. Similar to neurons, glial cells exhibit an extraordinary level of diversity in their size and shape throughout the brain. Their relative number and distribution is the subject of ongoing debate [51]. Thus, an accurate assessment of their influence on the detected diffusion MR signal is infeasible at this time. It is known, however, that the glial cells tend to be smaller than neurons, they lack axons, and many of them are star-shaped [52]. These structural features tend to disfavor the emergence of a q −1 decay. However, recent studies have suggested that the glial cells are significantly more prevalent in cerebral white-matter compared to gray-matter [53]. Thus, it can be argued that their contribution to the overall MR signal must be rather limited. We note that the vast variation in neuronal and glial morphology along with the reported regional differences justify future studies performed at high spatial resolutions, for understanding the influence of compositional variations on the diffusion MR signal.

# 3.4. Immobile Water

The detected orientationally-averaged signal would typically include contributions from many different compartments besides the neural projections, including extracellular matrix, cell bodies, and molecules trapped within very small regions (e.g., within certain organelles, between myelin layers, etc.). Among these, molecules diffusing relatively freely (for instance, between the neural processes) are expected to yield a decay rate faster than q −1 so that much of their contribution is expected to disappear at larger q values. Conversely, there would be no significant loss of signal for truly restricted particles. Presence of a substantial portion of signal originating from such compartments as well as noise-induced bias associated with employing magnitude-valued data would make the decay appear slower [54, 55] than q −1 . The decay exponent c ≈ 1.1 ± 0.1 reported by McKinnon et al. [1] suggests that contributions from such immobile spins could be negligible in their acquisitions. This can be attributed [23] to the relatively short transverse relaxation times those spins are expected to have, along with the long echo times employed at larger q-value acquisitions via clinical scanners.

# 3.5. A Possible Error

One may be tempted to employ Equation (1) for the compartmental signal, as it always has a Gaussian form for sufficiently small q values, and then to attribute the emergence of the q −1 behavior of the orientationally averaged signal through Equation (2) to the rank of the decay tensor **V** being 1. However, this may be permissible for large q only in regime C. This attribution would therefore be erroneous for large q in regimes A and B. Specifically, in regimes A and B, the q range in which Equation (1) applies is similar to the q range in which the orientational average (Equation 2) exhibits a Gaussian tail no matter the rank of the tensor **V** (see Appendix A in the Supplementary Material). Thus, the q −1 behavior is not a consequence of the signal decay tensor having rank 1. Importantly, such behavior naturally emerges at larger q-values in regimes A and B for reasonable shapes of neural projections.

# 3.6. Debye-Porod Law

The powder averaged signal may be envisioned as having originated from a porous specimen which is macroscopically isotropic, in which case one expects an asymptote of the form q −4 due to the Debye-Porod law [27] (see Appendix C in the Supplementary Material), while a steeper decay is predicted for longer pulses as the process approaches a Gaussian (Regime C). Since this behavior arises from quite general considerations, the observation of a tail of q −1 appears questionable. Indeed, the observed power is most likely valid only in an intermediate range, as opposed to the strict q → ∞ limit. The asymptotic expansion, whose leading term gives the Debye-Porod law, contains terms that decay faster than any pure power (e.g., exponential). These terms may very well exhibit decays slower than q −4 over a certain stretch of q-values, before the Debye-Porod asymptote takes over, which happens when q begins to compete with the inverse length of the smallest dimension of the pore geometry. Alas, the q values used in McKinnon et al.'s experiments are of the order 0.1µm−1 which is far from large enough compared to the inverse length scale of 1µm−1 afforded by the axon diameter.

# 3.7. Unaccounted Factors

One of the hallmarks of neuronal morphology is the axonal and dendritic branchings [56], which are not accounted for in our treatment of diffusion on curves. Although an accurate assessment of the influence of branchings could be achieved via careful numerical studies [57], the insight gained from our description above can be employed to some extent for making some predictions. For long structures, each branch can be considered a separate segment along which diffusion takes place. Thus, the presence of branchings would not impact the formulation for projections in regime A. However, since the total contour length can be considerably increased in the presence of branchings, it may be more difficult to satisfy the long pulse duration condition of regime C. When this condition is met, however, the rank of the signal decay tensor would almost certainly be greater than one. In other words, detecting q −1 decay in regime C would be nearly impossible for smaller cells unless most arborizations run parallel within a narrow cylindrical region.

The detected MR signal is known to be dependent on factors other than those accounted for in this work. Among these, spatial heterogeneity of magnetic susceptibility within the tissue have been reported to influence the diffusion decay [58, 59] as well. In fact, suppressing [60] or taking advantage [61, 62] of effects due to susceptibility variations is an active area of research. We note that the presence of internal gradients could be investigated as yet another mechanism that could explain the reported features in the orientationally averaged signal; doing so would require an extension of existing studies relating the diffusion MR signal to microscopic perturbations in susceptibility [63, 64].

# 4. CONCLUSION

In an attempt to interpret new experimental findings, we studied the influence of diffusion along parameterized curves on orientationally-averaged diffusion MR signal. We examined the problem in three distinct temporal regimes of the Stejskal-Tanner experiment and investigated the appearance of a slow decay. We have found that for smaller cells, the q −1 decay of the orientationally-averaged signal is predicted only for straight fibers. This decay is more general for cells with longer projections, while it fades away for curvy structures as the pulse duration

# REFERENCES

of the gradient sequence increases. Finally, we stressed that the q −1 decay could represent an intermediate range as the true asymptotic behavior is governed by a steeper attenuation. The findings of this paper are expected to provide insight into the link between the diffusion weighted MR acquisitions and geometry of the neural cells.

# AUTHOR CONTRIBUTIONS

EÖ conceptualized the problem, and performed the main derivations. EÖ and CY developed and refined the theory, performed the numerical simulations, and wrote the manuscript. MH and HK provided inputs to the theory while C-FW provided guidelines. All authors collaborated in bringing the manuscript to its final state.

# FUNDING

This study was supported by the Swedish Foundation for Strategic Research AM13-0090, the Swedish Research Council CADICS Linneaus research environment, the Swedish Research Council 2015-05356 and 2016-04482, Linköping University Center for Industrial Information Technology (CENIIT), VINNOVA/ITEA3 13031 BENEFIT, and National Institutes of Health P41EB015902, R01MH074794, P41EB015898.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy. 2018.00017/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Özarslan, Yolcu, Herberthson, Knutsson and Westin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# On the Vanishing of the t-term in the Short-Time Expansion of the Diffusion Coefficient for Oscillating Gradients in Diffusion NMR

Frederik B. Laun1 , Kerstin Demberg2 , Armin M. Nagel 1 , Micheal Uder 1 and Tristan A. Kuder 2 \*

1 Institute of Radiology, University Hospital Erlangen, Erlangen, Germany, 2 Medical Physics in Radiology, German Cancer Research Center, Heidelberg, Germany

Nuclear magnetic resonance (NMR) diffusion measurements can be used to probe porous structures or biological tissues by means of the random motion of water molecules. The short-time expansion of the diffusion coefficient in powers of t 1/2 , where t is the diffusion time related to the duration of the diffusion-weighting magnetic field gradient profile, is universally connected to structural parameters of the boundaries restricting the diffusive motion. The t 1/2 -term is proportional to the surface to volume ratio. The t-term is related to permeability and curvature. The short time expansion can be measured with two approaches in NMR-based diffusion experiments: First, by the use of diffusion encodings of short total duration and, second, by application of oscillating gradients of long total duration. For oscillating gradients, the inverse of the oscillation frequency becomes the relevant time scale. The purpose of this manuscript is to show that the oscillating gradient approach is blind to the t-term. On the one hand, this prevents fitting of permeability and curvature measures from this term. On the other hand, the t-term does not bias the determination of the t 1/2 -term in experiments.

Keywords: magnetic resonance imaging, diffusion, short-time limit, surface-to-volume ratio, gradient profile, oscillating gradients

# INTRODUCTION

This article builds on and extends a previous article [1], which investigated the term linear in time of the short-time expansion of the diffusion coefficient [2–5] which is given by:

$$D\left(t\right) = D\_0 \left(1 - \frac{4}{3d\sqrt{\pi}} \frac{\mathcal{S}}{V} \sqrt{D\_0 t} + \left(\frac{1}{d} \frac{\mathcal{S}}{V} \kappa + \frac{1}{2d} \frac{\mathcal{S}}{V} \varrho - \frac{1}{12} \frac{\mathcal{S}}{V} D\_0 \overline{\mathcal{R}^{-1}}\right) t + O\left(t^{3/2}\right)\right), \tag{1}$$

where D0 is the free diffusion coefficient, S/V is the surface-to-volume ratio, κ is the surface permeability, ̺ is the surface relaxivity, R−1 is a mean curvature term, d is the spatial dimension, and t is the observation time. This universal expansion is valuable, since it connects a measurable quantity, i.e., D (t), to structural parameters of barriers restricting the diffusive motion.

Using magnetic resonance diffusion experiments [6–9], information about the diffusive motion of spin-bearing particles can be encoded into the signal by using diffusion-weighting magnetic

#### Edited by:

Sune Nørhøj Jespersen, Aarhus University, Denmark

#### Reviewed by:

Valerij G. Kiselev, Universitätsklinikum Freiburg, Germany Junzhong Xu, Vanderbilt University, United States

\*Correspondence:

Tristan A. Kuder<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 03 August 2017 Accepted: 24 October 2017 Published: 14 November 2017

#### Citation:

Laun FB, Demberg K, Nagel AM, Uder M and Kuder TA (2017) On the Vanishing of the t-term in the Short-Time Expansion of the Diffusion Coefficient for Oscillating Gradients in Diffusion NMR. Front. Phys. 5:56. doi: 10.3389/fphy.2017.00056 field gradient pulses. Regarding the diffusion time t linked to the total duration of the diffusion-weighting gradient profile, one often considers the long-time and short-time limit. In the first case, the limit of long diffusion time, detailed information about the porous structure of the investigated material can be obtained [10–12] such as actual pore shapes [13–15]. On the other hand, D (t) can be measured in the short-time limit to obtain the structural parameters in Equation (1). For this purpose, a pair of bipolar gradient pulses is applied to achieve diffusion encoding [16]. In the short gradient pulse approximation [6, 17], the measured diffusion coefficient in such experiments is D (t) (Equation 1). If the gradient pulses cannot be considered to be short, Equation (1) must be modified to take into account the effect of the gradient pulse duration and of the temporal evolution of the gradients G(t) [18–22]:

$$D\_{\rm app}(t) = D\_0 \left( 1 - c\_1 \frac{4}{3d\sqrt{\pi}} \frac{\mathcal{S}}{V} \sqrt{D\_0 t} + c\_2 \left( \frac{1}{d} \frac{\mathcal{S}}{V} \kappa \right. \right. $$

$$- \frac{1}{12} \frac{\mathcal{S}}{V} D\_0 \overline{R^{-1}} \left( t + O\left( t^{3/2} \right) \right) \tag{2}$$

Here, the influence of the temporal gradient profile is expressed solely by the constants c1 and c2, which can be computed from G(t), so that an elegant decoupling takes place. Note that surface relaxation is neglected in Equation (2) and in the remainder of the manuscript, thus avoiding the difficulties in the mathematical treatment [23], and that t is the total duration of the diffusion gradients in Equation (2).

In Laun et al. [1], it was shown that c2 can be tuned to values between 0 and 1. Tuning c2 to zero can be advantageous, for example, if the aim of the experiment is to measure the √ t-term without bias from the t-term.

A striking result [24–30] is that the short-time expansion is moreover valid for the diffusion spectrum D (ω), or Dapp (τ ): = D 2π τ , that can be measured by the use of oscillating gradients. Note that Dapp (τ ) and Dapp (t) are different functions as outlined in more detail below. Then Equation (1) can be cast in a similar form for the diffusion spectrum, where t is replaced by the duration τ of one gradient oscillation:

$$\mathcal{D}\_{\text{app}}\left(\mathbf{r}\right) = D\_0 \left(1 - C\_1 \frac{4}{3d\sqrt{\pi}} \frac{\mathcal{S}}{V} \sqrt{D\_0 \mathbf{r}} + C\_2 \left(\frac{1}{d} \frac{\mathcal{S}}{V} \kappa\right) \right. \tag{3}$$

$$-\frac{1}{12} \frac{\mathcal{S}}{V} D\_0 \overline{R^{-1}} \left(\mathbf{r}\right) + O\left(\mathbf{r}^{3/2}\right) \tag{3}$$

The constants Cn in Equation (3) are printed in capital letters because they differ, in general, from the constants cn of Equation (2) as will be described below.

Equations (2) and (3) were successfully applied in experiments to obtain information about the first order term and thus about the surface-to-volume ratio [27, 28, 31–37]. The thereto necessary constants C1 were derived for some gradient waveforms such as the bipolar waveform and oscillating gradients [18, 27, 29, 38].

The aim of the work at hand is to investigate the constant C2. For this purpose, Equation (3) is derived starting from Equation (1).

# MATERIALS AND METHODS

Numerical simulations were performed as in Laun et al. [1] 1 . The diffusion coefficient Dapp (t) and the diffusion spectrum Dapp (τ ) [using the signal decrease as recalled in Equation (A13) in Appendix A (Supplementary Material)] were computed using the multiple correlation function (MCF) approach [22, 39–45] (using Equation 114 in Grebenkov [42]). The MCF approach decomposes the magnetization into the eigenfunctions of the Laplace operator. One important parameter is the number Nλ of employed eigenfunctions, which should be sufficiently large to ensure numerical accuracy. In the presented results, the accuracy was verified by increasing Nλ and checking whether numerical results remained unchanged. A detailed description of the MCF approach is beyond the scope of this article, but can be found in Grebenkov and Grebenkov [46, 47], for example.

The following closed domains were considered: slab, cylinder, sphere, "bi-slab" (see **Figure 1**). The bi-slab domain consists of three parallel planes. The inner plane is permeable, while the two outer ones are impermeable. Particles only reside within the volume between the two impermeable slabs. The radii of cylinder and sphere were 5µm, the separation of the slabs was 10µm, and the separation of the planes of the bi-slab domain was 10µm (thus the bi-slab domain was in total 20µm wide). The free diffusion coefficient D0 was set to 1µm2 /ms. The boundaries were fully reflecting except for the inner wall of the bi-slab domain, which had a permeability of 50 µm/s. Nλ was 100 for the bi-slab domain, 500 for slab domain and cylinder, and 200 for the sphere. Oscillating cosine gradients were simulated with a total duration Tcos of 0.05, 0.1, and 0.5 s. The number of oscillations N was varied in twenty steps. For bipolar gradients, δ was set to 10−3 ms and t was varied between 0.1 and 15 ms.

Additionally, the difference between simulated diffusion coefficients and first order short time expansion was calculated. This difference is labeled 1D in the plots and represents Dapp,simulated (t) − D0 − M1c1t 1/2 or Dapp,simulated (τ ) − D0 − M1C1τ 1/2 .

# RESULTS

# Derivation of the t-Term for Oscillating Gradients

First, a shorthand-notation for Equation (2) is introduced:

$$D\left(t\right) = \sum\_{n=0} M\_n c\_n t^{n/2} \tag{4}$$

with the coefficients

$$\begin{aligned} M\_0 &= D\_0\\ M\_1 &= -\frac{4}{3d\sqrt{\pi}} \frac{\mathcal{S}}{V} D\_0^{3/2} \\ M\_2 &= \frac{1}{d} \frac{\mathcal{S}}{V} \kappa D\_0 - \frac{1}{12} \frac{\mathcal{S}}{V} \overline{R^{-1}} D\_0^2 \end{aligned}$$

1A considerable overlap of the Methods sections with the corresponding sections of the earlier article is present.

FIGURE 1 | Impact of the t-term on the short-time expansion of the diffusion coefficient. Markers indicate values obtained by numerical simulations. Solid lines represent the short-time expansion to order √ t and dotted lines represent it to order t. Red square markers represent bipolar gradients (i.e., Dapp (t)) and black circular markers represent oscillating cosine gradients (i.e., Dapp (τ )). (Continued)

#### FIGURE 1 | Continued

(A) Slab domain. The t-term is zero because curvature and permeability of the sample are zero. (B–D) Cylinder, sphere, and bi-slab. In case of oscillating cosine gradients, the t-term is zero, because C2 is zero. For this reason, the markers stay close to the solid line in contrast to the markers indicating the bipolar gradients, which stay close to the dotted line. Tcos was 500 ms.

and so on (with c0 = 1).

As outlined in Appendix A (Supplementary Material), the short-time expansion for the position correlation function that generates an experimentally detectable signal attenuation reads:

$$\left<\mathbf{x}(t\_2)\mathbf{x}(t\_1)\right> = -D\left(|t\_{21}|\right)\cdot \left|t\_{21}\right| = -\sum\_{n=0} M\_n \left|t\_{21}\right|^{1+n/2} \tag{5}$$

with t21 = t2−t1, where the brackets h. . .i denote the expectation value. Note that the terms hx(t2) 2 i and hx(t1) 2 i where neglected in Equation (5) because they do not contribute to the signal attenuation. Equation (5) can be related to the diffusion spectrum D (ω) (see Appendix A in Supplementary Material) via:

$$\begin{split} \frac{\mathfrak{D}\left(\omega\right)}{\omega^{2}} &= \frac{1}{2} \int\_{-\infty}^{\infty} \langle \mathbf{x}\left(t\_{2}\right) \mathbf{x}\left(t\_{1}\right) \rangle \ e^{-i\omega t\_{21}} dt\_{21} \\ &= -\frac{1}{2} \sum\_{n=0} \int\_{-\infty}^{\infty} M\_{n} |t\_{21}|^{1+n/2} e^{-i\omega t\_{21}} dt\_{21} . \end{split} \tag{6}$$

This Fourier integral exists (see Appendix B in Supplementary Material):

$$\begin{aligned} \frac{1}{2} \int\_{-\infty}^{\infty} |t\_{21}|^{1+n/2} e^{-i\alpha t\_{21}} dt\_{21} &= -\omega^{-2-\frac{n}{2}} \cos\left(\frac{n\pi}{4}\right) \Gamma\left(2+\frac{n}{2}\right),\\ \omega &> 0, \ n \ge 0 \end{aligned} \tag{7}$$

and thus by inserting Equation (7) in Equation (6), one finds:

$$\mathfrak{D}(\omega) = D\_0 + M\_1 \sqrt{\frac{\pi}{2}} \frac{3}{4\alpha^{1/2}} + 0 \cdot M\_2 + 0 \left(\omega^{-3/2}\right). \tag{8}$$

Note that the gamma function Γ makes the constants cn increase swiftly at larger n. Defining the time parameter τ = t/n, entailing ω = 2π/ τ , one finds:

$$\mathfrak{D}\_{\rm app}(\mathfrak{r}) := \mathfrak{D}\left(\frac{2\pi}{\mathfrak{r}}\right) = D\_0 + M\_1 \frac{3}{8} \mathfrak{r}^{1/2} + \mathfrak{d} \cdot M\_2 + O\left(\mathfrak{r}^{3/2}\right) \\ \text{(9)}$$

with Dapp (τ ) being identical to D (ω) except for taking a different argument. Dapp (τ ) has exactly the form of Equation (3) as desired and one can read off the coefficients Cn directly: C1 = 3/8 and C2 = 0. Note that the value of C1 was reported previously (e.g., in Novikov and Kiselev [29]). The vanishing of C2 has not been reported so far to our knowledge.

Using the expression for M1, one finds:

$$\mathcal{D}\left(\omega\right) = D\_0 - \frac{D\_0}{d\sqrt{2}} \frac{\mathcal{S}}{V} \sqrt{\frac{D\_0}{\omega}} + O\left(\omega^{-3/2}\right),\tag{10}$$

$$\mathfrak{D}\_{\text{app}}\left(\tau\right) := \mathfrak{D}\left(\frac{2\pi}{\tau}\right) = D\_0 - \frac{D\_0}{2d\sqrt{\pi}} \frac{\mathcal{S}}{V} \sqrt{D\_0} \tau^{1/2} + \mathcal{O}\left(\tau^{3/2}\right) \cdot \tag{11}$$

It is interesting to calculate the coefficients c1 and c2 for a short-time cosine gradient with one oscillation [with methods as described, e.g., in Laun et al. [1] and references therein]. We find c1 = 3 · 4π FresnelC(2) + 3 FresnelS(2) /16/π ≈ 0.428 and c2 = 0. These values bear great similarity to C1 ≈ 0.375 and C2 = 0. It should be noted that c2 of oscillating cosine gradients with any number of oscillations equals zero because they are "flow-compensated," i.e., because their first moment vanishes [1].

# Validation with Simulations

**Figure 1** displays Dapp (t) and Dapp (τ ). Markers indicate simulation results using the MCF approach and lines represent the short-time expansion. For Dapp (t), solid lines equal M0 + M1c1t 1/2 and dotted lines equal M0 + M1c1t 1/2 + M2c2t. For Dapp (τ ), solid lines equal M0 + M1C1τ 1/2 and dotted lines equal M0 + M1C1τ 1/2 + M2C 2 1 τ . The term M2C 2 1 τ shall represent a reasonable "guess" for the t-term with an effective diffusion time τeff = C 2 1 τ , where the coefficient C2 was set to one. The intention is to visualize a line with C2 6= 0, although this term does not occur in reality. Some remarks on effective diffusion times can be found in Appendix C (Supplementary Material).

For the slab domain (**Figure 1A**), Mn>1 = 0 holds true (see [42]). Hence, **Figure 1A** does not display a dotted line and markers stay close to the solid lines.

In **Figures 1B,C** (cylinder, sphere) and **Figure 1D** (bi-slab), it is clearly visible that the markers for the bipolar gradients (with c2 = 1 6= 0) stay close to the dotted lines indicating the importance of the t-term. The markers of the oscillating cosine gradients stay close to the solid line indicating that the t-term does not influence Dapp (τ ). Owing to higher order terms, deviations between the short-time expansion and markers are present at larger t.

**Figure 2** shows 1D, i.e., the difference between simulated diffusion coefficients and first order short time expansion. The dotted line represents the t-term, i.e., M2c2t for the bipolar gradients. For the oscillating gradients, the black dotted line shall represent an educated guess for the t-term, i.e., M2C 2 1 τ , as in **Figure 1**.

First, the bipolar gradients displayed in **Figure 2** are discussed (displayed in red color). For the slab domain (**Figure 2A**), the dotted line is flat because M2 = 0. However, deviations of 1D from zero are well visible for t > 10 ms. This does not result from the influence of higher order terms because all higher order terms are zero. It rather indicates the breakdown of the short-time expansion. For cylinder, sphere, and bi-slab (**Figures 2B–D**), the slope of 1D is identical to that of the dotted line at t = 0, but starts deviating already roughly at t = 2 ms indicating that either higher order terms are needed or, again, that the short-time expansion breaks down. This deviation is more pronounced for cylinder and sphere than for the bi-slab.

Next, the oscillating gradients in **Figure 2** are discussed (displayed in black color). For cylinder, sphere, and bi-slab, 1D does have zero slope at t = 0 and does not follow the dotted line for any of these domains, which supports the finding that C2 = 0. This holds true for Tcos = 500 ms, but also for reduced total duration of the oscillating gradients, i.e., for smaller Tcos. The difference of 1D between Tcos = 500 ms and Tcos = 50 ms is smaller than 0.015 µm2 /ms for all domains at τ = 10 ms, which is roughly equally large as the guessed t-term, but an order of magnitude smaller than the √ t-term. Thus, for the considered domains, Tcos = 50 ms still appears to be well suited for investigations of the √ t-term, even with as few as five oscillations.

# DISCUSSION

The main result of this work is that oscillating cosine gradients are blind with respect to the t-term of the short-time expansion of the apparent diffusion coefficient.

Oscillating gradients and extensions [48–51] have been used in several research studies [27, 28, 32, 35, 38, 52–65], among them applications to human brains in vivo [66]. Comparing oscillating gradients to pulsed gradients, the advantage of the oscillating gradients is that the obtainable b-value is higher allowing the assessment of shorter times. This is particularly useful if strong gradient amplitudes are not available or if the structure of interest is too small. The disadvantage is the need for longer echo times entailing decreased signal-to-noise ratio due to transversal relaxation, which also entails a longer acquisition time.

As oscillating gradients are blind to the t-term, estimates of S/V as in Reynaud et al. [36] are not biased by this term, but, obviously, the membrane permeability, for example, cannot be estimated using the t-term. This is in line with the findings by Li et al. [67], who reported that the membrane permeability has little effect on oscillating gradient derived diffusion coefficients at high frequencies. This is presumably not a major limitation given the smallness of the t-term that is visible in **Figures 1**, **2**, which makes a fit challenging. The permeability information must have some influence on Dapp (τ ) at long τ ; otherwise diffusion in the bi-slab would have to be identical to that of a single slab domain of double size. Therefore, the estimation of membrane permeability using oscillating gradients might in principle be possible.

As different versions of Equations (10) and (11) can be found in the literature, a comparison is worthwhile. Equation (10) is identical to Equation 10 of the article by Novikov and Kiselev [29]. Except for a small deviation, which may be due to numerics, Equation (10) is also identical to Equation (3) of the article by Xu et al. [35], but, to our understanding, not to the respective equations in an earlier article [27]. In general, care must be taken concerning the definition of τ . For example, Zielinski et al. [38] use the definition τZielinski = τ/2, which is closer to the classical timing definitions of CPMG echo trains than our definition. Considering this difference in definitions,

# REFERENCES

their respective coefficient C1 for the CPMG condition as stated in their Equation (6) is almost identical to 3/8, which is in agreement with the finding that the difference between C1 of CPMG and cosine gradients should be almost negligible as stated in section 3.3 of Novikov and Kiselev [29]. Further, we found our coefficient C1 to be a factor of six smaller than the one stated in Equation 14 in the article by Stepišnik et al. [28]. This difference was noted by the authors themselves and in Novikov and Kiselev [29].

Interestingly, the disappearance of the t-term in the Mitra expansion of Equation (3) using oscillating gradients is due to its disappearance in D (ω), or Dapp (2π/τ ), respectively. Thus, optimizing oscillating gradient profiles instead of using, for example, just cosine gradients, which was a successful approach in other regards [68, 69], does not help to make the t-term reappear in the signal attenuation.

In practice, diffusion measurements use spin echoes and hence two gradients at both sides of the refocusing pulse (as in Baron and Beaulieu [66]). This effectively introduces an extra variable, i.e., the separation of two gradients, which can affect the spectrum of diffusion gradients. When interpreting oscillating gradient experiments, this effect must be taken into account.

A limitation of the presented simulations is that they cannot prove the disappearance of the t-term. In principle, a very small t-term might be present and go unnoticed.

In conclusion, oscillating gradients are blind to the t-term and hence no bias in fits of the surface-to-volume ratio arises from the t-term.

# AUTHOR CONTRIBUTIONS

FL performed the simulations and initial computations. TK and MU were involved in the design of the evaluations. TK, AN, KD, and FL were involved in implementation and testing of the MCF code and of the mathematical derivations.

# FUNDING

Financial support by the DFG (grant numbers SFB TRR 125/2 R01 and KU 3362/1-1) is gratefully acknowledged.

# SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy. 2017.00056/full#supplementary-material

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with the authors TK and FL.

Copyright © 2017 Laun, Demberg, Nagel, Uder and Kuder. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Effective Potential for Magnetic Resonance Measurements of Restricted Diffusion

Evren Özarslan1 \*, Cem Yolcu1 , Magnus Herberthson2 , Carl-Fredrik Westin1, 3 and Hans Knutsson1

1 Department of Biomedical Engineering, Linköping University, Linköping, Sweden, 2 Division of Mathematics and Applied Mathematics, Department of Mathematics, Linköping University, Linköping, Sweden, 3 Laboratory for Mathematics in Imaging, Department of Radiology, Harvard Medical School, Brigham and Women's Hospital, Boston, MA, United States

The signature of diffusive motion on the NMR signal has been exploited to characterize the mesoscopic structure of specimens in numerous applications. For compartmentalized specimens comprising isolated subdomains, a representation of individual pores is necessary for describing restricted diffusion within them. When gradient waveforms with long pulse durations are employed, a quadratic potential profile is identified as an effective energy landscape for restricted diffusion. The dependence of the stochastic effective force on the center-of-mass position is indeed found to be approximately linear (Hookean) for restricted diffusion even when the walls are sticky. We outline the theoretical basis and practical advantages of our picture involving effective potentials.

#### Edited by:

Itamar Ronen, Leiden University, Netherlands

#### Reviewed by:

Silvia Capuani, Consiglio Nazionale Delle Ricerche (CNR), Italy Christian Herbert Ziener, German Cancer Research Center, Germany

#### \*Correspondence:

Evren Özarslan<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 29 September 2017 Accepted: 08 December 2017 Published: 19 December 2017

#### Citation:

Özarslan E, Yolcu C, Herberthson M, Westin C-F and Knutsson H (2017) Effective Potential for Magnetic Resonance Measurements of Restricted Diffusion. Front. Phys. 5:68. doi: 10.3389/fphy.2017.00068 Keywords: diffusion, magnetic resonance, anisotropy, microscopic, local, confinement, tensor, parabolic

# 1. INTRODUCTION

Recovering the mesoscopic structure of porous media and biological tissues via diffusion sensitized NMR methods has been an active area of research since 1960s [1, 2]. As a molecule diffuses within the medium, its magnetic moment acquires a certain phase depending upon the particle's trajectory and the impressed magnetic field gradient profile. The total signal generated by all molecules can be expressed mathematically as

$$E = \left\langle e^{-i\boldsymbol{\gamma}\cdot\int\_0^t \mathbf{ds} \cdot \mathbf{G}(s) \cdot \mathbf{x}(s)} \right\rangle,\tag{1}$$

where γ is the gyromagnetic ratio while **G** and **x** denote the time-dependent magnetic field gradient and the position of the particle, respectively. The average is taken over all particle trajectories.

Conventional techniques for relating the NMR signal above to microstructural features of the medium vary from "localized" models in which the aggregate signal is envisioned to arise from isolated (e.g., restricted) compartments [3, 4] to more "global" models, which attempt to capture the medium's disorder [5, 6]. Here, we propose to approach the problem of relating the NMR signal to microstructural features of the medium with an alternative paradigm wherein diffusion is thought to be taking place within a potential energy landscape.

When the gradient waveform, **G**(s), comprises two well-separated pulses of infinitesimal durations, the relationship between the NMR signal and the compartment shape is the same as that in scattering experiments. For closed pores, the signal intensity is given by E(**q**) = | ˜ρ(**q**)| 2 , i.e.,the power spectrum of the equilibrium distribution of spins ρ(**r**). Here, the wavevector **q** is taken to be the time integral of each gradient pulse multiplied by γ . As the pulse duration is prolonged, the oscillatory diffraction pattern of the signal stretches toward larger q-values (q = |**q**|), and eventually becomes unobservable [7, 8]. This is a manifestation of the loss of information in experiments featuring long pulses; such information loss is key to our endeavor to identify an effective potential for the diffusion process.

In this Perspective, we argue that the theory of diffusion under a Hookean restoring force [7, 9–13] can be regarded as the effective theory of restricted diffusion for a wide class of highly relevant NMR signal acquisition scenarios. This is evident at the macroscopic scale when one considers the dependence of the signal on experimental parameters for sequences featuring long gradient pulses or plots the average effective force experienced by the particles against their mean position (center of mass of their trajectories) during the application of a long pulse. For NMR examinations of microscopic diffusion anisotropy [14, 15], the quadratic potential profile is ideally suited to represent the observable properties of small pores making up a compartmentalized medium.

# 2. EFFECTIVE POTENTIAL FOR RESTRICTED DIFFUSION

Here we demonstrate that, under certain experimental conditions, the influence of restricted diffusion is essentially the same as that for the Hookean potential model, which was studied in-depth recently [13]. After the theoretical grounds for the effective theory are established, we proceed with presenting simulation results that provide additional justification.

## 2.1. Long-Pulse Train

Consider a pulse sequence which consists of N rectangular gradient pulses of varying (vector) magnitudes **G**n and durations δn. With the definition

$$\mathbf{q}\_n = \boldsymbol{\chi} \mathbf{G}\_n \boldsymbol{\delta}\_n,\tag{2}$$

the NMR signal (Equation 1) is given by

$$E = \left\langle \exp\left\{-i\sum\_{n=1}^{N} \frac{1}{\delta\_n} \int\_{t\_n}^{t\_n + \delta\_n} \mathbf{q}\_n \cdot \mathbf{x}(t) \mathbf{d}t \right\} \right\rangle,\tag{3}$$

where tn denotes the leading edge of the nth pulse. Angular brackets indicate averaging over the possible trajectories **x**(t).

Further, introducing the stochastic center-of-mass of the particle trajectory [7] during the nth pulse through

$$\mathfrak{F}\_n = \frac{1}{\delta\_n} \int\_{t\_n}^{t\_n + \delta\_n} \mathfrak{x}(t) \,\mathrm{d}t,\tag{4}$$

the stochastic signal can be rewritten as

$$E = \left\langle e^{-i\left(q\_1 \cdot \mathfrak{k}\_1 + q\_2 \cdot \mathfrak{k}\_2 + ... + q\_N \cdot \mathfrak{k}\_N\right)} \right\rangle. \tag{5}$$

In words, the NMR signal elicited by spatially constant gradient pulses of finite duration is sensitive to the center of mass (average position) of Brownian trajectories, rather than instantaneous coordinates. The duration of the pulses therefore serve to smear out fine spatial features.

Let pcm(ξ , δ) denote the distribution of the center of mass during a time interval of duration δ. In the long duration regime, δn → ∞, the dependence on the pulse separation disappears and the signal intensity (Equation 5) factorizes, leading to

$$E = \tilde{p}\_{\rm cm}(\mathfrak{q}\_1, \delta\_1) \tilde{p}\_{\rm cm}(\mathfrak{q}\_2, \delta\_2) \dots \tilde{p}\_{\rm cm}(\mathfrak{q}\_N, \delta\_N), \tag{6}$$

where

$$\tilde{p}\_{\rm cm}(q,\delta) = \int p\_{\rm cm}(\xi,\delta) e^{-iq\cdot\xi} \,\mathrm{d}^3\xi \tag{7}$$

is the Fourier transform of the center of mass distribution. Due to its construction (Equation 4), the random variable ξ has a distribution pcm(ξ , δ) that approaches a Gaussian in the limit of long duration1 . Therefore, both its Fourier transform (Equation 7) and hence the signal (Equation 6) approach Gaussians as well.

The Gaussian distribution is determined simply by its variance matrix (with its mean set to zero for convenience). The relevance of this fact here is that when a train of long pulses is employed, the signal has no means of encoding for fine features of the microscopic environment where diffusion takes place2 . What survives in the signal is merely that component of the microscopic forces that provide a stable equilibrium (since the particles do not escape to infinity) for the center of mass and no more. More specifically, a Hookean force.

While restricted diffusion becomes a difficult problem to tackle in higher dimensions for all but a few special geometries, diffusion under a Hookean force is much more tractable, and it behooves one to adopt the latter model when its features are all that can be observed, as argued above. Hence we consider the case of diffusing particles subject to a (dimensionless) parabolic confining potential U(**x**) = (1/2)**x** ⊺**Cx**, where **C** is the confinement tensor [13]. Under this potential, the magnetization density evolves according to the Bloch-Torrey-Smoluchowski equation [13, 17, 18]. For this process, the center of mass distribution is Gaussian for any duration, and its Fourier transform is given by

$$
\tilde{p}\_{\rm cm}(q,\delta) = e^{-\frac{1}{2}q^{\tau}Vq},\tag{8}
$$

with the variance matrix

$$\mathbf{V} = 2(D\delta)^{-2}\mathbf{C}^{-3}(\mathbf{D}\mathbf{C}\delta + e^{-DC\delta} - 1),\tag{9}$$

which is just a straightforward generalization of an expression in Mitra and Halperin [7] to higher dimensions. Via Equation (6), the signal encodes the confinement tensor, **C**, for the restricted region.

1A rigorous proof of this can be found in the mathematics literature [16].

2 Indicated also by the absence of powers of "momentum" **q** higher than 2 in its cumulant expansion.

# 2.2. Insights from the One-Dimensional Problem

The eigenvalues of the inverse of the confinement tensor **C** can be associated with the square of three lengths pertaining to the overall size of a 3-dimensional region. Of course, these three lengths cannot characterize the irregular boundary of a general 3-dimensional region, as alluded to above. However, for the illustration of the size-confinement correspondence, we momentarily take up the simple 1-dimensional example of diffusion restricted by two parallel plates separated by a distance L. The variance of the center of mass position is provided by Mitra and Halperin [7] to be

$$\left< \xi(\delta)^2 \right>\_{\text{res}} = \frac{16L^4}{\pi^6 D \delta} \sum\_{n=1,3,5,\dots}^{\infty} \frac{1}{n^6} \left[ 1 - \frac{1 - e^{-\pi^2 n^2 D \delta/L^2}}{\pi^2 n^2 D \delta/L^2} \right],\tag{10}$$

where the subscript stands for restricted diffusion3 .

In the long pulse duration regime wherein the statistics of the restricted problem approach those of the confinement problem, the variances of the center of mass position for the two problems are found by taking the long time limit of Equations (9) and (10), respectively,

$$\left<\xi(\delta)^2\right>\_{\text{con}} \sim \frac{2}{C^2 D \delta}, \quad \delta \to \infty,\tag{11}$$

$$\left<\xi(\delta)^2\right>\_{\text{res}} \sim \frac{L^4}{60D\delta}, \quad \delta \to \infty. \tag{12}$$

Demanding that the two variances above agree results in the relation

$$C = \frac{\sqrt{120}}{L^2},\tag{13}$$

meaning that the confinement model with parameter C becomes the appropriate effective theory for the restriction model (the "true" theory) with length L, yielding the same observed quantities under the coarse-graining furnished by the long averaging pulse.

While in the above it was the long pulse duration that ensured the Gaussianity of the signal (through the center of mass distribution) for restricted diffusion, an alternative situation is one where the gradient strength is so small that qL is much smaller than unity. Because the observed signal is Gaussian for such q values, the confinement model would perfectly represent the data, yet Equation (13) would not apply since the pulse duration is not necessarily long. If one adopts a matching scenario as above, a correspondence can be obtained numerically for general pulse durations. As shown in **Figure 1A**, the effective confinement value obtained by matching the variances of the center-of-mass distributions for the two problems, which is equivalent to matching the signals in an experiment employing well-separated pulses, varies within a very narrow range, namely between C = 12/L 2 and C = √ 120/L 2 . Therefore, employing the asymptotic correspondence (Equation 13) in lieu of the numerical one (blue curve in **Figure 1A**) would still provide an acceptable estimate of L.

FIGURE 1 | (A) Effective dimensionless confinement value (gray circles) computed by matching the variances of the center of mass distributions in the two problems plotted vs. the dimensionless time Dδ/L 2 . The variation in CeffL 2 value is well-captured (error less than 0.12%) by the expression y = 12 − (12 − √ 120)(αx) kγ [1 + (αx) k ] −γ , with α = 9.495, γ = 1.210, and k = 1.266 (blue line). The dashed line indicates the asymptotic value √ 120. (B) The absolute value of the difference in the normalized signals implied by the two problems (restricted diffusion with separation L, and the Hookean potential with a confinement value taken to be C = √ 120L −2 ) for a traditional Stejskal-Tanner sequence [32] with two gradient pulses whose leading edges are separated by 1 with D1/L 2 = 100. The restricted diffusion signal was computed using the multiple correlation function method [33, 34].

Contrary to the two cases above, the signal for restricted diffusion is not Gaussian for other values of q or δ. Thus, the signal profile afforded by the confinement model can only be an approximation to the true signal when the pulses are not long. To investigate the discrepancy, it is insightful to compare the signal for the restricted diffusion problem (for a slab geometry of separation L) to that for the Hookean model featuring the longtime asymptotic value (Equation 13) of the effective confinement. In **Figure 1B**, we show that the error (defined as the difference in the normalized signal intensities implied by the restriction and confinement pictures) is below 0.01 for a very large portion of the relevant parameter space that could be probed. Such differences are not detectable when the signal-to-noise ratio is 100 or less, thus justifying the confinement description as a

3Coordinates are chosen such that hξ i = hxi = 0.

reliable substitute in studies involving restricted diffusion. In fact, the most significant exception to this statement is at the peaks of the diffraction pattern, which occur roughly at qL = 3/2, 5/2, . . . for very short pulses (i.e., along the left edge of this figure) and at slighly larger qL values for longer durations. For long pulse durations, the diffraction pattern does not prevail, and the confinement model eventually provides an accurate description of the detectable signal.

Consider, e.g., a 4 µm pore—roughly the size of a yeast cell. Assuming a bulk diffusivity value of 2 µm2 /ms, and for a pulse duration of 20 ms, which is typical for diffusion measurements via clinical scanners, the Dδ/L 2 value is 2.5. Thus, the signal response of the confinement picture would be indistinguishable from a much more elaborate theory involving restricted diffusion. The situation is even more favorable for smaller sized compartments. In laboratory spectrometers, however, pulse durations an order of magnitude shorter are feasible; in such scenarios, detectable differences between the restricted diffusion and confinement pictures can be encountered unless the pores are smaller. However, let us remark once more that it is only in this simplified one-dimensional scenario that there exists a single unambiguous size of the region, allowing a one-to-one correspondence to be drawn between C and L. In higher dimensions, one's measurements have access simply to the tensor **C**; the actual dimensions of the region can only be speculated upon to the extent that its geometry is known.

# 2.3. Boundary Force: Simulations

The confinement model is based on a Hookean force assumption, i.e., the presence of a restoring force whose magnitude increases linearly as the particles move away from an attractive center. We shall define an effective force, Feff, based on the impulse the particles experience during a time interval of duration δ, i.e.,

$$F\_{\rm eff} \delta = \int\_0^\delta F(t) \,\mathrm{d}t.\tag{14}$$

Because the time-dependent force F(t) is proportional to −Cx(t) for particles subject to a quadratic potential, the effective force is proportional to −Cξ (δ), implying a perfectly linear dependence on the center-of-mass of each and every trajectory. On the other hand for restricted diffusion, the effective force is proportional to −1N/δ, where 1N denotes the number of collisions the particle makes with the wall on the left subtracted from that with the wall on the right.

By performing random walk simulations, we investigated whether a similar (Hookean) force model could emerge for the restricted diffusion process as well. In the top panel of **Figure 2A**, we plot the histograms of mean particle position for different durations. As expected, the center of mass distribution approaches a Gaussian as the time is prolonged. Note that the net force is a stochastic variable just like the mean position of the trajectories. Thus, in the bottom panel of **Figure 2A**, we plot the average of the net effective force against the mean trajectory position. Remarkably, the effective force is linearly related to the mean position for an overwhelming portion of the particles even for moderate times. Moreover, by repeating the simulations for different L values (results not shown), we found that the slope of the curves at long durations, i.e., the effective force constant, is proportional to L −2 just like in Equation (13). These findings further support the idea of employing a Hookean force model as a substitute for restricted diffusion at long durations. Importantly, the Hookean model is valid even for more complicated problems, which are difficult to treat analytically. We illustrate this point by introducing stickiness to the walls of the same restricted geometry. As shown in **Figure 2B**, the same conclusions can be drawn for this scenario albeit for a different value of the effective confinement.

# 3. DISCUSSION AND CONCLUSION

While characterizing a compartmentalized specimen or biological tissue via diffusion NMR, one is faced with the problem of determining a reliable representation of the local compartment. An accurate mathematical description of the pore shape would typically necessitate numerous parameters. However, such parameters are simply unavailable in diffusion NMR measurements featuring long gradients as a result of coarse-graining associated with the diffusion process taking place during the application of the gradients. The confinement tensor model [13] contains just the right number of parameters as dictated by the central limit theorem.

Similar to the diffusion tensor, the confinement tensor is real and symmetric, thus is described by 6 independent numbers. However, the confinement tensor model has at least one extra parameter—the bulk diffusivity, which can be scalar (isotropic) or tensorial. For the latter case, the diffusion tensor model emerges simply in the C → 0 limit [13]. It is common to express the detected NMR signal for complex specimens as the sum of contributions from Gaussian and restricted compartments [3]. Our work implies that, when appropriate acquisition parameters are employed, both these compartments could be represented by the confinement tensor model though with different confinement values.

Diffusion tensors have also been employed for representing diffusion in microscopic compartments [4, 19, 20]. Since the two models (employing diffusion and confinement tensors) have the same q dependence for traditional Stejskal-Tanner measurements performed with varying gradient vectors, their predictions of the signal would be the same, and one model cannot be preferred over the other; this was observed when the confinement tensor is employed to represent the voxel-level signal [21] similar to what is done in diffusion tensor imaging [22]. However, the diffusion tensor model is based on free diffusion, and exhibits a time-dependence which is paradoxical for a compartmentalized structure. Manifestations of the difference in time-dependencies would be observable, e.g., when the timing parameters of the measurements are varied [23]. In fact, the influence of the gradient waveform on the free-diffusion signal is fully captured [24] by a measurement tensor, referred to as the b-matrix [25]. Thus, any set of acquisitions featuring the same b-matrix would yield the same signal value if diffusion tensor model is employed. However, the quadratic potential profile would lead to qualitatively different signal profiles, approximating the response

Different curves represent different time intervals, which can be associated with durations of the pulses in NMR measurements. The distributions of the positions of the center of mass are shown at the top. The average effective force exerted on the diffusing particles plotted vs. the trajectories' center of mass (bottom). (B) The simulations were repeated by introducing surface adsorption effect so that the particles spend, on average, 25% of the time stuck at the walls.

of a restricted diffusion process, which is more consistent with the underlying compartmentalized structure.

From a practical point of view, treating restricted diffusion is quite difficult even for simple geometries. For example, if we consider an ellipsoidal pore—the simplest geometry with the same number of parameters as a general confinement tensor the problem would be very difficult to solve analytically. We argue that doing so would also be unnecessary if long pulses are employed. The separability of the confinement model enables one to write [13] the expected signal simply in terms of the solution of the one-dimensional Hookean problem, provided that we work in a frame of reference whose axes coincide with those of the ellipsoid; in this frame, **C** is diagonal. Such simplicity due to separation of variables in the confinement model is valid even for more complicated shapes since the statistics of the center of mass converge to Gaussian as long as the durations of the gradient pulses are long.

We would like to note that, as shown by Bauer et al. [26] accumulation of Gaussian phase, if true, implies a propagator (not just stationary distribution) that is Gaussian. In fact, the Ornstein-Uhlenbeck propagator [17, 27, 28] associated with the quadratic potential profile is the general form following the basic properties of the process: diffusive, Markovian, stationary, with Gaussian transition probabilities [26].

In a very recent study [29], some of us considered a general potential energy landscape, U(**x**), in which all disturbances to otherwise free diffusion are envisioned to originate from the variations in this landscape. For this scenario, the equilibrium density of the particles is governed by Maxwell-Boltzmann statistics, i.e., ρ(**x**) ∝ e −U(**x**) , which could account for inhomogeneities in the bulk, and incorporate the effects of boundaries if needed, in a mathematically wieldy manner. We demonstrated that pulse sequences involving very short and long pulses [30, 31] could be employed to map this underlying landscape. Here, we have extended this finding to effective potential energy landscapes associated with restricted diffusion, which can be mapped without narrow pulses.

In conclusion, we have presented a new perspective in which findings of common and relevant NMR signal acquisition scenarios can be interpreted. At the heart of this perspective lies modeling the diffusion as taking place in an effective quadratic potential landscape instead of a restricted domain. We have argued that when probed via waveforms featuring long pulses, the two models become indistinguishable, and the signal should rather be taken to reflect the parameters of such an effective model. Simulations suggest that the stochastic effective force has a linear (Hookean) dependence on the average particle position. The signal for quadratic potential indeed provides a very good approximation to that for restricted diffusion in small (micron-scale) pores when examined via commonly available hardware (see **Figure 1B**). This observation can be generalized in a straightforward way to higher dimensions, thus making the confinement tensor model [13] relevant for applications such as clinical and pre-clinical MRI.

# REFERENCES

# AUTHOR CONTRIBUTIONS

EÖ conceptualized the problem. EÖ, CY, MH, and HK developed and refined the perspective. HK and EÖ performed the numerical simulations. CY and EÖ wrote the manuscript. C-FW provided guidelines. All authors collaborated in bringing the manuscript to its final state.

# FUNDING

This study was supported by the Swedish Foundation for Strategic Research AM13-0090, the Swedish Research Council CADICS Linneaus research environment, the Swedish Research Council 2015-05356 and 2016-04482, Linköping University Center for Industrial Information Technology (CENIIT), VINNOVA/ITEA3 13031 BENEFIT, and National Institutes of Health P41EB015902, R01MH074794, P41EB015898.

34. Grebenkov DS. NMR survey of reflected Brownian motion. Rev Mod Phys. (2007) **79**:1077–137. doi: 10.1103/RevModPhys.79.1077

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2017 Özarslan, Yolcu, Herberthson, Westin and Knutsson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Improving the Realism of White Matter Numerical Phantoms: A Step toward a Better Understanding of the Influence of Structural Disorders in Diffusion MRI

Kévin Ginsburger 1,2 \*, Fabrice Poupon2,3, Justine Beaujoin1,2, Delphine Estournet 1,2 , Felix Matuschke4 , Jean-François Mangin2,3,5, Markus Axer 4 and Cyril Poupon1,2

1 CEA DRF/ISVFJ/Neurospin/UNIRS, Gif-sur-Yvette, France, 2 Université Paris-Saclay, Orsay, France, 3 CEA DRF/ISVFJ/Neurospin/UNATI, Gif-sur-Yvette, France, 4 Research Centre Jülich, Institute of Neuroscience and Medicine, Jülich, Germany, 5 CATI, Multicenter Neuroimaging Platform, Orsay, France

#### Edited by:

Julien Valette, Commissariat à l'Energie Atomique et aux Energies Alternatives (CEA), France

#### Reviewed by:

Olivier Reynaud, École Polytechnique Fédérale de Lausanne, Switzerland Marco Palombo, University College London, United Kingdom

> \*Correspondence: Kévin Ginsburger<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biomedical Physics, a section of the journal Frontiers in Physics

Received: 29 November 2017 Accepted: 02 February 2018 Published: 20 February 2018

#### Citation:

Ginsburger K, Poupon F, Beaujoin J, Estournet D, Matuschke F, Mangin J-F, Axer M and Poupon C (2018) Improving the Realism of White Matter Numerical Phantoms: A Step toward a Better Understanding of the Influence of Structural Disorders in Diffusion MRI. Front. Phys. 6:12. doi: 10.3389/fphy.2018.00012 White matter is composed of irregularly packed axons leading to a structural disorder in the extra-axonal space. Diffusion MRI experiments using oscillating gradient spin echo sequences have shown that the diffusivity transverse to axons in this extra-axonal space is dependent on the frequency of the employed sequence. In this study, we observe the same frequency-dependence using 3D simulations of the diffusion process in disordered media. We design a novel white matter numerical phantom generation algorithm which constructs biomimicking geometric configurations with few design parameters, and enables to control the level of disorder of the generated phantoms. The influence of various geometrical parameters present in white matter, such as global angular dispersion, tortuosity, presence of Ranvier nodes, beading, on the extra-cellular perpendicular diffusivity frequency dependence was investigated by simulating the diffusion process in numerical phantoms of increasing complexity and fitting the resulting simulated diffusion MR signal attenuation with an adequate analytical model designed for trapezoidal OGSE sequences.This work suggests that angular dispersion and especially beading have non-negligible effects on this extracellular diffusion metrics that may be measured using standard OGSE DW-MRI clinical protocols.

Keywords: diffusion time-dependence, white matter microstructure, trapezoidal OGSE sequences, axonal diameter, Monte-Carlo simulations, biomimicking numerical phantoms

# 1. INTRODUCTION

Diffusion-weighted magnetic resonance imaging (dMRI), sensitized to the diffusive motion of water along the direction of an applied magnetic field gradient, has become a well-established technique to non-invasively probe the cellular organization of tissues in vivo. Diffusion NMR measurements embed some information about the inhibition of particles motion due to the presence of barriers in the local environment, and can therefore be exploited to map some specific microstructural features characterizing the brain white matter ultrastructure at cellular scales. While Pulsed Gradient Spin Echo (PGSE) sequences [1] are still widely used in clinical routine, alternative schemes such as Oscillating Gradient Spin Echo (OGSE) sequences seem a promising approach since they enable to explore the diffusion pattern in the frequency domain dual to the diffusion time domain and are able to probe shorter diffusion times compared to conventional PGSE. Some OGSE studies have reported that, at frequencies below 400 Hz, the OG-measured extra-axonal diffusivity transverse to axons in white matter is linearly dependent on the frequency of the employed OGSE sequence [2], whereas state-of-the art multi-compartment models of white matter relying on PGSE or OGSE sequences usually assume a Gaussian diffusion in the extra-axonal space [3–6]. A theoretical explanation was given in Burcaw et al. [7], where the observed frequency-dependence is interpreted as resulting from the extra-axonal 2D short-range disorder of axonal packings in the plane transverse to white matter fibers.

The extra-axonal perpendicular diffusivity transverse to axons was thus written as

$$D(t) = D\_{\infty} + A. \frac{\ln \left( t/t\_c \right)}{t}, \quad t\_c \ll t \tag{1}$$

and, equivalently, in the frequency domain:

$$\operatorname{Re}\mathcal{D}(\omega) \sim D\_{\infty} + A. \frac{\pi}{2}|\omega|, \quad |\omega| t\_{\mathcal{C}} \ll 1 \tag{2}$$

where t is the diffusion time, ω is the dual frequency and tc represents the time to diffuse across the correlation length lc of the packing geometry (lc closely follows the mean external radius rext of the axon packing [7, 8]).

Similar to recent models accounting for extra-axonal time-dependence in the case of single diffusion encoding sequences [9], a multi-compartment model for cosine OGSE sequences was proposed in Ginsburger et al. [10] which added a frequency-dependent term in the extra-axonal diffusion tensor perpendicular diffusivity based on Equation (2), showing a significant improvement of the model fit quality.

The first contribution of this article is to show the relevance of such a frequency-dependent correction using clinically feasible cosine trapezoidal OGSE sequences. The theoretically predicted linear relationship between the extracellular perpendicular diffusivity and the OGSE frequency was observed using 3D simulations of the diffusion process with different values of signal-to-noise ratios (SNR).

Having introduced a physically plausible frequencydependent correction in our model, the next step is to study the dependence of its scaling coefficient A (Equation 2) to the geometrical features of the extracellular space. Indeed, to our knowledge, this dependence is very little known. An empirical law A ∼ l 2 c relating A to the correlation length lc was given in Burcaw et al. [7] and Fieremans et al. [8] but is not sufficient to catch the complexity of the scaling coefficient A. The value of A is a measure of the strength of the structural disorder [7, 11], thus related to the geometrical properties characterizing the spatial organization of white matter at various scales. A possible approach to decipher the complex relationship between A and white matter features is to perform Monte-Carlo simulations of the diffusion process in diffusing media with increasing level of structural disorder.

The simulation of the diffusion process in state-of-the-art Monte-Carlo simulators such as CAMINO [12] and DMS [13] is decomposed into three main steps: (1) the generation of a 3D numerical phantom representing the diffusion medium (2) the Monte-Carlo simulation of the Brownian motion of spins (3) the synthesis of a DW-NMR signal. Simulators like CAMINO or DMS are able to extract triangle meshes from histological samples in step (1) in order to simulate diffusion in ultrarealistic media. However, this approach does not allow to have access to the plethora of possible geometries but only to the limited set of configurations provided by the reduced set of histological samples. Beside the possibility to use geometries extracted from histological samples in step (1), state-of-the-art simulators can only generate a limited number of geometries which might not represent white matter sufficiently well. For instance, the CAMINO and DMS simulators are able to simulate the diffusion process in any triangle mesh, but the algorithm used to construct simulation meshes from input geometrical parameters only generates substrates with straight cylinders of various diameters (including crossing between two populations of fibers). Other simulation tools like Fiberfox [14] rely on analytical models associated to each particular cell geometry, including various combinations of sticks, tensors, zeppelins, balls, dots, and astrosticks. They are inherently limited by the realism of the used geometries and the employed analytical models do not account for the presence of structural disorder in the extracellular space. There is thus a real need to propose alternative generative algorithms able to create more complex geometries while controlling the parameters driving the various sources of geometrical disorder to explore more extensively the vast domain of possible geometries.

The main contribution of this article is therefore the development of a novel algorithm to produce a wide variety of biomimicking numerical phantoms representing more realistic white matter tissue configurations from a reduced set of control parameters. Embedded in the Diffusion Microscopist Simulator (DMS) [13], this algorithm enables to control the degree of complexity of the generated geometrical configurations (induction of global angular dispersion and local tortuosity, presence of Ranvier nodes along the axonal membrane, presence of beading) with few design parameters, and without the necessity of any input histological sample. Such numerical white matter phantoms are then used to perform Monte-Carlo simulations of the diffusion process from which simulated diffusion-weighted NMR signals can be synthesized using trapezoidal OGSE sequences at different values of SNR. The obtained signal is fed into our analytical model to explore the evolution of the structural disorder coefficient A for various well-characterized geometrical configurations of the extracellular space.

# 2. MATERIALS AND METHODS

# 2.1. White Matter Numerical Phantoms

Actual simulation tools rely either on membrane surfaces extracted from manual or automatic segmentations of microscopic imaging data stemming from histological tissue samples or on algorithms generating designs of axonal membranes built from a distribution of rectilinear cylinders (possibly with angular dispersion) which diameters d follow a Gamma distribution h. The function h can be tuned using a α shape parameter and a β inverse scale parameter such that:

$$h(d; \alpha, \beta) = \frac{\beta^{\alpha} d^{\alpha - 1} e^{-\beta d}}{\Gamma(\alpha)} \tag{3}$$

where Ŵ is a complete Gamma function. The first approach is hard to achieve and generally yields a limited collection of membrane configurations, whereas the second approach provides over-simplistic representations of white matter. Here, we propose a novel framework to design more realistic membrane geometries better mimicking white matter structure. One of the challenging issue of such an approach is to design a tool enabling to cover a wide range of actual geometries while using a limited number of parameters. Several observations can be taken into account to improve the development of more realistic white matter geometries:

$$\frac{d}{D} = k \text{g} \sqrt{\log\left(\frac{1}{\mathcal{g}}\right)}\tag{4}$$

where k is a constant, D is the external diameter of the axon (including the myelin sheath) and g is the g-ratio defined as the ratio between the axonal membrane and the external myelin sheath outer membrane diameters.

• It is not clear whether the axon diameter and myelin sheath thickness remain constant along the axon; several studies have assumed this absence of variation [18, 19] whereas there is no clear assessment of such a property; in particular, it is known that membrane injury can induce axonal beading for instance due to cytoskeletal damage. According to Budde and Frank [20], beading-induced changes in cellmembrane morphology are sufficient to significantly hinder water mobility and thereby decrease the apparent diffusion coefficient; it is therefore recommended to account for this and allow axon diameter variation.

Accounting for all these observations, we propose an algorithm relying on a six-fold strategy to design white matter mimicking numerical phantoms, that do not present the actual limitations of existing phantom design tools. Such phantoms will allow to go further into the study of the impact of both intra- and extra-axonal compartments on the diffusion signal. In particular, current achievable numerical phantoms do not provide the possibility to induce local fiber tortuosity, which might have a significant impact on the signal stemming from the extraaxonal compartment by modifying the parallel and transverse diffusivities along and perpendicular to the fibers.

The phantom generation algorithm takes a maximal number of 14NPopulations + 2 parameters to generate complex axonal geometries, where NPopulations is the number of fiber populations (NPopulations > 1 in the case of crossing fibers). The list of control parameters used to design a fiber population is summarized in **Table 1**.