Diffusion Tensor Imaging Group Analysis Using Tract Profiling and Directional Statistics

Group analysis in diffusion tensor imaging is challenging. Comparisons of tensor morphology across groups have typically been performed on scalar measures of diffusivity, such as fractional anisotropy (FA), disregarding the complex three-dimensional morphologies of diffusion tensors. Scalar measures consider only the magnitude of the diffusion but not directions. In the present study, we have introduced a new approach based on directional statistics to use directional information of diffusion tensors in statistical group analysis based on Bingham distribution. We have investigated different directional statistical models to find the best fit. During the experiments, we confirmed that carrying out directional statistical analysis along the tract is much more effective than voxel- or skeleton-guided directional statistics. Hence, we propose a new method called tract profiling and directional statistics (TPDS) applicable to fiber bundles. As a case study, the method has been applied to identify connectivity differences of patients with major depressive disorder. The results obtained with the directional statistic-based analysis are consistent with those of NBS, but additionally, we found significant changes in the right hemisphere striatum, ACC, and prefrontal, parietal, temporal, and occipital connections as well as left hemispheric differences in the limbic areas such as the thalamus, amygdala, and hippocampus. The results are also evaluated with respect to fiber lengths. Comparison with the output of the network-based statistical toolbox indicated that the benefit of the proposed method becomes much more distinctive as the tract length increases. The likelihood of finding clusters of voxels that differ in long tracts is higher in TPDS, while that relationship is not clearly established in NBS.

Introduction

Diffusion tensor imaging (DTI) can reveal complicated structural differences in patient groups by using the orientation and integrity of white matter tracts to identify white matter abnormalities. The diffusion tensor is the covariance matrix of diffusion coefficients calculated from gradient directions for each voxel. Although DTI is by nature a nonscalar image which provides directional information for the neural tracts, group-based DTI analyses are mainly conducted using scalar descriptors such as fractional anisotropy (FA) (Basser, 1995), relative anisotropy (RA) (Basser and Pierpaoli, 2011), axial diffusivity (AD), and radial diffusivity (RD) (Song et al., 2002). Such scalar metrics do not describe the full tensor shape or distribution and do not capture all of the information available in the data. By developing advanced metrics for connectivity analysis between groups of subjects in a nonscalar fashion, findings regarding abnormalities can be improved.

The principal diffusion direction (PDD), which is the eigenvector that corresponds to the largest eigenvalue of the tensor, captures the estimation of the fiber direction within the voxel. PDD has been used mainly in directionally encoded color (DEC) maps (Pajevic and Pierpaoli, 1999) which facilitate visual comparison but not quantitative group analysis. In order to evaluate PDD, which is a vector, statistical methods that analyze vector and tensor data are needed.

Directional statistics is conducted on vectors and directions based on observations on compact Riemannian manifolds (Pennec, 2006). Hence, it can encapsulate much more information than scalar metrics about the diffusion. Without the limitation of scalar statistics, one can evaluate dispersion and coherence values among the populations, fit directional model to the data, and perform hypothesis testing for group-based studies.

In the literature, directional statistics have been used to characterize fiber orientation distribution functions, to estimate fiber dispersion quantitatively via fanning and bending fiber geometries throughout the brain (Sotiropoulos et al., 2012; Tariq et al., 2016). In addition, directional statistics have also been utilized to extract bundle-specific metrics from crossing fiber models (Riffert et al., 2014) and fiber tractography (Parker et al., 2003). However, Watson distribution, which has been used in previous directional statistics in group analysis, contains limited parameters (Schwartzman et al., 2005; Hutchinson et al., 2012). Watson distribution is a bimodal probability distribution on a two-dimensional unit sphere S2 in R3 which is symmetrical around mean direction, where each direction and its negative have the same probability. In our previous study (Metin and Gökçay, 2014), it has been shown that Bingham distribution better fits into PDD distributions for white matter tracts and improves the depiction of variability among subjects in anisotropic tensors areas, such as fiber crossings. This is because Bingham distribution is a generalization of Watson distribution: it is bimodal and elliptic around mean direction.

Group analysis methods on DTI or DWI data can be classified into three: (1) region of interest (ROI)-based methods, (2) voxel-based analysis, and (3) fiber tract-based analysis. ROI-based methods are very labor intensive plus error-prone. On the other hand, voxel-wise comparison is open to misalignment of voxels because during registration of individual subject’s data to a common space, topological variabilities may not be thoroughly resolved (Jones and Cercignani, 2010) for each fine structure. The amount of smoothing can greatly affect the final results, but there is no principled way of deciding how much smoothing is “correct” (Jones et al., 2005). For instance, tract-based spatial statistics (TBSS) tackles the alignment and smoothing problem for voxel-wise statistics by combining strengths of VBM-style analyses and tractography-based approaches (Smith et al., 2006). In short, analyses that involve fiber tracts are contingent upon computation of quantitative parameters of interest along the tracts (Goodlett et al., 2009) within diffusion tensor images. The properties of the fiber tract can be scalar values derived from tensors such as MD, FA, or trace, as well as shape information such as curvature and torsion of the specific tract (Mandl et al., 2010).

In this study, we propose a new tract-based framework using directional information in diffusion tensors to improve statistical group analysis, named as track profiling and directional statistics (TPDS). For this purpose, we have (1) generated a new data structure called tract profile by clustering fibers across subjects and (2) developed a method based on directional statistics to compare white matter (WM) differences of different groups across each tract profile. Overall, this new DTI group analysis method is called TPDS.

In order to demonstrate the superiority of the proposed framework, we compared the tract profiling method with two widely used techniques: TBSS (Smith et al., 2006) and voxel-based analysis (VBA) (Hecke et al., 2009). Furthermore, we ran a third comparison with the network-based statistic (NBS) toolbox (Zalesky et al., 2010) which utilizes nonparametric statistical testing to identify the components of an N × N undirected connectivity matrix that differ significantly between two distinct populations.

As a proof of concept, we demonstrated the strength of TPDS in the identification of differences of structural connectivity in major depressive disorder in a small data set (n = 30). Although depression has traditionally been viewed as an affective disorder, the last few decades of research have shown that MDD is also associated with considerable disturbances in cognitive functioning, including executive functions, attention, memory, and psychomotor speed (Castaneda et al., 2008; McClintock et al., 2010). In MDD, multidimensional, systems-level differences are reported in discrete, but functionally integrated pathways (Mayberg, 2003). Therefore, differences in MDD can be expected to cover a wide range of WM tracts. So far, especially white matter disturbances and connectivity differences have been analyzed using DTI-based analysis in MDD (Seminowicz et al., 2004; Zou et al., 2008; Cullen et al., 2010; Kieseppä et al., 2010; McClintock et al., 2010; Helm et al., 2018). Most of these studies state that loss of integrity occurs in the WM fiber tracts of the frontal, temporal, and cingulate cortex of MDD patients. White matter integrity can be described as biophysical white matter changes as a result of microstructural characteristic in both intra- and extra-axonal environments of WM such as axonal water fraction (AWF), intra-axonal diffusivity, and extra-axonal axial and radial diffusivities. More specifically, reported abnormalities in the connectivity of the DLPFC and ACC circuits (Helm et al., 2018), as well as subcortical regions, complement other findings specified in affective disorders (Sexton et al., 2009).

Materials and Methods

Data Acquisition

In order to demonstrate the benefits of TPDS, we used T1-weighted, T2-weighted, and DTI MR data obtained from healthy subjects and patients with MDD.

Subjects

The control group consisted of 14 healthy subjects (8 female and 6 male) with age 31.71 ± 7.62, who had no history of neurological disease and also are not taking any medication. The depression group consisted of medication-naïve 16 subjects (8 female and 8 male, age: 31.12 ± 8.95)¹. The data was collected as part of a local institutional project funded by METU (BAP-07-04-2012). Project management and subject recruitment were handled by a larger project² for which the results will be published elsewhere.

MRI Parameters

Whole-brain MRI scans were collected using the Siemens MAGNETOM 3 T scanner situated at the Bilkent University UMRAM center. T1-weighted [repetition (TR): 2,500 ms, echo time (TE): 3 ms, inversion time (TI): 1,000 ms, flip angle (FA): 8°, sagittal plane 1 mm isotropic resolution], T2-weighted (TR: 5,900 ms, TE: 108 ms, FA: 120°, spacing: 2.2, slice thickness 2 mm), and DWII scans (TR: 8,270 ms, TE: 83 ms, FA: 90°, spacing: 2.2, seven images with b-factor = 0 s/mm², 45 directions b-factor = 700 s/mm²) are collected from the participants in a single session.

Data Processing

Pre Processing

We have implemented a fully automated pipeline to perform preprocessing as illustrated in Figure 1. The overall pipeline has been designed using the Connectome Mapper (Daducci et al., 2012). At the individual subject level, preprocessing steps are performed using several software toolkits. The first step is intrasubject registration of T1, T2, and DWI images using FSL’s FLIRT as described in Jenkinson and Smith (2001) and Jenkinson et al. (2012). The registration is first done between the T2-weighted image and DWI B0 images, and then the high-resolution T1-weighted image is registered to the T2-weighted image. To eliminate the problem of transforming diffusion tensors, all of the images are registered to the DWI B0 image. This way, all image operations are performed on the diffusion image. For segmentation and parcellation of ROIs, FreeSurfer (Fischl et al., 2002) has been used. These steps transform the subject’s MRI to uniform space and segment white and gray matter as well as cortical and subcortical structures based on the underlying atlas. The parcellation algorithm (Fischl et al., 2004) reveals 83 distinct cortical and subcortical structures of the brain using the Desikan–Killiany atlas (Desikan et al., 2006). All of these steps constitute the top row of Figure 1.

FIGURE 1

Figure 1. Components of the pre-processing pipeline before TPDS is performed.

DTI processing begins with motion and eddy current artifact correction in FSL. Tensor estimation is done by Diffusion Toolkit (DTK) (Wang et al., 2007). For tractography (Parker et al., 2003; Cook et al., 2005), streamline fiber-tracking algorithm in Camino has been used. Each voxel in the parcellated image is selected as seeds. Eighty-three distinct cortical and subcortical areas are masked, and the generated binary image is used as the seed file of the algorithm for particular ROIs. For tracking, the fourth-order Runge–Kutta method has been chosen to propagate the tracks using a constant step size. Nearest-neighbor interpolation is applied around local voxel data. A minimum length criterion, 10 mm, is enforced to eliminate premature tract termination due to low SNR and low pathway anisotropy (Behrman-Lay et al., 2015). Each fiber bundle is pruned so that it only contains fibers connecting relevant regions. The number of streamlines depends on the size of the ROI. No additional elimination technique has been applied other than minimum length. These steps are illustrated in the second row of Figure 1.

Using the Connectome Mapper (Daducci et al., 2012), a connection matrix is generated to calculate the connectivity of the areas via the fiber tracts obtained in the first and second rows of Figure 1. After this step, the fiber tracts that connect corresponding brain areas will be bundled to construct relevant fiber bundles. In order to perform group analysis, one last step is necessary: the corresponding bundles of all subjects must be aligned. Therefore, both control and patient images are registered to the ICBM DTI-81 atlas using affine registration. The transformation obtained during this registration is applied to the fiber bundles as seen in the last row of Figure 1.

Tract Profiling

Tract profiles are cross sections of the fiber tracts that connect the ROIs specified by the connection matrix generated in preprocessing. For the connections in each ROI pair, a fiber bundle is formed based on the intersections of cross-sectional areas of all subjects’ DWI. Then, the medial line of the fiber bundle is computed. Finally, a cross-sectional profile is generated along the medial line so that the distribution of PDDs along each cross section is aggregated separately for each subject group.

Overlapping Fiber Calculation

Overlapping fibers/voxels are calculated across all of the subjects. This is done for each fiber bundle by calculating its maximum overlap. During this process, some specific bundles might be left out as outliers. In Figure 2, the overlapping fiber bundle is shown between the two ROIs: thalamus (green) and rostral anterior cingulate (purple). The bundles shown with the yellow, cyan, green, and red colors are marked as outliers and left out of the overlapping area.

FIGURE 2

Figure 2. Tract profiling: Generation of an overlapping tract bundle between two ROIs (shown by green and blue) for all of the subjects regardless of the groups.

The voxel image can be represented as image 𝒫 where 𝒫 = (Z³, m, n, B) (Kong and Rosenfeld, 1989). Each element in Z³ is called a point of 𝒫 and each point in B⊆Z³ is called a black point and assigned 1. Each point in Z³\B is called a white point and assigned 0. m holds black points and n holds white points.

In order to be used in multisubject analysis, adaptation of this definition can be made as follows. For given ROI pairs (i,j), let 𝒫₀, 𝒫₁, …, 𝒫_K be a set where K is the number of subjects, and 𝒫_k is the fiber bundle image from subject k. A point in P is assigned as black point if and only if it is also black point for all sets in 𝒫₀, P₁, …, 𝒫_K for a given ROI(i,j).

Medial Line Generation

The skeleton of the overlapping bundles is calculated. The curve skeleton is a one-dimensional set which runs through the center of the overlapping bundles in such a way that it preserves the topological properties of the overlapping area. Connectivity conditions are defined as follows. The sequence of points (x₀, x₁, …, x_n) is a j-path of length n ≥ 0 from the point x₀ to point x_n in a nonempty set of points X if each point of the sequence is in X and x_i is j-adjacent to x_{i – 1} for each 1 ≤ i ≤ n. The adjacency can be defined as N_j(p) the set of points j-adjacent, to the point p, where j = 6, 18, 26. Connectivity can be defined as j-connected if there is a j-path between them in X.

In order to construct the aforementioned skeleton, first of all, curve thinning (Blum, 1967; Kong and Rosenfeld, 1989) is used on P. The medial line of the fiber bundles was generated as depicted in Palágyi et al. (2001). As such, in each iteration, border points of P were deleted until no more deletion was possible. The algorithm is implemented as sequential iterations where each step checks for six subroutines for each of the six-directions that are immediate neighbors of a black point in P. In each iteration, border points are deleted upon satisfying a condition called simple point condition². In this way, the object is shrunk uniformly in each direction. The operation is continued until no more shrinking is possible for each direction. By adding connectivity conditions, the skeleton ends up with the medial line in the near center of the object. In Figure 2, the example medial line for the fiber bundle is shown with dark blue.

Finally, the resulting medial line is smoothed by generating a b-spline representation as follows. In order to generate b-spline representation of the medial line, the voxel coordinates on the medial line are represented as data points {P_k}, k ∈ MedialLine. A b-spline curve that fits the data is parameterized by t ∈ [0,1], where $X\left(t\right)=\sum _{i=0}^n{U}_{i,d}(t)Q_i$, the control points Q_i are unknown quantities that have been evaluated using the least-squares fitting method described below:

For n control points $\widehat{Q}=\left[\begin{array}{c}\hfill Q_0\hfill \\ \hfill Q_1\hfill \\ \hfill \begin{array}{c}\hfill \mathrm{⋮}\hfill \\ \hfill Q_n\hfill \end{array}\hfill \end{array}\right]$, and m sample points $\widehat{P}=\left[\begin{array}{c}\hfill P_0\hfill \\ \hfill P_1\hfill \\ \hfill \begin{array}{c}\hfill \mathrm{⋮}\hfill \\ \hfill P_m\hfill \end{array}\hfill \end{array}\right]$, the least-square error function between the b-spline curve and the sample points is the scalar valued function:

$$E(\widehat{Q)}=\frac{1}{2}\sum _{k=0}^m|\sum _{j=0}^n{U}_{j,d}\left(t_k\right)Q_j-P{|}^2$$

To minimize the error function, E, where it is quadratic in the components of $\widehat{Q}$, it is a graph of a paraboloid, so it has global minimum that can be found when all its first-order derivatives are 0. The first-order partial derivatives can be written as control points, Q_i

$$\frac{\partial E}{\partial Q_i}=\sum _{k=0}^m\left(\sum _{j=0}^n{U}_{j,d}\left(t_k\right)Q_j-P_k\right)U_{j,d}\left(t_k\right)$$

$$\frac{\partial E}{\partial Q_i}=\sum _{k=0}^m\sum _{j=0}^n{U}_{i,d}\left(t_k\right)U_{j,d}\left(t_k\right)Q_j-\sum _{k=0}^m{U}_{i,d}\left(t_k\right)P_k$$

It can be written as $\sum _{k=0}^m\sum _{j=0}^n{a}_{k,i}{a}_{k,j}{Q}_j-\sum _{k=0}^m{a}_{k,i}{P}_k$, where a_k,i=U_k,d(t_k) for 0 ≤ i≤ n, by setting the partial derivatives to zero vector, and it leads to the system of equations:

$$0=\sum _{k=0}^m\sum _{j=0}^n{a}_{k,i}{a}_{k,j}{Q}_j-\sum _{k=0}^m{a}_{k,i}{P}_k=A^{T}A\widehat{Q}-A^T\widehat{P}$$

Where A = [a_rc] is a matrix with m + 1 rows and n + 1 columns.

$$\widehat{Q}={(A^{T}A)}^{-1}{A}^T\widehat{P}=\left[{(A^{T}A)}^{-1}{A}^T\right]\widehat{P}=X\widehat{P}$$

Since A is tridiagonal where it has a contiguous set of upper bands and lower bands, the equation can be solved with the Cholesky decomposition and the vector of control points $\widehat{Q}$ can be found.

Since derivative of spline is 1 less order of yet another b-spline where new control points are defined as $Q_i=\frac{p}{u_{i+1+1}-u_{i+1}}[cpsbreak](P_{i+1}-P_i)$ from the surface tangent, a normal vector has been computed and cross-sectional areas have been extracted.

Calculation of Tract Cross Sections

The skeleton is sliced with 2-mm regular intervals so that cross-sectional areas that are perpendicular to the b-spline are obtained using normal vectors computed from the surface tangents in Figure 3. For each voxel in P that intersects with these cross-sectional areas, PDDs that represent individual subjects are added as tract profiles representing that slice. Hence, for a tract with J slices, there are J tract profiles that contain PDDs which are representative of the subject group. An example tract profile (i.e., a slice with PDDs) from a single subject is shown in Figure 4. The PDDs from the subjects for a specific group are aggregated as follows. At each slice, there are fixed number of voxels, and at each voxel, there can be multiple PDDs, each coming from a different subject, depending on whether the subject’s tract goes through that voxel or not.

FIGURE 3

Figure 3. Tract profiling: Representation of the medial line of the overlapping bundle with b-splines and generation tract profiles.

FIGURE 4

Figure 4. Tract profiling: Illustration of the PDDs from a single subject in a sample tract profile.

Directional Statistics

Statistical analysis is executed exclusively on areas that are defined by tract profiles eliminates voxel-wise comparison. Hence, misalignment problems no longer exist. Hypothesis testing is conducted only at cross-sectional tract profiles that are separated by 2-mm regular intervals. For the set of PDDs embodied in each tract profile j, a parametric directional statistic distribution is fitted. Through such parametrization, the PDDs of all subjects that fit into the tract profile j are projected onto a sphere.

Watson distribution in Figure 5 is bimodal and symmetrical around mean direction. Watson distribution assumes that diametrically opposite points have the same probability. Also, the probability density function of axial distributions process antipodal symmetry [i.e., f(−l,−m,−n) = g(l,m,n)]. The probability distribution of random vectors that belong to the Watson’s family is spherical on a sphere. Directional statistics have been used in the analysis of DTI previously (Schwartzman et al., 2005; Hutchinson et al., 2012), and it has been shown that DTI principal direction analysis using directional statistics can better identify the differences in anatomic structure between populations compared with statistical tests of scalar values such as FA. Both of these studies used Watson distribution to analyze principal directions. On the other hand, Bingham distribution (Figure 6) is bimodal and elliptical (Fisher et al., 1993; Cheng et al., 2014). Bingham distribution is free from symmetrical constrains; hence, it provides more advanced distribution fitting options in comparison with Watson distribution.

FIGURE 5

Figure 5. Watson Distribution.

FIGURE 6

Figure 6. Bingham Distribution.

Watson distribution is defined as follows (Mardia and Jupp, 1999):

    $\mathrm{Watson}\mathrm{Distribution}{W}_p(x;\mathrm{\mu },\mathrm{\kappa })=c_p\left(\mathrm{\kappa }\right)e^{\mathrm{\kappa }{\left({\mathrm{\mu }}^{T}x\right)}^2};$

    $c_p\left(\mathrm{\kappa }\right)={\displaystyle \frac{\text{Γ}\left(\frac{p}{2}\right)}{2{\mathrm{\pi }}^{p/2}M(\frac{1}{2},\frac{p}{2},\mathrm{\kappa })}}$

where x is the unit random vector, μ is the mean vector, is the concentration value, M is Kummer’s confluent hypergeometric function, Γ is a gamma function, and p is the dimension of the distribution. To estimate maximum likelihood of this function, we take logarithm. Hence, the log-likelihood function is

$$l(\mathrm{\mu },\mathrm{\kappa }\pm x_1,\mathrm{\dots },\pm x_n)=\mathrm{\kappa }\sum _{i=1}^n{\left(x_i^T\mathrm{\mu }\right)}^2-n\text{log}M(\frac{1}{2},\frac{p}{2},\mathrm{\kappa })$$

$$=n\left\{\mathrm{\kappa }{\mathrm{\mu }}^T\overline{T}\mathrm{\mu }-\text{log}M(\frac{1}{2},\frac{p}{2},\mathrm{\kappa })\right\}$$

where $\overline{T}$ is the scatter matrix of the given data. Differentiation with respect to κ gives

$$D_p\left(\mathrm{\kappa }\right)={\widehat{\mathrm{\mu }}}^T\overline{T}\widehat{\mathrm{\mu }};\text{for}p=3;=\frac{M(1.5,3.5,\mathrm{\kappa })}{3^{*}M(0.5,1.5,\mathrm{\kappa })}.$$

And to find its maximum likelihood estimate, we need a derivative of D_p(κ) for p = 3

$$D_3^{\prime }=\frac{M(2.5,3.5,\mathrm{\kappa })}{5M(0.5,1.5,\mathrm{\kappa })}-\frac{1}{9}*{\left(\frac{M(1.5,2.5,\mathrm{\kappa })}{M(0.5,1.5,\mathrm{\kappa })}\right)}^2.$$

The Newton–Raphson method can be used tfind maximum values for D_p(κ) and the biggest eigenvalue of scatter matrix, t₁, for a bipolar distribution or t₃ for a girdle distribution.

Bingham distribution is defined as a trivariatnormal distribution on a unit sphere. Different from Watson distribution, it has three orthogonal directions as μ₁, μ₂, μ₃ and concentration values (κ_n)for each orientation vector (Watson and Williams, 1956).

Concentration values define the dispersion of the distribution, where

(1) κ₁=κ₂ = 0 results in a spherical distribution of axes.

(2) κ₁=κ₂≪0 results in a symmetric bipolar distribution.

(3) κ₁ < κ₂≪0 results in an asymmetric bipolar distribution.

(4) κ₁≪κ₂ < 0 results in an asymmetric girdle distribution.

(5) If κ₁≪0andκ₂ = 0, then Watson distribution is obtained.

The probability distribution function of Bingham distribution is defined as follows (Bingham, 1974):

    $\mathrm{Bingham}\mathrm{Distribution}{B}_p(x;K)=c_p\left(K\right)e^{x^{T}Kx};$

    $c_p\left(K\right)={\displaystyle \frac{\text{Γ}\left(\frac{p}{2}\right)}{2{\mathrm{\pi }}^{p/2}F(\frac{1}{2},\frac{p}{2},K)}}$

where x is the unit random vector, K is the 3 × 3 orthogonal orientation matrix with concentration values, F denotes the confluent hypergeometric function of matrix argument, Γ is the gamma function, and p is the dimension of the distribution. For a given random sample ±x₁,…,±x_n, the log-likelihood function can be written as:

$$l(K;\pm x_1,\mathrm{\dots },\pm x_n)=n\left\{\log tr\left(A\overline{T}\right)-\text{log}F(\frac{1}{2},\frac{p}{2},K)\right\}$$

We can write K and $\overline{T}$ in polar form as = UKU^T, $\overline{T}=Vt{V}^T$ with U and V being orthogonal. K = diag(κ₁,…,κ_p) and $t=({\overline{t}}_1,\mathrm{\dots },{\overline{t}}_p)$, where κ₁≥…≥κ_p and ${\overline{t}}_1\ge \mathrm{\dots }\ge {\overline{t}}_p$. As suggested by Bingham himself, the following approximations can be used.

For the bipolar case:

    $d={\overline{t}}_2-{\overline{t}}_3,s={\overline{t}}_1+{\overline{t}}_2,{\mathrm{\kappa }}_0=-D_3^{-1}({\overline{t}}_1)$

    ${\mathrm{\kappa }}_1\approx 0,{\mathrm{\kappa }}_2\approx {\mathrm{\kappa }}_0+\mathrm{\delta },{\mathrm{\kappa }}_3={\mathrm{\kappa }}_0-\mathrm{\delta }$

For the girdle case:

    $d={\overline{t}}_1-{\overline{t}}_2,s={\overline{t}}_1+{\overline{t}}_2,{\mathrm{\kappa }}_0=-D_3^{-1}({\overline{t}}_3)$

    ${\mathrm{\kappa }}_1\approx 0,{\mathrm{\kappa }}_2=-2\mathrm{\delta },{\mathrm{\kappa }}_3={\mathrm{\kappa }}_0-\mathrm{\delta }$

where $\mathrm{\delta }=\frac{2d{\mathrm{\kappa }}_0}{s\left({\mathrm{\kappa }}_0-1.5\right)+1}$

After parametric representation through either Watson or Bingham distribution, two group of subjects can be compared by using an eclipse of confidence defined by the p value. For fitting a single group’s data, the mean direction vector of the group is computed. If it lies inside the eclipse of confidence of the targeted distribution, then the null hypothesis is likely, justifying a reasonable fit to the associated directional distribution. On the other hand, if the confidence ellipse around the mean direction does not overlap for a given confidence level, then the null hypothesis is unlikely, rejecting the fit. For two groups, the case with different means is indicated by separated cones of confidence, which in turn indicates significant differences. On the other hand, overlapping cones of confidence indicate insignificant differences, hence acceptance of the null hypothesis. An example distribution is provided in Figure 7, for two groups of subjects, for representation with the Bingham distribution.

FIGURE 7

Figure 7. PDD projections modeled by the Bingham Distribution. (A) Separated (Left) versus overlapping (Right) vector projections of PDDs on unit sphere for two different subject groups shown with blue and red. (B) Statistically significant (Left) versus insignificant (Right) differencs between populations.

Details of the eclipse of confidence can be given as follows. The maximum likelihood estimates of concentration parameters κ₁,κ₂ can be obtained from maximizing the log-likelihood function, where w_n are the eigenvalues of the principal eigenvector of the orientation matrix:

$$F=-N\text{log}\left(4\mathrm{\pi }\right)-N\log d({\mathrm{\kappa }}_1,{\mathrm{\kappa }}_2)+{\mathrm{\kappa }}_1{w}_1{\mathrm{\kappa }}_2{w}_2$$

Maximum likelihood estimators of κ₁,κ₂ in the Bingham distribution for given eigenvalues w₁,w₂ can be estimated as calculated by Mardia and Zemroch (1977).

The confidence ellipse around the mean direction within the specified percentage (%) of the estimated concentration values of distribution as

$$e_{\%}^{mn}=\sqrt{\left[\frac{X_{\%}^2}{2N({\mathrm{\Delta }}_{mn})}\right]}\text{for}{\mathrm{\Delta }}_{mn}=({\mathrm{\kappa }}_m-{\mathrm{\kappa }}_n)(w_m-w_n)\text{and}{X}_{\%}^2$$

is the chi-squared value for two degrees of freedom and % is p value for confidence interval.

For p = 0.01 and having κ₃ = 0 (Fisher et al., 1993) ends up with the semi-axes of the confidence eclipse about the mean direction associated with w₃ as below:

$$e_{32}=-1.517\frac{1}{k_{2}N(w_3-w_{2)}}\text{and}{e}_{31}=-1.517\frac{1}{k_{1}N(w_3-w_{1)}}$$

Performance Analysis

We have conducted two performance tests to analyze the effectiveness of the proposed method. First, we have analyzed the models generated by TPDS in comparison with VBA and TBSS. Hereby, we have adapted directional statistics to TPDS, TBSS, and VBA to compare their overall efficiency in representing vector-based statistical models. This test aimed to show the efficiency of tract profiling over voxel-based and skeleton-based analysis. Second, we have applied the full TPDS algorithm to the two subject populations (i.e., MDD versus healthy controls) and compared the results with NBS. This test aimed to show the efficiency of combining tract profiles with directional statistics over conventional methods. In this test, the effects of fiber length in estimating group differences were also evaluated.

Analysis of the Strengths of VBA, TBSS, and TPDS in Tract Modeling

In this test, we used a single group (i.e., healthy subjects). The statistics were derived using three different methods, VBA, TBSS, and TPDS, only on white matter areas—not using GM ROIs. As seen in Figure 8, the white matter areas that have been segmented using FreeSurfer are mapped to ICBM DTI-81 atlas (Mori et al., 2008) to allow for intersubject data aggregation. For VBA analysis, the atlas-based white matter areas are overlayed for all subjects for further processing. For TPDS analysis, tract profiles are generated from the atlas mappings of all subjects. In TBSS, before performing atlas mapping, skelotonized areas are generated from individual subject tracts. The rest of the data processing pipeline is the same for all three methods. At the first step, for each WM ROI, based on which method is used for defining the tract, PDDs are generated. Then these PDDs are parametrically modeled by two separate directional distributions, namely Bingham and Watson. Finally, in the last step, several PDDs are generated to represent the entire group using the newly developed parametrical models, and goodness of fit is computed to evaluate how good the chosen model is.

FIGURE 8

Figure 8. Comparison of VBA/TBSS with TPDS data processing pipeline.

PDD Generation

For each subject, primary diffusion directions are extracted for each voxel inside the given WM area using the primary eigenvector of the diffusion tensor. The WM area differs based on the chosen representation. In VBA, the WM area is extracted based on segmentation of the specific WM ROI. In TBSS, it is based on the skeleton of the tract in the WM ROI. In TPDS, it is embodied within each tract profile that composes the entire tract in the WM ROI. Aggregated data from all subjects compose the data to be fitted for each WM area.

Distribution Fitting

Watson and Bingham distributions were fitted to model each tract using the maximum likelihood method. For each tract, the parameters of the theoretical model were estimated from the pdf at hand. Then this theoretical probability density function was evaluated iteratively using synthetic random vector data for a total of 700 vectors that were almost uniformly distributed along a sphere. Finally, the difference between the estimated pdf and the random pdf is tested for null hypothesis.

Goodness of Fit Testing

Pearson’s chi-square tests have been used for goodness of fit tests to evaluate whether the observed frequency distribution differs from the theoretical distribution. Comparison of distributions is done using ANOVA and the chi-square test statistics was also used for each ROI.

In order to apply Pearson’s chi-square tests to check whether the observed frequency distribution differs from a theoretical distribution, the following steps are applied on the original data and synthetic random vector data and the respective models.

(1) For Watson distribution, the sample mean direction, $\overline{R}$, has been evaluated as a regular vector sum of the vectors under a population of vectors. The mean direction is a unit vector that is in the same direction with R: $\overline{x}=\frac{{\sum }_i{x}_i}{R}$, $\overline{y}=\frac{{\sum }_i{y}_i}{R}$, $\overline{z}=\frac{{\sum }_i{z}_i}{}R$.

(2) For Bingham distribution, the axis of moment of inertia of sample, $\overline{t}$, has been evaluated using the scatter matrix of distribution S. For the bipolar case, it is the biggest eigenvector, and for girdle case, it is the smallest eigenvector.

(3) The transformations $\overline{\mathrm{\theta }},\overline{\mathrm{\varphi }}$ have been evaluated in order to shift either $\overline{R}$ or $\overline{t}$. to positive z-axis.

(4) The transformation has been applied to original and synthetic data.

(5) The angle θ has been calculated as the angle between the positive x-axis and the projected vector on the x–y plane: 0 < θ < 2π.

(6) The observed frequencies and the excted frequencies were θ₁ < θ < θ₂, where the number frequency bins is 50 and θ₂ – θ₁ ≈7.2^o.

Analysis of the Group Difference Maps Generated by NBS and TPDS

In this part, the proposed framework will be applied to test for differences of fiber tract profiles between MDD patients and control subjects. Based on the same fiber tracts and connectivity matrix for healthy volunteers, comparisons will be made with the results of the network-based statistics. For this purpose, we used the 83 × 83 connectivity matrix generated at the end of data preprocessing by the Connectome Mapper (Figure 1).

In NBS, for each group, each pairwise association (i,j) between ROI i and ROI j is treated separately. First, Fisher’s r-to-z transform has been applied to ensure normality. Then, the test statistic of interest—which is the normalized number of fiber bundles—is compared between the groups using t-statistic. In order to correct for multiple comparisons, permutation testing was used to select the p value controlled for the FWE for each connected component. For each permutation, the same threshold is applied to define a set of suprathreshold links of connected components. Suprathreshold and the number of permutations were set according to the default parameter settings of NBS with corrected p < 0.005.

In TPDS, the following procedure is repeated for each possible connection between distinct ROI pairs (i.e., 83 × 83 times divided by 2). Tract profiles between each ROI i and ROI j are extracted for the healthy and MDD groups. Then for each slice in the tract profiles, significance is tested with a threshold value of p < 0.005. If there are n contiguous slices that satisfy this, it is indicated that the connection between ROIs i and j is significantly different between the control and patient groups. It is possible that there are multiple clusters of n contiguous slices that satisfy this condition. In order to reflect this information, we prepared a new 83 × 83 connectivity matrix, which contained the number of significantly different clusters between the two groups that are compared. Therefore, the difference map that is achieved through TPDS reflects a weighted graph, weight being the number of significantly different clusters between the two groups for that particular i to j connection. The more the number of significantly different n contiguous slices, the more the weight of the difference map.

Selection of n must be done according to a criterion related to the plausible tract lengths. In order to eliminate premature tract termination that result from low SNR and low pathway anisotropy (Behrman-Lay et al., 2015), 10 mm is the shortest tract length to be considered. Since DTI image has 2.2 mm spacing, choosing n as 4 satisfies this constraint. In other words, at least four consecutive cross-sectional areas must be found within a fiber bundle where the PDD of each cross-sectional area belongs to significantly different Bingham distributions for the control and MDD groups.

Results

The results of the performance tests that we performed to investigate the effectiveness of TPDS are as follows.

Comparison of VBA and TBSS With TPDS Using Directional Statistics

As can be seen in Table 1, among VBA, TBSS, and TPDS, the best fitted distribution is more representative in TPDS because the goodness of fit scores are better according to p values. In addition, based on the results of TPDS, the Bingham distribution is reported to be more favorable than the Watson distribution because only 2 out of 48 white matter tracts are represented better with Watson. Obviously, it is evident that TPDS is a better alternative to represent tracts in comparison with VBA and TBSS, because it favors a more parametrical fit to the entire set of fiber tracts.

TABLE 1

Table 1. Comparison of voxel-based analysis (VBA) and tract-based spatial statistics (TBSS) with tract profiling and directional statistics (TPDS) (VBA and TBSS have been adapted to run directional statistics).

A close inspection of Table 1 reveals that in terms of representing a given WM tract parametrically, TBSS is superior to VBA, and TPDS is superior to TBSS. It is evident that VBA contains more noise than TBSS and TPDS, because it contains the entire WM area from all subjects. Due to high noise, VBA fails to represent some of the tracts parametrically. On the other hand, TBSS is better than VBA, because it removes the areas—hence the noise associated in these parts—that lie outside the fiber bundles which constitute the skeleton. However, TBSS is not better than TPDS, because it smooths out the tracts while forming the skeleton and loses specificity. Overall, the tract profiles computed in TPDS are selective in choosing representative samples of the DWIs that are more informative, because outliers are removed while computing the medial line. Since the data points all belong to the same tract and on the same cross section over the medial line, very similar diffusion properties are expected for each analysis point. This tends to eliminate all negative effects of misalignment of images and partial volume effect. Due to this property, the computational effectiveness of TPDS is higher than other methods, because the model can be decided with much less number of data points.

The advantage of the Bingham distribution might be explained through the ease of fitting a girdle distribution in comparison with fitting a homogeneous mean direction distribution. The girdle distribution allows for more parameters; hence, it makes the development of a more general model possible. Furthermore, the computational accuracy of the Bingham distribution is better because the tracts represented with this distribution fit to the PDD of the actual tracts with a smaller p value.

Comparison of the Group Differences in Connectivity Maps Using Network-Based Statistics and TPDS

In NBS, with corrected p < 0.005, seven regions and eight connections have been observed to contain lower FA in MDD. Particularly, the connections in the right hemisphere and between the superior frontal cortex and rostral/caudal components of the anterior cingulate cortex, caudate, and inferior parietal cortices had lower FA in MDD. These connections are shown in Figure 9 as green lines.

FIGURE 9

Figure 9. Map of ROIS with statistically different connectivity between control and patient groups. Green lines represent the common connections that are found different between the groups using NBS. Red lines represent the significantly different connections detected by the directional statistics using tract profiling.

In TPDS, significantly different connections between the healthy and MDD groups are seen in Figure 9 as red lines. The thickness of the lines reflects the weights or in other words the number of cross-sectional areas above the threshold n (e.g., A weight value of 1 indicates that there exists only one slice cluster with significantly different n contiguous tract profiles, whereas a weight value of 6 indicates that there exist 6 disjoint clusters of n contiguous tract profiles that are significantly different). The right hemisphere differences reported by NBS, namely the frontal (superior frontal and rostral middle frontal) and medial (caudal and rostral anterior cingulate), are also detected by our method. But, additionally, TPDS revealed differences between the healthy and MDD populations in limbic, temporal cortex, occipital cortex, and hippocampal connections, as well as a few left hemisphere areas such as the amygdala, hippocampus, and thalamus.

The strength of the tract profile structure lies in the reduction of the misalignment problem. Furthermore, observations of the directional changes become more specific because contributions of the local changes can be reported along the tract not by the contribution of isolated voxels but by several slices across the two ROIs. Therefore, the proposed directional statistics comparison is expected to be a superior differentiator for especially long tracts.

In order to verify this, the following analysis has been done. Using TPDS, for each tract connecting 83 different regions, the z-score of each length is plotted against the z-score of the number of significantly different profile slices. For this purpose, the maximum overlapping shape (skeleton) is used. When regression lines are fitted to investigate the relationship with tract length and the number of different clusters, it is seen that the likelihood of finding clusters of voxels that differ in long tracts increase with respect to path length. This has been also tested using a linear regression model, where it has been found that the z-score of tract length significantly correlated with the z-score of the number of significantly different profile slices (p < 0.05, adjusted R²: 0.00162) as seen in Figure 10. Although the effect size is small, we can indicate that TPDS is a powerful method to find differences in two populations, especially as the tract lengths get longer.

FIGURE 10

Figure 10. Scatter diagram of z-score of tract lengths versus significantly different clusters.

Discussion

In this study, we proposed a novel framework for WM fiber connectivity analysis using TPDS. In contrast with other group studies (Goodlett et al., 2009) that are based on FA values, directional statistics deals with compact Riemannian manifolds, which allow observations regarding local diversities of principal diffusion directions of voxels in different groups of subjects.

Comparison of TPDS With Other Techniques Used in Analysis of Groups of DWI

Our pipeline implementation can be regarded as quantitative tractography. We analyze diffusion properties on the exact tracts and derive the statistics over sample points taking neighborhood cells into consideration. A similar method has been offered by Corouge et al. (2006) where diffusion properties along the fiber tracts, called fiber property profiles, are extracted. In that study, fiber tract parameterization was based on arc length parameter, starting from each fiber’s intersection with an “origin” plane. Goodlett et al. also proposed a similar tract profiling approach, where diffusion properties are calculated along the tract for each fiber bundle (Goodlett et al., 2009). Our method introduces three main improvements to these quantitative tractography methods. First, we are not just limiting the method with known anatomical fiber bundles but can derive statistics from any pair of connected gray matter areas. Second, we have introduced skeletonization and pruning to allow for applying statistics only within common areas across the groups. Third, we introduced vector analysis using directional statistics over scalar analyses such as FA, MD, etc.

There exist other methods which use directional statistics in DTI (Schwartzman et al., 2005; Hutchinson et al., 2012). However, these methods analyze group differences based on ROIs, not fiber tracts, ignoring the underlying connectivity. We have devised the tract profiling algorithm to operate on relevant voxels among the fibers that connect each ROI obtained from fully automatic brain segmentation and parcellation. Local registration errors are reduced after calculating cross-sectional area of the fibers and finding medial lines (i.e., profiles) to continue tract analysis. Afterwards, Bingham distribution, which is the most general form of directional distribution, is used for tract-based directional analysis, ensuring minimum parametric assumptions about the dataset. To the best of our knowledge, this approach has not been implemented in group analysis of DTI before.

Neurite orientation dispersion and density imaging (NODDI) is a novel neurite imaging and analysis framework and provides sensible neurite density and orientation dispersion estimates. Unlike FA, NODDI analyzes density and orientation dispersion separately. NODDI uses orientation distribution function (ODF), defined as Watson distribution which constrains the dispersion about the dominant orientation (Zhang et al., 2012). However, Bingham distribution fits better to diffusion properties, in comparison with Watson. Bingham-NODDI extends the NODDI method by generalizing it with Bingham distribution to cover anisotropic orientation dispersions of neurites (Tariq et al., 2016). Similarly, in our study, we found that modeling ODF using Bingham distribution explains the data better regardless of the tract identification, be it through VBA, TBSS, or our method, TPDS. A major difference between our approach and NODDI is in the estimation of the dispersion. The modeling we used to implement the Bingham distribution estimates dispersion in the vicinity of the dominant orientation, separately for the primary and secondary dispersion orientations. This eliminates the key limitation of NODDI, failing to model complex neurite configurations such as those arising from fanning and bending axons. On another front, just like ours, orientation dispersion (ODI) generated by the NODDI method can be also used with the TBSS method instead of FA metric (Timmers et al., 2016; Taoka et al., 2020). In this aspect, the main difference between our method and NODDI is in extending fiber dispersion along the tracts that connect two ROIs. By allowing such extension, our method enables using fiber dispersions as track characteristics and analyzing disease-related effects on connectivity of the tracks rather than the voxel.

We have demonstrated that in addition to scalar diffusibility changes, analyzing principal diffusion directions along a tract detects local changes better than scalar values. The strength of the directional statistics-based analysis we proposed lies in its applicability to TBSS and VBA as well; it is not limited to tract profiles.

Voxel-based analysis needs to register the subject’s images to a common coordinate frame. However, the fiber tracts do not accurately align during this process due to variation in tract size and shape. Especially, long-range fiber tracts contain more shape variation across subjects (Wassermann et al., 2011), so they are more prone to such misalignment. This problem is still valid for TBSS because even the voxel skeletons do not ensure that all relevant voxels correspond to the same tract (de Groot et al., 2013).

In directional statistics, the misalignment problem though the tract becomes more critical compared with scalar statistics like FA. As seen from the results of the first set of performance tests, tract profiles are superior structures for resolving the shape differences in comparison with VBA and TBSS, because tract profiles are better in terms of fitting a model to PDD vectors. We investigated the goodness of fit characteristics of VBA and TBSS, respectively, on all WM areas and on skeletonized WM areas using directional statistics. We found that several tracts in VBA and TBSS are rejected to fit to the most general Bingham distribution which contains minimum assumptions about the data. In comparison when tract profiling is used, most tracts could be fit parametrically, except a few. A parametrical model is advantageous in data processing, since it facilitates population-based comparisons.

The aforementioned tests also show how directional statistics can be adapted to the widely used analysis methods such as TBSS or VBA. Instead of FA values, PDD vectors can be used over each voxel within the skeleton. FA metric uses eigenvalues of the underlying diffusion characteristics of the voxel and defines only the amount of diffusion asymmetry where PDD uses the first eigenvector of the diffusion characteristic. The FA metric is sensitive to the underlying fiber architecture and correlates with PDD changes in disease conditions. However, FA does not have direction property. Different orientations might result in the same FA value simply because orientational changes of the diffusion property of the voxel might not end with FA changes, when there is a difference in eigenvector orientation but not its value. So, the FA metric is not as sensitive as PDD in detecting diffusion characteristic differences along the fiber track. As can be seen in Table 1, Bingham distribution fits better to describe the differences in the majority of white matter tracks. Further studies should be conducted to ease adaptation of directional statistics to TBSS skeletons and also to resolve issues related to the multiple comparison problem.

PDD analysis using directional statistics is not a summary statistics of each track but a measurement of diffusional properties of the fiber bundle connecting a pair of ROIs. The statistics of each voxel along the fiber track are summarized by many points using directional statistics along the fiber bundle. Fiber bundle skeletonization and normalization of PDD over tract cross sections allows for error correction and noise cancelation that might arise from tractography artifacts or misalignment. This should also be valid for trajectory changes of tracts under disease-related conditions, as long as a prominent disfiguration or an abnormal morphological change caused by a tumor deviation does not severely divert the alignment of the fiber bundles. In such a case, a lot of false positives may affect the model along the fiber bundles, hindering the correct estimation of PDDs along the actual but diverted tract.

During the second set of performance tests, the results of TPDS and NBS are compared to see whether these methods report the differences between the healthy and MDD populations consistently. We found that most of the right hemisphere-specific connectivity differences reported earlier in MDD have been detected by both of these approaches. The results are much more consistent among the shorter tracts such as frontal connections of the anterior cingulate. However, TPDS reveals additional connectivity differences mainly among longer tracts such as those between temporal and occipital cortex as well as those that contain areas with low FA values and higher crossing fibers such as the amygdala, hippocampus, and thalamus. Another strength of TPDS is due to its revelation about weights, which indicate the amount of difference between the subject populations along the tracts.

These findings are also consistent with MDD models proposed by Drevets et al. (2008) and Mayberg (2003) where MDD can be defined through a limbic–cortical dysregulation model. In this model, the limbic–thalamo–cortical (LTC) circuits, involving the amygdala, thalamus, and orbital and medial PFC, and the limbic–cortical–striatal–pallidal–thalamic (LCSPT) circuits are mainly the affected areas. These connections are found be affected both using NBS and TPDS. Additionally, TPDS revealed temporal, parietal, and occipital cortex connections that are different in MDD. Mainly the differences on inferior fronto-occipital tracts can be also supported by other DTI studies that report significantly decreased FA values among MDD patients (Cheng et al., 2014).

The ROIs reported to have statistically significant connectivity differences in MDD versus healthy participants are consistent with the two well-known lateralization models of emotion. According to the right hemisphere hypothesis, the right hemisphere is dominant in processing emotions (Alves et al., 2008). On the other hand, the valence hypothesis posits that the left hemisphere processes positive (or approach-related) information, but the right hemisphere processes negative (or avoidance-related) information (Alves et al., 2008). Within the context of MDD, hypoactivity in the left hemisphere fronto-striatal loops indicates the lack of downregulation of the subcortical areas. In Figure 9, TPDS—but not NBS—reported differences in the connectivity of the left hemisphere, amygdala, thalamus, and hippocampus, consistent with the valence hypothesis. However, the abundant presentation of right hemisphere ROIs in Figure 9 supports the right hemisphere hypothesis indicating that the connectivity within the right hemisphere may be a biomarker for MDD. TPDS revealed a larger right hemisphere network which was sidestepped by NBS. This network is predominantly composed of the basal temporal lobe structures as well as occipital ROIs such as precuneus and pericalcarine. The difference in the temporal and parietal functionality in MDD is reported less in comparison with those in front striatal structures; however, there is a growing body of literature that focuses on the hypoactivity of the right hemisphere temporal areas in MDD (Bruder et al., 2017). The detection of such ROIs by TPDS is supportive of these studies reported in Bruder et al. (2017). Finally, several rsfMRI biomarkers of MDD are reported in Drysdale et al. (2017). After clustering these biomarkers through machine learning techniques, four different subtypes of MDD can be derived, based on four different clusters of ROIs. Unfortunately, the temporal areas of the brain are excluded in this study, due to a lack of data collection from several participating research sites. However, the ROI network reported by both NBS and TPDS in Figure 9 is also reported in Drysdale et al. (2017), verifying our results in a much larger sample size.

In their meta-analysis of over 231 patients with MDD and 261 comparison participants, Yi Liao et al. found four consistent locations of decreased FA: white matter in the right frontal lobe, right fusiform gyrus, left frontal lobe, and right occipital lobe. Mainly, the right inferior longitudinal fasciculus, right inferior fronto-occipital fasciculus, and right posterior thalamic radiation were involved in such changes (Liao et al., 2013). This covers most of the connection pairs we have found in Figure 9, especially the right fusiform gyrus connections with R. Inferior temporal, parahipppocampal, and temporal gray matter are important because the NBS method failed to reveal all of these areas consistent with the meta-analysis.

In another meta-analysis (Wen et al., 2014), reduced FA is reported in the DLPFC and UF of patients with late-life depression (Wen et al., 2014). Those regions are part of frontostriatal and limbic networks consistent with our findings in Figure 9. This is also consistent with NBS analysis, especially the connections colored in green.

Another recent meta-analysis study has analyzed WM anisotropy and diffusivity in 1,305 MDD patients and 1,602 healthy controls (age range 12–88 years) from 20 samples worldwide (van Velzen et al., 2020). On adults, lower FA was observed in 16 of the 25 ROIs. The largest changes have been found mainly in the anterior corona radiata (ACR), corona radiata (CR), corpus callosum (CC), genu of the corpus callosum (GCC), body of the corpus callosum (BCC), and anterior limb of the internal capsule (ALIC). Significantly lower FA was also observed in the superior fronto-occipital fasciculus (SFO), sagittal stratum (SS), internal capsule (IC), posterior corona radiata (PCR), superior corona radiata (SCR), inferior fronto-occipital fasciculus (IFO), fornix/stria terminalis (FXST), external capsule (EC), and cingulate gyrus of the cingulum bundle (CGC). It is quite important to note that most of these regions are better fitted by TPDS in comparison with TBSS and VBA as revealed by our first test on these methods. The superior fronto-occipital fasciculus (left–right), sagittal stratum (left–right), superior corona radiata (left–right), posterior corona radiata (left–right), superior fronto-occipital fasciculus (left–right), inferior fronto-occipital fasciculus (left–right), external capsule (left–right), fornix (cres)/stria terminalis (left–right), and cingulum (left–right) are all better modeled using TPDS. This is also true for the anterior and superior corona radiata where only the right anterior corona radiata is modeled better with TBSS skeleton. The parts of the corpus callosum are on the other hand fitted better as the genu of the corpus callosum for TPDS, the body of the corpus callosum for VBA, and the splenium of the corpus callosum for TBSS. Overall, the benefit of TPDS is demonstrated in two different ways: 1. By fitting the underlying structural connections to an analytical model in a better way 2. By capturing wider network connectivity differences especially along longer tracts.

Conclusion

To conclude, we have shown that by analyzing PDDs using directional statistics, more insight is gained about fiber tracts regarding differences between populations. While other connectivity-based analysis methods may disregard the differences between longer fibers, TPDS becomes more robust as fiber tract length increases. In areas with low FA values, the distribution of PDDs among the fiber tracts can differentiate connectivity-based dysfunctions better, due to the power of directional statistics. The directional statistics analysis suggested here can also be applied by augmenting the existing methods, namely TBSS and VBA. Such an addition to the existing methods is valuable because it opens up the possibility to use parametric fitting along with directional statistics. The proposed method could be extended considering second and third directions of the diffusion tensor. In a future study, this can be modeled separately, fitting different distribution models for each direction and analyzing the statistical changes of each direction in disease conditions.

When we implemented TPDS in two subject populations, one healthy and the other with MDD, we found several WM tract differences that are not reported in other methods such as NBS and TBSS. It is imperative to use TPDS on other subject populations and with more subjects to justify its strength in comparison with other methods that perform WM tract-based group analysis.

Data Availability Statement

The data analyzed in this study is subject to the following licenses/restrictions: The data was collected from medication naive patients for a larger project. The entire project data will be released when the main research is published. Requests to access these datasets should be directed to DG, ZGlkZW1nb2tjYXlAZ21haWwuY29t.

Ethics Statement

The studies involving human participants were reviewed and approved by Ankara university faculty of medicine clinical research ethics committee. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

MM designed the model, computational framework, and carried out the implementation. DG directed the study on the MDD population and designed the data collection pipeline. Both authors analysed the data, discussed the results, and contributed to the final manuscript.

Funding

The authors appreciate METU (Middle East Technical University) and TUBITAK (Turkish National Science Foundation) for partially funding this project (project numbers BAP-07-04-2012 and 109E081).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to express their gratitude to Prof. Bora Baskak from Ankara University, Ergin Atalay from Bilkent University, and Zeynep Başgöze for diagnosis of the MDD patients, data acquisition at the UMRAM center, and data collection. Special thanks to Gözde Ünal and Vilda Purutçuoglu for their comments and suggestions during the doctoral studies of MM through which the current study was accomplished.

Footnotes

1. ^ TUBITAK 1001, no: 109E081, Ethical board approval: Ankara University Medical College.
2. ^ A black point is called a border point if it is six-adjacent to at least one white point. A black point is called an end point if it has exactly one black 26-neighbor. Black point p is simple in (Z³, 26, 6,B) if and only if all the following conditions hold (Palágyi et al., 2001):
   (1) The set N₂₆(p)∩(B\{p}) is not empty (p is not an isolated point).
   (2) The set N₂₆(p)∩(B\{p}). is 26-ected
   (3) The set (Z³\B)∩N₆(p) is not empty (p is a border point).
   (4) The set (Z³\B)∩N₆(p) is six-connected in the set (Z³\B)∩N₁₈(p).

References

Alves, N. T., Fukusima, S. S., and Aznar-Casanova, J. A. (2008). Models of brain asymmetry in emotional processing. Psychol. Neurosci. 1, 63–66. doi: 10.3922/j.psns.2008.1.010

CrossRef Full Text | Google Scholar

Basser, P. J. (1995). Inferring microstructural features and the physiological state of tissues from diffusion-weighted images. NMR Biomed. 8, 333–344.

Google Scholar

Basser, P. J., and Pierpaoli, C. (2011). Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. J. Magn. Reson. 213, 560–570. doi: 10.1016/j.jmr.2011.09.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Behrman-Lay, A. M., Usher, C., Conturo, T. E., Correia, S., Laidlaw, D. H., Lane, E. M., et al. (2015). Fiber bundle length and cognition: a length-based tractography MRI study. Brain Imaging Behav. 9, 765–775. doi: 10.1007/s11682-014-9334-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Bingham, C. (1974). An antipodally symmetric distribution on the sphere. Ann. Stat. 2, 1201–1225. doi: 10.1214/aos/1176342874

CrossRef Full Text | Google Scholar

Blum, H. (1967). “A transformation for extracting new descriptors of shape,” in Models for the Perception of Speech and Visual Form, Vol. 4. ed. W. Wathen-Dunn (MIT Press Cambridge: MIT Press), 153–171.

Google Scholar

Bruder, G. E., Stewart, J. W., and McGrath, P. J. (2017). Right brain, left brain in depressive disorders: clinical and theoretical implications of behavioral, electrophysiological and neuroimaging findings. Neurosci. Biobehav. Rev. 78, 178–191. doi: 10.1016/j.neubiorev.2017.04.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Castaneda, A. E., Tuulio-Henriksson, A., Marttunen, M., Suvisaari, J., and Lönnqvist, J. (2008). A review on cognitive impairments in depressive and anxiety disorders with a focus on young adults. J. Affect. Disord. 106, 1–27. doi: 10.1016/j.jad.2007.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Y., Xu, J., Yu, H., Nie, B., Li, N., Luo, C., et al. (2014). Delineation of early and later adult onset depression by diffusion tensor imaging. PLoS One 9:e0112307. doi: 10.1371/journal.pone.0112307

PubMed Abstract | CrossRef Full Text | Google Scholar

Cook, P. A., Bai, Y., Nedjati-Gilani, S., Seunarine, K. K., Hall, M. G., Parker, G. J., et al. (2005). Camino: open-source diffusion-MRI reconstruction and processing. 14th Scientific Meeting of the International Society for Magnetic Resonance in Medicine. 2759–2760. http://www.cs.ucl.ac.uk/research/medic/camino

Google Scholar

Corouge, I., Fletcher, P., Joshi, S., Gouttard, S., and Gerig, G. (2006). Fiber tract-oriented statistics for quantitative diffusion tensor MRI analysis. Med. Image Anal. 10, 786–798. doi: 10.1016/j.media.2006.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cullen, K. R., Klimes-Dougan, B., Muetzel, R., Mueller, B. A., Camchong, J., Houri, A., et al. (2010). Altered white matter microstructure in adolescents with major depression: a preliminary study. J. Am. Acad. Child Adolesc. Psychiatry 49, 173–183.e1.

Google Scholar

Daducci, A., Gerhard, S., Griffa, A., Lemkaddem, A., Cammoun, L., Gigandet, X., et al. (2012). The connectome mapper: an open-source processing pipeline to map connectomes with MRI. PLoS One 7:e48121. doi: 10.1371/journal.pone.0048121

PubMed Abstract | CrossRef Full Text | Google Scholar

de Groot, M., Vernooij, M. W., Klein, S., Ikram, M. A., Vos, F. M., Smith, S. M., et al. (2013). Improving alignment in tract-based spatial statistics: evaluation and optimization of image registration. NeuroImage 76, 400–411. doi: 10.1016/j.neuroimage.2013.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Desikan, R. S., Ségonne, F., Fischl, B., Quinn, B. T., Dickerson, B. C., Blacker, D., et al. (2006). An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage 31, 968–980. doi: 10.1016/j.neuroimage.2006.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Drevets, W. C., Price, J. L., and Furey, M. L. (2008). Brain structural and functional abnormalities in mood disorders: implications for neurocircuitry models of depression. Brain Struct. Funct. 213, 93–118. doi: 10.1007/s00429-008-0189-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Drysdale, A. T., Grosenick, L., Downar, J., Dunlop, K., Mansouri, F., Meng, Y., et al. (2017). Resting-state connectivity biomarkers define neurophysiological subtypes of depression. Nat. Med. 23, 28–38. doi: 10.1038/nm.4246

PubMed Abstract | CrossRef Full Text | Google Scholar

Fischl, B., Salat, D. H., Busa, E., Albert, M., Dieterich, M., Haselgrove, C., et al. (2002). Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain. Neuron 33, 341–355. doi: 10.1016/S0896-6273(02)00569-X

CrossRef Full Text | Google Scholar

Fischl, B., van der Kouwe, A., Destrieux, C., Halgren, E., Ségonne, F., Salat, D. H., et al. (2004). Automatically parcellating the human cerebral cortex. Cereb. Cortex 14, 11–22. doi: 10.1093/cercor/bhg087

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, N. I., Lewis, T., and Embleton, B. J. J. (1993). Statistical Analysis of Spherical Data. Cambridge: Cambridge University Press.

Google Scholar

Goodlett, C. B., Fletcher, P. T., Gilmore, J. H., and Gerig, G. (2009). Group analysis of DTI fiber tract statistics with application to neurodevelopment. NeuroImage 45, S133–S142. doi: 10.1016/j.neuroimage.2008.10.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Hecke, W. V., Leemans, A., de Backer, S., Jeurissen, B., Parizel, P. M., and Sijbers, J. (2009). Comparing isotropic and anisotropic smoothing for voxel-based DTI analyses: a simulation study. Hum. Brain Mapp. 31, 98–114. doi: 10.1002/hbm.20848

PubMed Abstract | CrossRef Full Text | Google Scholar

Helm, K., Viol, K., Weiger, T. M., Tass, P. A., Grefkes, C., del Monte, D., et al. (2018). Neuronal connectivity in major depressive disorder: a systematic review. Neuropsychiatr. Dis. Treat. 14, 2715–2737. doi: 10.2147/NDT.S170989

PubMed Abstract | CrossRef Full Text | Google Scholar

Hutchinson, E. B., Rutecki, P. A., Alexander, A. L., and Sutula, T. P. (2012). Fisher statistics for analysis of diffusion tensor directional information. J. Neurosci. Methods 206, 40–45. doi: 10.1016/j.jneumeth.2012.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Jenkinson, M., Beckmann, C. F., Behrens, T. E. J., Woolrich, M. W., and Smith, S. M. (2012). Review FSL. NeuroImage 62, 782–790. doi: 10.1016/j.neuroimage.2011.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Jenkinson, M., and Smith, S. (2001). A global optimisation method for robust affine registration of brain images. Med. Image Anal. 5, 143–156. doi: 10.1016/S1361-8415(01)00036-6

CrossRef Full Text | Google Scholar

Jones, D., Symms, M., Cercignani, M., and Howard, R. (2005). The effect of filter size on VBM analyses of DT-MRI data. Neuroimage 26, 546–554. doi: 10.1016/j.neuroimage.2005.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, D. K., and Cercignani, M. (2010). Twenty-five pitfalls in the analysis of diffusion MRI data. NMR Biomed. 23, 803–820. doi: 10.1002/nbm.1543

PubMed Abstract | CrossRef Full Text | Google Scholar

Kieseppä, T., Eerola, M., Mäntylä, R., Neuvonen, T., Poutanen, V. P., Luoma, K., et al. (2010). Major depressive disorder and white matter abnormalities: a diffusion tensor imaging study with tract-based spatial statistics. J. Affect. Disord. 120, 240–244. doi: 10.1016/j.jad.2009.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Kong, T. Y., and Rosenfeld, A. (1989). Digital topology: introduction and survey. Comput. Vis. Graph. Image Process. 48, 357–393. doi: 10.1016/0734-189X(89)90147-3

CrossRef Full Text | Google Scholar

Liao, Y., Huang, X., Wu, Q., Yang, C., Kuang, W., Du, M., et al. (2013). Is Depression a disconnection syndrome? Meta- analysis of diffusion tensor imaging studies in patients with MDD. J. Psychiatry Neurosci. 38, 49–56. doi: 10.1503/jpn.110180

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandl, R. W., Schnack, H. G., Luigjes, J., van den Heuvel, M. P., Cahn, W., Kahn, R. S., et al. (2010). Tract-based analysis of magnetization transfer ratio and diffusion tensor imaging of the frontal and frontotemporal connections in schizophrenia. Schizophrenia Bull. 36, 778–787. doi: 10.1093/schbul/sbn161

PubMed Abstract | CrossRef Full Text | Google Scholar

Mardia, K. V., and Jupp, PE (eds) (1999). “Directional statistics,” in Wiley Series in Probability and Statistics, (Hoboken, NJ: John Wiley & Sons, Inc), doi: 10.1002/9780470316979

CrossRef Full Text | Google Scholar

Mardia, K. V., and Zemroch, P. J. (1977). Table of maximum likelihood estimates for the bingham distribution. J. Stat. Comput. Simul. 6, 29–34. doi: 10.1080/00949657708810165

CrossRef Full Text | Google Scholar

Mayberg, H. S. (2003). Modulating dysfunctional limbic-cortical circuits in depression: towards development of brain-based algorithms for diagnosis and optimised treatment. Br. Med. Bull. 65, 193–207. doi: 10.1093/bmb/65.1.193

PubMed Abstract | CrossRef Full Text | Google Scholar

McClintock, S. M., Husain, M. M., Greer, T. L., and Cullum, C. M. (2010). Association between depression severity and neurocognitive function in major depressive disorder: a review and synthesis. Neuropsychology 24, 9–34. doi: 10.1037/a0017336

PubMed Abstract | CrossRef Full Text | Google Scholar

Metin, M. O., and Gökçay, D. (2014). Direction statistics methods evaluation for principal diffusion directions. Hum. Brain Mapp.

Google Scholar

Mori, S., Oishi, K., Jiang, H., Jiang, L., Li, X., Akhter, K., et al. (2008). Stereotaxic white matter atlas based on diffusion tensor imaging in an ICBM template. NeuroImage 40, 570–582. doi: 10.1016/j.neuroimage.2007.12.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Pajevic, S., and Pierpaoli, C. (1999). Color schemes to represent the orientation of anisotropic tissues from diffusion tensor data. Magn. Reson. Med. 42, 526–540.

Google Scholar

Palágyi, K., Balogh, E., Kuba, A., Halmai, C., Erdõhelyi, B., Sorantin, E., et al. (2001). “A sequential 3D thinning algorithm and its medical applications,” in Information Processing in Medical Imaging. IPMI 2001. Lecture Notes in Computer Science, eds M. F. Insana and R. M. Leahy (Berlin: Springer), 409–415. doi: 10.1007/3-540-45729-1_42

CrossRef Full Text | Google Scholar

Parker, G. J. M., Haroon, H. A., and Wheeler-Kingshott, C. A. M. (2003). A framework for a streamline-based probabilistic index of connectivity (PICo) using a structural interpretation of MRI diffusion measurements. J. Magn. Reson. Imag. 18, 242–254. doi: 10.1002/jmri.10350

PubMed Abstract | CrossRef Full Text | Google Scholar

Pennec, X. (2006). Intrinsic statistics on riemannian manifolds: basic tools for geometric measurements. J. Math. Imaging Vis. 25, 127–154. doi: 10.1007/s10851-006-6228-4

CrossRef Full Text | Google Scholar

Riffert, T. W., Schreiber, J., Anwander, A., and Knösche, T. R. (2014). Beyond fractional anisotropy: extraction of bundle-specific structural metrics from crossing fiber models. NeuroImage 100, 176–191. doi: 10.1016/j.neuroimage.2014.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartzman, A., Dougherty, R. F., and Taylor, J. E. (2005). Cross-subject comparison of principal diffusion direction maps. Magn. Reson. Med. 53, 1423–1431. doi: 10.1002/mrm.20503

PubMed Abstract | CrossRef Full Text | Google Scholar

Seminowicz, D. A., Mayberg, H. S., McIntosh, A. R., Goldapple, K., Kennedy, S., Segal, Z., et al. (2004). Limbic–frontal circuitry in major depression: a path modeling metanalysis. NeuroImage 22, 409–418. doi: 10.1016/j.neuroimage.2004.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sexton, C. E., Mackay, C. E., and Ebmeier, K. P. (2009). A systematic review of diffusion tensor imaging studies in affective disorders. Biol. Psychiatry 66, 814–823. doi: 10.1016/j.biopsych.2009.05.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, S. M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T. E., Mackay, C. E., et al. (2006). Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data. NeuroImage 31, 1487–1505. doi: 10.1016/j.neuroimage.2006.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, S. K., Sun, S. W., Ramsbottom, M. J., Chang, C., Russell, J., and Cross, A. H. (2002). Dysmyelination revealed through MRI as increased radial (but Unchanged Axial) diffusion of water. NeuroImage 17, 1429–1436.

Google Scholar

Sotiropoulos, S. N., Behrens, T. E. J., and Jbabdi, S. (2012). Ball and rackets: inferring fiber fanning from diffusion-weighted MRI. NeuroImage 60, 1412–1425. doi: 10.1016/j.neuroimage.2012.01.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Taoka, T., Aida, N., Fujii, Y., Ichikawa, K., Kawai, H., Nakane, T., et al. (2020). White matter microstructural changes in tuberous sclerosis: evaluation by neurite orientation dispersion and density imaging (NODDI) and diffusion tensor images. Sci. Rep. 10, 1–9. doi: 10.1038/s41598-019-57306-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Tariq, M., Schneider, T., Alexander, D. C., Wheeler-Kingshott, C. A. G., and Zhang, H. (2016). Bingham-NODDI: mapping anisotropic orientation dispersion of neurites using diffusion MRI. NeuroImage 133, 207–223. doi: 10.1016/j.neuroimage.2016.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Timmers, I., Roebroeck, A., Bastiani, M., Jansma, B., Rubio-Gozalbo, E., and Zhang, H. (2016). Assessing microstructural substrates of white matter abnormalities: a comparative study using DTI and NODDI. PLoS One 11:e167884. doi: 10.1371/journal.pone.0167884

PubMed Abstract | CrossRef Full Text | Google Scholar

van Velzen, L. S., Kelly, S., Isaev, D., Aleman, A., Aftanas, L. I., Bauer, J., et al. (2020). White matter disturbances in major depressive disorder: a coordinated analysis across 20 International cohorts in the ENIGMA MDD Working Group. Mol. Psychiatry 25, 1511–1525. doi: 10.1038/s41380-019-0477-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Benner, T., Sorensen, A., and Wedeen, V. J. (2007). Diffusion toolkit: a software package for diffusion imaging data processing and tractography. Proc Int Soc Mag Reson Med 15.

Google Scholar

Wassermann, D., Rathi, Y., Bouix, S., Kubicki, M., Kikinis, R., Shenton, M., et al. (2011). “White matter bundle registration and population analysis based on gaussian processes,” in Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 6801 LNCS, (Stapleton, NY: NIH Public Access), 320–332. doi: 10.1007/978-3-642-22092-0_27

CrossRef Full Text | Google Scholar

Watson, G. S., and Williams, E. J. (1956). On the construction of significance tests on the circle and the sphere. Biometrika 43:344. doi: 10.2307/2332913

CrossRef Full Text | Google Scholar

Wen, M. C., Steffens, D. C., Chen, M. K., and Zainal, N. H. (2014). Diffusion tensor imaging studies in late-life depression: systematic review and meta-analysis. Int. J. Geriatr. Psychiatry 29, 1173–1184. doi: 10.1002/gps.4129

PubMed Abstract | CrossRef Full Text | Google Scholar

Zalesky, A., Fornito, A., and Bullmore, E. T. (2010). Network-based statistic: identifying differences in brain networks. NeuroImage 53, 1197–1207. doi: 10.1016/j.neuroimage.2010.06.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Schneider, T., Wheeler-Kingshott, C. A., and Alexander, D. C. (2012). NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain. NeuroImage 61, 1000–1016. doi: 10.1016/j.neuroimage.2012.03.072

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, K., Huang, X., Li, T., Gong, Q., Li, Z., Ou-Yang, L., et al. (2008). Alterations of white matter integrity in adults with major depressive disorder: a magnetic resonance imaging study. J. Psychiatry Neurosci. 33, 525–530.

Google Scholar