AnatomiCuts: Hierarchical clustering of tractography streamlines based on anatomical similarity

Accepted Manuscript AnatomiCuts: Hierarchical clustering of tractography streamlines based on anatomical similarity Viviana Siless, Ken Chang, Bruce Fischl, Anastasia Yendiki PII:

S1053-8119(17)30892-3

DOI:

10.1016/j.neuroimage.2017.10.058

Reference:

YNIMG 14436

To appear in:

NeuroImage

Received Date: 22 March 2017 Revised Date:

18 October 2017

Accepted Date: 26 October 2017

Please cite this article as: Siless, V., Chang, K., Fischl, B., Yendiki, A., AnatomiCuts: Hierarchical clustering of tractography streamlines based on anatomical similarity, NeuroImage (2017), doi: 10.1016/ j.neuroimage.2017.10.058. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

RI PT

AnatomiCuts: Hierarchical Clustering of Tractography Streamlines Based on Anatomical Similarity Viviana Silessa , Ken Changa , Bruce Fischla,b , Anastasia Yendikia a

SC

Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA b Computer Science and Artificial Intelligence Laboratory, Massachusetts Institute of Technology, Cambridge, MA, USA

M AN U

Abstract

Diffusion MRI tractography produces massive sets of streamlines that contain a wealth of information on brain connections. The size of these datasets creates a need for automated clustering methods to group the streamlines into meaningful bundles. Conventional clustering techniques group streamlines

D

based on their spatial coordinates. Neuroanatomists, however, define white-

TE

matter bundles based on the anatomical structures that they go through or next to, rather than their spatial coordinates. Thus we propose a similarity

EP

measure for clustering streamlines based on their position relative to cortical and subcortical brain regions. We incorporate this measure into a hierarchical clustering algorithm and compare it to a measure that relies on Euclidean

AC C

distance, using data from the Human Connectome Project. We show that the anatomical similarity measure leads to a 20% improvement in the overlap

of clusters with manually labeled tracts. Importantly, this is achieved without introducing any prior information from a tract atlas into the clustering algorithm, therefore without imposing the existence of any named tracts. Keywords: hierarchical clustering, normalized cuts, tractography, diffusion

Preprint submitted to Neuroimage

October 27, 2017

ACCEPTED MANUSCRIPT

MRI 1. Introduction

RI PT

1

Diffusion MRI (dMRI) allows us to estimate the preferential direction of

3

water molecule diffusion at each voxel in the white matter (WM) (Le Bihan

4

et al., 1986). Tractography algorithms follow these directions to reconstruct

5

continuous paths of diffusion (Mori et al., 1999; Basser et al., 2000). The

6

most common approach to segmenting the WM from dMRI data is to use

7

every voxel in the brain as a seed for tractography and to group the resulting

8

streamlines into bundles. Recent advances in dMRI acquisition hardware and

9

software (Wright & Wald, 1997; Setsompop et al., 2013) have increased both

10

spatial and angular resolution, yielding large tractography datasets that are

11

difficult to parse manually. This creates a need for computational methods

12

that can group streamlines into anatomically meaningful bundles automati-

13

cally.

TE

D

M AN U

SC

2

Several unsupervised clustering algorithms have been proposed to assign

15

tractography streamlines into bundles. Each method is characterized by (i)

16

the measure that it uses to quantify the similarity between two streamlines

17

and (ii) the algorithm that it follows to group streamlines based on their

18

similarity.

20

21

AC C

19

EP

14

Typical similarity measures for unsupervised clustering of streamlines are

functions of the streamline coordinates in Euclidean space. These measures express either the distance between streamlines (Ding et al., 2003; Gerig

22

et al., 2004; Moberts et al., 2005; Tsai et al., 2007; Wassermann et al., 2010;

23

Visser et al., 2011; Siless et al., 2013) or the similarity of some shape de-

2

ACCEPTED MANUSCRIPT

scriptor computed from spatial coordinates (Gerig et al., 2004; Brun et al.,

25

2004; Ros et al., 2013; Zhang et al., 2014). A subset of these methods use the

26

distance metrics to project the streamlines onto an eigenspace and perform

27

the clustering in that space (Tsai et al., 2007; O’Donnell & Westin, 2007).

RI PT

24

However, bundling streamlines based on their proximity or their shape

29

is not consistent with the approach followed by neuroanatomists. In neu-

30

roanatomy, WM bundles are defined based on the brain structures that they

31

go through or next to, rather than their coordinates in a reference space. For

32

example, “The uncinate fasciculus extends between the anterior temporal

33

lobe and the frontal lobe passing though the floor of the insula in an area be-

34

tween the anterior claustrum and the amygdala” (Makris & Pandya, 2008).

35

Here we propose a similarity measure that adopts the neuroanatomist’s ap-

36

proach of comparing streamlines based on their anatomical neighborhood.

M AN U

SC

28

Previous attempts to incorporate anatomical information in unsupervised

38

streamline clustering used the structures that the streamlines terminate in,

39

either in the similarity measure itself (Tunc et al., 2014), in a post-hoc man-

40

ner (Wassermann et al., 2010), or for initialization (Guevara et al., 2011).

41

Alternatively, supervised clustering approaches introduced prior information

42

on WM anatomy from an atlas of predefined bundles, labeled by an ex-

43

pert (O’Donnell & Westin, 2007; Maddah et al., 2008; Ziyan et al., 2009;

45

46

TE

EP

AC C

44

D

37

Guevara et al., 2012; Wang et al., 2013; Ros et al., 2013; Jin et al., 2014). Another class of supervised methods used deterministic rules on the position of streamlines with respect to anatomically defined regions, thus going

47

directly from the whole-brain tractography data to a set of pre-determined

48

WM bundles (Xia et al., 2005; Li et al., 2010; Wassermann et al., 2016). In

3

ACCEPTED MANUSCRIPT

Xia et al. (2005) and Li et al. (2010) clustering based on Euclidean dis-

50

tance was performed as a subsequent step, to refine these bundles, but the

51

anatomical regions were only used in the form of deterministic rules before

52

that clustering step.

RI PT

49

In this work we propose an anatomical similarity measure for unsuper-

54

vised streamline clustering that uses not only the anatomical regions that

55

the streamlines terminate in, but all regions that form the anatomical neigh-

56

borhood of a streamline, everywhere along its trajectory. We have previously

57

used a similar description of anatomical neighborhoods to incorporate prior

58

information from a set of training subjects into the tractography step it-

59

self (Yendiki et al., 2011). However, that was a supervised approach, limited

60

to a set of predefined bundles from an atlas. Here we are interested in ex-

61

ploratory studies that do not rely on training data. Therefore, we do not

62

impose the existence of any named pathways based on a WM atlas; instead,

63

we cluster whole-brain tractography streamlines based on their position with

64

respect to a set of cortical and subcortical structures.

TE

D

M AN U

SC

53

The algorithms that have been used for unsupervised streamline cluster-

66

ing are either partitioning or hierarchical. Partitioning approaches include

67

K-means (Tsai et al., 2007; O’Donnell & Westin, 2007; Siless et al., 2013),

68

affinity propagation (Zhang et al., 2014), maximum likelihood (Maddah et al.,

70

71

AC C

69

EP

65

2008), and maximum a posteriori (Wang et al., 2011; Tunc et al., 2014). Hierarchical algorithms that have been applied to unsupervised streamline clustering are either agglomerative or divisive. Agglomerative algorithms fol-

72

low a bottom-up approach, first assigning each streamline to its own cluster

73

and then merging clusters iteratively (Ding et al., 2003; Gerig et al., 2004;

4

ACCEPTED MANUSCRIPT

Wassermann et al., 2010; Visser et al., 2011; Guevara et al., 2012; Ros et al.,

75

2013). Divisive algorithms follow a top-down approach, first assigning all

76

streamlines to a single cluster and then splitting clusters iteratively (Brun

77

et al., 2004; Ziyan et al., 2009; Wu et al., 2012). Hierarchical methods orga-

78

nize streamlines into a tree structure, where larger clusters are composed of

79

smaller clusters. This is particularly suitable for describing WM organization,

80

as anatomical tracing studies have revealed that the large WM pathways are

81

subdivided into multiple smaller bundles that share a large part of their tra-

82

jectory but originate or terminate in different brain regions (Lehman et al.,

83

2011). Hence we follow the hierarchical approach here.

M AN U

SC

RI PT

74

Our proposed method, AnatomiCuts, incorporates an anatomical similar-

85

ity measure into normalized cuts, an algorithm that has been used previously

86

for streamline clustering with conventional similarity measures (Brun et al.,

87

2004; O’Donnell et al., 2006; Ziyan et al., 2009). We follow an implementa-

88

tion of normalized cuts where a hierarchy of clusters is generated by recursive

89

bipartitions of the full data set (Shi & Malik, 2000; Brun et al., 2004). The

90

present work expands on a prior conference publication (Siless et al., 2016),

91

by adding much more extensive experimental results on the performance of

92

the proposed method and the effects of various tuning parameters that are

93

involved in its implementation.

95

96

TE

EP

AC C

94

D

84

Specifically, we compare the performance of our anatomical similarity

measure to one based on Euclidean distance between streamlines, using data from the MGH-USC Human Connectome Project. We show that cluster-

97

ing streamlines based on their anatomical neighborhood rather than their

98

spatial coordinates leads to a 20% improvement in the agreement of the clus-

5

ACCEPTED MANUSCRIPT

ters with WM bundles labeled manually by a human rater. Importantly, we

100

achieve this improvement without using prior information from clusters la-

101

beled manually in training subjects; we only use a subject’s own anatomical

102

segmentation. Our approach is unsupervised, allowing exploratory studies

103

of whole-brain tractography data, without being constrained to a predeter-

104

mined set of bundles.

105

2. Methods

106

2.1. Similarity measures

M AN U

SC

RI PT

99

107

Let fi be a tractography streamline, defined as a sequence of N points,

108

fi = [xi1 , . . . , xiN ], where xik ∈ R3 , k = 1, ..., N . A tractography dataset is a

109

set of M streamlines, F = {f1 , ..., fM }.

In the following, we assume that all streamlines have the same number of

111

points N . Although the similarity measures can be adapted to streamlines

112

with variable numbers of points, the approach of downsampling streamlines

113

to an equal number of points has been used previously by others to make

114

computation tractable (O’Donnell & Westin, 2007; Visser et al., 2011; Gary-

115

fallidis et al., 2012; Guevara et al., 2012; Wu et al., 2012; Siless et al., 2013).

116

2.1.1. Euclidean similarity

118

119

120

TE

EP

AC C

117

D

110

We define a similarity measure based on the mean Euclidean distance

between pairs of points on two streamlines fi and fj , as previously used in the literature (Visser et al., 2011; Garyfallidis et al., 2012; Guevara et al., 2012; Wu et al., 2012; Siless et al., 2013): N −1 1 X ωE (fi , fj ) , 1 + kxik − xjk k2 . N k=1



6

ACCEPTED MANUSCRIPT

Since the ordering of points is not consistent across streamlines, it is possible

122

for the first point of fi to be closer to the last point of fj and vice versa. We

123

account for this by also evaluating the similarity between fi and the reversed

124

fj . This leads to the following definition for the similarity measure:

RI PT

121

125

(1)

SC

wE (fi , fj ) , max{ωE (fi , fj ), ωE (fi , fjrev )}, where fjrev = [xjN , . . . , xj1 ].

Another option for converting a Euclidean distance into a similarity mea-

127

sure is by means of a Gaussian kernel. However that formulation involves the

128

choice of a kernel width, and the optimal width can vary greatly depending

129

on the input (Brun et al., 2004; O’Donnell & Westin, 2007). Here we opt

130

for a similarity measure that does not involve free parameters (Sammut &

131

Webb, 2017).

132

2.1.2. Anatomical similarity

D

M AN U

126

We introduce a new similarity measure that uses anatomical information

134

in the form of a cortical and subcortical segmentation, S(x), x ∈ R3 . This is

135

done by finding the segmentation labels that each streamline goes through

136

or next to, at each point along its trajectory. Specifically, each point x on a

137

streamline is associated with a set of segmentation labels, S(x + dl (x)vl ), l =

139

140

EP

AC C

138

TE

133

1, ..., P , where dl (x) is the minimum d > 0 such that S(x+dvl ) 6= S(x). That is, for each point x, we find the nearest neighboring segmentation labels in a set of directions vl , l = 1, ..., P . A neighborhood of P = 26 elements includes

141

neighboring labels in the directions vl = [e1 , e2 , e3 ], where e1,2,3 ∈ {−1, 0, 1}.

142

We also define v0 = [0, 0, 0] to represent the segmentation label that the

143

streamline passes through. Fig. 1 shows an illustration of this neighborhood. 7

ACCEPTED MANUSCRIPT

144

Fig. 2 shows an example of a point in the WM and some of its neighboring

SC

6 neighbors

RI PT

segmentation labels in different directions.

14 neighbors

M AN U

26 neighbors

Figure 1: Directions in which nearest neighboring segmentation labels are found, for a streamline point that lies in the center of the cube. Red directions belong to the 6, 14 and 26-element neighborhoods. Blue directions belong to the 14 and 26-element neighbor-

D

hoods. Green directions belong only to the 26-element neighborhood. 145

For each direction l = 0, ..., P , we compute a label histogram Hil ∈ RK ,

147

where K is the total number of labels in the anatomical segmentation. This

148

histogram represents the frequency with which different segmentation labels

149

are the l-th neighbor across all points on the i-th streamline. We introduce a

150

similarity measure between two streamlines fi and fj that expresses the joint

151

probability of their anatomical neighborhoods. We define this anatomical

EP

AC C

152

TE

146

similarity measure as: wA (fi , fj ) , |Li ∩ Lj |

P X

hHil , Hjl i,

(2)

l=0

153

where h·, ·i is the inner product, and Li , Lj are the sets of all labels found to

154

be neighbors of streamlines fi , fj . The normalization factor |Li ∩ Lj |, which 8

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

vl

Neighbor

0

[0, 0, 0]

Left unsegmented white matter

1

[1, 0, 0]

Left pallidum

2

[−1, 0, 0]

Left thalamus proper

3

[0, 0, 1]

Left paracentral white matter

4

[0, 0, −1]

Left ventral diencephalon

5

[1, 0, 1]

Left insular white matter

6

[1, 0, −1]

Left hippocampus

7

[−1, 0, 1]

Left caudate

EP

TE

D

l

8

[−1, 0, −1] Third ventricle

...

...

AC C

...

Figure 2: Example of the anatomical neighborhood of a point in the corticospinal tract. Neighboring FreeSurfer segmentation labels in the L-R and I-S directions are shown in the table.

9

ACCEPTED MANUSCRIPT

is the number of common neighbors between the two streamlines, penalizes

156

trivial streamlines with too few neighbors.

157

2.2. Normalized cuts

RI PT

155

We approach clustering as a graph partitioning problem, where the nodes

159

of the graph are the elements to be clustered and the weights of the edges

160

between nodes are the similarities between elements. These weights are used

161

to form a similarity matrix, and clusters are defined based on the eigenvectors

162

of this matrix. Given a connected graph G, the normalized cuts algorithm

163

searches for a graph cut that divides G into sets A and B (A ∩ B = ∅), by

164

minimizing similarity between A and B and maximizing similarity within

165

A and B (Shi & Malik, 2000). Clusters are split recursively, generating a

166

top-down hierarchical structure.

M AN U

SC

158

When clustering tractography data, each streamline represents a node on

168

a graph G, and the weight of each edge is the similarity between the nodes

169

it connects. In this case graph G is a clique, and thus each streamline is

170

connected to every other streamline. For the purpose of finding the optimal

171

cut of the graph into A and B, we need to quantify the similarity between A

174

175

176

TE

EP

173

and B. This is done by summing the weights of the edges that connect them: P s(A, B) = u∈A,v∈B w(u, v), where w(·, ·) is a similarity function between two

AC C

172

D

167

streamlines. To avoid trivial solutions where A or B is a single isolated node, P s is normalized by an association measure a(A) = u∈A,t∈G w(u, t). Thus the minimum cut is defined based on the following cost function: min snorm (A, B), A,B

177

snorm (A, B) =

s(A, B) s(A, B) + . a(A) a(B)

(3)

We can formulate this equation with matrix notation where W is the similar10

ACCEPTED MANUSCRIPT

P

178

ity matrix with Wij = w(i, j), D is a diagonal matrix with Dii =

179

and x the partitioning vector, where xi = 1 if the i-th element belongs to A

180

and xi = −1 is the i-th element belongs to B. The association measure will

181

be encoded in k =

183

:

(1 + x)T (D − W)(1 + x) (1 + x)T (D − W)(1 + x) + . k1T D1 (1 − k)1T D1

The problem is embedded in the real-valued domain and can be efficiently solved by the following equation (Shi & Malik, 2000): 1

M AN U

182

di

SC

snorm (A, B) =

xi >0 P di

1

D 2 (D − W)D 2 z = λz, 184

w(i, j)

RI PT

P

j

(4)

where z = Dy, and D − W is a Laplacian matrix.

The second eigenvector of the Laplacian matrix is a solution to the above

186

equation. The optimal cut is approximated by assigning the i-th node to A

187

if yi > 0, where yi ∈ R is the i-th element of y, and to B otherwise (Golub

188

& Van Loan, 1996). The algorithm is run iteratively by cutting the largest

189

remaining cluster, until a desired number of clusters is reached or a threshold

190

for minimum cluster size or dispersion is met.

191

2.3. Prototype streamlines

193

194

TE

EP

AC C

192

D

185

Building an affinity matrix requires computing the similarity between

every pair of streamlines, which is extremely time- and memory-intensive for the large, high-resolution data sets that are increasingly common. In

195

previous work, the size of the affinity matrix was reduced by processing left,

196

right, and inter-hemisphere streamlines separately (Guevara et al., 2012);

197

by using a random subset of streamlines to guide the clustering of the full 11

ACCEPTED MANUSCRIPT

set (O’Donnell & Westin, 2007; Siless et al., 2013; Ros et al., 2013; Zhang

199

et al., 2014); or by generating clustering solutions from multiple subsets of

200

streamlines and merging them (Visser et al., 2011).

RI PT

198

The approach of using a subset of the streamlines to compute the affinity

202

matrix has been advocated as a way to reduce sensitivity to outliers in hier-

203

archical clustering (Guha et al., 2001). When clustering tractography data,

204

each bundle contains hundreds to thousands of similar streamlines, resulting

205

in redundancy. Therefore it is reasonable to expect a subset of streamlines

206

to be sufficient for capturing the general shape of the clusters. In the divi-

207

sive algorithm that we use in this work, we evaluate the similarity matrix to

208

make a binary cut. In early stages of the algorithm, where computational de-

209

mands are highest because of the large number of streamlines to be divided,

210

these divisions are rather coarse (for example, separating the left hemisphere

211

from the right hemisphere, or the SLF from the arcuate fasciculus). This

212

means that we only need to select enough streamlines to ensure that we have

213

representive examples from these coarse sets.

TE

D

M AN U

SC

201

Here we select a random subset of the streamlines to be clustered at each

215

iteration. We generate the affinity matrix from this subset, which we refer to

216

as prototype streamlines. After finding the minimum cut for the prototype

217

streamlines, we assign each non-prototype streamline to the same cluster as

219

220

AC C

218

EP

214

its closest prototype streamline. As part of the experiments described below, we investigate the effect that the number of prototype streamlines has on the performance of our method.

12

ACCEPTED MANUSCRIPT

221

2.4. Data acquisition For the experimental evaluation of our unsupervised clustering, we use

223

dMRI and structural MRI (sMRI) data from 32 healthy subjects, avail-

224

able publicly through the MGH-USC Human Connectome Project (http:

225

//humanconnectomeproject.org). The data was acquired on the MGH

226

Siemens Connectom, a Skyra 3T MRI system with a custom gradient capa-

227

ble of maximum strength 300 mT/m and slew rate 200 T/m/s (Setsompop

228

et al., 2013). The sMRI data was acquired with MEMPRAGE (Mugler &

229

Brookeman, 1990; Van der Kouwe et al., 2008), TR=2530ms, TE=1.15ms,

230

TI=1100ms, 1mm isotropic resolution. The dMRI data was acquired with

231

2D EPI, TR=8800ms, TE=57ms, 1.5mm isotropic resolution, 512 gradient

232

directions, bmax = 10, 000s/mm2 . The data were pre-processed as described

233

in Fan et al. (2016).

234

2.5. Data analysis

D

M AN U

SC

RI PT

222

We reconstruct orientation distribution functions from the dMRI data us-

236

ing the generalized q-sampling imaging model (Yeh et al., 2010) and perform

237

whole-brain deterministic tractography using DSI Studio (Yeh et al., 2013).

238

We seed every voxel in the segmentation map computed by FreeSurfer. Previ-

239

ously, separate processing strategies have been proposed for longer vs. shorter

241

242

EP

AC C

240

TE

235

streamlines (Guevara et al., 2012; Zhang et al., 2014). Here we focus on long-range connections, as short-range connections, such as those representing U-fibers, could be clustered simply based on the anatomical regions that

243

they terminate in. To obtain a large number of long-range connections while

244

keeping computation tractable, we generate a total of 500,000 streamlines

245

per subject and then exclude any streamlines shorter than 55mm, leaving 13

ACCEPTED MANUSCRIPT

between 100,000 and 150,000 streamlines per subject. We downsample all

247

streamlines to N = 10 equispaced points.

RI PT

246

For comparison with our unsupervised clustering, we use a set of ma-

249

jor WM bundles labeled manually in each subject by a trained rater, fol-

250

lowing the protocol described in Wakana et al. (2007). This includes the

251

following bundles: corticospinal tract (CST), inferior longitudinal fascicu-

252

lus (ILF), uncinate fasciculus (UNC), anterior thalamic radiation (ATR),

253

cingulum - supracallosal bundle (CCG), cingulum - infracallosal (angular)

254

bundle (CAB), superior longitudinal fasciculus - parietal bundle (SLFP), su-

255

perior longitudinal fasciculus - temporal bundle (SLFT), corpus callosum -

256

forceps major (FMAJ), corpus callosum - forceps minor (FMIN). The intra-

257

hemispheric bundles, i.e., all except FMAJ and FMIN, are labeled on both

258

hemispheres (L and R), leading to a total of 18 manually labeled WM bun-

259

dles per subject. The rater interacts with the whole-brain tractography data

260

in TrackVis (http://trackvis.org), using the fractional anisotropy map and

261

the color-coded map of the primary eigenvector of the tensor for anatomical

262

guidance. For each of the bundles, at least two inclusion masks are hand-

263

drawn on slices of those maps (as specified in Wakana et al. (2007)) to isolate

264

the streamlines that belong to the bundle. Additional inclusion and exclu-

265

sion masks are drawn as needed to refine the labeling. Fig. 3 shows the 18

267

268

M AN U

D

TE

EP

AC C

266

SC

248

manually labeled WM bundles from an example subject. Anatomical segmentations are obtained by processing the sMRI data with

the automated cortical parcellation and subcortical segmentation tools in

269

FreeSurfer (Fischl et al., 2002, 2004). As part of this processing, subcortical

270

WM labels are defined by classifying each WM voxel that is within 5mm

14

(b)

SC

(a)

RI PT

ACCEPTED MANUSCRIPT

Figure 3: Sagittal (a) and axial (b) view of the 18 manually labeled bundles that we use for

M AN U

comparison to the unsupervised clustering results, shown in a randomly selected subject.

from the cortex based on its nearest cortical label. This results in a total

272

of K = 261 cortical and subcortical labels. Each subject’s dMRI and sMRI

273

data are aligned with a within-subject, boundary-based, affine registration

274

method, which optimizes the contrast of the b = 0 dMRI volumes across the

275

grey-white boundary surfaces reconstructed from the sMRI data (Greve &

276

Fischl, 2009). The subject’s anatomical segmentation is then mapped from

277

sMRI space to dMRI space by applying the inverse of that transformation.

TE

D

271

For each subject, we perform unsupervised clustering with normalized

279

cuts and each of the two similarity measures defined in section 2.1. As the

280

selection of the total number of clusters is arbitrary, we evaluate the effect

281

of this number throughout our experiments. We do this by running the

283

284

AC C

282

EP

278

clustering algorithm until a total of 200 clusters are generated, and then pruning the hierarchical clustering tree to keep the first 75, 100, 125, 150 or 200 clusters.

15

ACCEPTED MANUSCRIPT

2.6. Evaluation metrics

286

2.6.1. Comparison with manual labeling

RI PT

285

For each subject, we evaluate the extent to which the clusters obtained

288

with each of the two similarity measures resemble bundles labeled by a human

289

rater in the same subject. First, we compute the overlap, as quantified by the

290

Dice coefficient (Dice, 1945; Sørensen, 1948), between each manually labeled

291

bundle and the union of all clusters for which at least 5% of streamlines

292

belong to the manually labeled bundle. The goal of this evaluation is to

293

determine the extent to which streamlines that a human rater would classify

294

into a single anatomical structure become divided among disparate clusters.

295

Second, we compute the homogeneity and the completeness of each clus-

296

ter of streamlines with respect to the 18 manually labeled bundles in each

297

subject. These measures are defined in Rosenberg & Hirschberg (2007). Ho-

298

mogeneity is highest when each cluster contains elements from a single manu-

299

ally labeled bundle only and completeness is highest when elements from the

300

same manually labeled bundle are assigned to a single cluster. It is important

301

to evaluate both measures, as it is possible to achieve perfect homogeneity

302

by compromising completeness, i.e., by assigning each streamline to its own

303

cluster.

305

306

M AN U

D

TE

EP

AC C

304

SC

287

2.6.2. Anatomical and spatial consistency of clusters One might expect streamlines that are close to each other in Euclidean

space to also have common anatomical neighbors, and vice versa. Hence,

307

we evaluate the extent to which streamlines that are clustered based on

308

the anatomical similarity measure wA are also similar with respect to the

309

Euclidean similarity measure wE and vice versa. To this end, we find the 16

ACCEPTED MANUSCRIPT

centroid of each cluster, defined as the streamline closest to the average of all

311

streamlines in a cluster. We then evaluate wE and wA between every element

312

in a cluster and the cluster centroid.

RI PT

310

We also compute a distance metric that neither the wE -optimal nor the

314

wA -optimal clusters are optimized for. This is the mean of the closest-point

315

distance (Corouge et al., 2004; Moberts et al., 2005; O’Donnell et al., 2006;

316

Guevara et al., 2012) between two streamlines fi and fj , defined as:

317

N 1 X N min kxik − xjl k2 , N k=1 l=1

M AN U

ωCP (fi , fj ) ,

SC

313

where ωCP (fi , fj ) and ωCP (fj , fi ) are averaged to ensure symmetry: 1 wCP (fi , fj ) , (ωCP (fi , fj ) + ωCP (fj , fi )). 2 In our case, fj is the centroid of the cluster that fi belongs to.

319

2.6.3. Robustness to errors in the anatomical segmentation

D

318

(5)

As our proposed similarity measure relies on an anatomical segmentation,

321

we evaluate the robustness of the clustering to errors in the boundaries of the

322

anatomical segmentation labels. To this end, we perturb the boundaries of

323

the FreeSurfer cortical and subcortical segmentation labels in each subject,

324

and repeat the clustering using the perturbed segmentation. We perform

326

327

EP

AC C

325

TE

320

the perturbations by replacing the segmentation label of each voxel with a label selected randomly among the labels of its neighboring voxels. This is done on segmentations that have been mapped into dMRI space, thus the

328

perturbations are in the order of the dMRI voxel size (1.5mm). To establish

329

that this is comparable to the actual accuracy of the FreeSurfer anatomical

330

segmentation, we also show plots of the mean closest point distance between 17

ACCEPTED MANUSCRIPT

automatically segmented and manually labeled brain structures in 26 sub-

332

jects from the FreeSurfer atlas. We note that no manual editing has been

333

performed on the FreeSurfer reconstructions that we use in our experiments.

334

2.6.4. Effect of prototype streamlines

RI PT

331

For the evaluations described above, we used 500 randomly selected pro-

336

totype streamlines to compute the similarity matrices. We assess the impact

337

of this choice by repeating the clustering using 5 different sets of 500 proto-

338

type streamlines, and also by varying the number of prototype streamlines

339

between 50 and 3000.

340

2.6.5. Effect of anatomical neighborhood size

M AN U

SC

335

For the evaluations described above, we computed our anatomical similar-

342

ity measure using the 26-element neighborhood from Fig. 1. We investigate

343

the impact of this neighborhood size by repeating the clustering with the

344

6-element and 14-element anatomical neighborhoods.

345

2.6.6. Effect of streamline downsampling

TE

D

341

When computing the distance of tractography streamlines in Euclidean

347

space, it is common for clustering algorithms to downsample streamlines to

348

a constant number of points N (O’Donnell & Westin, 2007; Visser et al.,

350

351

AC C

349

EP

346

2011; Garyfallidis et al., 2012; Guevara et al., 2012; Wu et al., 2012; Siless et al., 2013). This reduces computation time for Euclidean distances dramatically, as it makes the task of finding corresponding points between two

352

streamlines trivial. Although our anatomical similarity measure does not re-

353

quire all streamlines to have an equal number of points, we did downsample

354

streamlines for the purpose of comparing to the Euclidean distance measure. 18

ACCEPTED MANUSCRIPT

We evaluate the effect of this downsampling on our anatomical similarity

356

measure by comparing results obtained with streamlines sampled to N = 10

357

and N = 50 points.

358

2.6.7. Comparison to an alternative segmentation

RI PT

355

Our method is unsupervised, i.e., it does not impose the existence of

360

a set of named WM bundles from an atlas. Nonetheless, we explore here

361

the use of labels from a WM atlas, as an alternative to the subjects’ indi-

362

vidual FreeSurfer cortical parcellation and subcortical segmentation labels.

363

Specifically, we use the well-known JHU-ICBM WM atlas (Wakana et al.,

364

2007; Hua et al., 2008; Oishi et al., 2010), as distributed with the FSL soft-

365

ware package (Jenkinson et al., 2012). We map this WM atlas to each for

366

our subjects’ individual dMRI space by a non-linear transformation that we

367

obtain from co-registering the fractional anisotropy maps of the atlas and

368

the individual using the FNIRT tool from FSL. We then cluster each sub-

369

ject’s streamlines by computing our anatomical similarity measure based on

370

the WM labels from this atlas, instead of the labels from the subject’s own

371

FreeSurfer cortical parcellation and subcortical segmentation.

372

2.6.8. Comparison to shape similarity

EP

TE

D

M AN U

SC

359

AC C

Finally, we compared our anatomical similarity to a similarity measure

376

377

of prototype streamlines.

373

374

375

based on streamline shape descriptors, as defined in Brun et al. (2004). Shape descriptors reduce the dimensionality of streamlines, so that the full set of streamlines can be used for clustering, as an alternative to using only a set

19

ACCEPTED MANUSCRIPT

378

2.6.9. Computation time The computational complexity of evaluating the elements of the similarity

380

matrix is O(Mp2 )Ow , where Mp the number of prototype streamlines and Ow

381

the computational complexity of the similarity measure for a pair of stream-

382

lines. The latter is Ow = O(N ) for the Euclidean similarity measure and

383

Ow = O(L) for the anatomical similarity measure, where L ≤ N the number

384

of different anatomical labels that are neighbors of a streamline in a single

385

direction, assuming that computations for neighbors in different directions

386

are parallelized. The Euclidean distance measure can benefit from a substan-

387

tial speed-up when all streamlines have the same number of points N , which

388

is why we chose to downsample streamlines to a constant N here. However,

389

the complexity of our anatomical similarity measure does not depend on

390

whether streamlines have equal lengths or not. Furthermore, our implemen-

391

tation saves histograms of anatomical neighbors using the unordered map

392

class of the standard C++ library, where the complexity of element search

393

is on average O(1).

394

3. Results

395

3.1. Comparison with manual labeling

399

AC C

EP

TE

D

M AN U

SC

RI PT

379

400

all subjects for each WM bundle, when the total number of clusters is 200.

401

The anatomical similarity measure performed 20% better than the Euclidean

396

397

398

Fig. 4(a) shows the Dice coefficient between the manually labeled WM

bundles and the streamline clusters obtained with each of the two similarity measures, averaged over all 18 bundles and 32 subjects, as a function of the total number of clusters. Fig. 4(b) shows the average Dice coefficient over

20

ACCEPTED MANUSCRIPT

similarity measure in terms of its agreement with bundles defined by a human

403

rater. A repeated measures analysis of variance, with factors of similarity

404

measure and number of clusters, showed a statistically significant effect of

405

similarity measure on the Dice coefficient (p < 0.0001).

M AN U

SC

RI PT

402

(a)

(b)

D

Figure 4: (a) Dice coefficient between the manually labeled WM bundles and the streamline clusters obtained with each of the two similarity measures, averaged over all 18 bundles

TE

and 32 subjects, as a function of the total number of clusters. (b) Average Dice coefficient over all subjects by tract, when the total number of clusters is 200.

Fig. 5 and Fig. 6 depict the clusters used in the above evaluation, aver-

407

aged across the 32 subjects in the space of the FMRIB58 FA template (Smith

408

et al., 2006) and displayed as isosurfaces. (Note that we mapped clusters to

410

411

AC C

409

EP

406

this template only for the purposes of this illustration. Otherwise, all computations were performed in the native dMRI space of each individual.) For reference, the figures also show population averages of the corresponding

412

manually labeled bundles. Although individual variability and finer differ-

413

ences between methods cannot be appreciated from these population aver-

414

ages, the figures suggest that the 5% threshold that we use to select the 21

ACCEPTED MANUSCRIPT

415

clusters included in this evaluation is reasonable. Fig. 7 shows plots of homogeneity (a) and completeness (b) for both

417

the anatomical and Euclidean similarity measure. Our anatomical similarity

418

measure outperforms the Euclidean similarity measure in both homogeneity

419

and completeness. A repeated measures analysis of variance, with factors of

420

similarity measure and number of clusters, showed a statistically significant

421

effect of similarity measure on both homogeneity and completeness (p <

422

0.0001).

SC

RI PT

416

We also evaluated a combined similarity measure, summing the anatom-

424

ical and Euclidean similarities after normalizing them to account for their

425

different variances. We did not find any statistically significant difference

426

between the performance of the combined similarity measure and the anatom-

427

ical similarity itself (p = 0.85 for Dice coefficient, p = 0.60 for homogeneity,

428

and p = 0.72 for completeness, based on a repeated measures analysis of

429

variance with factors of similarity measure and number of clusters). This

430

suggests that, once the anatomical neighborhood of the streamlines has been

431

taken into account, the spatial coordinates of the streamlines may not con-

432

tribute additional information to the clustering. Hence, no further results

433

are shown for the combination of similarity measures.

435

436

437

D

TE

EP

Finally, we evaluated clustering the streamlines after mapping them to

AC C

434

M AN U

423

the subject’s sMRI space, instead of mapping the cortical parcellation and subcortical segmentation labels to dMRI space. We did not find statistically significant differences between the performance of the anatomical similarity

438

measure in sMRI and dMRI space (p = 0.62 for Dice coefficient, p = 0.55 for

439

homogeneity, and p = 0.9 for completeness, based on a repeated measures

22

ACCEPTED MANUSCRIPT

L-ATR

L-CCG

M AN U

SC

Anatomical

AC C

EP

TE

D

Euclidean Euclidean Anatomical Manual

L-CAB

RI PT

FMIN

Manual

FMAJ

Figure 5: Population averages of the manually labeled bundles and the streamline clusters that contain more than 5% streamlines overlapping with the corresponding manually labeled bundle. Each color represents a different WM pathway. The population averages across all subjects are shown as isosurfaces in axial and sagittal views.

23

ACCEPTED MANUSCRIPT

L-ILF

L-SLFP

Manual

L-UNC

M AN U

SC

Anatomical

AC C

EP

TE

D

Euclidean Euclidean Anatomical Manual

L-SLFT

RI PT

L-CST

Figure 6: Population averages of the manually labeled bundles and the streamline clusters that contain more than 5% streamlines overlapping with the corresponding manually labeled bundle (continued from Fig. 5). Each color represents a different WM pathway. The population averages across all subjects are shown as isosurfaces in axial and sagittal views.

24

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

M AN U

(a)

Figure 7: Homogeneity (a) and completeness (b) of unsupervised clustering, averaged over 18 WM bundles and 32 subjects, as a function of the number of clusters.

analysis of variance with factors the number of clusters and the space that

441

the clustering was performed in). Throughout the paper, the anatomical

442

similarity measure is computed in dMRI space.

443

3.2. Anatomical and spatial consistency of clusters

TE

D

440

Fig. 8(a) shows plots of wE for clusters derived by optimizing either wE

445

or wA . Fig. 8(b) shows plots of wA for clusters derived by optimizing either

446

wE or wA . Finally, Fig. 8(c) shows plots of the mean closest-point distance,

447

wCP , for clusters derived by optimizing either wE or wA .

451

AC C

EP

444

452

of clusters (Fig. 8(a)). This is confirmed by a repeated measures analy-

453

sis of variance on these differences, with factors of similarity measure and

448

449

450

Of course we expect that the wE -optimal clusters have higher wE than the

wA -optimal clusters, and vice versa. This is confirmed in the plots. However, we find that the difference in wA between the wE -optimal and wA -optimal clusters (Fig. 8(b)) is greater than the difference in wE between the two sets

25

ACCEPTED MANUSCRIPT

number of clusters (p < 10−8 ). Furthermore, there is no difference in the

455

within-cluster mean closest-point distance between the two sets of clusters

456

(Fig. 8(c)). This is confirmed by a repeated measures analysis of variance

457

on wCP , with factors of similarity measure and number of clusters (p = 0.3).

458

These results imply that the anatomical similarity measure tends to group

459

together streamlines that are also close to each other in Euclidean space,

460

whereas the Euclidean distance similarity measure may not necessarily yield

461

clusters that are anatomically consistent.

M AN U

SC

RI PT

454

(b)

(c)

D

(a)

TE

Figure 8: Euclidean similarity measure wE (a), anatomical similarity measure wA (b), and mean closest-point distance wCP (c), for streamlines clustered together either by the Euclidean or the anatomical similarity measure (wE -optimal vs. wA -optimal). All measures

EP

are plotted as a function of the number of clusters, averaged over all subjects.

A visual illustration of this is shown in Fig. 9. For four pairs of anatomical

463

segmentation labels, we show the clusters for which at least 5% of streamlines

466

AC C

462

467

For example, streamlines that lie on opposite sides of the midline but are

468

close to each other in space may be erroneously clustered together by wE ,

464

465

pass through both labels, when the number of clusters is 200. As seen in the figure, the Euclidean distance similarity measure tends to produce noisier and less anatomically consistent clusters than the anatomical similarity measure.

26

ACCEPTED MANUSCRIPT

RI PT

Euclidean similarity

M AN U

SC

Anatomical similarity

(a)

(b)

(c)

(d)

Figure 9: Clusters from a single subject, obtained with each of the two similarity measures.

D

Clusters were selected so that at least 5% of the streamlines in a cluster pass through a pair of anatomical segmentation labels: (a) precentral and brainstem, (b) superior parietal and

cingulate.

TE

brainstem, (c) superior temporal and precentral, (d) isthmus cingulate and rostral anterior

but not by wA (see Fig. 9(d)).

470

3.3. Robustness to errors in the anatomical segmentation

474

AC C

EP

469

475

unsupervised clustering, when the original and perturbed segmentations are

476

used to compute the anatomical similarity measure. Based on repeated mea-

471

472

473

Fig. 10 shows examples of an original and perturbed segmentation that we

used to assess the robustness of the clustering to errors in the boundaries of the anatomical segmentation labels. Fig. 11 shows plots of the Dice coefficient with manually labeled bundles, the homogeneity, and the completeness of the

27

RI PT

ACCEPTED MANUSCRIPT

SC

Figure 10: Example of an original (left) and perturbed (right) indivudual anatomical

M AN U

segmentation used to evaluate robustness to segmentation errors.

(b)

(c)

D

(a)

Figure 11: Dice coefficient with respect to manually labeled bundles (a), homogeneity (b),

TE

and completeness (c) of unsupervised clustering with our anatomical similarity measure, when the subjects’ original and perturbed anatomical segmentations are used.

sures analyses of variance with factors of segmentation type and number of

478

clusters, the effect of the segmentation on the performance measures was not

479

statistically significant (p = 0.79 for Dice coefficient, p = 0.35 for homogene-

481

482

AC C

480

EP

477

ity, and p = 0.36 for completeness). This suggests that our results are robust

to errors in the individual anatomical segmentation. The perturbations that we performed for the above evaluation were in the

483

order of one dMRI voxel. For reference, Fig. 12 shows the accuracy of the au-

484

tomated FreeSurfer segmentation for a set of internal brain structures. The

28

ACCEPTED MANUSCRIPT

plots show the mean closest-point distance between the automatically seg-

486

mented and manually labeled structures. The mean error over all structures

487

is 0.9mm, which is less than the sMRI resolution, and much less than the

488

dMRI resolution. Hence, typical FreeSurfer segmentation errors are smaller

489

than the perturbations that we performed in our experiment.

490

3.4. Effect of prototype streamlines

SC

RI PT

485

Supplemental Fig. S1 shows how the number of prototype streamlines

492

affects the Dice coefficient with manually labeled bundles, the homogeneity,

493

the completeness, and the within-cluster anatomical similarity. As seen in

494

the plots, results obtained with the anatomical similarity measure are very

495

similar for the different numbers of prototype streamlines (50, 100, 250, 500,

496

1000, 3000). The anatomical similarity measure with 50 prototype stream-

497

lines still outperforms the Euclidean distance similarity measure with 500

498

prototype streamlines (p = 0.002 for Dice and p < 0.0001 for homogene-

499

ity and completeness, based on a repeated measures analysis of variance

500

with factors of similarity measure and number of clusters). Supplemental

501

Fig. S2 shows results using 5 different randomly selected sets of 500 pro-

502

totype streamlines to compute the anatomical similarity measure. For all

503

sets, the anatomical similarity measure outperforms the Euclidean distance

505

506

D

TE

EP

AC C

504

M AN U

491

similarity measure with 500 prototype streamlines (p < 0.0001 based on a repeated measures analysis of variance with factors of similarity measure and number of clusters). Moreover, the difference between the results obtained

507

with the 5 random sets of prototype streamlines used for the anatomical

508

similarity measure was not significant (p > 0.6 for Dice, homogeneity and

509

completeness, based on a repeated measures analysis of variance with factors 29

ACCEPTED MANUSCRIPT

of prototype set and number of clusters, for each pair of sets).

511

3.5. Effect of anatomical neighborhood size

RI PT

510

Supplemental Fig. S3 shows how the size of the anatomical neighborhood

513

used for clustering affects the Dice coefficient with manually labeled bundles,

514

the homogeneity, the completeness, and the within-cluster anatomical sim-

515

ilarity. Although the 6-element neighborhood leads to somewhat degraded

516

performance, all three anatomical neighborhood sizes perform better than

517

the Euclidean distance similarity measure (p < 0.0001 based on repeated

518

measures analysis of variance with factors of similarity measure and number

519

of clusters).

520

3.6. Effect of streamline downsampling

M AN U

SC

512

Supplemental Fig. S4 shows how downsampling streamlines to N = 10

522

vs. N = 50 points affects the Dice coefficient with manually labeled bundles,

523

the homogeneity, the completeness, and the within-cluster anatomical simi-

524

larity. Although the plots show a modest improvement in the Dice coefficient

525

and homogeneity when N = 50 points are used, the difference is not signif-

526

icant (p = 0.1 for Dice coefficient, p = 0.2 for the homogeneity, p = 0.7 for

527

the completeness and p = 0.1 for the anatomical similarity evaluated with

529

530

TE

EP

AC C

528

D

521

repeated measures analysis of variance with factors of similarity measure and number of clusters). 3.7. Comparison to an alternative segmentation

531

Supplemental Fig. S5 shows results from performing clustering with our

532

anatomical similarity measure, after replacing the subjects’ individual FreeSurfer

30

ACCEPTED MANUSCRIPT

cortical parcellation and subcortical segmentation with the JHU-ICBM WM

534

atlas. We find that, for the purposes of unsupervised clustering of whole-brain

535

tractography with our anatomical similarity measure, using the individual’s

536

cortical and subcortical labels from FreeSurfer outperforms using labels from

537

a WM atlas (p = 0.0007 for Dice coefficient, and p = 0.0001 for homogeneity

538

and completeness, based on a repeated measures analysis of variance with

539

factors of label set and number of clusters).

540

3.8. Comparison to shape similarity

M AN U

SC

RI PT

533

Results from the comparison of our anatomical similarity measure to a

542

shape similarity measure are shown in Supplemental Fig. S6. The anatomical

543

similarity measure outperformed the shape similarity measure (p = 0.0003 for

544

Dice coefficient, and p < 0.0001 for homogeneity and completeness, based on

545

a repeated measures analysis of variance with factors of similarity measure

546

and number of clusters). As seen in the figure, the anatomical similarity

547

measure yielded a 30% improvement of the Dice coefficient between clusters

548

and manually labeled WM bundles, when compared to the shape similarity

549

measure.

550

3.9. Computation time

554

AC C

EP

TE

D

541

555

putations of anatomical similarity for different neighborhood elements are

556

parallelized. Otherwise, the computation time for the anatomical similarity

551

552

553

Table 1 shows computation times for clustering a randomly selected sub-

ject’s data with the anatomical and Euclidean distance similarity measure, for six different numbers of prototype streamlines and two different numbers of points per streamline. These times are obtained by assuming that com-

31

ACCEPTED MANUSCRIPT

Number of prototype streamlines

50

100

250

500

1000

3000

9.09

21.37 60.07 243.46 1951.30

Euclidean similarity (N = 10) 2.45

4.16

10.14 41.26 282.30 1878.08

RI PT

Anatomical similarity (N = 10) 5.31

Anatomical similarity (N = 50) 8.46 12.57 27.44 74.16 260.20 2347.28 Euclidean similarity (N = 50) 9.20 15.26 36.03 93.29 321.49 2392.37

SC

Table 1: Computation times (in minutes) for clustering a randomly selected subject’s

data with the two similarity measures and different numbers of prototype streamlines or numbers of points N , assuming that computations of anatomical similarity for different

M AN U

neighborhood elements are parallelized.

measure would also depend on the size of the anatomical neighborhood that

558

is used. Computation is faster for the Euclidean similarity measure when the

559

number of sampled points is small, but it is faster for the anatomical similar-

560

ity measure when the number of sampled points is large. Note, however, that

561

it would be possible to accelerate computation for the Euclidean similarity

562

measure by parallelizing over the sampled points. Times are reported for a

563

quad-core Xeon 5472 with 3.0GHz CPUs and 7GB of RAM.

564

3.10. Visual illustration

EP

TE

D

557

Fig. 13 shows an example of the hierarchical tree generated with the

566

anatomical similarity measure for a randomly selected subject. The root of

567

568

569

AC C

565

the tree is the full set of streamlines and the leaves of the tree are the 200 streamline clusters. For a set of nodes that are marked (a)-(i) on the tree, we show images of the clusters at the corresponding iterations of the algorithm:

570

• Low-level nodes represent coarse anatomical divisions of the stream-

571

lines: (a) left hemisphere; (b) posterior components of the corpus cal32

ACCEPTED MANUSCRIPT

losum; (c) right hemisphere; (d) anterior and middle components of the

573

corpus callosum.

RI PT

572

• Mid-level nodes represent streamline groupings with greater anatomical

575

specificity: (e) left corticospinal tract and thalamic radiations; (f) left

576

superior longitudinal fasciculus and arcuate fasciculus; (g) left cingulum

577

bundle.

579

580

581

• High-level nodes are individual WM pathways: (h) the right corticospinal tract.

M AN U

578

SC

574

• Leaves (terminal-level nodes) are subdivisions of larger WM pathways: (i) sub-bundle of the left arcuate fasciculus.

In the supplemental material, we provide a video with an example of the

583

tree hierarchy for the anatomical similarity measure. The video shows the left

584

hemisphere clusters, including the subdivisions of the superior longitudinal

585

fasciculus and the arcuate fasciculus.

586

4. Discussion

EP

TE

D

582

In this work we propose a novel anatomical similarity measure for un-

588

supervised clustering of tractography streamlines. Our similarity measure

591

AC C

587

592

anatomical similarity measure into a divisive hierarchical clustering algo-

593

rithm, normalized cuts, and compare its performance to that of a conven-

589

590

does not use the spatial coordinates of the streamlines, relying instead on histograms of the anatomical segmentation labels that each streamline goes through or next to, at every point along its trajectory. We incorporate this

33

ACCEPTED MANUSCRIPT

tional similarity measure based on the Euclidean distance between stream-

595

lines.

RI PT

594

Our results show that streamline clusters obtained with our anatomical

597

similarity measure are in better agreement with WM bundles labeled manu-

598

ally by a trained rater. Specifically, there is a 20% improvement in the overlap

599

with the manually labeled bundles when streamlines where clustered based

600

on their anatomical similarity than their Euclidean distance (Fig. 4). Visual

601

inspection of average clusters across subjects suggests that, when streamlines

602

are clustered based on Euclidean distance, clusters that overlap with these

603

manually labeled bundles include streamlines that are close to the bundles

604

but that would not be assigned to these bundles by a human rater (Fig. 5,

605

Fig. 6). Furthermore, we find that the anatomical similarity measure leads

606

to an increase in both the homogeneity and the completeness of clusters with

607

respect to the manually labeled bundles, when compared to the Euclidean

608

distance-based similarity (Fig. 7). Importantly, the increased agreement be-

609

tween clusters and manually labeled bundles is achieved without introducing

610

any prior information from manual labeling. Our method is unsupervised, us-

611

ing only the subject’s own cortical and subcortical segmentation labels. The

612

improvement stems from the fact that, by not grouping streamlines based on

613

their proximity but based on their anatomical neighborhood, the algorithm

615

616

M AN U

D

TE

EP

AC C

614

SC

596

follows an approach that resembles the anatomical criteria used by human raters to identify WM bundles. We also evaluate the anatomical and spatial consistency of clusters ob-

617

tained by optimizing the anatomical or Euclidean distance similarity mea-

618

sure (Fig. 8). We find that the clusters obtained with the two methods differ

34

ACCEPTED MANUSCRIPT

significantly more in terms of within-cluster anatomical similarity than they

620

differ in terms of within-cluster Euclidean distance. Furthermore, the within-

621

cluster mean closest-point distance is not significantly different between the

622

clusters obtained with the two methods. This suggests that streamlines that

623

share common anatomical neighbors tend to also be close together in space,

624

but the reverse is not necessarily true (see also Fig.9). It has been pointed

625

out previously that clustering based on spatial coordinates cannot capture

626

all anatomical subdivisions, if the streamlines in the subdivisions do not vary

627

significantly in terms of shape (O’Donnell et al., 2006).

M AN U

SC

RI PT

619

We assess the robustness of our anatomical similarity measure to errors

629

in the anatomical segmentation by repeating the clustering after perturbing

630

the boundaries of the segmentation labels. We find that there is no signif-

631

icant deterioration in the performance of the anatomical similarity measure

632

when the perturbed segmentations are used instead of the original segmen-

633

tations (Fig. 11). We also show that errors in the automated anatomical

634

segmentation labels obtained by FreeSurfer are smaller than the range of the

635

perturbations performed in our experiment, and smaller than typical dMRI

636

voxel sizes (Fig. 12).

EP

TE

D

628

We have chosen to use a randomly selected set of prototype streamlines

638

to compute the similarity matrix at each iteration of the normalized cuts

639

640

641

AC C

637

algorithm. This makes computation tractable for the large data sets that result from high-resolution dMRI acquisitions such as the one used here. The number of prototype streamlines has a modest impact on the results,

642

and does not change the conclusions of our performance comparison between

643

similarity measures (Supplemental Fig. S1 and Supplemental Fig. S2). We

35

ACCEPTED MANUSCRIPT

conclude that a choice of Np = 500 prototype streamlines is sensible. Note

645

that the number of streamlines to be divided with normalized cuts decreases

646

at each iteration of the algorithm. At the final stages of the clustering (as

647

we get further away from the root of the hierarchical tree of Fig. 13, thus

648

the clusters become finer and greater accuracy becomes crucial), Np = 500

649

represents a large portion of the total number of streamlines to be divided.

SC

RI PT

644

The number of anatomical neighbors included in our proposed similarity

651

measure is another factor that affects computational complexity, although

652

histogram computations for different neighbors are fully parallelizable. We

653

find a modest improvement in performance when anatomical neighborhoods

654

of size greater than 6 are used (Supplemental Fig. S3), but the anatomical

655

similarity measure outperforms the Euclidean distance similarity measure

656

for any neighborhood size. Finally, for all the results shown here, only long-

657

range connections (streamlines longer than 55mm) are included to reduce

658

computation. Our approach can be applied to full tractography data sets

659

without a reduction in performance at the expense of longer computation

660

times for both similarity measures.

TE

D

M AN U

650

Our evaluation of the anatomical similarity measure after downsampling

662

streamlines to N = 10 vs. N = 50 equispaced points shows that increas-

663

ing this number may improve results somewhat (Supplemental Fig. S4),

665

666

AC C

664

EP

661

although this improvement is not statistically significant. Here we downsampled streamlines to a fixed N for a more balanced comparison to the Euclidean distance similarity measure, which relies on downsampling to re-

667

duce computation time. However, our proposed anatomical similarity mea-

668

sure does not require streamlines to have the same number of points, and its

36

ACCEPTED MANUSCRIPT

669

computational complexity is not impacted by this number. We have chosen to compute the anatomical similarity of streamlines with

671

respect to a set of cortical parcellation and subcortical segmentation labels,

672

instead of labels from a WM atlas, for two reasons. First, our goal was to

673

develop an unsupervised method, rather than impose the existence of certain

674

WM bundles through the use of an atlas. Our approach allows exploratory

675

analyses to be performed on tractography data from the whole brain, instead

676

of being limited to a set of pre-defined bundles. Second, our similarity mea-

677

sure relies not only on the structures that a streamline goes through, but also

678

on the structures that surround it, in a set of different directions. Therefore,

679

one might expect a segmentation that includes cortical and subcortical struc-

680

tures to perform well in that setting. Indeed, our results suggest that our

681

anatomical similarity measure works better when we compute it using labels

682

of cortical and subcortical structures, rather than labels from a WM atlas

683

(Supplemental Fig. S5). That is, although a WM atlas can be an extremely

684

valuable tool for supervised reconstruction of specific pathways, it does not

685

appear to be optimal for the purposes of whole-brain tractography clustering

686

based on anatomical similarity.

EP

TE

D

M AN U

SC

RI PT

670

The similarity measure that we have proposed involves the histograms

688

of the anatomical neighbors of tractography streamlines. These histograms

689

690

691

AC C

687

provide an “anatomical signature” that allows us to describe streamlines without relying on their spatial coordinates in a native or template space. We expect this signature to be useful for finding corresponding bundles, e.g.,

692

between subjects or between hemispheres, without the need for an exact

693

point-by-point transformation between them. This is an exciting prospect

37

ACCEPTED MANUSCRIPT

for large population studies, and it will be very important to evaluate our

695

method on larger data sets than the one used here. Early results from apply-

696

ing AnatomiCuts to a much larger set of data from the MGH-Harvard-USC

697

pilot Lifespan Human Connectome Project show promise for using the pro-

698

posed anatomical similarity measure to find corresponding clusters across

699

subjects aged 8–90 without co-registration (Siless et al., 2017). We intend to

700

investigate this further and to evaluate our method on other large data sets

701

from different disease populations in the future.

SC

RI PT

694

In previous work, we used a similar description of the anatomical neigh-

703

borhood of a path to introduce priors into a global probabilistic tractography

704

algorithm (Yendiki et al., 2011). In that work, histograms of the anatomical

705

neighbors of tractography streamlines were computed from training subjects,

706

where the bundles of interest had been labeled manually, and used as anatom-

707

ical priors when performing tractography in a novel subject. In the present

708

work, we do not rely on training data; we introduce anatomical informa-

709

tion only from the subject’s own anatomical segmentation, and only after

710

the tractography has been performed. On the one hand, this allows us to

711

perform exploratory, whole-brain tractography analyses, without being con-

712

strained to a set of bundles included in an atlas. On the other hand, because

713

the tractography is performed in an unconstrained way, it is prone to all the

715

716

D

TE

EP

AC C

714

M AN U

702

usual errors due to crossing WM bundles, such as streamlines ending prematurely or jumping from one bundle to another. Post hoc clustering methods, like the one presented here, cannot recover from such errors, as they do not

717

modify the tractography streamlines—they only assign the streamlines into

718

groups. It would be possible to use the two approaches in a complementary

38

ACCEPTED MANUSCRIPT

way. One could apply the unsupervised clustering presented here to high-

720

quality dMRI data (e.g., from the Connectom scanner) to generate clusters

721

from the whole brain, include the clusters that occur consistently across sub-

722

jects into an atlas to generate priors for the approach in Yendiki et al. (2011),

723

and then apply the latter to reconstruct the same pathways robustly in more

724

routine-quality dMRI data.

725

5. Conclusion

M AN U

SC

RI PT

719

We present AnatomiCuts, a method for unsupervised hierarchical cluster-

727

ing of dMRI tractography data based on anatomical similarity. We compare

728

this to the conventional approach of using a similarity based on Euclidean

729

distance. We find that the anatomical similarity yields results that are more

730

consistent with manual labeling relative to the Euclidean similarity. That

731

is, without introducing any training data from human raters, we are able to

732

obtain results that are in closer agreement with such a rater. We achieve this

733

simply by using a similarity measure that is better at replicating how a hu-

734

man with neuroanatomical expertise would segment WM tracts, i.e., based

735

on the anatomical structures that they either intersect or neighbor, every-

736

where along the tracts’ trajectory. This allows us to obtain anatomically

738

739

TE

EP

AC C

737

D

726

meaningful WM bundles without being limited to a set of tracts included in an atlas. Therefore we expect our approach to facilitate fully unsupervised analyses of whole-brain tractography data in large population studies.

39

ACCEPTED MANUSCRIPT

740

Acknowledgements Support for this research was provided in part by the National Institute

742

for Biomedical Imaging and Bioengineering (P41EB015896, 1R01EB023281,

743

R01EB006758, R21EB018907, R01EB019956), the National Institute for Men-

744

tal Health (U01-MH108168; Boston Adolescent Neuroimaging of Depres-

745

sion and Anxiety project), the National Institute on Aging (R01AG008122,

746

R01AG016495), the National Institute of Diabetes and Digestive and Kidney

747

Diseases (R21DK108277), the National Institute for Neurological Disorders

748

and Stroke (R01NS0525851, R21NS072652, R01NS070963, R01NS083534,

749

U01NS086625), and was made possible by the resources provided by Shared

750

Instrumentation Grants 1S10RR023401, 1S10RR019307, and 1S10RR023043.

751

Additional support was provided by the NIH Blueprint for Neuroscience Re-

752

search (U01-MH093765; part of the multi-institutional Human Connectome

753

Project, T90DA022759/R90DA023427). In addition, BF has a financial in-

754

terest in CorticoMetrics, a company whose medical pursuits focus on brain

755

imaging and measurement technologies. BF’s interests were reviewed and

756

are managed by Massachusetts General Hospital and Partners HealthCare in

757

accordance with their conflict of interest policies.

758

Basser, P.J., Pajevic, S., Pierpaoli, C., Duda, J., & Aldroubi, A. 2000. In vivo

760

761

762

SC

M AN U

D

TE

EP

AC C

759

RI PT

741

fiber tractography using DT-MRI data. Magnetic Resonance in Medicine, 44(4), 625–632.

Brun, A., Knutsson, H., Park, H.J., Shenton, M.E., & Westin, C.F. 2004. Clustering Fiber Traces Using Normalized Cuts. International Confer-

40

ACCEPTED MANUSCRIPT

ence on Medical Image Computing and Computer-Assisted Intervention,

764

3216/2004(3216), 368–375.

RI PT

763

Corouge, I., Gouttard, S., & Gerig, G. 2004 (April). Towards a shape model

766

of white matter fiber bundles using diffusion tensor MRI. Pages 344–

767

347 Vol. 1 of: 2004 2nd IEEE International Symposium on Biomedical

768

Imaging: Nano to Macro (IEEE Cat No. 04EX821).

770

Dice, L.R. 1945. Measures of the Amount of Ecologic Association Between Species. Ecology, 26(3), 297–302.

M AN U

769

SC

765

771

Ding, Z., Gore, J.C., & Anderson, A.W. 2003. Classification and quantifi-

772

cation of neuronal fiber pathways using diffusion tensor MRI. Magnetic

773

Resonance in Medicine, 49(4), 716–721.

Fan, Q., Witzel, T., Nummenmaa, A., Van Dijk, K.R.A., Van Horn, J.D.,

775

Drews, M.K., Somerville, L.H., Sheridan, M.A., Santillana, R.M., Snyder,

776

J., Hedden, T., Shaw, E.E., Hollinshead, M.O., Renvall, V., Zanzonico, R.,

777

Keil, B., Cauley, S., Polimeni, J.R., Tisdall, D., Buckner, R.L., Wedeen,

778

V.J., Wald, L.L., Toga, A.W., & Rosen, B.R. 2016. MGH-USC Human

779

Connectome Project datasets with ultra-high b-value diffusion MRI. Neu-

780

roImage, 124, 1108–1114.

782

783

TE

EP

AC C

781

D

774

Fischl, B., Salat, David H., Busa, E., Albert, M., Dieterich, M., Haselgrove, C., Van Der Kouwe, A., Killiany, R., Kennedy, D., Klaveness, S., Montillo, A., Makris, N., Rosen, B., & Dale, A.M. 2002. Whole brain segmentation:

784

Automated labeling of neuroanatomical structures in the human brain.

785

Neuron, 33(3), 341–355. 41

ACCEPTED MANUSCRIPT

Fischl, B., Van Der Kouwe, A., Destrieux, C., Halgren, E., S´egonne, F., Salat,

787

D.H., Busa, E., Seidman, L.J., Goldstein, J., Kennedy, D., Caviness, V.,

788

Makris, N., Rosen, B., & Dale, A.M. 2004. Automatically Parcellating the

789

Human Cerebral Cortex. Cerebral Cortex, 14(1), 11–22.

RI PT

786

Garyfallidis, E., Brett, M., Correia, M.M., Williams, G.B., & Nimmo-Smith,

791

I. 2012. QuickBundles, a Method for Tractography Simplification. Fron-

792

tiers in neuroscience, 6(December), 175.

SC

790

Gerig, G., Gouttard, S., & Corouge, I. 2004. Analysis of brain white matter

794

via fiber tract modeling. Pages 4421–4424 of: Engineering in Medicine and

795

Biology Society, 2004. IEMBS ’04. 26th Annual International Conference

796

of the IEEE, vol. 2.

M AN U

793

Golub, G.H., & Van Loan, C.F. 1996. Matrix Computations.

798

Greve, D.N., & Fischl, B. 2009. Accurate and robust brain image alignment using boundary-based registration. NeuroImage, 48(1), 63–72.

TE

799

D

797

Guevara, P., Poupon, C., Rivi`ere, D., Cointepas, Y., Descoteaux, M.,

801

Thirion, B., & Mangin, J.F. 2011. Robust clustering of massive trac-

802

tography datasets. NeuroImage, 54(3), 1975–1993.

806

AC C

EP

800

807

Guha, S., Rastogi, R., & Shim, K. 2001. CURE: An efficient clustering

803

804

805

808

Guevara, P., Duclap, D., Poupon, C., Marrakchi-Kacem, L., Fillard, P., Le Bihan, D., Leboyer, M., Houenou, J., & Mangin, J.F. 2012. Automatic fiber bundle segmentation in massive tractography datasets using a multisubject bundle atlas. NeuroImage, 61(4), 1083–1099.

algorithm for large databases. Information Systems, 26(1), 35–58. 42

ACCEPTED MANUSCRIPT

Hua, K., Zhang, J., Wakana, S., Jiang, H., Li, X., Reich, D.S., Calabresi,

810

P. A., Pekar, J.J., Van Zijl, P.C.M., & Mori, S. 2008. Tract probability

811

maps in stereotaxic spaces: Analyses of white matter anatomy and tract-

812

specific quantification. NeuroImage, 39(1), 336 – 347.

814

Jenkinson, M., Beckmann, C.F., Behrens, T.E.J., Woolrich, M.W., & Smith, S.M. 2012. FSL.

SC

813

RI PT

809

Jin, Y., Shi, Y., Zhan, L., Gutman, B.A, de Zubicaray, G.I., McMahon,

816

K.L., Wright, M.J., Toga, A.W., & Thompson, P.M. 2014. Automatic

817

clustering of white matter fibers in brain diffusion MRI with an application

818

to genetics. NeuroImage, 100, 75–90.

M AN U

815

Le Bihan, D., Breton, E., Lallemand, D., Grenier, P., Cabanis, E., & Laval-

820

Jeantet, M. 1986. MR imaging of intravoxel incoherent motions: appli-

821

cation to diffusion and perfusion in neurologic disorders. Radiology, 161,

822

401–407.

TE

D

819

Lehman, J.F., Greenberg, B.D., McIntyre, C.C., Rasmussen, S.A., & Haber,

824

S.N. 2011. Rules ventral prefrontal cortical axons use to reach their tar-

825

gets: implications for diffusion tensor imaging tractography and deep brain

826

stimulation for psychiatric illness. Journal of Neuroscience, 31(28), 10392–

828

829

830

AC C

827

EP

823

10402.

Li, H., Xue, Z., Guo, L., Liu, T., Hunter, J., & Wong, S.T.C. 2010. A hybrid approach to automatic clustering of white matter fibers. NeuroImage, 49(2), 1249 – 1258.

43

ACCEPTED MANUSCRIPT

Maddah, M., Z¨ollei, L., Grimson, W.E.L., Westin, C.F., & Wells, W.M. 2008.

832

A mathematical framework for incorporating anatomical knowledge in DT-

833

MRI analysis. Pages 105–108 of: 2008 5th IEEE International Symposium

834

on Biomedical Imaging: From Nano to Macro, Proceedings, ISBI.

RI PT

831

Makris, N., & Pandya, D.N. 2008. The extreme capsule in humans and

836

rethinking of the language circuitry. Brain Structure and Function, 213(3),

837

343–358.

SC

835

Moberts, B., Vilanova, A., & Van Wijk, J.J. 2005 (Oct). Evaluation of

839

fiber clustering methods for diffusion tensor imaging. Pages 65–72 of:

840

Visualization, 2005. VIS 05. IEEE.

M AN U

838

Mori, S., Crain, B.J., Chacko, V.P., & Zijl, P.V. 1999. Three-dimensional

842

tracking of axonal projections in the brain by magnetic resonance imaging.

843

Ann. Neurol., 45(2), 265–9.

TE

D

841

Mugler, J.P., & Brookeman, J.R. 1990. Three-dimensional magnetization-

845

prepared rapid gradient-echo imaging (3D MP RAGE). Magnetic Reso-

846

nance in Medicine, 15(1), 152–157.

848

849

850

O’Donnell, L., Kubicki, M., Shenton, M. E., Dreusicke, M., Grimson, W.E.L.,

AC C

847

EP

844

& Westin, C.F. 2006. A Method for Clustering White Matter Fiber Tracts. AJNR, 27(5), 1032–1036.

O’Donnell, L.J., & Westin, C.F. 2007. Automatic tractography segmenta-

851

tion using a high-dimensional white matter atlas. IEEE Transactions on

852

Medical Imaging, 26(11), 1562–1575.

44

ACCEPTED MANUSCRIPT

854

Oishi, K., Faria, A.V., Van Zijl, P.C.M., & Mori, S. 2010. MRI Atlas of Human White Matter. Academic Press.

RI PT

853

Ros, C., G¨ ullmar, D., Stenzel, M., Mentzel, H.J., & Reichenbach, J.R. 2013.

856

Atlas-guided cluster analysis of large tractography datasets. PLoS ONE,

857

8(12).

SC

855

Rosenberg, A., & Hirschberg, J. 2007. V-measure: A conditional entropy-

859

based external cluster evaluation measure. Pages 410–420 of: Proceedings

860

of the Joint Conference on Empirical Methods in Natural Language Pro-

861

cessing and Computational Natural Language (EMNLP-CoNLL’07), vol.

862

1.

864

Sammut, C., & Webb, G. I. (eds). 2017. Classification Algorithms. Boston, MA: Springer US. Pages 208–209.

D

863

M AN U

858

Setsompop, K., Kimmlingen, R., Eberlein, E., Witzel, T., Cohen-Adad, J.,

866

McNab, J.A., Keil, B., Tisdall, M.D., Hoecht, P., Dietz, P., Cauley, S.F.,

867

Tountcheva, V., Matschl, V., Lenz, V.H., Heberlein, K., Potthast, A.,

868

Thein, H., Van Horn, J., Toga, A., Schmitt, F., Lehne, D., Rosen, B.R.,

869

Wedeen, V., & Wald, L.L. 2013. Pushing the limits of in vivo diffusion

870

MRI for the Human Connectome Project. NeuroImage, 80, 220–233.

873

AC C

EP

TE

865

874

Siless, V., Medina, S., Varoquaux, G., & Thirion, B. 2013 (June). A Com-

871

872

Shi, J., & Malik, J. 2000. Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(8), 888– 905.

45

ACCEPTED MANUSCRIPT

parison of Metrics and Algorithms for Fiber Clustering. Pages 190–193 of:

876

2013 International Workshop on Pattern Recognition in Neuroimaging.

877

Siless, V., Chang, K., Fischl, B., & Yendiki, A. 2016. Hierarchical Clustering

878

of Tractography Streamlines Based on Anatomical Similarity. Pages 184–

879

191 of: Medical Image Computing and Computer-Assisted Intervention.

880

Springer International Publishing.

SC

RI PT

875

Siless, V., Davidow, J.Y., Nielsen, J., Fan, Q., Hedden, T., Hollinshead,

882

M., Bustamante, C.V., Drews, M.K., Van Dijk, K.R.A., Sheridan, M.A.,

883

Buckner, R. L., Fischl, B., Somerville, L., & Yendiki, A. 2017. Registration-

884

free analysis of diffusion MRI tractography data across subjects through

885

the human lifespan. In: Proc. Intl. Soc. Mag. Reson. Med., vol. 25.

M AN U

881

Smith, S.M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T.E.,

887

Mackay, C.E., Watkins, K.E., Ciccarelli, O., Cader, M.Z., Matthews, P.M.,

888

& Behrens, T.E.J. 2006. Tract-based spatial statistics: Voxelwise analysis

889

of multi-subject diffusion data. NeuroImage, 31(4), 1487–1505.

TE

D

886

Sørensen, T.J. 1948. A method of establishing groups of equal amplitude in

891

plant sociology based on similarity of species and its application to analyses

892

of the vegetation on Danish commons. Biol. Skr., 5, 1–34.

894

895

896

AC C

893

EP

890

Tsai, A., Westin, C.F., Hero, A.O., & Willsky, A.S. 2007. Fiber tract clustering on manifolds with dual rooted-graphs. In: Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition.

46

ACCEPTED MANUSCRIPT

Tunc, B., Parker, W.A., Ingalhalikar, M., & Verma, R. 2014. Automated

898

tract extraction via atlas based adaptive clustering. NeuroImage, 102(P2),

899

596–607.

RI PT

897

Van der Kouwe, A., Benner, T., Salat, D.H., & Fischl, B. 2008. Brain

901

morphometry with multiecho MPRAGE. NeuroImage, 40(2), 559–569.

902

Visser, E., Nijhuis, E.H.J., Buitelaar, J.K., & Zwiers, M.P. 2011. Partition-

903

based mass clustering of tractography streamlines. NeuroImage, 54(1),

904

303–312.

M AN U

SC

900

Wakana, S., Caprihan, A., Panzenboeck, M.M., Fallon, J.H., Perry, M., Gol-

906

lub, R.L., Hua, K., Zhang, J., Jiang, H., Dubey, P., Blitz, A., van Zijl,

907

P., & Mori, S. 2007. Reproducibility of quantitative tractography methods

908

applied to cerebral white matter. NeuroImage, 36(3), 630–644.

D

905

Wang, Q., Yap, P.T., Wu, G., & Shen, D. 2013. Application of neuroanatom-

910

ical features to tractography clustering. Human Brain Mapping, 34(9),

911

2089–2102.

EP

TE

909

Wang, X., Grimson, W.E.L, & Westin, C.F. 2011. Tractography segmenta-

913

tion using a hierarchical Dirichlet processes mixture model. NeuroImage,

914

915

916

AC C

912

54(1), 290–302.

Wassermann, D., Bloy, L., Kanterakis, E., Verma, R., & Deriche, R. 2010. Unsupervised white matter fiber clustering and tract probability map gen-

917

eration: Applications of a Gaussian process framework for white matter

918

fibers. NeuroImage, 51(1), 228–241.

47

ACCEPTED MANUSCRIPT

Wassermann, D., Makris, N., Rathi, Y., Shenton, M., Kikinis, R., Kubicki,

920

M., & Westin, C.F. 2016. The white matter query language: a novel

921

approach for describing human white matter anatomy. Brain Structure

922

and Function, 221(9), 4705–4721.

924

Wright, S.M., & Wald, L.L. 1997. Theory and application of array coils in MR spectroscopy. NMR in Biomedicine, 10(8), 394–410.

SC

923

RI PT

919

Wu, X., Xie, M., Zhou, J., Anderson, A.W., Gore, J.C., & Ding, Z. 2012.

926

Globally optimized fiber tracking and hierarchical clustering – a unified

927

framework. Magnetic resonance imaging, 30(4), 485–495.

M AN U

925

Xia, Y., Turken, U., Whitfield-Gabrieli, S.L., & Gabrieli, J.D. 2005.

929

Knowledge-Based Classification of Neuronal Fibers in Entire Brain. Berlin,

930

Heidelberg: Springer Berlin Heidelberg. Pages 205–212.

932

Yeh, F.C., Wedeen, V.J., & Tseng, W.Y.I. 2010. Generalized q-sampling

TE

931

D

928

imaging. IEEE Transactions on Medical Imaging, 29(9), 1626–1635. Yeh, F.C., Verstynen, T.D., Wang, Y., Fern´andez-Miranda, J.C., & Tseng,

934

W.Y.I. 2013. Deterministic diffusion fiber tracking improved by quantita-

935

tive anisotropy. PLoS ONE, 8(11).

937

938

AC C

936

EP

933

Yendiki, A., Panneck, P., Srinivasan, P., Stevens, A., Z¨ollei, L., Augustinack, J., Wang, R., Salat, D., Ehrlich, S., Behrens, T., Jbabdi, S., Gollub, R., & Fischl, B. 2011. Automated Probabilistic Reconstruction of White-

939

Matter Pathways in Health and Disease Using an Atlas of the Underlying

940

Anatomy. Frontiers in Neuroinformatics, 5, 23.

48

ACCEPTED MANUSCRIPT

Zhang, T., Chen, H., Guo, L., Li, K., Li, L., Zhang, S., Shen, D., Hu, X., &

942

Liu, T. 2014. Characterization of U-shape streamline fibers: Methods and

943

applications. Medical Image Analysis, 18(5), 795–807.

RI PT

941

Ziyan, U., Sabuncu, M.R., Grimson, W.E.L., & Westin, C.F. 2009. Consis-

945

tency clustering: A robust algorithm for group-wise registration, segmen-

946

tation and automatic atlas construction in diffusion MRI. Pages 279–290

947

of: International Journal of Computer Vision, vol. 85.

AC C

EP

TE

D

M AN U

SC

944

49

AC C

EP

TE

D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 12: Mean closest-point distance between automatically segmented and manually labeled brain structures that are included in the FreeSurfer subcortical segmentation.

50

D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

AC C

EP

TE

(a)

Figure 13: Top: Tree representation of the hierarchical cuts performed on a subject’s

51

whole-brain tractography data using our anatomical similarity measure. Bottom: Clusters associated with the tree nodes labeled (a)-(i).

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

TE

D

(a)

(d)

EP

(c)

Figure S1: Dice coefficient with respect to manually labeled bundles (a), homogeneity (b), completeness (c), and anatomical similarity (d) of unsupervised clustering, using different

AC C

numbers of prototype streamlines to compute the similarity matrix for normalized cuts. Results are shown for clustering with the anatomical similarity measure and 50, 100, 250, 500, 100 or 3000 prototype streamlines (A-50, A-100, A-250, A-500, A-1000, A-3000). Results for clustering with the Euclidean distance similarity measure and 500 prototype streamlines (E-500) are also shown for reference.

52

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

TE

D

(a)

(d)

EP

(c)

Figure S2: Dice coefficient with respect to manually labeled bundles (a), homogeneity

AC C

(b), completeness (c), and anatomical similarity (d) of unsupervised clustering, using randomly selected sets of 500 prototype streamlines to compute the similarity matrix for normalized cuts. Results are shown for clustering with the anatomical similarity measure and 5 different sets of 500 prototype streamlines (A-500). Results for clustering with the Euclidean distance similarity measure and a single set 500 prototype streamlines (E-500) are also shown for reference.

53

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

EP

TE

D

(a)

(a)

(b)

AC C

Figure S3: Dice coefficient with respect to manually labeled bundles (a), homogeneity (b), completeness (c), and anatomical similarity (d) of unsupervised clustering, using different sizes of anatomical neighborhoods. Results are shown for clustering with the anatomical similarity measure and a 6-, 14-, or 26-element neighborhood. Results for clustering with the Euclidean distance similarity measure are also shown for reference.

54

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

EP

TE

D

(a)

(c)

(d)

Figure S4: Dice coefficient with respect to manually labeled bundles (a), homogeneity

AC C

(b), completeness (c), and anatomical similarity (d) of unsupervised clustering, using the anatomical similarity measure on streamlines downsampled to N = 10 vs. N = 50 points. Results for clustering with the Euclidean distance similarity measure for streamlines downsampled to N = 10 points are also shown for reference.

55

TE

D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure S5: Dice coefficient between the manually labeled WM bundles and the streamline

EP

clusters obtained with the anatomical similarity measure and two different sets of labels: the subjects’ individual FreeSurfer cortical parcellation and subcortical segmentation la-

AC C

bels (FS) or the labels from the JHU-ICBM WM atlas, mapped into individual space (JHU). The plot shows Dice coefficients averaged over all subjects by tract, when the total number of clusters is 200.

56

D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

TE

Figure S6: Dice coefficient between the manually labeled WM bundles and the streamline clusters obtained with our proposed anatomical similarity and a similarity measure based on shape descriptors. The plot shows Dice coefficients averaged over all subjects by tract,

AC C

EP

when the total number of clusters is 200.

57