Robust assessment of the brain's sheet structure using normalized convolution

Chantal M.W. Tax^{1,2}, Carl-Fredrik Westin^{2}, Tom Dela Haije^{3}, Andrea Fuster^{3}, Max A. Viergever^{1}, Luc Florack^{3}, and Alexander Leemans^{1}

The hypothesis that
brain pathways cross nearly orthogonal forming 2D sheet-like structures
(similar to the “warp and weft of fabric”^{1}) layered in 3D space
“like pages of a book”^{1} has been a topic of debate in literature^{1,2,3}.
This theory was mainly supported by qualitative findings using diffusion MRI
(dMRI) tractography^{1}, and a comprehensive *quantitative *analysis is necessary to reliably assess the
occurrence of sheet structure in the brain and across individuals^{2,4}.
To this end, an approach was developed to robustly quantify the *normal component of the Lie bracket* from
dMRI data^{4,5}, which indicates the existence of sheet structure
when close to zero^{1}. This method is based on the “flows-and-limits”
definition, in which the Lie bracket $$$[V,W]$$$ at a
point $$$p$$$ is
computed as the vector deviation from $$$p$$$ (or
“closure gap”) when moving around in an infinitesimal loop along the integral
curves of two vector fields $$$V$$$ and $$$W$$$. These are assumed to coincide with the peak
directions of fiber ODFs (fODFs).

The flows-and-limits
approach relies on tractography of a large number of loops in a neighborhood, which is computationally expensive. In addition, it requires the clustering of fODF peaks along the paths and it is unclear how to
proceed in the case of missing peaks. Finally,
the used definition is an approximation for non-infinitesimal loops. Here we
present an alternative, robust, and elegant approach for Lie bracket
computation, based on a conceptually different but equivalent definition of $$$[V,W]$$$. This approach uses *normalized convolution*^{6}, which can cope with
missing or uncertain data.

**Coordinate definition**

$$$[V,W]$$$ can be computed with respect to a given coordinate system. Here, the vector fields can be represented by vector-valued functions (i.e. $$$V=(v^1,v^2,v^3)^T$$$) in standard Cartesian coordinates $$$\textbf{x}=(x^1,x^2,x^3)^T$$$, giving

$$[V,W]=J_W\cdot{V}-J_V\cdot{W}\,\,[1],$$

where $$$J$$$ is the Jacobian. The Lie bracket normal component is then $$$[V,W]^\bot=[V,W]\cdot(W\times{V})$$$.

Implementation: The fODFs in a neighborhood of size
$$$N$$$ are clustered
using a front-propagation approach^{7}. We estimate the Jacobian $$$J$$$ with normalized
convolution^{6}, where both $$$V$$$ and the operator filter basis $$$B$$$ are accompanied by a scalar
component ($$$c$$$ and $$$a$$$ respectively) representing their certainty:

$$\{aB\hat{\otimes}cV\}_N=\{aB\otimes{B}\hat{\cdot}c\}^{-1}\{aB\hat{\otimes}cV\}\,\,[2].$$

Here, general convolution is defined as $$$\{aB\hat{\otimes}cV\}(p)=\sum_\xi{a}(\xi)B(\xi)\otimes{c}(p-\xi)V(p-\xi)$$$, $$$\otimes$$$ denotes tensor product, $$$\cdot$$$ scalar product, $$$\hat{}$$$ convolution, $$$^{-1}$$$ inverse, and $$$\xi$$$ is the local spatial coordinate. We set $$$c$$$ at a particular location to zero when a vector is missing, and choose

$$a(r)=\begin{cases}\cos^\beta(\frac{\pi{r}}{2r_{max}})&\textrm{if}\;r<r_{max}\\0&\textrm{else}\end{cases} [3],$$

where $$$r$$$ is the distance to the neighborhood center, and $$$r_{max}$$$ and $$$\beta$$$ thus determine the “weight” of the neighborhood (Fig.1a). In practice, we set $$$r_{max}=0.5N\delta$$$, with $$$N$$$ and $$$\delta$$$ the kernel and voxel size . Basis functions for gradient estimation include linear ramps in three dimensions, i.e., we choose $$$B=(1,x^1,x^2,x^3)$$$:

$$\{aB\hat{\otimes}cV\}_N=\begin{bmatrix}\{a\hat{\cdot}c\}&\{ax^1\hat{\cdot}c\}&\{ax^2\hat{\cdot}c\}&\{ax^3\hat{\cdot}c\}\\\{ax^1\hat{\cdot}c\}&\{a(x^1)^2\hat{\cdot}c\}&\{ax^1x^2\hat{\cdot}c\}&\{ax^1x^3\hat{\cdot}c\}\\\{ax^2\hat{\cdot}c\}&\{ax^1x^2\hat{\cdot}c\}&\{a(x^2)^2\hat{\cdot}c\}&\{ax^2x^3\hat{\cdot}c\}\\\{ax^3\hat{\cdot}c\}&\{ax^1x^3\hat{\cdot}c\}&\{ax^2x^3\hat{\cdot}c\}&\{a(x^3)^2\hat{\cdot}c\} \end{bmatrix}^{-1} \begin{bmatrix}\{a\hat{\cdot}cv^1\}&\{a\hat{\cdot}cv^2\}&\{a\hat{\cdot}cv^3\}\\\{ax^1\hat{\cdot}cv^1\}&\{ax^1\hat{\cdot}cv^2\}&\{ax^1\hat{\cdot}cv^3\}\\\{ax^2\hat{\cdot}cv^1\}&\{ax^2\hat{\cdot}cv^2\}&\{ax^2\hat{\cdot}cv^3\}\\\{ax^3\hat{\cdot}cv^1\}&\{ax^3\hat{\cdot}cv^2\}&\{ax^3\hat{\cdot}cv^3\} \end{bmatrix}.$$

The result is a regularized vector field $$$\hat{V}$$$ and an estimate of the Jacobian $$$\hat{J_V}$$$ :

$$\{aB\hat{\otimes}cV\}_N^T=\begin{bmatrix}\hat{v}^1&\hat{\frac{\partial{v^1}}{\partial{x^1}}}&\hat{\frac{\partial{v^1}}{\partial{x^2}}}&\hat{\frac{\partial{v^1}}{\partial{x^3}}}\\\hat{v}^2&\hat{\frac{\partial{v^2}}{\partial{x^1}}}&\hat{\frac{\partial{v^2}}{\partial{x^2}}}&\hat{\frac{\partial{v^2}}{\partial{x^3}}}\\\hat{v}^3&\hat{\frac{\partial{v^3}}{\partial{x^1}}}&\hat{\frac{\partial{v^3}}{\partial{x^2}}}&\hat{\frac{\partial{v^3}}{\partial{x^3}}}\end{bmatrix}=\begin{bmatrix}\hat{V}&\hat{J_V}\end{bmatrix}\,\,[4].$$

**Vector field simulations**

We define three vector fields that are tangent to a sphere with radius $$$\rho$$$: $$$U$$$ and $$$V$$$ are tangent to the upper hemisphere, $$$W$$$ is tangent to the lower hemisphere (Fig.1b). $$$U$$$ and $$$V$$$ form a sheet such that $$$[U,V]_p^{\perp{}}=0\,\forall{}\,p\in{}S_{UV}=\left\{x\in{}R^3\vert{}{\left(x^1\right)}^2+{\left(x^2\right)}^2<{\rho{}}^2\right\}$$$ and $$$U$$$ and $$$W$$$ generally do not form a sheet at $$$p\in{}S_{UW}=S_{UV}\setminus\{\textbf{0}\}$$$:

$${\left[U,W\right]}_p^{\perp{}}=\frac{6x^1x^2{\left({\rho{}}^2-{\left(x^1\right)}^2-{\left(x^2\right)}^2\right)}^{\frac{3}{2}}}{{\rho{}}^2\left(\rho{}-{\left(x^1\right)}^2\right)\left(\rho{}-{\left(x^2\right)}^2\right)}\,\,[5].$$

We add noise using a Von-Mises-Fisher distribution with concentration parameter $$$k$$$ . These simulations allow us to test different aspects independently of the dMRI method.

**dMRI simulations**

dMRI signals (60
directions, $$$b=3000\,s/mm^2$$$
)
were simulated using a ZeppelinStickDot model^{8} with the fiber
direction defined by the previously used vector fields. fODFs and their peaks
were calculated using constrained spherical deconvolution^{9}.

**Vector field simulations**

The color plots in Fig.2 show the mean error of the estimated Lie bracket normal component $$$\widehat{{\left[\cdot{},\cdot{}\right]}}_p^{\perp{}}$$$, suggesting that high $$$h_{max}$$$ (maximum distance along one pathway) with low $$$\Delta{h}$$$ (step size), and high $$$N$$$ with low $$$\beta$$$ give the most accurate estimate for the corresponding implementation. Fig.3 shows the range and mean of $$$\widehat{{\left[\cdot{},\cdot{}\right]}}_p^{\perp{}}$$$ as function of the simulated $$${\left[\cdot{},\cdot{}\right]}_p^{\perp{}}$$$, indicating that the coordinate implementation is more precise. Fig.4 shows that the coordinate implementation is also more robust to missing peaks.

**dMRI simulations**

Fig.5 shows $$$\widehat{{\left[\cdot{},\cdot{}\right]}}_p^{\perp{}}$$$ as function of SNR, confirming that the coordinate implementation is more precise, and has a higher resolving power to distinguish sheet from non-sheet.

Fig. 1: (a) Operator applicability function $$$\alpha(r)$$$. (b) Vectorfields $$$U=(-\sin\phi_1,\cos\phi_1\cos\theta_2,\cos\phi_1\sin\theta_2),$$$ $$$V=(\cos\phi_2\cos\theta_1,-\sin\phi_2,\cos\phi_2\sin\theta_1),$$$
and $$$W=(\cos\phi_2\cos\theta_1,-\sin\phi_2,-\cos\phi_2\sin\theta_1),$$$
where $$$\theta_i=\tan^{-1}\frac{x^i}{\sqrt{\rho^2-(x^1)^2-(x^2)^2}},$$$ and $$$\phi_i=\cos^{-1}\frac{x^i}{\rho}$$$.

Fig. 2: Vector field
simulations showing the influence of parameter settings for both
implementations. Color plots show the error of the estimated Lie bracket normal
component.
Voxel size $$$1\,mm,\,\rho=26\,mm,\,p=(10,10,0),\,[U,V]^\perp_p=0$$$ (top row), $$$[U,W]^\perp_p=0.0278$$$ (bottom row), 50 noise iterations with $$$k=150$$$.

Fig. 3: Vector
field simulations showing $$$\widehat{[U,W]}^\perp_p$$$ for
different simulated $$$[U,W]^\perp_p$$$ in red
(obtained by changing $$$p$$$ at a
fixed $$$\rho=13\,mm$$$, see Eq. [5]). For reference, $$$\widehat{[U,W]}^\perp_p$$$ is
evaluated at the same $$$p$$$ and
plotted in green. Voxel
size 1 mm, 50
noise iterations with $$$k=150$$$. (a) $$$h_{max}=5,\,\Delta{h}=0.5$$$. (b) $$$N=11,\,\beta=1$$$.

Fig. 4: Vector field simulations showing $$$\widehat{[\cdot,\cdot]}^\perp_p$$$ as a function of the percentage artificial
dropouts of vectors (random throughout the whole dataset). Voxel size 1 mm, $$$\rho=26\,mm,\,p=(10,10,0)$$$, 50 noise iterations with $$$k=150$$$. (a) $$$h_{max}=5,\,\Delta{h}=0.5$$$. (b) $$$N=11,\,\beta=1$$$.

Fig. 5: dMRI simulations showing $$$\widehat{[\cdot,\cdot]}^\perp_p$$$ as a function of SNR. Voxel size 1 mm, $$$\rho=13\,mm,\,p=(3,3,0)$$$, 50 noise iterations . (a) $$$h_{max}=5,\,\Delta{h}=0.5$$$. (b) $$$N=11,\,\beta=1$$$.

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)

2073