# Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography

Authors: S. Jbabdi; P. Bellec; R. Toro; J. Daunizeau; M. Pélégrini-Issac; H. Benali
Journal: International Journal of Biomedical Imaging
 (2008)

---

## Abstract

Using geodesics for inferring white matter fibre tracts from diffusion-weighted MR data is an attractive method for at least two reasons: (i) the method optimises a global criterion, and hence is less sensitive to local perturbations such as noise or partial volume effects, and (ii) the method is fast, allowing to infer on a large number of connexions in a reasonable computational time. Here, we propose an improved fast marching algorithm to infer on geodesic paths. Specifically, this procedure is designed
to achieve accurate front propagation in an anisotropic elliptic medium, such as DTI data. We evaluate the numerical performance of this approach on simulated datasets, as well as its robustness to local perturbation induced by fiber crossing. On real data, we demonstrate the feasibility of extracting geodesics to connect an extended set of brain regions.


---

## Body

## 1. Introduction

For decades, dissection, lesion studies, or axonal
transport of tracers have been the only available techniques for studying the
brain's anatomical connections. It is not surprising that due to the
invasiveness of these methods, most of the data concerning the large-scale,
white matter tracts of the brain were collected on animals, for example, cats [1] or
monkeys [2], while structural data for the human brain were largely missing [3].
Diffusion weighted MR imaging now offers a propitious and unique framework to
explore noninvasively the organisation of white matter in the living human
brain [4, 5]. Despite the poor spatial resolution of this technique, already diffusion
data are beginning to inform us about human brain large-scale connections [6–8] and how
they relate to the functional role of cortical and subcortical networks [9, 10].Inferring on white matter architecture from diffusion
data relies on the properties of water diffusion in the tissues. Water
molecules diffuse more easily along the fibre tracts than across them, and this
anisotropy is captured by the diffusion-weighted MR signal. Inferring on
connexions given this local feature is challenging, since the observations
(diffusion properties) are indirectly related to the actual structure (axonal
orientations, size, and packing). The tractography algorithms use the
information of directionality contained in diffusion data to infer connectivity
between brain regions. Usually, information about the orientation of white
matter fibres is estimated locally, via models (e.g., diffusion tensor imaging
(DTI) [11], mixture models [12], or partial volume models [13, 14]) or in a
model-free manner (e.g., Q-ball imaging [15]). Fibre tracking consists then in
inferring connexions between distant brain regions, given this local
orientation. This can be done either in a deterministic way, by trusting the
local orientation information and following these directions until reaching a
target region (i.e., streamline tractography [16–19]), or in a probabilistic way,
by building distributions of connexions, using local probabilistic models for
fibre orientation distributions [13, 14, 20].In both cases, when tracking a fibre between two
regions of the brain, these algorithms start in one seed region, and try to
find the tracts, or distribution of tracts, that will end up in the target
region. In cases where the local orientation information present in the
diffusion data is consistent with the presence of this pathway, then these
tractography algorithms manage in general to recover the connexion between the
seed and the target. However, it often happens that in some parts of the
trajectory, the local diffusion information no longer supports the presence of
the pathway. This can either be due to a high level of noise compared to the
actual signal, or to the presence of a high number of crossing fibres
heterogeneous in their orientations. This issue is crucial in streamlining
algorithms, and is also met in probabilistic algorithms when a single
orientation per voxel is modelled [21]. The problem with those algorithms is
that when tracking from a seed, the algorithm has no information about the
region it will end up in.A possible solution to the problem of local
perturbations in the diffusion data may be provided by global tractography,
that is, optimising a global criterion while seeking for connexions. A global
tractography algorithm can potentially overcome errors in estimating local
structure, because its goal is to connect two given regions. In other words, if
we tell the algorithm which connexion we are looking at, that is, which pair of
regions is to be connected, it is better at finding it. Geodesic tractography
(GT), first proposed by Parker et al. [22], falls into this category. GT is
based on the hypothesis that brain fibers can be interpreted as minimal
distance paths (geodesics) for a metric derived from the water diffusion
profile. This distance criterion is global by definition.The basic idea for constructing a geodesic in a metric
space is to build a distance field from a seed region, the very same region one
would use as a seed for streamline tractography. This is done by solving the
so-called Eikonal equation, a partial differential equation (PDE) that
describes the time of arrival at each point of the space, as a function of the
local speed. In a constant speed field, this PDE can be easily integrated, and
the geodesics are simply straight lines. When the speed varies across the
space, the geodesics can curve, preferring high local speed locations to decrease
the arrival time. Finally, if the speed depends on the direction of travel
(e.g., along versus across a fibre tract), then the PDE   is said to be
anisotropic.Solving the Eikonal equation in a heterogeneous and
highly anisotropic medium, as is the human brain, is a technically challenging
problem [23]. This is especially true if one uses single-pass algorithms, which
is particularly important when dealing with data containing hundreds of
thousands of voxels. There have been a few attempts at solving this problem in
the context of diffusion-based tractography [22, 24–27].We describe a method for constructing geodesics in an
anisotropic medium, and apply it to the problem of DTI-based tractography. This
method relies on works in optimal path planning [28] and, more recently, vessel
extraction in 3D angiography images [29]. It has been shown to be very accurate
in anisotropic media [29],
and requires less computation than the exact method
proposed in Sethian and Vladimirsky [30] in a general framework for anisotropic optimal
path planning. The main contribution of this work is to show how this method
applies to the case of an elliptic medium, where the algorithm performs
extremely well both in terms of accuracy and efficiency, as shown in the
simulations. We also show the feasibility of applying such method to the
extraction of structural connectivity in an extended brain network using
diffusion data from a healthy subject.


## 2. Methods

In this section, we will give some theoretical
background on geodesics and the Eikonal equation, and describe a single-pass
algorithm for building geodesics.


## 2.1. Geodesics and the Eikonal Equation

A geodesic is a pathway minimising an integral of the
form(1)J(γ)=∫F(s,γ,γ′)ds,  where F(s,γ,γ′)=γ′(s)TM(γ(s))γ′(s) describes an
infinitesimal distance along a pathway γ, relative to a metric tensor M.Now, letu(x) be the arrival
time function starting from a location x0, that is, u(x) is equal to the
minimum value of the integral J(γ) along a
geodesic connecting x0 to x. Then, the arrival time function and the geodesics
satisfy these two fundamental equations:
(2a)∇uTM−1∇u=1,(2b)γ′∝M−1∇u,
where ∇u is the spatial
gradient of u. Equation (2a) is the anisotropic version of the so-called Eikonal
equation. In the isotropic case, this equation is usually written |∇u|=1/v, where v is the local
speed. Hence, this equation tells us two things: (i) it is a generalisation of
the speed equation, stating that the time of arrival is inversely proportional
to the speed, and (ii) changing the local metric tensor can be seen as changing
the local speed. Equation (2b) shows that the tangent of the geodesic lines is
parallel to the gradient of the time of arrival function with respect to the
inverse metric. This is very important because it gives us a convenient way to
reconstruct geodesics from any point in space, given the solution to the
Eikonal equation. Figure 1 shows example geodesics in an isotropic space
composed of two subsets with different local speeds.Figure 1
Example geodesics in a double isotropic space.
Black arrows show the local orientations of the geodesics. The speed in the
dark grey region is twice as high as that in the light grey one. Notice that in
each separate space, the geodesics are straight lines. Also, notice how one of
the geodesics (bold dashed lines) travels backward to the high speed part
before getting back to the low speed one.Proof.
Recall that the functionu(x) is the minimum
value of J along the
geodesic from point x0 to an arbitrary
point x: (3)u(x)=minγ∫x0xF(s,γ,γ′)ds.  A general
variation of (3) is given (see, e.g., [31]) as(4)δu=∂F∂γ′δγ+∫x0x(∂F∂γ−dds∂F∂γ′)ds. Since we have
integrated along a geodesic, the second term on the right-hand side of (4) equals zero (Euler condition). We obtain(5)∇u=∂u∂γ=∂F∂γ′=Mγ′(γ′TMγ′)1/2.  Equation (2b) directly follows. Finally, and using the symmetry of
the metric tensor M, we get the Eikonal equation: (6)∇uTM−1∇u=γ′TMTM−1Mγ′γ′TMγ′=1.Equations (2a) and (2b) summarise  the two steps for building geodesics: (i)
solve the Eikonal equation for u, given a metric tensor M and a starting
point x0; (ii) construct geodesics between any given point and
the starting point  x0 by following
the gradient of u with respect to
the inverse metric M−1.


## 2.2. Fast-Marching Algorithm

A few
algorithms  have been proposed in the literature for computing the functionu on a discrete
grid. The most popular are Tsitsiklis's method [28] and Sethian's method [32],
which are based on the construction of the time of arrival function u(x) using front
propagation. These methods are also referred to as fast marching methods
because they construct the function u in a
single-pass through the grid nodes. Tsitsiklis's method relies on (1) while Sethian's method uses the Eikonal equation (2a). Both methods are suitable in the case of isotropic
media, that is, where the metric M is proportional
to the identity matrix, but they fail in anisotropic media [23]. An exact
scheme to deal with anisotropy has been proposed by Sethian and Vladimirsky [30],
but while remaining a single-pass algorithm, it still requires a computational
effort that is growing with the amount of anisotropy. A variant of the initial
fast-marching algorithm of Tsitsiklis [28] has been proposed to deal with
anisotropic media [29], which is more computationally efficient than the exact
scheme of Sethian [30]. Yet, it relies on a generic optimisation procedure that
was undocumented for the special case of the elliptical media we face with DTI
tractography. We extended this method by deriving a solution to the
optimisation procedure in this case.The general idea of the fast-marching algorithm was
borrowed from the graph theory. It is a direct extension of Dijkstra's
algorithm for finding minimal paths in a graph [33]. The algorithm relies on a
very simple observation: suppose that the time of arrival is known inside a
close set of grid nodes (a set we will refer to as the known set). Then, the
first nodes that will be encountered by the propagating front are the nodes on
the edge of the known set (this narrow band of grid nodes will be called
the trial set). Secondly, the first node that will be encountered by the
propagating front is the closest one to known (in terms of geodesic distance),
and crucially, there will be no other way to make this distance smaller after
propagating the front further. This means that the arrival time at this voxel
will not change, and can be frozen. In other words, the value of the
time of arrival u can be
calculated, starting from x0, in a single-pass through the voxels, only by considering,
at each iteration, the neighbouring voxels of the propagating front. The other
voxels (the far set) are not examined. Figure 2(a) schematises this front propagation scheme. The
fast-marching algorithm is summarised in the appendices.(a) Grid representation of the different sets
involved during the fast-marching algorithm. (b) Position of the optimal point
on a simplex such as to minimise the geodesic distance tox. (c) Geometry of the 48 simplexes surrounding a voxel
(central grey dot). The little red dots represent the centres of the 26
neighbouring voxels.
(a)(b)(c)The crucial step in this front propagation is the
computation of the distance between the front and the neighbouring voxels in
thetrial set. In our case, this distance is anisotropic, and we cannot use the
standard methods, because they rely on the assumption that the gradients of u are parallel to
its geodesic lines (see [23] for further details). To account for the
anisotropy, we consider a set of simplexes (triangles) that cover the whole
neighbourhood around a voxel of the narrow band [29], and minimise the distance
function between the simplexes and that voxel (see Figures  2(b) and 2(c)). The introduction of these simplexes allows to
describe the trajectories on a continuous rather than a discrete grid. The
definition of a simplex neighbouring a point x is simply a set
of three points (x1,x2,x3) that are
26 neighbours of x, defining a triangle that we denote x1x2x3¯. There are 48 such triangles around x for the
26 connexities (Figure 2(c)). The procedure for computing the anisotropic
distance between the propagating front and the voxels in the trial set is given
in the appendices.During the updating procedure, the time of arrival at
a voxelxm of the trial
set is calculated from its neighbours on a simplex using an approximation
(strictly speaking, two approximations!). Normally, if the geodesic passing by xm comes from
simplex x1x2x3¯, then the time of arrival is given by (7)u(xm)ming∈x1x2x3¯{u(g)+∫gxmF(s,γ,γ′)ds}.  We use a
parametric approximation to this formula, given by the minimisation of the
following function: (8)f(α)=∑i=13αiu(xi)︸(I)+∥x−∑i=13αixi∥M︸(II),  where
∥⋅∥M is the
quadratic norm with respect to the metric M and α=(α1,α2,α3). Equation (8) follows the approximations of Tsitsiklis [28]. Term
(I) approximates the distance from the starting point x0 to the simplex
centre of mass g as a weighted
sum of the distances to the nodes of the simplex. Term (II) approximates the
remaining distance by considering the local metric as being constant, equal to
its value at xm.Minimisingf in the simplex
can be written as a constrained optimisation problem that can be solved
explicitly, since f and the simplex
are convex. The analytical solution is detailed in the appendices.


## 2.3. How to Choose the Metric?

In the GT
framework, we make the hypothesis that white matter fibres are geodesics with
respect to a metric tensor. But so far, we have not specified which metric
tensor we mean. In DTI, the inverse tensor (M=D−1) seems to be
the natural choice. Intuitively, water molecules diffusion is faster along the
tract than across them. When inverting the diffusion tensor, the highest
eigenvalues become the lowest, and the shortest distance is parallel to the
fibres. One can also notice that the inverse tensor defines a metric in a
Riemannian space that induces a Laplace-Beltrami operator (generalisation of
the Laplace operator) which is encountered in the diffusion equation [25, 34].However, the inverse tensor is not suitable in all
circumstances. Consider the situation described in Figure3 were a circular tract of radius r connects points
A and B, with diffusion tensors tangent to the tract having the same shape.
Suppose the rest of the space is isotropic, with the same mean diffusion as
along the tract. If one considers the inverse tensor metric M=D−1, the distance between A and B through the circular
path is(9)∫𝒞dxTD−1dx=πrλ1, where λ1 is the largest
eigenvalue of the tensors along the circular pathway. On the other hand, the
straight line distance between A and B is equal to 23r/trace(D). Hence, a necessary condition for the circular tract
to be a geodesic is that its length is smaller than a straight line, that is, (10)πrλ1≤23rtrace(D),  which leads to λ1≥π2trace(D)/12, that is, a condition on the tensor shape to be peaky
enough. Of course, one can imagine that even if this condition is satisfied, a
geodesic path might certainly lie somewhere in between a straight line and the
circular line, as shown in Figure 4. Which metric to choose is hence still debatable.
Nonetheless, in our simulations and real data applications, we will use the
inverse diffusion tensor as a metric for defining geodesics.Figure 3
Comparison between a
straight line and a geodesic.(a) Simulated circular
tensor field. (b)–(f) Increasing the anisotropy of the
circular tensor makes the geodesic path (red line) closer to a circle.
(a)(b)(c)(d)(e)(f)


## 3. Applications

## 3.1. Simulations

We have
evaluated the GT method on simulated data. The purpose of these simulations is
twofold. First, they show how the anisotropic fast-marching algorithm performs
on elliptic media, in both homogeneous field (where the analytical solution is
available) and a heterogeneous field. Second, they allow to compare GT with
streamlining in cases where the data present local perturbations (crossing
fibres).In a homogeneous medium, where the data support the
same diffusion tensorD in every voxel,
the analytic solution to the Eikonal equation is given by (11)u(x)=(x−x0)TD−1(x−x0).  It is easy to
check that in this case, u(x0)=0 and ∇uTD∇u=1. We generated a tensor where the two smaller
eigenvalues are equal, and gradually increased the anisotropy. Figure 5 shows the level curves of the analytic versus the
numerical solution to the Eikonal equation. The two solutions are very close
even for a large anisotropy, corresponding to a ratio of 50 between the largest
and the lowest tensor eigenvalues. Table 1 summarises the mean and standard deviations of the
relative error for different values of the anisotropy, which is expressed both
in terms of the ratio between the largest and the lowest tensor eigenvalue, or
in terms of the more widely used fractional anisotropy (FA, see, e.g., [35]).Table 1
Summary of the simulation results with an
increasing ratio between the largest and the lowest tensor eigenvalue (the
corresponding FA value is shown on the second row). Top: mean and standard
deviations of the relative error between numerical and analytic solutions for
the Eikonal equation in a homogeneous medium. Bottom: mean and standard
deviations of the value of∇uTD∇u in a circular
tensor field.
ratio1251050FA00.170.590.790.96mean (%)0.790.931.251.542.16SD (%)0.620.861.532.163.71ratio5102050100FA0.590.790.900.960.98mean0.9950.9930.9890.9971.059SD0.0680.0860.1120.2130.634Contour plots of the numerical
solution (top) and the analytic solution (bottom) to the Eikonal equation in a
homogeneous medium. Anisotropy levels are increasing from left (isotropic) to
right (ratio of 50 between the extreme tensor eigenvalues).
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)In a heterogeneous medium, such an analytical solution
does not exist. However, we can verify that the Eikonal equation is satisfied,
that is,∇uTD∇u is equal to
one. We used the same circular tensor field as shown in Figure 4. In Table 1, we show the mean and standard deviations of ∇uTD∇u for different
anisotropies. Notice that these are close to one, but with a higher deviation
from one with increasing anisotropy.Finally, we show results of GT in the case of local
perturbations. We generated a tensor field simulating a crossing fibre
situation. The zone where the two fibres cross has a diffusion tensor that is
the average of the two crossing fibres' tensors. We increased the crossing
fibre area and compared the behaviour of GT to streamlining tractography
(Figure6). As expected, because the streamlining simply
follows the direction of highest diffusion given by the tensor, the fibre
trajectory was deviated. In the case of GT, there was little, if any, deviation
from the straight line.Comparison between streamline
(top) and geodesic (bottom) tractography in the presence of a crossing fibre
bundle, the width of which increases from zero (left) to twice the width of the
principal bundle (right). Note how streamlining gets deviated from the straight
line because of partial volume effect.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


## 3.2. Real Data

Acquisition
Data from a single healthy subject were acquired atService de Neuroradiologie (CHNO des
Quinze-Vingts, Paris). Six gradient weighted and one T2-weighted images were acquired
on a 1.5 Tesla MR Scanner (GE Signa) using the following scan parameters:
128×128 image matrix, 2.03mm in-plane pixel size; 3.5mm slice thickness;
b=1000; (TR; TE) = (5000; 91.8) milliseconds; Number of averages = 8. Thirty-six contiguous
slices covering the whole brain were acquired. The total scanning time was
approximately 14 minutes.Regions of Interest
Five hundred and sixty-seven (N=567) regions
covering the whole cortex were manually selected in the DTI space. Each region
was represented by a single voxel. The anatomical localization of these regions
is shown in Figure 7. We performed a front propagation from each region,
which provided the distance functions (ui)i=1N. Then back propagation allowed us to construct the N(N−1)/2=160,461 geodesics
connecting the whole set of voxel pairs. We computed a heuristic connectivity
index consisting of the mean diffusivity along each geodesic, multiplied by the
mean FA along the pathways.
In order to better visualize this anatomical
connectivity index in a matrix form, the set of brain regions were grouped with
respect to their localization. The regions were divided into five groups,
including the frontal lobe (left: 99 voxels, right: 101 voxels), the limbic
cortex (left: 31, right: 30), the occipital lobe (left: 56; right: 54), the
parietal lobe (left: 64; right: 62), and the temporal lobe (left: 34; right:
36). This classification was based on an automatic labelling of the voxels
locations given by the Talairach Daemon (http://ric.uthscsa.edu/projects/tdc),
after registering the DTI data into the MNI standard space, and subsequent
correction from MNI to Talairach space (see, e.g., [36]). Figure 8 shows the distribution of the connectivity index, in
the matrix form, between any two regions, arranged by group and by hemisphere.
The matrix shown in Figure8  reveals an organization of the connectivity index
that follows the anatomical organization of the brain regions regarding their
locations. Since the connectivity index encompasses the anisotropy factor, its
value highly depends on which regions we are connecting, which means which
global pathways the geodesics are close to.
First, the diagonal blocks of the matrix show clearly
a lower level of connectivity than the extradiagonal blocks. This seems to
indicate that the connectivity index penalises short fibers, and inversely
favors long fibers, especially interhemispheric fibers. Secondly, the blocks
that show the highest connectivity index are the blocks that connect the right
and left occipital lobes.
This result is not surprising since the fiber tracts
that connect right and left occipital lobes follow a trajectory through the
splenium of the corpus callosum (forceps major), which is a highly anisotropic
area.Figure 7
Localisation of the regions of
interest on the cortex. 3D fronto-sagittal view.Figure 8
Anatomical connectivity matrix
rearranged into anatomical groups: F (frontal lobe), L (limbic), O (Occipital),
P (parietal), T (temporal). In each group, the left and right hemispheres are
also separated.Geodesics
We further
investigated which of the constructed geodesics may represent actual fiber
trajectories. To approach this question, we thresholded the connectivity matrix
in order to emphasize the geodesics with the largest connectivity indices.
Specifically, we considered the10% geodesics with the highest connectivity
indices for each interhemispheric block connecting symmetrical groups, taken
independently. Figure 9 represents each group of geodesics in different
colors. The most probable geodesics paths follow the principal long association
fasciculi. The frontal lobe is connected to the occipital lobe via the
fronto-occipital fasciculus. The temporal lobe is connected to the occipital
via the inferior longitudinal fasciculus, and to the frontal lobe via the
uncinate fasciculus. All major long association tracts are represented by these
geodesics.(a)10% most probable intrahemispheric
geodesics shown in the left hemisphere. Blue paths connect the occipital lobe
to the temporal lobe. Purple paths connect the frontal to the occipital lobe.
Green paths connect the frontal lobe to the temporal lobe. (b) 10% most
probable interhemispheric geodesics connecting symmetrical regions. Green:
frontal lobe, red: limbic lobe, brown: occipital lobe, blue: parietal lobe,
yellow: temporal lobe.
(a)(b)Geodesics versus Streamlining
Finally, in
order to compare the results of our method to a conventional fiber tracking
method, we performed a streamline tractography from theN seed voxels,
with four tracts per voxel. As a stopping criterion, we chose a maximum step
angle of 60∘, and an anisotropy threshold of 0.1 [19]. To
compare the results to GT, we selected the four geodesics, having the highest
probability index, for each voxel in the set of seed voxels. This way, we have
the same number of tracts using both methods (4×N tracts). Figure 10
shows the results of these two procedures. The
streamline method produces many incomplete tracts, especially association
tracts, while the proposed GT method succeeded in reconstructing the major
association and commissural tracts, including the uncinate, the inferior
fronto-temporal, and the callosal fibers. Note that the fronto-occipital tract
is not present at this level of threshold (we only considered four geodesics
per voxel).(a) Results of the streamline tractography
algorithm applied to the set of brain voxels. Four streams per voxel are
computed. The stopping criteria are60∘ for the maximal
angle step, and 0.1 for the minimal
anisotropy value. (b) Geodesics computed by the GT method. For each brain voxel
of the set, four geodesics with the highest probability index are shown.
(a)(b)


## 4. Discussion

Global optimisation is a valuable strategy in the context of path planning. When one
has the information of where to start and where to go, this information is used
to overcome local poor optimality. In the context of white matter diffusion-based
tractography, where we often have strong hypotheses about the localisation of
the regions in the brain, global optimisation can overcome some serious
weaknesses of the process. Mainly, uncertainty about local fibre orientation,
reflecting partial volume effects caused by crossing fibres, or local low
signal to noise, can be handled efficiently using GT.We have presented here a method to perform such
global-based path planning in an anisotropic medium. The method is very robust
to high anisotropy, and provides an extremely accurate numerical solution to
the Eikonal equation.On real-data experiments, the reconstructed geodesics
that have a high connectivity index correspond to known fiber tract fasciculi
connecting the cortex. These fasciculi can all be retrieved by other
tractography methods that use DTI data, providing priors on their location
using one or more regions of interest  [37, 38], especially intermediate regions
located in white matter. GT automatically depicted these fasciculi with no
prior.However, the U-shaped fibers, that is, the short
association tracts, are not favored by our connectivity index. This can be
easily seen by looking at the diagonal blocks of the matrix in Figure8. The long association tracts, as well as the
commissural fibers, are more present with a higher connectivity index.GT also allows one to construct interhemispheric
tracts between each pair of regions located in different hemispheres. These
tracts include homotopic and heterotopic connexions, that is, tracts connecting,
respectively, symmetrical and asymmetrical regions lying in different
hemispheres. It is worth noting that standard tractography methods usually fail
to recover most callosal connexions, apart from the medial ones. This is a good
illustration of the problem of crossing fibres, as those connexions cross the
superior longitudinal fasciculus. However, recent probabilistic tractography
with more complex local models has successfully traced those types of
connexions [14, 20, 21].There is an intuitive relationship between geodesic,
for the inverse tensor metric, and probabilistic tractographies. Probabilistic
tractography consists of constructing a distribution of connexions, by sampling
tracts using local orientation distributions. In the basic case where this
local probability model for fibre orientations is defined using the tensor
model (i.e., a Gaussian local model with a covariance matrix proportional to the
diffusion tensorD), the
probability of a tract following an orientation given by dx at a location x writes(12)p(x+dx∣x)=𝒩(x,D),  then, for some
pathway γ connecting x0 to x1, and for some discretisation of this pathway, the
probability of moving along γ is the product
of the infinitesimal step probabilities:(13)p(x0→x1)=∏k=1np(x0+kdx∣x0+(k−1)dx)∝∏k=1nexp{−12dxTD−1dx}=exp{−12∑k=1ndxTD−1dx}→exp{−12∫γdxTD−1dx}≤exp{−12(∫γdxTD−1dx)2}.   Maximising this
probability could then be related to minimising the geodesic distance, relative
to the inverse tensor metric. While the probabilistic method gives a
distribution of connexions, GT gives the mode of this distribution, that is, the
path with highest probability. Note also that the probabilistic model given by
(12) can be improved to fit the data more accurately
(e.g., multiple tensors, etc.), which can be seen as a change in the metric
tensor in GT.Using GT, it is possible to study the organisation of
large brain networks in terms of their anatomical connexions. Such networks
have been studied in terms of structural invariants in a graph theoretical
framework by several authors [39–41]. These works have been conducted for studying the
structural organisation of the cat or primate brain, as well as for the human
functional brain organisation, but have never been applied to large human
anatomical networks, because no method has been proposed to construct such
networks. GT could provide this structural information, via a graph that has
been thresholded or not, since the connectivity index in itself contains
information about the connectional structure.There are two major issues when using geodesics for
the tractography. First, choosing a metric for which geodesics represent fibre
pathway trajectories is not straightforward. The correct metric might show more
anisotropy than the diffusion tensor, as discussed earlier. Also, the choice of
the metric might depend on the white matter fibres under investigation. The
second issue is that, for any pair of regions in the brain, there exists a
geodesic between those regions. However, this is not true for white matter
fibres. One then has to decide when a geodesic is a fibre trajectory, for example, by
defining indices and performing statistical thresholding under some null
hypothesis. This problem of thresholding tractography results is not specific
to GT, but is met by any other tractography method. It is though a bigger
problem in the case of GT because every pair of regions is potentially
connected. Another problem with GT is that, in the presence of two separate
connexions between two regions, we are only able to detect one of them (the
shortest one in terms of geodesic distance).One way to validate GT results would be by comparison
with another measure of connectivity. For example, measures of functional
connectivity using functional magnetic resonance imaging (fMRI) by means of
correlations [42] or partial correlations [43] are thought to be closely linked
to the anatomical structure sustaining the brain regions, seen as graph nodes.
The GT technique provides a unique tool for performing a comparison between
anatomical and functional connectivity, since it can apply to large networks,
and provide a measure of anatomical connectivity between each pair of nodes of
the brain network. It can readily be used to compare the architectures of brain
networks that have been studied in humans from the functional perspective (e.g.,
Salvador et al. [44] used partial correlations of fMRI data on a set of 100
regions), or using voxel-based morphometry to correlate cortical thickness
between different cortical areas (e.g., He et al. [45] used this technique to
study 100 cortical areas in humans). Such investigations have considerable
possible applications, both cognitive and clinical. On the one hand, this
method could serve as a basis for comparing anatomical and functional
connectivities, as said earlier, and could help to understand how the brain works
as an evolving network. On the other hand, the structure of restricted networks
has already helped to distinguish between healthy subjects and patients, for
example, Alzheimer disease in the case of functional connectivity [46], and
Schizophrenia in the case of white matter morphology [47]. The GT method could
serve for the characterisation of the structural organisation of those brain
networks in terms of their connectional fingerprints.



---

*Source: 320195-2007-12-06.xml*