When exploring long-range cortical connectivities using tractography, streamlines are biased to terminate in the gyral crowns rather than walls or sulci^1,2. This is partly due to the often sharp curves fibers make when crossing the white/gray matter (WM/GM) boundary^1,3, which is unresolved at typical diffusion MRI (dMRI) resolutions. We propose a geometric model of the fiber orientations across the WM/GM boundary based on histology data, which is able to describe this sharp curvature and show that we can estimate the model parameters from lower resolution dMRI data.

To model the transition in fiber orientations from WM to cortical GM, we define a radial orientation in every voxel by linearly interpolating the surface normals between the gyral walls (in the WM) and between the WM/GM boundary and pial surface (in the cortex). The gyral walls were identified as the vertices closest in Euclidean space, where the connecting line passes through the center of the voxel (Figure 1).

We model the angle $$$\theta$$$ between the radial orientation and the fiber orientation to smoothly transition across the WM/GM boundary using a sigmoid function:

$$\theta(d) = \pi / 2 - \frac{\pi / 2}{1 - \exp([d - o] / w)}, \tag{1}$$

where $$$d$$$ is the signed distance from the WM/GM boundary (taken as negative within the white matter), $$$o$$$ is the offset of the transition from the WM/GM matter boundary and $$$w$$$ is the width of the transition band. Figure 2 illustrates the tangential-to-radial transition of equation 1 in 2D. To define the fiber orientation in 3D we add a free parameter $$$\phi$$$, which determines the fiber orientation within the tangential plane.

We model the attenuation of the dMRI signal continuously at every point in space as a multi-tensor model with one isotropic and one or more anisotropic tensors. Each anisotropic tensor can have a different tangential orientation $$$\phi$$$, but the deviation from the tangential plane is determined by its distance from the transition band (Eq. 1). This continuous attenuation function is convolved with a Gaussian point spread function to obtain the discrete voxelwise attenuation.

We parcellate the cortical surface into random patches using k-means and the model variables are simultaneously fitted to all voxels corresponding to a single surface patch (in both the cortex and the superficial white matter). Therefore, given a set of dMRI measurements and a WM/GM boundary surface, we can for every surface patch estimate features of the transition band (w and o) and therefore the orientations within this band (Eq.1).

To assess the model's ability to resolve sub-voxel features of the WM/GM transition, we use in-vivo dMRI data acquired with HCP-like⁵ acquisition parameters of a human subject at three different resolutions (1.35, 2, and 2.5 mm). The dMRI data were pre-processed and the WM/GM boundary was extracted using the HCP pipeline⁶.

The fibers approximate a radial orientation to the WM/GM boundary throughout the cortex and a tangential orientation in the superficial white matter for both high-resolution 2D histology data (Figure 3)⁷ and lower resolution 3D dMRI data (Figure 4). Because radial orientations are defined relative to the gyral walls and not to the nearest surface element, this holds true even for voxels close to the gyral crown. Thus, our assumption in equation 1 of tangential orientations in the white matter and radial orientations in the cortex is well founded.

The model fitted to the dMRI data shows that consistently across all three spatial resolutions for the gyral walls and crowns the transition from tangential orientations starts at the WM/GM boundary (Figure 5a) and only become fully radial about 1.5 mm into the cortex (Figure 5b). This is in rough agreement with the histology data (i.e. Figures 2, 3), where the fibers remain tangential until the edge of the heavily myelinated region (white line in Figure 3) and then smoothly transition to radial orientations in the lower cortical layers. Figure 5c, d shows high reproducibility of the estimated parameters across different resolutions for the gyral walls and crowns, despite the sub-voxel width of the transition band and the small patch sizes used in the fit (covering only 0.15% of the surface). This reproducibility suggests that the transition band estimates capture the sub-voxel orientation patterns instead of being driven by voxel size.

We have proposed a geometric model to represent the transition of fiber orientations across the WM/GM boundary. This model is able to consistently reproduce the transition location and width from dMRI data acquired at different spatial resolutions and shows qualitative agreement with the high-resolution fiber orientations estimated from histology.