Diffusion-weighted (DW) magnetic resonance imaging (MRI) has empowered the analysis of the brain’s white matter (WM) anatomy and its relationship with function and pathologies in recent years. DWI has provided great tools to advance the study of neuroscience and neuropathology. However, the relationship between the measures obtained from the DWI signal and their underlying biological process is still unclear. In this work, we build upon current advances in DWI acquisition and signal modelling to develop a new technology to trace WM fascicles while simultaneously characterizing the distribution of axonal calibers (i.e. diameters) within the fascicle. In doing this, we expect to provide better tools to characterize the WM. Examples of these possible applications could be pathology-oriented, for instance through the detection of axonal swelling and quantification of the connectivity between two cortical areas; or neuroscience oriented, as axonal caliber has been shown to be closely related to the efficiency of electrical signal propagation within the axons [12].

Despite recent advances in axonal caliber quantification from DWI, to the best of our knowledge, none of them has been used to improve tractography at the moment of the tract tracing process. Two approaches to combine microstructure information and tractography have been proposed [9, 14]. To solve complex WM areas (e.g. branching, crossing, merging, and bottlenecks), these approaches reject tracts from a full brain tractography based on microstructure information. This is sensitive to the choice of the tractography algorithm used, since this requires all plausible configurations of tracts inside the complex region. Our novel tractography algorithm, AxTract, addresses the complex configuration problem differently. AxTract incorporates the axonal caliber estimation in the tractography algorithm and uses it during the tracing process. This produces tracts with embedded microstructure information and enables the possibility of solving the tracing through WM areas using axonal calibers information.

In developing AxTract to be useful for reasonably long DWI acquisitions, we base our model solely on the ensemble average propagator (EAP) [11] at a fixed gradient separation time. This relaxes the requirements of axonal caliber estimation techniques such as AxCaliber [3] and ActiveAx [1] which focus on the signal attenuation and need a sampling over different gradient separation times. Our approach has the main advantage of simultaneously modelling, through the Fourier slice theorem, all measurements on the perpendicular plane to the cylinder population. To prove the soundness of our model, we developed the first contribution of this work: a generalized return-to-axis probability (RTAP) measure [11] showing that even at gradient separation times used in current clinical protocols, axonal caliber can be quantified. This enabled our second contribution: the AxTract algorithm which estimates axonal caliber during the tractography process using it as a prior for the traced tract.

The DWI signal within a voxel measures diffusion of water particles within different compartments, such as axons and astrocytes. It is possible to characterize the diffusion process as the displacement probability density of water particles within these compartments, the ensemble average propagator (EAP) [11]:
where are the particle’s start and end positions; the diffusion time; C is the set of compartments; is the probability of being inside compartment c; the probability of the initial position in c; and the compartment-specific propagator. The EAP is related to the attenuation of the DWI signal by the Fourier transform:
Within the study of the human brain WM, the axons present a specific interest. Within a voxel, axons can be modelled as cylindrical segments, for which specific formulations of P and E exist [6]. In compartments where diffusion takes place within a set of cylinders oriented along direction with negligible tortuosity and permeability, the displacement is decomposable in the parallel and perpendicular directions to the cylinder [4]. We write this as leading to separable formulations of P and E [4]:
This decomposition enables the use of theoretical models for  [6]. If the propagator is measured at a cylinder of cross-sectional area A and it’s filled with water with diffusion coefficient D, we derive the following expressions [6]:

where is the indicator function for , and are the respective angles of and when expressed in polar coordinates, is the n-th cylindrical Bessel function and the k-th the root of its derivative: .

Having a specific model for cylindrical compartments, i.e. axons, we can derive the theory to estimate the axonal cross-sectional area. Different techniques to capitalize the theoretical model in Eq. (4) for cylinder compartments and measure axonal radii have been proposed. The main exponents of these are AxCaliber [3]; ActiveAx [1] and the Return-to-Axis-Probability (RTAP) [11]. However, the applicability of these techniques is limited, even in current state-of-the-art whole human brain acquisitions such as the HCP project: AxCaliber relies on a relatively dense sampling along the q and dimensions; ActiveAx estimates a single-parameter which experimentally correlates with the mean caliber without an explicit formal relationship to it; and RTAP which needs very large diffusion times () for the perpendicular EAP to become [11]
and converge to the reciprocal of the cross-sectional area of the axonal population
The first contribution of our work is to prove that, even at small values, the propagator along the cylinder has a specific relationship with the distribution of cross-sectional areas in a cylinder population. We base our model on the EAP as opposed to the AxCaliber and RTAP approaches which focus on the signal attenuation. This has the main advantage of simultaneously modelling, through the Fourier slice theorem, all measurements on the perpendicular plane to the cylinder population. We start our model in the style of AxCaliber and “infinite ” RTAP, we attach a density to the cross-sectional area of a cylinder population. Our density is based on three hypotheses given by Özarslan et al. [11]. First, each particular cylinder’s contribution to the overall signal is proportional to the ratio of water particles in it, which is in direct relationship with the cylinder’s cross-sectional area. Second, the cylinder population is Gamma-distributed [3]. This leads to specific EAP formulation, Eq. (1), for N cylinders:
where each is an independent and identically distributed random variable with Gamma distribution, of shape and rate , of the cross-sectional area. Finally, our third hypothesis assumes that the population is large enough to be approximated by an infinite number of cylinders. Combining Equations 1 and 5
where the integral over A takes in account all possible cross-sectional areas; is the probability density function of a Gamma distribution with shape and rate ; and is the average cross-sectional area under the distribution . By using the separability of the EAP (see Eq. 3) and assuming a uniform probability of finding a water particle within the cylinder population, we marginalize Eq. 6 for the return-to-axis probability, i.e. ,