The hypothesis that brain pathways cross nearly orthogonal forming 2D sheet-like structures (similar to the “warp and weft of fabric”¹) layered in 3D space “like pages of a book”¹ has been a topic of debate in literature^1,2,3. This theory was mainly supported by qualitative findings using diffusion MRI (dMRI) tractography¹, 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¹. 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⁶, 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⁷. We estimate the Jacobian $$$J$$$ with normalized convolution⁶, 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⁸ with the fiber direction defined by the previously used vector fields. fODFs and their peaks were calculated using constrained spherical deconvolution⁹.

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.

Our simulations show that the coordinate Lie bracket implementation is more accurate, precise, and robust to missing peaks than the flows-and-limits implementation proposed in previous work^1,4. In future work we will apply this method in real data to assess the occurrence of sheet structure, where we expect the clustering of genuine peaks in a neighborhood into distinct vector fields to be the biggest challenge.