††thanks: Co-last authors. These senior authors contributed equally to this
work.††thanks: Co-last authors. These senior authors contributed equally to
this work.

# Magic DIAMOND: Multi-Fascicle Diffusion Compartment Imaging with Tensor
Distribution Modeling and Tensor-Valued Diffusion Encoding

A. Reymbaut1,2<EMAIL_ADDRESS>A. Valcourt Caron1, G. Gilbert2,
F. Szczepankiewicz3,4, M. Nilsson3, S. K. Warfield5 M. Descoteaux1 B.
Scherrer5 1Université de Sherbrooke, Sherbrooke, QC J1K 2R1, Canada
2MR Clinical Science, Philips Healthcare Canada, Markham, ON L6C 2S3, Canada
3Department of Clinical Sciences, Lund University, 22184, Lund, Sweden
4Random Walk Imaging AB, 22224, Lund, Sweden
5Department of Radiology, Boston Children’s Hospital, Boston, MA 02115, United
States

(August 28, 2024)

###### Abstract

Diffusion tensor imaging provides increased sensitivity to microstructural
tissue changes compared to conventional anatomical imaging but also presents
limited specificity. To tackle this problem, the DIAMOND model subdivides the
voxel content into diffusion compartments and draws from diffusion-weighted
data to estimate compartmental non-central matrix-variate Gamma distribution
of diffusion tensors, thereby resolving crossing fascicles while accounting
for their respective heterogeneity. Alternatively, tensor-valued diffusion
encoding defines new acquisition schemes tagging specific features of the
intra-voxel diffusion tensor distribution directly from the outcome of the
measurement. However, the impact of such schemes on estimating brain
microstructural features has only been studied in a handful of parametric
single-fascicle models. In this work, we derive a general Laplace transform
for the non-central matrix-variate Gamma distribution, which enables the
extension of DIAMOND to tensor-valued encoded data. We then evaluate this
"Magic DIAMOND" model in silico and in vivo on various combinations of tensor-
valued encoded data. Assessing uncertainty on parameter estimation via
stratified bootstrap, we investigate both voxel-based and fixel-based metrics
by carrying out multi-peak tractography. We show that our estimated metrics
can be mapped along tracks robustly across regions of fiber crossing, which
opens new perspectives for tractometry and microstructure mapping along
specific white-matter tracts.

Abbreviations used: MRI, magnetic resonance imaging; DWI, diffusion-weighted
MRI; DW, diffusion-weighted; DTI, diffusion tensor imaging; DIAMOND,
distribution of anisotropic microstructural environments in diffusion
compartment imaging; CHARMED, composite hindered and restricted model of
diffusion; NODDI, neurite orientation dispersion and density imaging; ADC,
apparent diffusion coefficient; FW, free-water; fMD, fascicle mean
diffusivity; fAD, fascicle axial diffusivity; fRD, fascicle radial
diffusivity; fFA, fascicle fractional anisotropy; PGSE, pulsed gradient spin
echo; AIC, Akaike information criterion; SNR, signal-to-noise ratio; IQR,
interquartile range; LL, linear-linear diffusion encoding combination; LP,
linear-planar diffusion encoding combination; LS, linear-spherical diffusion
encoding combination.

##  I Introduction

Measuring water diffusion with diffusion-weighted MRI (DWI) enables the non-
invasive characterization of biological tissues in vivo. In particular,
diffusion tensor imaging (DTI) (Basser:1994, ), a very common DWI technique,
has proven sensitive to microstructural tissue changes but only provides a
voxel-scale average of the intra-voxel diffusion profile. Indeed, DTI
considers the DW signal $\mathcal{S}(b,\mathbf{n})$ acquired during a typical
Stejskal-Tanner sequence (Stejskal_Tanner:1965, ) with b-value $b$ and unit
orientation $\mathbf{n}$ for the diffusion-probing magnetic field gradient,
and interprets it as the following monoexponential:

$\frac{\mathcal{S}(b,\mathbf{n})}{\mathcal{S}_{0}}=\exp(-b\,\mathbf{n}^{\mathrm{T}}\cdot\langle\mathbf{D}\rangle\cdot\mathbf{n})=\exp(-bD_{\mathbf{n}})\,,$
(1)

where $\mathcal{S}_{0}$ is the non-DW signal, $\langle\mathbf{D}\rangle$ is
the voxel-scale averaged diffusion tensor, $D_{\mathbf{n}}$ is the effective
diffusivity along orientation $\mathbf{n}$, "$\cdot$" denotes the
vector/tensor multiplication, and the superscript "$\mathrm{T}$" indicates
vector/tensor transposition. As a result, DTI shows poor specificity in
depicting the precise nature of microstructural tissue changes
(Alexander:2001, ; Alexander:2007, ; Jones_Cercignani:2010, ; Jones:2012, ),
especially in voxels of crossing fascicles (Pierpaoli:1996, ; Douaud:2011, ),
which make up 60 to 90% of the brain (Jeurissen:2013, ). Fig. 1 illustrates a
few archetypal voxels where DTI’s lack of specificity may arise.

Figure 1: Illustration of archetypal voxel contents in terms of colored
microscopic diffusion tensors. While the voxel contents typically associated
to a white-matter (WM) fiber (a) and to WM crossing fibers (b) would for
instance correspond to the healthy corpus callosum and its crossing with the
cingulum, respectively, the contents typically associated to demyelination (c)
and WM inflammation (d) usually correspond to unhealthy white matter. DTI
would typically yield comparable voxel-averaged diffusion tensors for the
three voxels on the right.

DTI’s lack of specificity stems from the inherent heterogeneity of typical DWI
voxels, whose volume of a few cubic millimeters encompasses multiple cell
types, sizes, geometries and orientations, and the extra-cellular space
(Norris:2001, ; Sehy:2002, ; Minati:2007, ; Mulkern:2009, ). In other words,
the measured DW signal is a combination of signals arising from a variety of
microstructural environments. Approximating that diffusion is a Gaussian
process within these environments and that no exchange occurs between them
(Johnson:1993, ), which usually holds in healthy tissues for typical
acquisition times (NilssonThesis:2011, ; Nilsson:2013a, ; Nilsson:2013b, ;
Lampinen:2017, ), intra-voxel heterogeneity can be accounted for by
considering a weighted sum of all microstructural signals, thus obtaining the
following DW signal:

$\frac{\mathcal{S}(b,\mathbf{n})}{\mathcal{S}_{0}}=\int\mathcal{P}(\mathbf{D})\,\exp(-b\,\mathbf{n}^{\mathrm{T}}\cdot\mathbf{D}\cdot\mathbf{n})\,\mathrm{d}\mathbf{D}\,,$
(2)

where the intra-voxel diffusion tensor distribution $\mathcal{P}(\mathbf{D})$
weighs the different microstructural contributions (Jian:2007, ) and the
integral covers the space $\mathrm{Sym}^{+}(3)$ of 3$\times$3 symmetric
positive-definite tensors. To better capture intra-voxel heterogeneity and
thus overcome DTI’s limitations, one can either attempt to improve the
interpretation of the diffusion signal post-acquisition, via an appropriate
signal representation or model, or to enhance the specificity of the set of
acquired signals itself.

In this work, we combine two techniques that respectively target these two
solutions: tensor-valued diffusion encoding (Cory:1990, ; Mori:1995, ;
Wong:1995, ; Eriksson:2013, ; Westin:2014, ; Eriksson:2015, ; Westin:2016, ;
Topgaard:2017, ) and the Distribution of anisotropic microstructural
environments in diffusion compartment imaging (DIAMOND) model
(Scherrer_DIAMOND:2016, ; Scherrer_aDIAMOND:2017, ), that has shown great
promise in mapping the Meyers loop (Chamberland:2017, ), imaging the preterm
human cortex (Eaton-Rosen:2017, ), characterizing mild traumatic brain
injuries (Scherrer_ISMRM:2017, ) and studying infantile autism (Scherrer:2019,
). After introducing these techniques in sections II.1 and II.2, we extend
DIAMOND to tensor-valued encoded data in section II.3 and obtain the "Magic
DIAMOND" model, named in homage to the magic-angle spinning DW acquisitions of
Ref. Eriksson:2013,  and its physical chemistry inspiration (Andrew:1959, ).
This new approach, first to draw from tensor-valued diffusion encoding to
estimate compartmental distributions of diffusion tensors, is designed to
tackle various problems such as edema, neuroinflammation or axonal degradation
in voxels of crossing fascicles, all relevant to the majority of the brain and
to the study and understanding of neurodegenerative diseases. In section III,
we describe the methods used to evaluate Magic DIAMOND on in silico signals
and in vivo human data. In particular, we address how to select the
appropriate number of intra-voxel fascicles, estimate parameters, acquire the
in silico and in vivo datasets, and compare our model with the original
DIAMOND model on these datasets using stratified bootstrap (in silico and in
vivo), cost function computation (in vivo) and tractography (in vivo). Results
reporting on various combinations of diffusion encodings are presented in
section IV and discussed in section V. While the in silico signals enable
preliminary quantification of the accuracy and precision of Magic DIAMOND’s
estimations, the in vivo data offers a proof of principle of Magic DIAMOND’s
potential.

##  II Theory

###  II.1 The DIAMOND model

Mathematical models enable a direct translation from the diffusion signal’s
features to metrics that should describe the biological properties of the
voxel content, as discussed in (Yablonskiy:2010, ; Jelescu:2017, ;
Novikov:2018, ). Among them, diffusion compartment imaging (DCI) reflects the
presence of intra-voxel tissue compartments and relates compartmental features
to microstructural tissue properties at sub-voxel resolution, effectively
providing insight into microstructural tissue changes with critically improved
sensitivity and specificity. A brief review of major DCI techniques includes
the Composite hindered and restricted model of diffusion (CHARMED)
(Assaf_CHARMED:2004, ; Assaf_CHARMED:2005, ), its extension AxCaliber
(Assaf:2008, ), and the Neurite orientation dispersion and density imaging
(NODDI) model (Zhang_NODDI:2012, ). However, these methods consider key
assumptions that are inconsistent with the known tissue microstructure
(Alexander:2010, ; Aboitiz:1992, ; Lampinen_CODIVIDE:2017, ; Nilsson:2017, ).
As a solution, fitting either compartment-specific distributions of apparent
diffusion coefficients (ADCs) (Yablonskiy:2003, ) or one-dimensional Gamma
distributions of ADCs in the ball-and-stick model (Jbabdi:2012, ) was
proposed. However, while the former could not characterize the 3D anisotropy
of diffusion observed in the brain (Moseley:1990, ), the latter only
implemented a Gamma distribution with same shape and scale parameters for all
compartments. Finally, none of the previous models estimates crossing
anisotropic compartments.

A recent DCI model has been proposed to circumvent these limitations: the
Distribution of anisotropic microstructural environments in diffusion
compartment imaging (DIAMOND) model (Scherrer_DIAMOND:2016, ;
Scherrer_aDIAMOND:2017, ). Illustrated in Fig. 2, DIAMOND considers a free-
water (FW) compartment and a series of anisotropic compartments per voxel. In
the brain, these anisotropic compartments are usually referred to as
"fascicles," as anisotropic diffusion is conventionally recognized to arise
from the presence of white matter fascicles at the mesoscale. The signal
arising from each fascicle is described by a distinct unimodal non-central
matrix-variate Gamma distribution
$\mathcal{P}_{\Gamma}(\mathbf{D},\kappa,\mathbf{\Psi},\mathbf{\Theta})$ of
3$\times$3 symmetric positive-definite diffusion tensors
(Gupta_Nagar_Book:2000, ):

$\displaystyle\mathcal{P}_{\Gamma}(\mathbf{D},\kappa,\mathbf{\Psi},\mathbf{\Theta})=$
$\displaystyle\;\frac{\mathrm{Det}(\mathbf{D})^{\kappa-2}}{\mathrm{Det}(\mathbf{\Psi})^{\kappa}\,\Gamma_{3}(\kappa)}\,\exp(-\mathrm{Tr}(\mathbf{\Theta}+\mathbf{\Psi}^{-1}\cdot\mathbf{D}))\times\mathcal{F}_{0,1}(\kappa,\mathbf{\Theta}\cdot\mathbf{\Psi}^{-1}\cdot\mathbf{D})\;,$
(3)

where $\Gamma_{3}(\kappa)=\pi^{3/2}\prod_{m=1}^{3}\Gamma(\kappa-(m-1)/2)$ is
the multivariate Gamma function, $\mathcal{F}_{0,1}$ is the hypergeometric
(Bessel) function of matrix argument of order $(0,1)$, $\kappa>1$ is the shape
parameter, $\mathbf{\Psi}\in\mathrm{Sym}^{+}(3)$ is the scale tensor and
$\mathbf{\Theta}\in\mathrm{Sym}(3)$ is the noncentrality parameter
($\mathrm{Sym}(3)$ denotes the space of 3$\times$3 symmetric tensors).
Consequently, the overall modeled signal writes

$\frac{\tilde{\mathcal{S}}(b,\mathbf{n})}{\mathcal{S}_{0}}=f_{\mathrm{FW}}\,\exp(-bD_{\mathrm{FW}})+\sum_{j}f_{j}\,\frac{\tilde{\mathcal{S}}_{\mathrm{f},j}(b,\mathbf{n})}{\mathcal{S}_{0}}\,,$
(4)

with

$\frac{\tilde{\mathcal{S}}_{\mathrm{f},j}(b,\mathbf{n})}{\mathcal{S}_{0}}=\int\mathcal{P}_{\Gamma}(\mathbf{D},\kappa_{j},\mathbf{\Psi}_{j},\mathbf{\Theta}_{j})\,\exp(-b\,\mathbf{n}^{\mathrm{T}}\cdot\mathbf{D}\cdot\mathbf{n})\,\mathrm{d}\mathbf{D}\,,$
(5)

where $f_{\mathrm{FW}},f_{j}\in[0,1]$ are compartmental signal fractions,
normalized so that $f_{\mathrm{FW}}+\sum_{j}f_{j}=1$, and $D_{\mathrm{FW}}=3$
µ$\mathrm{m}^{2}/\mathrm{ms}$.

The mean diffusion tensor of the distribution is given as a function of the
distribution parameters by

$\langle\mathbf{D}\rangle=\mathbf{\Psi}\cdot(\kappa\mathbf{I}_{3}+\mathbf{\Theta})\,.$
(6)

The "width" of the peak-shaped tensor distribution is interpreted in terms of
"tissue heterogeneity", since a very heterogeneous compartment should be
described using a variety of distinct diffusion tensors. While the shape
parameter $\kappa$ relates to the isotropic component of tissue heterogeneity,
the noncentrality parameter $\mathbf{\Theta}$ allows to squeeze the
distribution in order to capture anisotropic tissue heterogeneity. Indeed,
this effect of $\mathbf{\Theta}$ on any fit of Eq. (4) is ensured by the
choice of parametrization made for $\mathbf{\Theta}$ in Ref.
Scherrer_aDIAMOND:2017, : if one expresses $\langle\mathbf{D}\rangle$ and
$\mathbf{\Theta}$ using the transfer matrix $\mathbf{V}$ from the laboratory
frame of reference to the compartment’s eigenbasis and to insert this result
into Eq. (6), one then obtains

$\langle\mathbf{D}\rangle=\mathbf{V}\cdot\mathrm{Diag}(\lambda^{\perp},\lambda^{\perp},\lambda^{\parallel})\cdot\mathbf{V}^{\mathrm{T}}\,,$
(7)

where $\lambda^{\parallel}\geq\lambda^{\perp}$, and

$\mathbf{\Theta}=\mathbf{V}\cdot\mathrm{Diag}(0,0,\kappa^{\prime})\cdot\mathbf{V}^{\mathrm{T}}\,.$
(8)

Note that the FW signal in Eq. (4) corresponds to a non-central matrix-variate
Gamma distribution with $\kappa\to+\infty$, $\kappa^{\prime}\to+\infty$ and
isotropic $\langle\mathbf{D}\rangle$.

DIAMOND thus effectively allows the evaluation of compartment-specific
diffusion features, such as fascicle mean diffusivity (fMD), fascicle axial
diffusivity (fAD), fascicle radial diffusivity (fRD) and fascicle fractional
anisotropy (fFA), while also accounting for intra-compartment heterogeneity,
especially in voxels of crossing fascicles. However, extracting compartment-
specific features from conventional DWI data is known to be a difficult
inversion problem, involving functionals with numerous local minima. Besides,
even if intra-compartment heterogeneity is accounted for by the use of a
compartmental diffusion tensor distribution - a mathematical achievement in
itself - properly quantifying this heterogeneity remains a challenge to this
day. Indeed, analytically tractable definitions for such heterogeneity metrics
have, to our knowledge, not yet been derived from parametric tensor
distributions.

Figure 2: Illustration of the way a complex voxel content (crossing fibers and
free water, left panel) is depicted withing the DIAMOND model (right panel).
Each anisotropic diffusion compartment is described in terms of a non-central
matrix-variate Gamma distribution
$\mathcal{P}_{\Gamma}(\mathbf{D},\kappa,\mathbf{\Psi},\mathbf{\Theta})$ Eq.
(3) whose mean diffusion tensor is associated to various fascicle metrics,
such as the fascicle mean diffusivity (fMD), the fascicle axial diffusivity
(fAD), the fascicle radial diffusivity (fRD) and the fascicle fractional
anisotropy (fFA). Each overall distribution is also given a total signal
fraction $f$.

###  II.2 Tensor-valued diffusion encoding

An alternative way to enhance the specificity of the information provided by
diffusion MRI is to obtain complementary pieces of diffusion information at
the acquisition stage. The most commonly used diffusion sequence, the pulsed
gradient spin echo (PGSE) sequence (Stejskal_Tanner:1965, ), corresponds to
applying a linear (unidirectional) diffusion gradient
$\mathbf{g}(t)=g(t)\,\mathbf{n}$. Such a sequence is characterized by a spin-
dephasing vector (or q-vector)
$\mathbf{q}(t)=\gamma\int_{0}^{t}\mathbf{g}(t^{\prime})\,\mathrm{d}t^{\prime}=q(t)\,\mathbf{n}$
that grows and vanishes through time along a unique direction given by the
unit vector $\mathbf{n}$. The time-dependent norm of this q-vector yields the
b-value $b=\int_{0}^{\tau}q^{2}(t)\,\mathrm{d}t$, where $\tau$ is the total
encoding time. Instead, tensor-valued diffusion encoding (or b-tensor
encoding) (Westin:2014, ; Eriksson:2015, ; Westin:2016, ; Topgaard:2017, ) has
introduced q-trajectories in which the q-vector
$\mathbf{q}(t)=\gamma\int_{0}^{t}\mathbf{g}(t^{\prime})\,\mathrm{d}t^{\prime}=q(t)\,\mathbf{n}(t)$
follows a continuous and non-trivial trajectory through time (Eriksson:2013, ;
Westin:2016, ). The name "b-tensor encoding" comes from the fact that the
time-dependency of the q-vector’s orientation $\mathbf{n}(t)$ does not encode
diffusion weighting in a vector, but rather in a 3$\times$3 symmetric
positive-semidefinite tensor,

$\mathbf{b}=\int_{0}^{\tau}q^{2}(t)\,\mathbf{n}(t)\cdot\mathbf{n}^{\mathrm{T}}(t)\,\mathrm{d}t\in\mathrm{Sym}^{(+)}(3)\,,$
(9)

where $\mathrm{Sym}^{(+)}(3)$ denotes the space of symmetric positive-
semidefinite tensors. Note that this b-tensor was introduced in earlier works
as a b-matrix accounting for cross terms between diffusion gradient pulses
(Basser:1994, ; Mattiello:1997, ).

While the shape of the b-tensor’s glyph gives its name to the corresponding
diffusion encoding (linear or spherical encoding for instance), its trace
defines the associated b-value $b=\mathrm{Tr}(\mathbf{b})$. In this general
b-tensor formalism, the specific tagging of a given diffusion pattern is
achieved through the expression of the overall DW signal (Westin:2014, ;
Westin:2016, ):

$\frac{\mathcal{S}(\mathbf{b})}{\mathcal{S}_{0}}=\int\mathcal{P}(\mathbf{D})\,\exp(-\mathbf{b}:\mathbf{D})\,\mathrm{d}\mathbf{D}\;,$
(10)

where

$\mathbf{b}:\mathbf{D}=\sum_{ij}b_{ij}D_{ij}=\int_{0}^{\tau}q^{2}(t)\;\mathbf{n}^{\mathrm{T}}(t)\cdot\mathbf{D}\cdot\mathbf{n}(t)\,\mathrm{d}t$
(11)

is the Frobenius inner product, a generalized scalar product that "projects"
the diffusion tensor onto the b-tensor, thereby effectively selecting a
certain diffusion pattern to probe. For instance, the signal Eq. (2),
associated to the PGSE q-vector $\mathbf{q}(t)=q(t)\,\mathbf{n}$, actually
maps back to the linear (stick-like) b-tensor
$\mathbf{b}=b\;\mathbf{n}\cdot\mathbf{n}^{\mathrm{T}}$ that was already
introduced in (Jian:2007, ). Also, the spherical encoding yielded by a
spherical b-tensor allows for the specific encoding of isotropic diffusion
during the acquisition process (Mori:1995, ; Wong:1995, ; Eriksson:2013, ).
Combinations of these two encodings, coupled to the fitting of a Gamma
distribution of ADCs on powder-averaged data, have already enabled the
extraction of novel diffusion metrics (Lasic:2014, ; Eriksson:2015, ) that may
be used to measure the microscopic diffusion anisotropy of healthy brain
tissues (Szczepankiewicz:2015, ), as well as to probe diffusional variance in
tumors (Szczepankiewicz:2016, ).

Even though the opportunity to measure non-conventional pieces of diffusion
information dates back almost 30 years, with the double diffusion encoding of
Ref. Cory:1990,  and the triple diffusion encoding of Refs. Mori:1995, ;
Wong:1995, , it is only within the last decade that non-trivial diffusion
encoding has drawn much attention (Ozarslan:2009, ; Shemesh:2010, ;
Lawrenz:2010, ; Valette:2012, ; Jespersen:2012, ; Jespersen:2013, ;
Lawrenz:2016, ; Ianus:2016, ; Lawrenz:2018, ; Ianus:2018, ;
Jensen_Helpern:2018, ; Szczepankiewicz_Jeurissen_ISMRM:2018, ;
Herberthson:2019, ). From a modeling standpoint, the main interest lies in
assessing whether the estimation of microstructural features benefits from
these non-trivial diffusion encodings or not. The impact of planar and
spherical encodings has already been investigated in ball-and-stick-like
models such as the white matter Standard Model (Novikov_WMSM:2018, ) and NODDI
(Zhang_NODDI:2012, ), focusing mainly on accuracy (Lampinen_CODIVIDE:2017, ),
precision (Fieremans_ISMRM:2018, ) and degeneracy (Reisert:2019, ;
Coelho:2019, ; Coelho_arxiv:2019, ) in parameter estimation. Also, a
generalized cumulant approach implicitly assuming a normal distribution of
diffusion tensors has already been implemented in the framework of b-tensor
encoding (Westin:2016, ). However, these models do not separately characterize
crossing fascicles. One exception lies in the Monte-Carlo inversion of Refs.
deAlmeidaMartins:2016, ; deAlmeidaMartins:2018, ; Topgaard:2019,  that shows
impressive results in teasing apart intra-voxel populations in the diffusion-
relaxation space but can unfortunately be significantly time-consuming and
noise-sensitive.

###  II.3 The Magic DIAMOND model

Extending the DIAMOND model to tensor-valued encoded diffusion data, a
procedure illustrated in Fig. 3, consists in merging Eqs. (4) and (10) to
yield the following Magic DIAMOND signal:

$\displaystyle\frac{\tilde{\mathcal{S}}(\mathbf{b})}{\mathcal{S}_{0}}=$
$\displaystyle\;f_{\mathrm{FW}}\,\exp(-bD_{\mathrm{FW}})+\sum_{j}f_{j}\underbrace{\int\mathcal{P}_{\Gamma}(\mathbf{D},\kappa_{j},\mathbf{\Psi}_{j},\mathbf{\Theta}_{j})\,\exp(-\mathbf{b}:\mathbf{D})\,\mathrm{d}\mathbf{D}}_{\tilde{\mathcal{S}}_{\mathrm{f},j}(\mathbf{b})/\mathcal{S}_{0}}\,,$
(12)

whose fascicle part
$\tilde{\mathcal{S}}_{\mathrm{f},j}(\mathbf{b})/\mathcal{S}_{0}$ remains to be
analytically computed. This fascicle part is the general Laplace transform of
the non-central matrix-variate Gamma distribution Eq. (3).

Figure 3: Diffusion compartement imaging (DCI) techniques subdivide the voxel
content into diffusion compartments (left panel). In particular, the DIAMOND
model estimates a non-central matrix-variate Gamma distribution of diffusion
tensors $\mathcal{P}_{\Gamma}(\mathbf{D})$ for each compartment (top center
panel). Alternatively, tensor-valued encoding (bottom center panel)
establishes novel sequences that yield for instance spherical encoding, which
specifically tags isotropic diffusion patterns ($D_{\mathrm{iso}}$, in purple)
of any diffusion tensor (grey glyph), contrary to the conventional linear
encoding that only probes diffusion along single directions $\mathbf{n}$
($D_{\mathbf{n}}$, in orange). The spherically encoded sequence used
throughout this paper (gradient waveform $\mathbf{g}(t)$, spin-dephasing
vector $\mathbf{q}(t)$ and q-trajectory) is presented as an example in the
right panel, where color codes for orientation. In this work, we incorporate
tensor-valued encoding in the DIAMOND fitting procedure at two key steps
(dashed arrow lines): 1) during the acquisition process, 2) in the estimation
of $\mathcal{P}_{\Gamma}(\mathbf{D})$ for each compartment.

####  II.3.1 General Laplace transform of the non-central matrix-variate Gamma
distribution

Let us define the general Laplace transform of an arbitrary diffusion tensor
distribution $\mathcal{P}(\mathbf{D})$ from Eq. (10) as

$\mathcal{T}_{\mathcal{P}}(\mathbf{b})=\int\mathcal{P}(\mathbf{D})\,\exp(-\mathbf{b}:\mathbf{D})\,\mathrm{d}\mathbf{D}\,,$
(13)

where $\mathcal{T}$ denotes the Laplace transformation. To compute this
Laplace transform, let us introduce the moment-generating function for any
distribution $\mathcal{P}$ of symmetric tensors (Gupta_Nagar_Book:2000, ):

$M_{\mathcal{P}}(\mathbf{Z})=\int\\!\mathcal{P}(\mathbf{D})\,\exp\\!\left(\mathrm{Tr}[\mathbf{Z}\cdot\mathbf{D}]\right)\mathrm{d}\mathbf{D}\,,$
(14)

where $\mathbf{Z}\in\mathrm{Sym}(3)$. By choosing the symmetric tensor

$\mathbf{Z}\equiv-\mathbf{b}=-\int_{0}^{\tau}\\!q^{2}(t)\,\mathbf{n}(t)\cdot\mathbf{n}^{\mathrm{T}}(t)\,\mathrm{d}t$
(15)

from Eq. (9) and by using the trace property
$\mathrm{Tr}(\mathbf{u}\cdot\mathbf{v}^{\mathrm{T}})=\mathbf{v}^{\mathrm{T}}\cdot\mathbf{u}$
and the definition of the Frobenius inner product Eq. (11), we prove that the
Laplace transform Eq. (13) equals

$\mathcal{T}_{\mathcal{P}}(\mathbf{b})=M_{\mathcal{P}}(-\mathbf{b})\,.$ (16)

In the particular case of the non-central matrix-variate Gamma distribution
$\mathcal{P}_{\Gamma}$ Eq. (3), it can be shown (Gupta_Nagar_Book:2000, ) that
its moment-generating function is given by

$M_{\mathcal{P}_{\Gamma}}(\mathbf{Z})=[\mathrm{Det}(\mathbf{I}_{3}-\mathbf{Z}\cdot\mathbf{\Psi})]^{-\kappa}\;\mathrm{exp}\\!\left[\mathrm{Tr}([(\mathbf{I}_{3}-\mathbf{Z}\cdot\mathbf{\Psi})^{-1}-\mathbf{I}_{3}]\cdot\mathbf{\Theta})\right],$
(17)

with $\mathbf{Z}$ such that
$(\mathbf{I}_{3}-\mathbf{Z}\cdot\mathbf{\Psi})\in\mathrm{Sym}^{+}(3)$. Since
the choice $\mathbf{Z}\equiv-\mathbf{b}$ Eq. (15) satisfies this last
condition, the general Laplace transform of the non-central matrix-variate
Gamma distribution writes

$\displaystyle\mathcal{T}_{\mathcal{P}_{\Gamma}}(\mathbf{b})$
$\displaystyle=\int\mathcal{P}_{\Gamma}(\mathbf{D},\kappa,\mathbf{\Psi},\mathbf{\Theta})\,\exp(-\mathbf{b}:\mathbf{D})\,\mathrm{d}\mathbf{D}$
$\displaystyle=M_{\mathcal{P}_{\Gamma}}(-\mathbf{b})$
$\displaystyle=[\mathrm{Det}(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})]^{-\kappa}\;\mathrm{exp}\\!\left[\mathrm{Tr}([(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})^{-1}-\mathbf{I}_{3}]\cdot\mathbf{\Theta})\right].$
(18)

This expression for $\mathcal{T}_{\mathcal{P}_{\Gamma}}(\mathbf{b})$ can be
rewritten by considering the Woodbury matrix identity to simplify the previous
trace: given a square invertible $n\times n$ matrix $\mathbf{A}$, an $n\times
m$ matrix $\mathbf{U}$ and an $m\times n$ matrix $\mathbf{V}$, assuming that
$(\mathbf{I}_{m}+\mathbf{V}\cdot\mathbf{A}^{-1}\cdot\mathbf{U})$ is
invertible, one has
$(\mathbf{A}+\mathbf{U}\cdot\mathbf{V})^{-1}=\mathbf{A}^{-1}-\mathbf{A}^{-1}\cdot\mathbf{U}\cdot(\mathbf{I}_{m}+\mathbf{V}\cdot\mathbf{A}^{-1}\cdot\mathbf{U})^{-1}\cdot\mathbf{V}\cdot\mathbf{A}^{-1}$.
Taking $\mathbf{A}\equiv\mathbf{I}_{3}$, $\mathbf{U}\equiv\mathbf{b}$ and
$\mathbf{V}\equiv\mathbf{\Psi}$, we yield

$\mathcal{T}_{\mathcal{P}_{\Gamma}}(\mathbf{b})=[\mathrm{Det}(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})]^{-\kappa}\;\mathrm{exp}\\!\left[-\mathbf{b}:[(\mathbf{I}_{3}+\mathbf{\Psi}\cdot\mathbf{b})^{-1}\cdot\mathbf{\Psi}\cdot\mathbf{\Theta}]\right],$
(19)

which is valid for any arbitrary b-tensor.

####  II.3.2 Diffusion signal within the Magic DIAMOND model

Returning to Eq. (12), using Eqs. (18) and (19), and omitting the
compartmental index $j$, we obtain the Magic DIAMOND fascicle signal

$\displaystyle\frac{\tilde{\mathcal{S}}_{\mathrm{f}}(\mathbf{b})}{\mathcal{S}_{0}}$
$\displaystyle=[\mathrm{Det}(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})]^{-\kappa}\;\mathrm{exp}\\!\left[\mathrm{Tr}([(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})^{-1}-\mathbf{I}_{3}]\cdot\mathbf{\Theta})\right]$
$\displaystyle=[\mathrm{Det}(\mathbf{I}_{3}+\mathbf{b}\cdot\mathbf{\Psi})]^{-\kappa}\;\mathrm{exp}\\!\left[-\mathbf{b}:[(\mathbf{I}_{3}+\mathbf{\Psi}\cdot\mathbf{b})^{-1}\cdot\mathbf{\Psi}\cdot\mathbf{\Theta}]\right]$
(20)

for any arbitrary b-tensor. Let us render that signal expression more explicit
in the case of an axisymmetric b-tensor. Using the previous expressions Eqs.
(6) and (8) for $\mathbf{\Psi}$ and $\mathbf{\Theta}$, the Magic DIAMOND
fascicle signal writes

$\displaystyle\frac{\tilde{\mathcal{S}}_{\mathrm{f}}(b_{\mathrm{L}},b_{\mathrm{S}},\beta)}{\mathcal{S}_{0}}=\left[f_{\perp}(b_{\mathrm{S}})\,F(b_{\mathrm{L}},b_{\mathrm{S}},\beta)\right]^{-\kappa}\mathrm{exp}\\!\left(-\kappa^{\prime}\left[1+\frac{f(b_{\mathrm{L}},b_{\mathrm{S}},\beta)\,[1-f_{\perp}(b_{\mathrm{S}})]-1}{F(b_{\mathrm{L}},b_{\mathrm{S}},\beta)}\right]\right),$
(21)

with the dimensionless quantities

$\displaystyle f_{\parallel,\perp}(b_{\mathrm{S}})=$
$\displaystyle\;1+\frac{b_{\mathrm{S}}}{3}\,\frac{\lambda^{\parallel,\perp}}{\kappa^{\parallel,\perp}}\,,$
(22) $\displaystyle f(b_{\mathrm{L}},b_{\mathrm{S}},\beta)=$
$\displaystyle\;1+\frac{3b_{\mathrm{L}}}{b_{\mathrm{S}}}\,\sin^{2}\beta\,,$
(23) $\displaystyle F(b_{\mathrm{L}},b_{\mathrm{S}},\beta)=$
$\displaystyle\;f_{\parallel}(b_{\mathrm{S}})\,f_{\perp}(b_{\mathrm{S}})+b_{\mathrm{L}}\left(\frac{\lambda^{\parallel}}{\kappa^{\parallel}}\cos^{2}\beta+\frac{\lambda^{\perp}}{\kappa^{\perp}}\sin^{2}\beta+\frac{b_{\mathrm{S}}}{3}\,\frac{\lambda^{\parallel}\lambda^{\perp}}{\kappa^{\parallel}\kappa^{\perp}}\right).$
(24)

Here, $\kappa^{\parallel}=\kappa+\kappa^{\prime}$ and $\kappa^{\perp}=\kappa$
are the axial and radial shape parameters, respectively, $\beta$ is the angle
separating the axes of revolution of the diffusion and b- tensors, and
$b_{\mathrm{L}}=bb_{\Delta}$ and $b_{\mathrm{S}}=b(1-b_{\Delta})$ are two
shape parameters drawn from Ref. Westin:2002,  for the b-tensor. These
parameters both depends on the normalized anisotropy parameter
$b_{\Delta}\in[-0.5,1]$ introduced in Ref. Haeberlen:1976, , so that
$\tilde{\mathcal{S}}_{\mathrm{f}}(b_{\mathrm{L}},b_{\mathrm{S}},\beta)\equiv\tilde{\mathcal{S}}_{\mathrm{f}}(b,b_{\Delta},\beta)$.
In particular, $b_{\Delta}=1$ yields the linear signal

$\frac{\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(lin.)}}(b,\beta)}{\mathcal{S}_{0}}=\left[1+b\left(\frac{\lambda^{\parallel}}{\kappa^{\parallel}}\cos^{2}\beta+\frac{\lambda^{\perp}}{\kappa^{\perp}}\sin^{2}\beta\right)\right]^{-\kappa^{\perp}}\mathrm{exp}\\!\left[\frac{\displaystyle-b\,(\kappa^{\parallel}-\kappa^{\perp})\,\lambda^{\parallel}\cos^{2}\beta}{\displaystyle\kappa^{\parallel}+b\left(\lambda^{\parallel}\cos^{2}\beta+\frac{\kappa^{\parallel}}{\kappa^{\perp}}\,\lambda^{\perp}\sin^{2}\beta\right)}\right],$
(25)

$b_{\Delta}=0$ yields the spherical signal

$\frac{\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(sph.)}}(b)}{\mathcal{S}_{0}}=\left[\left(1+\frac{b}{3}\,\frac{\lambda^{\perp}}{\kappa^{\perp}}\right)^{2}\left(1+\frac{b}{3}\,\frac{\lambda^{\parallel}}{\kappa^{\parallel}}\right)\right]^{-\kappa^{\perp}}\mathrm{exp}\\!\left[\frac{-b\,(\kappa^{\parallel}-\kappa^{\perp})\,\lambda^{\parallel}}{3\kappa^{\parallel}+b\,\lambda^{\parallel}}\right],$
(26)

and $b_{\Delta}=-0.5$ yields the planar signal

$\frac{\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(plan.)}}(b,\beta)}{\mathcal{S}_{0}}=\left(1+\frac{b}{2}\,\frac{\lambda^{\perp}}{\kappa^{\perp}}\right)^{-\kappa^{\perp}}\left[1+\frac{b}{2}\left(\frac{\lambda^{\parallel}}{\kappa^{\parallel}}\sin^{2}\beta+\frac{\lambda^{\perp}}{\kappa^{\perp}}\cos^{2}\beta\right)\right]^{-\kappa^{\perp}}\mathrm{exp}\\!\left[\frac{\displaystyle-b\,(\kappa^{\parallel}-\kappa^{\perp})\,\lambda^{\parallel}\sin^{2}\beta}{\displaystyle
2\kappa^{\parallel}+b\left(\lambda^{\parallel}\sin^{2}\beta+\frac{\kappa^{\parallel}}{\kappa^{\perp}}\,\lambda^{\perp}\cos^{2}\beta\right)}\right].$
(27)

This new set of equations satisfies two key conditions. First, the linear
signal Eq. (25) directly translates back to the original DIAMOND signal of
Ref. Scherrer_aDIAMOND:2017, . Second, this set of equations obeys the
relationship linking linear, planar and spherical diffusion encodings for a
purely homogeneous fascicle ($\kappa\to+\infty$ and
$\kappa^{\prime}\to+\infty$):

$\lim\limits_{\begin{subarray}{c}\kappa^{\prime}\to+\infty\\\
\kappa\to+\infty\end{subarray}}\left[\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(lin.)}}\left(\frac{b}{2},\beta\right)\,\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(plan.)}}\left(b,\beta\right)\right]=\mathcal{S}_{0}\lim\limits_{\begin{subarray}{c}\kappa^{\prime}\to+\infty\\\
\kappa\to+\infty\end{subarray}}\left[\tilde{\mathcal{S}}_{\mathrm{f}}^{\mathrm{(sph.)}}\left(\frac{3b}{2}\right)\right].$
(28)

##  III Material and methods

###  III.1 Model selection

As shown in Eq. (12), Magic DIAMOND considers at each voxel a single free
water (FW) compartment and a series of $N_{\mathrm{f}}$ fascicles (anisotropic
compartments):

$\frac{\tilde{\mathcal{S}}(\mathbf{b})}{\mathcal{S}_{0}}=f_{\mathrm{FW}}\,\exp(-bD_{\mathrm{FW}})+\sum_{j}^{N_{\mathrm{f}}}f_{j}\,\frac{\tilde{\mathcal{S}}_{\mathrm{f},j}(\mathbf{b})}{\mathcal{S}_{0}}\,,$
(29)

with $\tilde{\mathcal{S}}_{\mathrm{f},j}(\mathbf{b})/\mathcal{S}_{0}$ given by
Eq. (21). Notice that with this choice, intra-axonal restricted diffusion and
extra-axonal hindered diffusion arising from a fascicle are both modeled using
one tensor distribution. In order to set the number of intra-voxel fascicles
($N_{\mathrm{f}}\in\\{0,1,2,3\\}$ in this work), model selection was performed
using the Akaike information criterion (AIC) (Akaike:1974, ), promoting the
model of lowest index

$\mathrm{AIC}=2n_{\mathrm{param}}-2\ln(\mathcal{L})\,,$ (30)

with $n_{\mathrm{param}}$ the number of estimated model parameters and
$\mathcal{L}$ the maximum value of the model likelihood function. This
criterion favors the goodness of fit while discouraging overfitting, but is
not entirely impervious to noise in the data (Novikov:2018, ). To ensure a
fast and robust model selection, the Akaike model selection was performed on
ball-and-stick models.

###  III.2 Parameter estimation

The parameters of any evaluated diffusion tensor distribution were estimated
within each voxel using a maximum a posteriori approach where no prior
correlation between the expectations $\langle\mathbf{D}\rangle$, fractions
$\mathbf{f}$ and shape parameters $\bm{\kappa}$ is assumed. To ensure
positive-definiteness of the diffusion tensors, we parametrized them using
$\mathbf{L}\equiv\ln(\langle\mathbf{D}\rangle)$. The optimal map of parameters
was then obtained through the Bayes maximization

$\left\\{\mathbf{L}_{\mathrm{opt}},\bm{\kappa}_{\mathrm{opt}},\mathbf{f}_{\mathrm{opt}}\right\\}=\underset{\mathbf{L},\bm{\kappa},\mathbf{f}}{\mathrm{argmax}}\left[P(\mathbf{I}|\mathbf{L},\bm{\kappa},\mathbf{f})\,P(\mathbf{f}|\mathbf{L},\bm{\kappa})\,P(\bm{\kappa}|\mathbf{L})\,P(\mathbf{L})\right],$
(31)

where $\mathbf{I}$ denotes the set of acquired images. While the probability
$P(\mathbf{I}|\mathbf{L},\bm{\kappa},\mathbf{f})$ was assumed to obey a normal
distribution, $P(\mathbf{f}|\mathbf{L},\bm{\kappa})$ and
$P(\bm{\kappa}|\mathbf{L})$ were taken as uniform densities and
$P(\mathbf{L})$ was submitted to spatial regularization (Besag:1986, ;
Scherrer_DIAMOND:2016, ), merely reflecting the fact that the brain is a
continuous diffusion medium that does not contain actual voxels.

###  III.3 In vivo acquisitions

MRI acquisitions were performed on a clinical 3T system with 45 mT/m maximum
gradient amplitude (Ingenia, Philips Healthcare, Best, the Netherlands) using
a 32-channel head coil. Imaging was performed on two healthy young male
volunteers using a prototype diffusion-weighted spin-echo EPI sequence with
numerically optimized (Sjolund:2015, ) Maxwell-compensated
(Szczepankiewicz_Maxwell:2019, ) spherical, planar and linear encoding
waveforms. We also attempted to match their respective frequency contents
(Lundell:2019, ), but did not sacrifice encoding efficiency to do so.

Acquisition parameters were: TR=6500 ms, TE=121 ms, spatial
resolution=2.5$\times$2.5$\times$2.5 mm3, 48 slices, in-plane acceleration
factor=1.9, multiband factor (SENSE)=2, FoV=240$\times$240 mm2, echo
spacing=0.82 ms, multi-shell scheme of 45 signal samples at 1$\times b=0$,
6$\times b=0.1$, 6$\times b=0.7$, 12$\times b=1.4$ and 20$\times b=2$
(ms/µ$\mathrm{m}^{2}$) following the scheme suggested in Ref.
Szczepankiewicz_DIVIDE:2019, , and encoding times of 41 ms (pre), 41 ms
(post), 16 ms (pause) for linear gradient waveforms, and 45.33 ms (pre), 41.33
ms (post), 16 ms (pause) for planar and spherical gradient waveforms. All
waveforms (including the spherical ones) were rotated via a rotation scheme
similar to the one used in Ref. Szczepankiewicz_data:2019, . The acquired
images were corrected for motion and Eddy currents by applying an in-house
topup-eddy procedure based on a common averaged $b=0$ image for the linear,
planar and spherical parts of the dataset independently. They were then
resampled to 2$\times$2$\times$2 mm3 using linear interpolation. Other topup-
eddy strategies for tensor-valued diffusion encoded data can be found in the
literature (Nilsson:2015, ). The signal-to-noise ratio (SNR) of the in vivo
dataset was estimated to around 40 in the corona radiata using the method
described in the Supplemental Material of Ref. Szczepankiewicz_DIVIDE:2019,
on our $b=0.7$ ms/µ$\mathrm{m}^{2}$ spherical data. Finally, in order to
evaluate the impact of introducing tensor-valued diffusion encoding within
Magic DIAMOND, the acquired volumes were combined so as to create the three
following encoding combinations:

  * $\bullet$

LL: 45 linear $+$ 45 linear b-tensors,

  * $\bullet$

LP: 45 linear $+$ 45 planar b-tensors,

  * $\bullet$

LS: 45 linear $+$ 45 spherical b-tensors.

###  III.4 Numerical simulations

Numerical simulations were performed using a modified version of Fiberfox
(Neher:2014, ) that allows for tensor-valued diffusion encoding. We generated
a medium composed of three partially crossing fascicles surrounded by water,
described by the free isotropic diffusivity $D_{\mathrm{FW}}=3$
µ$\mathrm{m}^{2}/\mathrm{ms}$. Two fascicles span along the $x$\- and $y$\-
axes, respectively, and the third fascicle crosses with them at a $45^{\circ}$
angle between the $y$\- and $z$\- axes. Each fascicle is composed of
$900\,000$ fibers characterized by the axial and radial diffusivities
$D_{\parallel}=1.7$ µ$\mathrm{m}^{2}/\mathrm{ms}$ and $D_{\perp}=0.4$
µ$\mathrm{m}^{2}/\mathrm{ms}$, respectively (Pierpaoli:1996, ;
Descoteaux:2006, ). The radial distribution of these fibers follows a normal
distribution and the remaining extra-fiber space is filled up with free water.
Our Fiberfox phantom is illustrated in Fig. 4.

Datasets were simulated for linear, planar and spherical diffusion encodings
at two SNRs, namely the $\mathrm{SNR}=40$ of our in vivo dataset and the ideal
infinite SNR. The Rician noise of the $\mathrm{SNR}=40$ dataset was generated
as follows:

$\mathcal{S}_{\mathrm{noisy}}(\mathbf{b})=\mathcal{S}_{0,\mathrm{true}}\sqrt{\left(\frac{\mathcal{S}_{\mathrm{true}}(\mathbf{b})}{\mathcal{S}_{0,\mathrm{true}}}+\frac{\nu}{\mathrm{SNR}}\right)^{2}+\left(\frac{\nu^{\prime}}{\mathrm{SNR}}\right)^{2}}\,,$
(32)

where $\mathcal{S}_{\mathrm{true}}$ and
$\mathcal{S}_{0,\mathrm{true}}=\mathcal{S}_{\mathrm{true}}(\mathbf{b}=\mathbf{0})$
come from the infinite-SNR (noise-free) dataset, and $\nu$ and $\nu^{\prime}$
are drawn from a normal distribution with zero mean and unit standard
deviation. We set the same acquisition parameters and considered the same
encoding combinations as those described in section III.3. To assess whether
or not two encoding combinations yield distinct estimated values for a given
metric, we used non-parametric Mann-Whitney U-tests (McKnight_Najab:2010, ),
equivalent to two-sided Wilcoxon rank sum tests.

Figure 4: Illustration of the Fiberfox numerical phantom described in section
III.4 and used to evaluate Magic DIAMOND in section IV.1. The phantom is made
out of three fascicles: one along the $x$-axis (red), one along the $y$-axis
(green), and one crossing with them at a 45∘ angle between the $y$\- and $z$\-
axes (blue). Each fascicle is composed of 900 000 fibers characterized by the
axial and radial diffusivities $D_{\parallel}=$ 1.7 µm2/ms and $D_{\perp}=$
0.4 µm2/ms, respectively. The radial distribution of these fibers follows a
normal distribution and the remaining extra-fiber space is filled up with free
water of isotropic diffusivity $D_{\mathrm{FW}}=$ 3 µm2/ms.

###  III.5 Precision on parameter estimation with stratified bootstrap

The evaluation of diffusion models with in vivo data is challenging because
the underlying ground truth is not known. We compared the various encoding
strategies by evaluating their precision on parameter estimation. Intuitively,
the precision is high when multiple acquisitions of the same object leads to
similar results and does not depend on noise realizations. Assessing
estimation uncertainty allowed us to test whether the additional information
provided by the tensor-valued diffusion encoding better constrains the
estimation problem, while keeping a constant number of DW images.

The typical way to compute estimation uncertainty is to repeat the measurement
of the same object, estimate parameters for each and assess the variance in
parameters. This would unfortunately lead to unrealistic scan time in DW-MRI.
We instead used stratified bootstrap, which amounts to 1) acquire two
repetitions of the same DW-MRI experiment, and 2) create a large number of
"virtual" acquisitions by randomly choosing, for each diffusion gradient
direction, one of the two repetitions. Unlike wild bootstrap (Whitcher:2008, )
and repetition bootstrap (Chung:2006, ), stratified bootstrap allows for the
creation of a large number of noise realizations without any explicit noise
modeling. Specifically, we acquired twice each of the datasets described in
section III.3, namely four sets of 45 linear b-tensors (L1, L2, L3, L4), two
sets of 45 planar b-tensors (P1, P2), and two sets of 45 spherical b-tensors
(S1, S2), giving the following combinations:

  * $\bullet$

LL: [L1 $\longleftrightarrow$ L2] $+$ [L3 $\longleftrightarrow$ L4],

  * $\bullet$

LP: [L1 $\longleftrightarrow$ L2] $+$ [P1 $\longleftrightarrow$ P2],

  * $\bullet$

LS: [L1 $\longleftrightarrow$ L2] $+$ [S1 $\longleftrightarrow$ S2],

where $\longleftrightarrow$ denotes random choice between two datasets on each
acquired direction and $+$ denotes dataset concatenation. We generated 100
realizations of each of these combinations. Medians and interquartile ranges
(IQRs) of estimated parameters across these realizations were then computed
within a white-matter mask to quantify their average value and variance,
respectively. We compared the results from only linear b-tensors (LL,
equivalent to the original DIAMOND formulation) to those arising from
combinations of linear, planar and spherical b-tensors (LP and LS).

###  III.6 Visualization of the optimization landscape with in vivo data

The estimation of a DCI model is known to be a challenging inversion problem
involving multidimensional functionals with numerous local minima. We
evaluated the optimization landscape by investigating the shape of the cost
function used for parameter estimation (Scherrer:2012, ) in the parameter
space. This investigation was carried out in two voxels of interest, a typical
one-fiber voxel in the corpus callosum and a typical three-fiber voxel in the
centrum semiovale, after evaluating Magic DIAMOND on the following sets of 180
b-tensors (covering all our acquired signals, see section III.5):

  * $\bullet$

LL: L1 $+$ L2 $+$ L3 $+$ L4,

  * $\bullet$

LP: L1 $+$ L2 $+$ P1 $+$ P2,

  * $\bullet$

LS: L1 $+$ L2 $+$ S1 $+$ S2,

where $+$ denotes dataset concatenation. Finally, we compared the different
cost function shapes extracted from the LL, LP and LS acquisition schemes.

###  III.7 Tractography and fascicle-metric mapping

While computing whole-brain maps is the most common approach with single-
fascicle models, the lack of one-to-one correspondence between neighbouring
voxels with multi-fascicle models makes them more challenging to evaluate. For
instance, there may be fascicles that are not present in all voxels (one-to-
zero correspondence) or fascicles that are represented by different number of
compartments in different voxels (one-to-many correspondence). We therefore
evaluated the impact of tensor-valued diffusion encoding using a fixel-based
(fascicle-element based) analysis (Raffelt:2015, ), by delineating tract
streamlines with tractography and by extracting quantities related to the most
aligned fascicle compartment along the streamlines.

Multi-peak tractography (Chamberland:2014, ) was performed on the fFA-weighted
peaks extracted from Magic-DIAMOND evaluation on the datasets defined in
section III.6. To obtain dense tractograms, we upsampled our peaks to
1$\times$1$\times$1 mm3 using nearest-neighbor interpolation, used ten seeds
per voxel and generated an ensemble tractogram (Takemura:2016, ) for each
dataset with three separate tractography runs with acceptance angles of
$45^{\circ}$, $60^{\circ}$ and $75^{\circ}$. The output tracks were finally
filtered using the recent "Recobundles" algorithm (Garyfallidis:2018, ) that
relies on shape analysis and anatomical priors to recognize and extract
bundles of interest, namely the corpus callosum, the arcuate fasciculus and
the corticospinal tract, three tracts that intersect in the centrum semiovale,
making them difficult to reconstruct in their full extent (Rheault:2019, ).
Each bundle tractogram was then used as a skeleton onto which metrics
associated to the Magic-DIAMOND fascicle most aligned with any given local
track orientation (sub-voxel tract segment) can be color-mapped. Performing
this process across all bootstrap realizations described in section III.5
enabled color-mapping of the median and interquartile range of any metric
along streamlines.

##  IV Results

###  IV.1 In silico data

Fig. 5 compares physical quantities generated within a mask corresponding to
the fascicle aligned with $x$ in our Fiberfox phantom (see Fig. 4) with their
Magic-DIAMOND estimations for the encoding combinations LL, LP and LS at
infinite SNR and at the finite $\mathrm{SNR}=40$ characterizing the in vivo
dataset described in section III.3. These physical quantities comprise the
fascicle axial diffusivity (fAD), the fascicle radial diffusivity (fRD), the
FW signal fraction $f_{\mathrm{FW}}$, and the angular deviation
$\Delta\Theta$. This last quantity is defined voxel-wise as the average of the
angles between each Magic-DIAMOND estimated fascicle and its closest Fiberfox
fiber bundle.

Figure 5: In silico evaluation of Magic DIAMOND within a mask corresponding to
the fascicle aligned with $x$ in the Fiberfox phantom illustrated in Fig. 4
for LL (red), LP (green) and LS (blue) at infinite SNR and SNR=40. The
estimated metrics are the fascicle axial diffusivity (fAD, panel (a)), the
fascicle radial diffusivity (fRD, panel (b)), the FW signal fraction
$f_{\mathrm{FW}}$ (panel (c)), and the angular deviation $\Delta\Theta$ (panel
(d)). $\Delta\Theta$ is defined voxel-wise as the average of the angles
between each Magic-DIAMOND estimated fascicle and its closest Fiberfox fiber
bundle. Results of non-parametric Mann-Whitney U-tests applied between metric
distributions are reported in order to assess whether or not two Magic-DIAMOND
estimations yield statistically distinct metrics. Such tests consider the null
hypothesis $\mathcal{H}_{0}$ that data in two Magic-DIAMOND evaluations are
sampled from identically shaped non-median-shifted continuous distributions,
against the alternative that they are not. The $p$-values resulting from these
tests inform on the acceptance or rejection of $\mathcal{H}_{0}$ at a certain
significance level: $*\equiv 0.01\leq p<0.05$, $**\equiv 0.05\leq p<0.1$,
$***\equiv 0.1\leq p<0.15$, $****\equiv 0.15\leq p<0.2$, and $\star\equiv
0.6\leq p\leq 0.7$.

###  IV.2 Precision on parameter estimation by stratified bootstrap on in vivo
data

A typical map of local orientations in a coronal slice is presented in Fig. 6,
featuring orientations that are consistent with the known anatomy, from the
one-fiber corpus callosum to the three-way crossing of the centrum semiovale.
To assess the uncertainty on estimation for these orientations within white
matter, Fig. 7 presents the distribution of the median angular deviation
$\mathrm{Med}(\Delta\Theta(\mathrm{Max[fFA]}))$ for the orientation of the
intra-voxel fascicle with maximal fFA, chosen as a proxy to evaluate the same
fascicle orientation, across 100 stratified bootstrap runs of Magic DIAMOND
for LL, LP and LS.

Figure 6: Local orientations of the fascicles estimated by Magic DIAMOND for
our LP in vivo dataset in a white-matter masked coronal slice. These
orientations, consistent with the expected anatomy, were computed after one
stratified bootstrap signal realization. Inset: Zoom on the three-way crossing
of the centrum semiovale. Figure 7: Median angular deviation
$\mathrm{Med}(\Delta\Theta(\mathrm{Max[fFA]}))$ (in degrees) across 100
stratified bootstrap runs of Magic DIAMOND, computed for the orientation of
the intra-voxel fascicle with maximal fFA. For a given stratified bootstrap
run, $\Delta\Theta(\mathrm{Max[fFA]})$ was computed voxel-wise as the angular
difference between the orientation of the Max[fFA] fascicle within this run
and the median orientation of this fascicle across stratified bootstrap runs.
$\mathrm{Med}(\Delta\Theta(\mathrm{Max[fFA]}))$ is here presented in two ways.
Panels (a), (b) and (c): White-matter masked coronal views of
$\mathrm{Med}(\Delta\Theta(\mathrm{Max[fFA]}))$ (greyscale maps) for LL, LP
and LS, respectively, with superimposed median orientations of the Max[fFA]
fascicles across stratified bootstrap runs (orientation-colored peaks). Panel
(d): Distribution of the values of
$\mathrm{Med}(\Delta\Theta(\mathrm{Max[fFA]}))$ within our whole-brain white-
matter mask for LL (red), LP (green) and LS (blue). The main features of these
distributions, i.e. medians and interquartile ranges (IQRs), are quantified as
boxplots.

Figs. 8-9-10 report on the white-matter distributions of three scalar metrics
estimated throughout the stratified bootstrap experiment: the FW signal
fraction $f_{\mathrm{FW}}$, the maximal fascicle axial diffusivity Max[fAD]
and the maximal fascicle radial diffusivity Max[fRD], respectively. The
Max[$\,\cdot\,$] operator was used for fAD and fRD as a proxy to identify
between each bootstrap realization the same fascicle and to show a single
whole-brain map of Max[fAD] and Max[fRD] in Figs. 9-10.

Figure 8: Median value (Med) and interquartile range (IQR) of the free water
signal fraction $f_{\mathrm{FW}}$ over 100 stratified bootstrap runs of Magic
DIAMOND. Panels (a), (b) and (c): White-matter masked coronal views of
$\mathrm{Med}(f_{\mathrm{FW}})$ (greyscale maps) for LL, LP and LS,
respectively. Panels (d) and (e): Distributions of the values of
$\mathrm{Med}(f_{\mathrm{FW}})$ and $\mathrm{IQR}(f_{\mathrm{FW}})$,
respectively, within our whole-brain white-matter mask for LL (red), LP
(green) and LS (blue). The main features of these distributions, i.e. medians
and interquartile ranges, are quantified as boxplots. Figure 9: Median value
(Med) and interquartile range (IQR) of the maximal fascicle axial diffusivity
($\mathrm{Max[fAD]}$, in µm2/ms) over 100 stratified bootstrap runs of Magic
DIAMOND. Layout conventions are identical to those of Fig. 8. Figure 10:
Median value (Med) and interquartile range (IQR) of the maximal fascicle
radial diffusivity ($\mathrm{Max[fRD]}$, in µm2/ms) over 100 stratified
bootstrap runs of Magic DIAMOND. Layout conventions are identical to those of
Fig. 8.

###  IV.3 Fascicle metric mapping on tractograms

Fig. 11 shows tractograms where we color-mapped the median and interquartile
range of fFA across all bootstrap realizations along the corpus callosum, the
arcuate fasciculus and the corticospinal tract.

###  IV.4 Visualization of the optimization landscape with in vivo data

Figs. 12-13-14-15 provide visualizations of the minimization space with LL, LP
and LS. Specifically, they show the three-dimensional subspace of parameters
covering various axial diffusivities $\lambda^{\parallel}$, radial
diffusivities $\lambda^{\perp}$ and signal fractions $f$ for the intra-voxel
Max[fFA] fascicle (Figs. 12-13), and two-dimensional cuts of this subspace
(Figs. 14-15), in two different voxels of interest: one in the corpus callosum
(CC, single fiber), the other in the centrum semiovale (CS, three-way
crossing). Note that the Magic-DIAMOND minimization was in fact performed in a
six-dimensional space for the corpus callosum and in a sixteen-dimensional
space for the centrum semiovale (one parameter for the free-water compartment,
five parameters for each of the three fascicles).

Figure 11: Color-mapping of Magic DIAMOND metrics onto bundles of interest.
Panel (a): Coronal view of the tractograms associated to the corpus callosum,
the arcuate fasciculus and the corticospinal tract. Panel (b): Sagittal view
of the tractograms associated to the left arcuate fasciculus. The median (Med)
and interquartile range (IQR) across all stratified bootstrap realizations of
the fFA associated to the Magic-DIAMOND fascicle most aligned with the local
tract orientation were color-mapped along these bundles for each encoding
combination described in sections III.6-III.7. Figure 12: Visualization of the
minimization landscape in a corpus-callosum voxel for the LL, LP and LS signal
combinations described in section III.6. The cost function was computed on
sampled grid points in a three-dimensional parameter subspace encapsulating
the axial diffusivity $\lambda^{\parallel}$, the radial diffusivity
$\lambda^{\perp}$ and the signal fraction $f$ of the single fascicle found in
this corpus-callosum voxel (which has maximal fFA by default). Color encodes
the value of the cost function, from its minimal value (blue) to its maximal
value (red). Opacity decreases with increasing cost function to focus on
regions of the parameter subspace that are relevant for minimization.
Diffusivities are expressed in µm2/ms. While the filled blue point is
associated to the minimal cost function found by Magic DIAMOND, the filled red
point corresponds to the global minimum in the sampled parameter subspace.
Empty circles denote projections of the filled points. Figure 13:
Visualization of the minimization landscape in a centrum-semiovale voxel for
the LL, LP and LS signal combinations described in section III.6. The cost
function was computed on sampled grid points in a three-dimensional parameter
subspace encapsulating the axial diffusivity $\lambda^{\parallel}$, the radial
diffusivity $\lambda^{\perp}$ and the signal fraction $f$ of the fascicle with
maximal fFA in this centrum-semiovale voxel that contained three intra-voxel
fascicles. Layout conventions are identical to those of Fig. 12. Figure 14:
Visualization of two-dimensional cuts of interest in the three-dimensional
minimization landscape shown in Fig. 12 for the LL, LP and LS signal
combinations described in section III.6. The color maps indicate the value of
the cost function, from its minimal value (blue) to its maximal value (red),
within two-dimensional projections intersecting the Magic DIAMOND minimum
(blue point, top) and the global minimum in the aforementioned three-
dimensional subspace (red point, bottom). White areas correspond to unphysical
fiber configurations where $\lambda^{\perp}>\lambda^{\parallel}$.
Diffusivities are expressed in µm2/ms. Figure 15: Visualization of two-
dimensional cuts of interest in the three-dimensional minimization landscape
shown in Fig. 13 for the LL, LP and LS signal combinations described in
section III.6. Layout conventions are identical to those of Fig. 14.

##  V Discussion

The in silico results presented at finite SNR=40 and infinite SNR in Fig. 5
report on the accuracy and precision of Magic DIAMOND’s estimations for
various encoding combinations (LL, LP and LS, described in section III.3)
within a numerical phantom that matches DIAMOND’s compartmental assumptions
(see section III.4 and Fig. 4). The LL encoding combination yields the most
accurate and precise Magic-DIAMOND estimations, both at infinite and finite
SNRs. At constant SNR, precision decreases when going from LL to LP, and from
LP to LS. In terms of accuracy, LP and LS present infinite-SNR biases that are
absent in LL. While LS’s accuracy improves at finite SNR when estimating non-
orientational information, highlighting an acute effect of noise in the data,
LP’s accuracy exhibits a resilience to noise comparable to that of LL’s. Our
non-parametric Mann-Whitney U-tests indicate that there is no statistical
difference between the estimations yielded by LL and LS, in contrast with
those yielded by LP for most metrics.

In vivo, uncertainty on parameter estimation is evaluated for Magic DIAMOND
via the stratified bootstrap procedure detailed in section III.5. From the
standpoint of diffusion orientation information alone, while Fig. 6 features
estimated local orientations that are consistent with the known anatomy, Fig.
7 shows that angular uncertainty increases when going from LL to LP, and from
LP to LS, in agreement with the in silico results of Fig. 5. Spatially,
angular uncertainty tends to maximize in expected areas, i.e. cortical grey
matter, deep grey nuclei, and crossing-fascicle regions.

From the standpoint of non-orientational metrics, Figs. 8-9-10 demonstrate
that while LP presents similar uncertainties on estimation of these parameters
compared to LL, as quantified by the median values of the white-matter
distributions of $\mathrm{IQR}(f_{\mathrm{FW}})$,
$\mathrm{IQR}(\mathrm{Max[fAD]})$ and $\mathrm{IQR}(\mathrm{Max[fRD]})$
(panels (e) of Figs. 8-9-10), LS exhibits a rather substantial increase in
these uncertainties. The median values of the white-matter distributions of
$\mathrm{Med}(f_{\mathrm{FW}})$, $\mathrm{Med}(\mathrm{Max[fAD]})$ and
$\mathrm{Med}(\mathrm{Max[fRD]})$ (panels (d) of Figs. 8-9-10) underline the
overall changes in Magic DIAMOND’s estimations between diffusion-encoding
strategies. Both LP and LS tend to yield lower values of $f_{\mathrm{FW}}$ and
higher values of Max[fAD] and Max[fRD] than LL. Although not shown here, this
increase in Max[fAD] and Max[fRD] leads to an increase in maximal fascicle
mean diffusivity (Max(fMD)), from $\mathrm{Med(Max[fMD])}\simeq 0.65$
µ$\mathrm{m}^{2}/\mathrm{ms}$ (LL) to $\mathrm{Med(Max[fMD])}\simeq 0.8$
µ$\mathrm{m}^{2}/\mathrm{ms}$ (LS), but rather constant maximal fascicle
fractional anisotropy $0.65\leq\mathrm{Med(Max[fFA])}\leq 0.7$. We acknowledge
that a few discrepancies exist between the in silico results of Fig. 5 and the
in vivo results of Figs. 8-9-10, which may originate from an insufficient
match between our simulated phantom and the microstructural complexity of
actual in vivo tissues.

The fascicle-specific metrics estimated by Magic DIAMOND can be combined with
tractography using their associated local orientations, as shown in Fig. 11
with the fascicle fractional anisotropy fFA mapped onto tracks associated to
the corpus callosum, the corticospinal tract and the arcuate fasciculus.
Outputted by our stratified bootstrap procedure, the median fFA along bundles
is rather constant within the bulk of these bundles, except close to cortical
grey matter and deep grey nuclei, within fascicle-crossing regions, and in
high curvature areas (e.g. high curvature area in the arcuate fasciculus, see
panel (b) of Fig. 11). The lower Med[fFA] in these regions is explained by the
higher axonal dispersion which, as shown in Ref. Scherrer_aDIAMOND:2017, , is
captured by an increased heterogeneity but also by a decreased fAD and
increased fRD. These particular areas also appear bright in the IQR-colored
bundles, reflecting a loss in precision on parameter estimation, in agreement
with Fig. 7.

Finally, Figs. 14 and 15 show that the centrum-semiovale (CS) optimization
landscape is flatter than the corpus-callosum (CC) optimization landscape in
the $f$-$\lambda^{\parallel}$ and $f$-$\lambda^{\perp}$ projections, and is
more peaked than the CC optimization landscape in the
$\lambda^{\parallel}$-$\lambda^{\perp}$ projections. But we note that the
optimization landscape exhibits no consistent change upon introduction of non-
conventional b-tensors (LS and LP). This fact is reinforced by the three-
dimensional illustration of the optimization landscape presented in Figs. 12
and 13: LL, LP and LS feature comparable separations between the minimum found
by Magic DIAMOND and the global minimum of the sampled subspace.

##  VI Conclusions

In this work, we derive for the first time a general Laplace transform of the
non-central matrix-variate Gamma distribution. This Laplace transform, valid
for any arbitrary b-tensor, enables the design of the Magic DIAMOND model, i.e
the combination of the DIAMOND model, capable of separately characterizing
crossing fascicles, with tensor-valued diffusion encoding. Our new theoretical
approach aims at drawing from the complementary pieces of diffusion
information yielded by tensor-valued diffusion encoding in order to refine the
orientation and diffusional features of each intra-voxel fascicle estimated by
DIAMOND. By design, it therefore offers new possibilities to quantify the
optimization of tensor-valued encoded acquisition schemes (Bates_ISMRM:2019,
).

Our in vivo evaluations suggest that introducing non-conventional b-tensors
(either planar or spherical) leads to a lower estimated free-water signal
fraction in the white matter (Fig. 8). This change is in agreement with the
understanding that there is none or very little freely diffusing water in the
densely packed white matter when using typical clinically feasible diffusion
times (Dhital:2018, ), which in turn supports the fact that non-conventional
b-tensors sample new information about the underlying tissue microstructure
and allow improved separation of isotropic and anisotropic diffusion
components. Importantly, our statistical evaluation via stratified bootstrap
shows that spherical encoding consistently lowers the precision on parameter
estimation. This was verified in both our voxel-based evaluation throughout
white matter (Figs. 7-8-9-10) and in our fixel-based evaluation in which we
colored each tract streamline with the IQR of the corresponding estimated
compartment (Fig. 11). First, this is likely explained by the intrinsic lack
of information linking the fascicles’ orientations to their diffusional
features when using spherical encoding (see Eq. (26)). Second, the inherently
low signal-to-noise ratio (SNR) of spherical encoding may also be involved in
lowering precision. In a conventional linear acquisition, the diffusion-
attenuated signal is low when the diffusion gradient orientation is aligned
with a fascicle orientation. Because spherical encoding samples all
orientations, a decrease of spherical signal occurs for any fascicle
orientation, leading to an SNR that is always lower than that with linear
encoding. However, the SNR depends on the acquisition scheme and scanning
hardware, so that the difference in SNR between distinct tensor-valued
acquisitions may be mitigated in other acquisition setups. As a result, we
found that at constant scan duration time, the isotropic information brought
by spherical encoding has too low SNR and not enough information about the
fascicles themselves to improve the estimation of fascicle-specific
diffusivities and orientation of a multi-fascicle model without increasing
uncertainty. In contrast, the introduction of planar encoding did not lower
precision on parameter estimation. However, we surprisingly did not observe a
clear advantage of our linear-planar dataset (LP) over our purely linear
dataset (LL). Indeed, not only did it not reduce angular uncertainty on the
fascicles’ orientations (Fig. 7), unlike what could have been expected
(Cottaar:2019, ), it also did not provide a more constraining optimization
landscape (Figs. 12-13-14-15).

Nonetheless, combining non-trivial diffusion encodings with the multi-fascicle
DIAMOND model should not be discarded right away. Indeed, potential
improvements upon introducing tensor-valued diffusion encoding may be yielded
by the following families of solutions. First, from a data-acquisition
standpoint, one could for instance try to obtain higher SNRs and/or higher
b-values, by using higher gradient amplitudes and shortening TE as in Ref.
Szczepankiewicz_DIVIDE:2019, , in order to further differentiate linear and
planar pieces of diffusion information (Jespersen:2012, ). Alternatively, the
acquisition protocol could be further optimized by maximizing correlations
between acquisition dimensions so that to bring tighter constraints to the
Magic-DIAMOND minimization. Combining at least two diffusion encodings might
also simply make Magic DIAMOND’s parameter estimation more robust to data
undersampling, an effect that has already been observed in the past
(Szczepankiewicz_Jeurissen_ISMRM:2018, ). Second, from a data-processing point
of view, phase correction strategies (Eichner:2015, ; Pizzolato:2020, ) may be
practical given the inherently lower SNR of planar and spherical encodings
compared to linear encoding. Third, the current design of the DIAMOND model
may be too constrained to take advantage of the power of tensor-valued
diffusion encoding. Such limitation was previously noticed in Ref.
Lampinen_CODIVIDE:2017,  when using spherical encoding within the original
NODDI model (Zhang_NODDI:2012, ). The importance of DIAMOND’s constraints
could be investigated by opening up the model so that it becomes degenerate
for LL but not for LP or LS. As the mathematics of the Magic DIAMOND model are
now laid down, we invite the diffusion MRI community to use them and challenge
or reproduce our conclusions.

## Acknowledgments

The authors would like to thank the "Axe Imagerie Médicale du Centre de
recherche du CHUS" and the "Centre d’imagerie médicale de l’Université de
Sherbrooke" for sponsored image acquisition time. The study was also supported
by the Neuroinformatics Chaire of the Université de Sherbrooke. We also thank
the Quebec bio-imaging network (QBIN) for partial funding of first author A.
Reymbaut.

## References

  * (1) P. J. Basser, J. Mattiello, and D. LeBihan,  Biophys J 66, 259 (1994)http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1275686/.
  * (2) E. O. Stejskal and J. E. Tanner,  The Journal of Chemical Physics 42, 288 (1965)http://dx.doi.org/10.1063/1.1695690.
  * (3) A. L. Alexander, K. M. Hasan, M. Lazar, J. S. Tsuruda, and D. L. Parker,  Magnetic Resonance in Medicine 45, 770 (2001)http://dx.doi.org/10.1002/mrm.1105.
  * (4) A. L. Alexander, J. E. Lee, M. Lazar, and A. S. Field,  Neurotherapeutics 4, 316 (2007)https://doi.org/10.1016/j.nurt.2007.05.011.
  * (5) D. K. Jones and M. Cercignani,  NMR in Biomedicine 23, 803 (2010)http://dx.doi.org/10.1002/nbm.1543.
  * (6) D. K. Jones, T. R. Knösche, and R. Turner,  NeuroImage 73, 239 (2013)http://www.sciencedirect.com/science/article/pii/S1053811912007306.
  * (7) C. Pierpaoli, P. Jezzard, P. J. Basser, A. Barnett, and G. D. Chiro,  Radiology 201, 637 (1996)https://doi.org/10.1148/radiology.201.3.8939209,  PMID: 8939209.
  * (8) G. Douaud et al.,  NeuroImage 55, 880 (2011)http://www.sciencedirect.com/science/article/pii/S1053811910015909.
  * (9) B. Jeurissen, A. Leemans, J.-D. Tournier, D. K. Jones, and J. Sijbers,  Human Brain Mapping 34, 2747 (2013)http://dx.doi.org/10.1002/hbm.22099.
  * (10) D. G. Norris,  NMR in Biomedicine 14, 77 (2001)http://dx.doi.org/10.1002/nbm.682.
  * (11) J. V. Sehy, J. J. Ackerman, and J. J. Neil,  Magnetic Resonance in Medicine 48, 765 (2002)http://dx.doi.org/10.1002/mrm.10301.
  * (12) L. Minati and W. P. Węglarz,  Concepts in Magnetic Resonance Part A 30A, 278 (2007)http://dx.doi.org/10.1002/cmr.a.20094.
  * (13) R. V. Mulkern, S. J. Haker, and S. E. Maier,  Magn Reson Imaging 27, 1151 (2009)http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2894527/,  19520535[pmid].
  * (14) C. Johnson,  Journal of Magnetic Resonance, Series A 102, 214 (1993)http://www.sciencedirect.com/science/article/pii/S1064185883710934.
  * (15) M. Nilsson,  Biophysical modelling in diffusion MRI: The role of tissue microstructure and water exchange,  PhD thesis, Medical Radiation Physics, Lund, 2011.
  * (16) M. Nilsson et al.,  Magnetic Resonance in Medicine 69, 1572 (2013)http://dx.doi.org/10.1002/mrm.24395.
  * (17) M. Nilsson, D. van Westen, F. Ståhlberg, P. C. Sundgren, and J. Lätt,  Magnetic Resonance Materials in Physics, Biology and Medicine 26, 345 (2013)https://doi.org/10.1007/s10334-013-0371-x.
  * (18) B. Lampinen et al.,  Magnetic Resonance in Medicine 77, 1104 (2017)http://dx.doi.org/10.1002/mrm.26195.
  * (19) B. Jian, B. C. Vemuri, E. Özarslan, P. R. Carney, and T. H. Mareci,  Neuroimage 37, 164 (2007)http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2576290/,  17570683[pmid].
  * (20) D. G. Cory, A. N. Garroway, and J. B. Miller,  Polym Prepr 31, 149 (1990).
  * (21) S. Mori and P. C. M. Van Zijl,  Magnetic Resonance in Medicine 33, 41 (1995)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.1910330107.
  * (22) E. C. Wong, R. W. Cox, and A. W. Song,  Magnetic Resonance in Medicine 34, 139 (1995)http://dx.doi.org/10.1002/mrm.1910340202.
  * (23) S. Eriksson, S. Lasic, and D. Topgaard,  Journal of Magnetic Resonance 226, 13 (2013)http://www.sciencedirect.com/science/article/pii/S1090780712003370.
  * (24) C.-F. Westin et al.,  Medical image computing and computer-assisted intervention : MICCAI … International Conference on Medical Image Computing and Computer-Assisted Intervention 17, 209 (2014)https://www.ncbi.nlm.nih.gov/pubmed/25320801,  25320801[pmid].
  * (25) S. Eriksson, S. Lasič, M. Nilsson, C.-F. Westin, and D. Topgaard,  The Journal of Chemical Physics 142, 104201 (2015)http://dx.doi.org/10.1063/1.4913502.
  * (26) C.-F. Westin et al.,  NeuroImage 135, 345 (2016)http://www.sciencedirect.com/science/article/pii/S1053811916001488.
  * (27) D. Topgaard,  Journal of Magnetic Resonance 275, 98 (2017)http://www.sciencedirect.com/science/article/pii/S1090780716302701.
  * (28) B. Scherrer et al.,  Magnetic Resonance in Medicine 76, 963 (2016)http://dx.doi.org/10.1002/mrm.25912.
  * (29) B. Scherrer et al.,  Decoupling axial and radial tissue heterogeneity in diffusion compartment imaging,  in Information Processing in Medical Imaging (IPMI), pp. 440–452, 2017.
  * (30) M. Chamberland et al.,  Human Brain Mapping 38, 509 (2017)http://dx.doi.org/10.1002/hbm.23399.
  * (31) Z. Eaton-Rosen et al.,  NeuroImage 162, 65 (2017)http://www.sciencedirect.com/science/article/pii/S1053811917306596.
  * (32) B. Scherrer et al.,  Diffusion compartment imaging reveals microstructural injuries in a mouse model of mild traumatic brain injury,  in Proc Intl Soc Mag Reson Med, 2017.
  * (33) B. Scherrer et al.,  Cerebral Cortex , in press (2019).
  * (34) E. R. Andrew, A. Bradbury, and R. G. Eades,  Nature 183, 1802 (1959)https://doi.org/10.1038/1831802a0.
  * (35) D. A. Yablonskiy and A. L. Sukstanskii,  NMR in Biomedicine 23, 661 (2010)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.1520.
  * (36) I. O. Jelescu and M. D. Budde,  Frontiers in Physics 5, 61 (2017)https://www.frontiersin.org/article/10.3389/fphy.2017.00061.
  * (37) D. S. Novikov, V. G. Kiselev, and S. N. Jespersen,  Magnetic Resonance in Medicine 79, 3172 (2018)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.27101.
  * (38) Y. Assaf, R. Z. Freidlin, G. K. Rohde, and P. J. Basser,  Magnetic Resonance in Medicine 52, 965https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.20274.
  * (39) Y. Assaf and P. J. Basser,  NeuroImage 27, 48 (2005)http://www.sciencedirect.com/science/article/pii/S1053811905002259.
  * (40) Y. Assaf, T. Blumenfeld-Katzir, Y. Yovel, and P. J. Basser,  Magnetic Resonance in Medicine 59, 1347 (2008)http://dx.doi.org/10.1002/mrm.21577.
  * (41) H. Zhang, T. Schneider, C. A. Wheeler-Kingshott, and D. C. Alexander,  NeuroImage 61, 1000 (2012)http://www.sciencedirect.com/science/article/pii/S1053811912003539.
  * (42) D. C. Alexander et al.,  NeuroImage 52, 1374 (2010)http://www.sciencedirect.com/science/article/pii/S1053811910007755.
  * (43) F. Aboitiz, A. B. Scheibel, R. S. Fisher, and E. Zaidel,  Brain Research 598, 143 (1992)http://www.sciencedirect.com/science/article/pii/000689939290178C.
  * (44) B. Lampinen et al.,  NeuroImage 147, 517 (2017)http://www.sciencedirect.com/science/article/pii/S105381191630670X.
  * (45) M. Nilsson, S. Lasič, I. Drobnjak, D. Topgaard, and C.-F. Westin,  NMR in Biomedicine 30, e3711 (2017)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.3711, https://onlinelibrary.wiley.com/doi/pdf/10.1002/nbm.3711,  e3711 nbm.3711.
  * (46) D. A. Yablonskiy, G. L. Bretthorst, and J. J. Ackerman,  Magnetic Resonance in Medicine 50, 664 (2003)http://dx.doi.org/10.1002/mrm.10578.
  * (47) S. Jbabdi, S. N. Sotiropoulos, A. M. Savio, M. Graña, and T. E. J. Behrens,  Magnetic Resonance in Medicine 68, 1846 (2012)http://dx.doi.org/10.1002/mrm.24204.
  * (48) M. E. Moseley et al.,  Radiology 176, 439 (1990)https://doi.org/10.1148/radiology.176.2.2367658.
  * (49) A. K. Guptar and D. K. Nagar,  Matrix Variate Distributions (Chapman and Hall/CRC, Boca Raton, Florida, 2000).
  * (50) J. Mattiello, P. J. Basser, and D. Le Bihan,  Magnetic Resonance in Medicine 37, 292 (1997)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.1910370226.
  * (51) S. Lasič, F. Szczepankiewicz, S. Eriksson, M. Nilsson, and D. Topgaard,  Frontiers in Physics 2, 11 (2014)http://journal.frontiersin.org/article/10.3389/fphy.2014.00011.
  * (52) F. Szczepankiewicz et al.,  NeuroImage 104, 241 (2015)http://www.sciencedirect.com/science/article/pii/S105381191400799X.
  * (53) F. Szczepankiewicz et al.,  NeuroImage 142, 522 (2016)http://www.sciencedirect.com/science/article/pii/S1053811916303457.
  * (54) E. Özarslan, N. Shemesh, and P. J. Basser,  J Chem Phys 130, 104702 (2009)https://www.ncbi.nlm.nih.gov/pubmed/19292544,  19292544[pmid].
  * (55) N. Shemesh, E. Özarslan, M. E. Komlosh, P. J. Basser, and Y. Cohen,  NMR in Biomedicine 23, 757 (2010)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.1550.
  * (56) M. Lawrenz, M. A. Koch, and J. Finsterbusch,  Journal of Magnetic Resonance 202, 43 (2010)http://www.sciencedirect.com/science/article/pii/S1090780709002717.
  * (57) J. Valette et al.,  Magnetic Resonance in Medicine 68, 1705 (2012)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.24193.
  * (58) S. N. Jespersen,  NMR in Biomedicine 25, 813 (2012)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.1808.
  * (59) S. N. Jespersen, H. Lundell, C. K. Sønderby, and T. B. Dyrby,  NMR in Biomedicine 26, 1647 (2013)http://dx.doi.org/10.1002/nbm.2999.
  * (60) M. Lawrenz, S. Brassen, and J. Finsterbusch,  NeuroImage 141, 313 (2016)http://www.sciencedirect.com/science/article/pii/S105381191630338X.
  * (61) A. Ianuş, I. Drobnjak, and D. C. Alexander,  NMR in Biomedicine 29, 672 (2016)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.3496.
  * (62) M. Lawrenz and J. Finsterbusch,  Magnetic Resonance in Medicine 0 (2018)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.27451.
  * (63) A. Ianuş et al.,  NeuroImage 183, 934 (2018)http://www.sciencedirect.com/science/article/pii/S1053811918307328.
  * (64) J. H. Jensen and J. A. Helpern,  NMR in Biomedicine 31, e3930 (2018)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.3930,  e3930 NBM-17-0275.R2.
  * (65) B. Jeurissen and F. Szczepankiewicz,  Spherical deconvolution of diffusion MRI data with tensor-valued encodings,  in Proc Intl Soc Mag Reson Med, 2018.
  * (66) M. Herberthson, C. Yolcu, H. Knutsson, C.-F. Westin, and E. Özarslan,  Scientific Reports 9, 4899 (2019)https://doi.org/10.1038/s41598-019-41317-8.
  * (67) D. S. Novikov, E. Fieremans, S. N. Jespersen, and V. G. Kiselev,  NMR in Biomedicine 0, e3998 (2018)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.3998.
  * (68) E. Fieremans et al.,  Effect of combining linear with spherical tensor encoding on estimating brain microstructural parameters,  in Proc Intl Soc Mag Reson Med, 2018.
  * (69) M. Reisert, V. G. Kiselev, and B. Dhital,  Magnetic Resonance in Medicine 81, 3819 (2019)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.27685.
  * (70) S. Coelho, J. M. Pozo, S. N. Jespersen, D. K. Jones, and A. F. Frangi,  Magnetic Resonance in Medicine 82, 395 (2019)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.27714.
  * (71) S. Coelho, J. M. Pozo, S. N. Jespersen, and A. F. Frangi,  arXiv e-prints , arXiv:1907.06139 (2019)https://ui.adsabs.harvard.edu/abs/2019arXiv190706139C, 1907.06139.
  * (72) J. P. de Almeida Martins and D. Topgaard,  Phys. Rev. Lett. 116, 087601 (2016)https://link.aps.org/doi/10.1103/PhysRevLett.116.087601.
  * (73) J. P. de Almeida Martins and D. Topgaard,  Scientific Reports 8, 2488 (2018)https://doi.org/10.1038/s41598-018-19826-9.
  * (74) D. Topgaard,  NMR in Biomedicine 32, e4066 (2019)https://onlinelibrary.wiley.com/doi/abs/10.1002/nbm.4066.
  * (75) C.-F. Westin et al.,  Medical Image Analysis 6, 93 (2002)http://www.sciencedirect.com/science/article/pii/S1361841502000531.
  * (76) U. Haeberlen,  High resolution NMR in solids : selective averaging (Academic Press, New York, 1976).
  * (77) H. Akaike,  IEEE Transactions on Automatic Control 19, 716 (1974)http://ieeexplore.ieee.org/document/1100705/.
  * (78) J. Besag,  Journal of the Royal Statistical Society B 48, 48 (1986)http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.399.7501.
  * (79) J. Sjölund et al.,  Journal of Magnetic Resonance 261, 157 (2015)http://www.sciencedirect.com/science/article/pii/S1090780715002451.
  * (80) F. Szczepankiewicz, C.-F. Westin, and M. Nilsson,  Magnetic Resonance in Medicine 82, 1424 (2019)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.27828.
  * (81) H. Lundell et al.,  Scientific Reports 9, 9026 (2019)https://doi.org/10.1038/s41598-019-45235-7.
  * (82) F. Szczepankiewicz, J. Sjölund, F. Ståhlberg, J. Lätt, and M. Nilsson,  PLOS ONE 14, 1 (2019)https://doi.org/10.1371/journal.pone.0214238.
  * (83) F. Szczepankiewicz, S. Hoge, and C.-F. Westin,  Data in Brief 25, 104208 (2019)http://www.sciencedirect.com/science/article/pii/S2352340919305621.
  * (84) M. Nilsson, F. Szczepankiewicz, D. van Westen, and O. Hansson,  PLOS ONE 10, 1 (2015)https://doi.org/10.1371/journal.pone.0141825.
  * (85) P. F. Neher, F. B. Laun, B. Stieltjes, and K. H. Maier-Hein,  Magnetic Resonance in Medicine 72, 1460 (2014)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.25045.
  * (86) M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche,  Magnetic Resonance in Medicine 56, 395 (2006)https://onlinelibrary.wiley.com/doi/abs/10.1002/mrm.20948.
  * (87) P. E. McKnight and J. Najab,  Mann-Whitney U Test (American Cancer Society, 2010), pp. 1–1, https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470479216.corpsy0524.
  * (88) B. Whitcher, D. S. Tuch, J. J. Wisco, A. G. Sorensen, and L. Wang,  Human Brain Mapping 29, 346 (2008)https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.20395.
  * (89) S. Chung, Y. Lu, and R. G. Henry,  NeuroImage 33, 531 (2006)http://www.sciencedirect.com/science/article/pii/S1053811906007403.
  * (90) B. Scherrer and S. K. Warfield,  PLOS ONE 7, 1 (2012)https://doi.org/10.1371/journal.pone.0048232.
  * (91) D. A. Raffelt et al.,  NeuroImage 117, 40 (2015)http://www.sciencedirect.com/science/article/pii/S1053811915004218.
  * (92) M. Chamberland, K. Whittingstall, D. Fortin, D. Mathieu, and M. Descoteaux,  Frontiers in Neuroinformatics 8, 59 (2014)https://www.frontiersin.org/article/10.3389/fninf.2014.00059.
  * (93) H. Takemura, C. F. Caiafa, B. A. Wandell, and F. Pestilli,  PLOS Computational Biology 12, 1 (2016)https://doi.org/10.1371/journal.pcbi.1004692.
  * (94) E. Garyfallidis et al.,  NeuroImage 170, 283 (2018)http://www.sciencedirect.com/science/article/pii/S1053811917305839,  Segmenting the Brain.
  * (95) F. Rheault et al.,  NeuroImage 186, 382 (2019)http://www.sciencedirect.com/science/article/pii/S1053811918320883.
  * (96) A. Bates, A. Daducci, and E. Caruyer,  Multi-dimensional diffusion mri sampling scheme: B-tensor design and accurate signal reconstruction,  Proc Intl Soc Mag Reson Med, 2019.
  * (97) B. Dhital, E. Kellner, V. G. Kiselev, and M. Reisert,  NeuroImage 182, 398 (2018)http://www.sciencedirect.com/science/article/pii/S1053811917308856,  Microstructural Imaging.
  * (98) M. Cottaar et al.,  Improved fibre dispersion estimation using b-tensor encoding, 2019, 1901.05820.
  * (99) C. Eichner et al.,  NeuroImage 122, 373 (2015)http://www.sciencedirect.com/science/article/pii/S1053811915006941.
  * (100) M. Pizzolato, G. Gilbert, J.-P. Thiran, M. Descoteaux, and R. Deriche,  NeuroImage 206, 116274 (2020)http://www.sciencedirect.com/science/article/pii/S1053811919308651.


* Bacciagaluppi, (2005) Bacciagaluppi, G. (2005).  A conceptual introduction to Nelson’s mechanics.  In Rosolino Buccheri, A. E. and Saniga, M., editors, Endophysics, Time, Quantum And The Subjective: Proceedings of the ZiF Interdisciplinary Research Workshop, Bielefeld, Germany, 17-22 January 2005, pages 367–388. World Scientific.
  * Beebee, (2000) Beebee, H. (2000).  The non-governing conception of laws of nature.  Philosophy and Phenomenological Research, 61:571–94.
  * Ben-Menahem, (2018) Ben-Menahem, Y. (2018).  Causation in science.  Princeton University Press.
  * Bhogal, (2017) Bhogal, H. (2017).  Minimal anti-Humeanism.  Australasian Journal of Philosophy, 95(3):447–460.
  * Bhogal, (2020) Bhogal, H. (2020).  Humeanism about laws of nature.  Philosophy Compass, 15(8):e12696.
  * Bird, (2007) Bird, A. (2007).  Nature’s Metaphysics: Laws and Properties.  Oxford University Press, Oxford.
  * Blanchard, (2016) Blanchard, T. (2016).  Physics and causation.  Philosophy Compass, 11(5):256–266.
  * Callender, (2004) Callender, C. (2004).  Measures, explanations and the past: Should ‘special’ initial conditions be explained?  The British Journal for the Philosophy of Science, 55(2):195–217.
  * Callender, (2017) Callender, C. (2017).  What makes time special? Oxford University Press.
  * Carroll, (1994) Carroll, J. W. (1994).  Laws of nature.  Cambridge University Press.
  * Carroll, (2018) Carroll, J. W. (2018).  Becoming Humean.  In Ott, W. and Patton, L., editors, Laws of nature, pages 122–138. Oxford University Press.
  * Carroll, (2020) Carroll, J. W. (2020).  Laws of Nature.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Winter 2020 edition.
  * Cartwright, (1983) Cartwright, N. (1983).  How the Laws of Physics Lie.  Oxford University Press.
  * (18) Cartwright, N. (1994a).  Fundamentalism vs. the patchwork of laws.  In Proceedings of the Aristotelian Society, volume 94, pages 279–292. JSTOR.
  * (19) Cartwright, N. (1994b).  Nature’s Capacities and their Measurement.  Oxford University Press.
  * Chen, (2020) Chen, E. K. (2020).  The Past Hypothesis and the nature of physical laws.  In Loewer, Barry, W. E. and Weslake, B., editors, Time’s Arrows and the Probability Structure of the World. Harvard University Press, forthcoming.
  * Chen, (2022) Chen, E. K. (2022).  Fundamental nomic vagueness.  The Philosophical Review, 131(1).
  * Choi and Fara, (2021) Choi, S. and Fara, M. (2021).  Dispositions.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Spring 2021 edition.
  * Cohen and Callender, (2009) Cohen, J. and Callender, C. (2009).  A better best system account of lawhood.  Philosophical Studies, 145(1):1–34.
  * Deckert, (2010) Deckert, D.-A. (2010).  Electrodynamic absorber theory: a mathematical study.  PhD thesis, Ludwig Maximilian University of Munich.
  * Demarest, (2017) Demarest, H. (2017).  Powerful properties, powerless laws.  In Jacobs, J. D., editor, Causal powers, pages 38–53. Oxford University Press Oxford, United Kingdom.
  * Demarest, (2019) Demarest, H. (2019).  Mentaculus laws and metaphysics.  Principia: An International Journal of Epistemology, 23(3):387–399.
  * Demarest, (2021) Demarest, H. (2021).  Powers, best systems, and explanation.  Manuscript.
  * Dorst, (2021) Dorst, C. (2021).  Productive laws in relativistic spacetimes.  Manuscript.
  * Dretske, (1977) Dretske, F. (1977).  Laws of nature.  Philosophy of Science, 44:248–68.
  * Dürr et al., (1992) Dürr, D., Goldstein, S., and Zanghì, N. (1992).  Quantum equilibrium and the origin of absolute uncertainty.  Journal of Statistical Physics, 67(5-6):843–907.
  * Ellis, (2001) Ellis, B. (2001).  Scientific essentialism.  Cambridge University Press.
  * Ellis, (2014) Ellis, B. (2014).  The philosophy of nature: A guide to the new essentialism.  Routledge.
  * Fodor, (1974) Fodor, J. A. (1974).  Special sciences (or: The disunity of science as a working hypothesis).  Synthese, 28:97–115.
  * Friederich and Evans, (2019) Friederich, S. and Evans, P. W. (2019).  Retrocausality in quantum mechanics.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, summer 2019 edition.
  * Ghirardi, (2018) Ghirardi, G. (2018).  Collapse theories.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, fall 2018 edition.
  * Ghirardi et al., (1986) Ghirardi, G., Rimini, A., and Weber, T. (1986).  Unified dynamics for microscopic and macroscopic systems.  Physical Review D, 34(2):470.
  * Goldstein, (2001) Goldstein, S. (2001).  Boltzmann’s approach to statistical mechanics.  In Bricmont, J., Dürr, D., Galavotti, M. C., Ghirardi, G., Petruccione, F., and Zanghì, N., editors, Chance in Physics, pages 39–54. Berlin: Springer.
  * Goldstein, (2012) Goldstein, S. (2012).  Typicality and notions of probability in physics.  In Probability in physics, pages 59–71. Springer.
  * Goldstein, (2017) Goldstein, S. (2017).  Bohmian mechanics.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, summer 2017 edition.
  * Goldstein and Tumulka, (2003) Goldstein, S. and Tumulka, R. (2003).  Opposite arrows of time can reconcile relativity and nonlocality.  Classical and Quantum Gravity, 20(3):557.
  * Henderson, (2020) Henderson, L. (2020).  The Problem of Induction.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Spring 2020 edition.
  * Hicks and Schaffer, (2017) Hicks, M. T. and Schaffer, J. (2017).  Derivative properties in fundamental laws.  The British Journal for the Philosophy of Science, 68(2):411–450.
  * Hildebrand, (2013) Hildebrand, T. (2013).  Can primitive laws explain?  Philosophers’ Imprint.
  * Hildebrand, (2020) Hildebrand, T. (2020).  Non-Humean theories of natural necessity.  Philosophy Compass, 15(5):e12662.
  * Hildebrand and Metcalf, (2021) Hildebrand, T. and Metcalf, T. (2021).  The nomological argument for the existence of God.  Noûs.
  * Hoefer, (2019) Hoefer, C. (2019).  Chance in the world: A Humean guide to objective chance.  Oxford University Press.
  * Khlentzos, (2021) Khlentzos, D. (2021).  Challenges to Metaphysical Realism.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Spring 2021 edition.
  * Kimpton-Nye, (2017) Kimpton-Nye, S. (2017).  Humean laws in an unHumean world.  Journal of the American Philosophical Association, 3(2):129–147.
  * Kraut, (2017) Kraut, R. (2017).  Plato.  In Zalta, E. N., editor, The Stanford Encyclopedia of Philosophy. Metaphysics Research Lab, Stanford University, Fall 2017 edition.
  * Lange, (2009) Lange, M. (2009).  Laws and lawmakers: Science, metaphysics, and the laws of nature.  Oxford University Press.
  * Lange, (2016) Lange, M. (2016).  Because Without Cause: Non-Casual Explanations In Science and Mathematics.  Oxford University Press.
  * Laudisa, (2015) Laudisa, F. (2015).  Laws are not descriptions.  International Studies in the Philosophy of Science, 29(3):251–270.
  * Lazarovici, (2018) Lazarovici, D. (2018).  Against fields.  European Journal for Philosophy of Science, 8(2):145–170.
  * Lewis, (1980) Lewis, D. (1980).  A subjectivist’s guide to objective chance.  In Jeffrey, R. C., editor, Studies in Inductive Logic and Probability, volume 2, pages 263–93. University of California Press, Berkeley.
  * Lewis, (1983) Lewis, D. (1983).  New work for a theory of universals.  Australasian Journal of Philosophy, 61:343–77.
  * Lewis, (1986) Lewis, D. (1986).  Philosophical papers II.  Oxford: Oxford University Press.
  * Lewis, (1994) Lewis, D. (1994).  Humean supervenience debugged.  Mind, 103:473–90.
  * Loewer, (2001) Loewer, B. (2001).  Determinism and chance.  Studies in History and Philosophy of Science Part B: Studies in History and Philosophy of Modern Physics, 32(4):609–620.
  * Loewer, (2012) Loewer, B. (2012).  Two accounts of laws and time.  Philosophical Studies, 160(1):115–137.
  * Loewer, (2020) Loewer, B. (2020).  The package deal account of laws and properties (PDA).  Synthese, pages 1–25.
  * Loewer, (2021) Loewer, B. (2021).  What breathes fire into the equations.  Manuscript.
  * Maudlin, (2007) Maudlin, T. (2007).  The Metaphysics Within Physics.  Oxford University Press, New York.
  * Mumford, (2004) Mumford, S. (2004).  Laws in nature.  Routledge.
  * Nelson, (1966) Nelson, E. (1966).  Derivation of the Schrödinger equation from Newtonian mechanics.  Physical Review, 150(4):1079.
  * Penrose, (1974) Penrose, R. (1974).  The role of aesthetics in pure and applied mathematical research.  Bull. Inst. Math. Appl., 10:266–271.
  * Penrose, (1979) Penrose, R. (1979).  Singularities and time-asymmetry.  In Hawking, S. and Israel, W., editors, General relativity, pages 581–638. Cambridge: Cambridge University Press.
  * Penrose, (1989) Penrose, R. (1989).  The Emperor’s New Mind: Concerning Computers, Minds, and the Laws of physics.  Oxford: Oxford University Press.
  * Roberts, (2008) Roberts, J. T. (2008).  The law-governed universe.  Oxford University Press.
  * Schaffer, (2016) Schaffer, J. (2016).  It is the business of laws to govern.  dialectica, 70(4):577–588.
  * Strawson, (2015) Strawson, G. (2015).  ‘Humeanism’.  Journal of the American Philosophical Association, 1(1):96.
  * Sutherland, (2008) Sutherland, R. I. (2008).  Causally symmetric Bohm model.  Studies in History and Philosophy of Science Part B: Studies in History and Philosophy of Modern Physics, 39(4):782–805.
  * Tooley, (1977) Tooley, M. (1977).  The nature of laws.  Canadian Journal of Philosophy, 7(4):667–698.
  * Tooley, (1997) Tooley, M. (1997).  Time, Tense and Causation.  Oxford University Press, Clarendon Press.
  * Van Fraassen, (1989) Van Fraassen, B. C. (1989).  Laws and symmetry.  Oxford University Press.
  * Vetter, (2015) Vetter, B. (2015).  Potentiality: From dispositions to modality.  Oxford University Press.
  * (76) Weatherall, J. O. (2019a).  Part 1: Theoretical equivalence in physics.  Philosophy Compass, 14(5):e12592.
  * (77) Weatherall, J. O. (2019b).  Part 2: Theoretical equivalence in physics.  Philosophy Compass, 14(5):e12591.
  * Weaver, (2019) Weaver, C. G. (2019).  Fundamental causation: physics, metaphysics, and the deep structure of the world.  Routledge.
  * Weinberg, (1992) Weinberg, S. (1992).  Dreams of a Final Theory: The Search for the Fundamental Laws of Nature.  New York: Pantheon.
  * Wheeler and Feynman, (1945) Wheeler, J. A. and Feynman, R. P. (1945).  Interaction with the absorber as the mechanism of radiation.  Reviews of modern physics, 17(2-3):157.
  * Wheeler and Feynman, (1949) Wheeler, J. A. and Feynman, R. P. (1949).  Classical electrodynamics in terms of direct interparticle action.  Reviews of modern physics, 21(3):425.
  * Wigner, (1985) Wigner, E. (1985).  Events, laws of nature, and invariance principles.  In Zichichi, A., editor, How Far Are We from the Gauge Forces–Proceedings of the 21st Course of the International School of Subnuclear Physics, Aug 3-14, 1983, pages 699–708. Plenum.
  * Wigner, (1964) Wigner, E. P. (1964).  Symmetry and conservation laws.  Proceedings of the National Academy of Sciences of the United States of America, 51(5):956–965.
  * Wilson, (1987) Wilson, M. (1987).  What is a law of nature? (book review).  The Philosophical Review, 96(3):435–441.