YNIMG-12389; No. of pages: 14; 4C: 2, 5, 8 NeuroImage xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

3Q2

Robert E. Smith a,⁎, Jacques-Donald Tournier a,b,c,d, Fernando Calamante a,b,e, Alan Connelly a,b,e

O

F

2

SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography

4 5 6 7 8

a

9

a r t i c l e

10 11 12 13

Article history: Received 9 February 2015 Accepted 30 June 2015 Available online xxxx

14 15 16 17 18 19 20

Keywords: Magnetic resonance imaging Streamlines Tractography Diffusion SIFT Structural connectivity

Imaging Division, Florey Institute of Neuroscience and Mental Health, Heidelberg, Victoria, Australia Department of Medicine, Austin Health and Northern Health, University of Melbourne, Melbourne, Victoria, Australia Centre for the Developing Brain, King’s College London, London, United Kingdom d Department of Biomedical Engineering, Division of Imaging Sciences and Biomedical Engineering, King’s College London, London, United Kingdom e Florey Department of Neuroscience and Mental Health, University of Melbourne, Melbourne, Victoria, Australia b

a b s t r a c t

P

i n f o

R O

c

T

E

D

Diffusion MRI streamlines tractography allows for the investigation of the brain white matter pathways noninvasively. However a fundamental limitation of this technology is its non-quantitative nature, i.e. the density of reconstructed connections is not reflective of the density of underlying white matter fibres. As a solution to this problem, we have previously published the “spherical-deconvolution informed filtering of tractograms (SIFT)” method, which determines a subset of the streamlines reconstruction such that the streamlines densities throughout the white matter are as close as possible to fibre densities estimated using the spherical deconvolution diffusion model; this permits the use of streamline count as a valid biological marker of connection density. Particular aspects of its performance may have however limited its uptake in the diffusion MRI research community. Here we present an alternative to this method, entitled SIFT2, which provides a more logically direct and computationally efficient solution to the streamlines connectivity quantification problem: by determining an appropriate cross-sectional area multiplier for each streamline rather than removing streamlines altogether, biologically accurate measures of fibre connectivity are obtained whilst making use of the complete streamlines reconstruction. © 2015 Published by Elsevier Inc.

C

1Q1

E

34 38 36 35

R

37

Introduction

40

Diffusion MRI can be used to probe the tissue microstructure of the brain in vivo. By using mathematical models of the measured diffusion signal, followed by digital reconstruction based on estimates of the white matter fibre orientations, this imaging technique is capable of providing estimates of the axon pathways in the white matter of the brain non-invasively (Tournier et al., 2011). This modality therefore has the potential to interrogate, characterize, and longitudinally follow, both the structural connectivity of the healthy brain and the interruption of this connectivity in various disease processes.

47 48

N C O

45 46

U

43 44

R

39

41 42

21 22 23 24 25 26 27 28 29 30 31 32 33

Abbreviations: aTV, asymmetrical total variation (a form of regularisation); ACT, anatomically-constrained tractography; CSD, constrained spherical deconvolution; FOD, fibre orientation distribution; FoV, field of view; GMWMI, grey matter–white matter interface (in reference to streamline seeding); SIFT, spherical-deconvolution informed filtering of tractograms; TDI, track density imaging; Tik, Tikhonov (a form of regularisation); WM, white matter (in reference to streamline seeding). ⁎ Corresponding author at: Florey Institute of Neuroscience and Mental Health, Melbourne Brain Centre, 245 Burgundy Street, Heidelberg, Victoria 3084, Australia. Fax: +61 3 9035 7301. E-mail address: robert.smith@florey.edu.au (R.E. Smith).

In order to provide robust and interpretable measures of connectivity, such a reconstruction should ideally provide metrics that are biologically relevant; possibly the most direct of which is the density of white matter connections between two regions of the brain grey matter. Such measures of connectivity are a prerequisite for the emerging field of connectomics (Sporns et al., 2005). However, an inherent limitation of the streamlines tractography reconstruction process is that the density of reconstructed connections is not indicative of the density of underlying axon connections (Jones et al., 2013). Therefore, many studies are limited to drawing measures from the diffusion model estimated in each voxel as surrogate markers of ‘integrity’ or ‘connectivity’; for instance, the mean fractional anisotropy (FA) of those voxels along the connecting pathway (e.g. Cercignani, 2010). Such measures are typically indirect, and may have limitations on biological interpretability based on the diffusion model used to estimate them. We have previously presented a solution to this quantification problem, as the “spherical-deconvolution informed filtering of tractograms (SIFT)” method (Smith et al., 2013). Through selective removal of individual streamlines from the reconstruction, an estimate of whole-brain structural connectivity is provided where the density of reconstructed connections is proportional to the fibre density within each image element as estimated by the diffusion model; the number of streamlines connecting two regions of grey matter therefore provides an estimate

http://dx.doi.org/10.1016/j.neuroimage.2015.06.092 1053-8119/© 2015 Published by Elsevier Inc.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71

2

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

1. The vast majority (at least 90%) of streamlines generated must be discarded to provide an adequate fit with the fibre density estimates, requiring a long execution time in order to generate the necessary initial dense reconstruction. 2. Due to the need to generate a very large number of streamlines initially in order to accommodate discarding a large proportion of them, the filtering algorithm itself may require dedicated highperformance computing hardware with a large amount of volatile memory, depending on the density of the initial reconstruction. 3. At completion of filtering, the number of remaining reconstructed connections may be too few, such that any subsequent analysis demonstrates poor test–retest reproducibility, or quantitative estimates of connectivity between regions are noisy due to the quantised nature of the streamlines reconstruction itself.

115

Methods

116

Synthetic data

117 118

121

For demonstrating and comparing the proposed method to the original SIFT method on data with a well-defined reconstruction bias, the synthetic phantom and streamlines reconstruction presented in Smith et al. (2013) was used; a slice of the phantom is shown in Fig. 1, with relevant details in the figure legend.

122

In vivo data acquisition

123 124

In vivo image data were acquired from a healthy volunteer on a Siemens Tim Trio 3 T MRI system (Erlangen, Germany). A diffusionweighted image series was acquired using a single-shot, twicerefocused spin-echo (Reese et al., 2003) echo-planar imaging sequence: field of view (FoV) 240 mm, matrix size 96 × 96, phase partial Fourier 6/ 8, parallel acceleration factor 2, 60 slices, 2.5 mm isotropic resolution, 60 diffusion-sensitisation directions at b = 3000 s mm− 2 with 7 b = 0 volumes, TE/TR = 110/8400 ms, 10 min acquisition time. An additional pair of b = 0 images (one with reversed phase encoding) was acquired

110 111 112 113

119 120

125 126 127 128 129 130 131

C

E

R

108 109

R

106 107

O

105

C

103 104

N

101 102

U

99 100

Fig. 1. The synthetic phantom used to demonstrate the operation of the SIFT and SIFT2 algorithms. A fibre orientation distribution (FOD) glyph is displayed at the centre of each voxel. Each of the bundles has a cross-section of 3 × 3 voxels. Importantly, the length of the upper bundle is three times that of the lower bundle (12 vs. 4 voxels). One thousand streamlines were generated by seeding randomly throughout the image.

for estimation of the inhomogeneity field (Andersson et al., 2003). A T1weighted anatomical contrast image for the purposes of tissue segmentation was acquired using the magnetisation-prepared rapid acquisition gradient echo (MPRAGE) sequence (Mugler and Brookeman, 1990) (FoV 230 × 230 mm, matrix size 256 × 256, 192 slices, 0.9 mm isotropic resolution, TE/TI/TR = 2.6/900/1900 ms, flip angle 9°, 8 min acquisition time).

132 133

Data pre-processing

139

T

114

In addressing these concerns, we have devised an alternative solution to the streamlines connectivity quantization problem, entitled SIFT2 (Smith et al., 2014). Instead of removing streamlines to converge the reconstructed streamlines densities toward the estimated fibre densities (as is done in the SIFT method), we instead determine an effective cross-sectional area uniquely for each streamline, optimizing these parameters such that the reconstructed fibre volumes match those estimated directly from the diffusion signal throughout the brain white matter. Importantly, because a streamlines reconstruction is given quantitative properties without necessitating the removal of streamlines, the user does not need to make a compromise between the total execution time (including both streamlines tractography and filtering) and the density of the modified reconstruction: SIFT2 makes use of all generated streamlines, ensuring that the underlying fibre orientation field is still maximally sampled once the algorithm has completed. In this article, we describe the operation of the SIFT2 algorithm, compare its performance to the original SIFT method on both synthetic and in vivo data, and evaluate the effects of regularization on the operation of the algorithm.

97 98

F

81

O

79 80

R O

78

P

76 77

D

74 75

of the cross-sectional area of the white matter axons connecting those regions: a truly quantitative and biologically-relevant measure of ‘structural connectivity’. We have also performed a limited validation of this method by comparing whole-brain streamlines reconstructions with and without SIFT to quantitative and qualitative measures of white matter estimated from post mortem brain dissection (Smith et al., 2015). Although this method is effective at achieving its goal, a significant limitation to its uptake in the field is the computational expense required to provide such reconstructions:

E

72 73

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

The in vivo data were prepared as follows:

134 135 136 137 138

140

• Diffusion image data: correction for B0 field inhomogeneities/eddy currents/inter-volume motion using topup and eddy tools in FSL5 (Andersson et al., 2003; Smith et al., 2004); B1 field inhomogeneity correction using ITK-N4 (Tustison et al., 2010); estimation of singlefibre response function using an iterative optimisation method akin to (Tax et al., 2014); FODs estimated using constrained spherical deconvolution (CSD) (Tournier et al., 2007) as provided in the MRtrix software package (Tournier et al., 2012), with default parameters and a maximum spherical harmonic order of 8. • Anatomical image data: rigid body registration to diffusion image data using SPM8 (Frackowiak et al., 2003); tissue segmentation using FSL tools (Zhang et al., 2001; Smith, 2002; Patenaude et al., 2011) as described in Smith et al. (2012); default Freesurfer segmentation and parcellation pipeline (Dale et al., 1999; Desikan et al., 2006). • Whole-brain streamlines tractography: 10 million streamlines generated using the Second-order Integration over Fibre Orientation Distributions (iFOD2) algorithm (Tournier et al., 2010), incorporating the anatomically-constrained tractography (ACT) framework (Smith et al., 2012). Streamline seeding was performed using one of three mechanisms:

141 142

○ randomly throughout the white matter (as defined from segmentation of the T1 image), denoted ‘WM’; ○ randomly from the grey matter–white matter interface (Smith et al., 2012), denoted ‘GMWMI’; ○ a novel seeding mechanism, where the SIFT model is imposed during whole-brain streamlines tractography, and the probability of seeding from a particular location in the white matter is constantly updated based on the relative difference between the fibre density as estimated from the diffusion model and the current reconstructed

162 163

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

143 144 Q3 145 146 147 148 149 150 151 152 Q4 153 154 155 156 157 158 159 160 161

164 165 166 167 168 169 170

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

streamlines density; this is referred to as ‘dynamic’ seeding, and is described in detail in Appendix 1.

Each streamline s that traverses lobe l contributes streamline density according to the product of its length through the voxel |sl|, and weight eFs. We refer to Fs as the weighting coefficient for streamline s, and eFs as the weighting factor. The value of the proportionality coefficient μ, however, remains fixed throughout the SIFT2 algorithm:

173

179 180 181 182 183 184 185

• Targeted tracking: Following anatomical segmentation from FreeSurfer, the parcellation node corresponding to the left precentral gyrus was extracted, dilated by 2 voxels, and used as a GMWMI seeding mask, generating 10 million streamlines. This was then concatenated with the whole-brain dataset generated using dynamic seeding, to produce a tractogram with 20 million streamlines produced from a combination of whole-brain and targeted seeding.

μ¼

The original SIFT model/method

188

In the original SIFT method (Smith et al., 2013), we first perform FOD segmentation, and assign streamlines to the FOD lobes they traverse, based on both the voxels they pass through and their tangent through each voxel. All streamlines contribute track density to the FOD lobes to which they are assigned, based on their length through the relevant voxel only:

191 192 193

X

TDl ¼

f ¼

l¼1

ð2Þ

R

200

215

SIFT2

216

In SIFT2, we extend the definition of streamlines density in Eq. (1) to permit variable contribution weights from individual streamlines as follows:

207 208 209 210 211 212

217 218

U

205 206

N C O

213 214

Here index l loops over the total number of FOD lobes L in the image; FDl and TDl are respectively the fibre density (taken as the FOD lobe integral), and the track density (sum of streamline segment lengths within the relevant voxel) attributed to that lobe; and PMl is the value of the processing mask in the voxel where lobe l is located (this decreases the contribution of partial-volume-contaminated voxels toward the model fit; see Smith et al., 2013). The expression is simplified by considering all FOD lobes independently, regardless of the number of lobes present in any particular voxel; note however that all FOD lobes within a voxel will possess the same value of the processing mask. In the original SIFT method, the streamlines density in each FOD lobe (Eq. (1)) is altered dynamically to reflect the decrease in streamlines density as streamlines are removed from the reconstruction; this concomitantly causes an increase in the proportionality coefficient μ according to Eq. (2) as filtering proceeds.

203 204

X   TDl ¼ jsl j  e F s 220

s:jsl j N 0

ð4Þ

i

L h X

N h i i X PM l  ðμ  TDl −FDl Þ2 þ A  λreg  f reg ðsÞ :

227 228 229 230 231 232 233 234 235

ð5Þ

s¼1

ð3Þ

This is similar to the cost function in SIFT; the differences are the weighted calculations of TD l as in Eq. (3) (not explicitly shown in Eq. (5)), and the additional regularisation term freg . A is a scaling parameter designed such that for a given value of user-controllable parameter λreg , the effect of the regularisation is comparable for different imaging and reconstruction parameters:

E

T

C

:

R

l¼1 L X ½PMl  TDl 

E

L X ½PMl  FDl 

201 202

224 225

237

(TDl is the track density in FOD lobe l; |sl| is the length of streamline s attributed to FOD lobe l). We then define a proportionality coefficient ‘μ’ that determines an appropriate global scaling between streamlines density and FOD lobe integrals:

μ¼

PM l 

l¼1

ð1Þ

jsl j

195

197 198

TD0l

i.e. unlike in Eq. (2), the streamlines density in each FOD lobe in this particular expression does not change during optimization; the ‘0’ superscript in TDl0 indicates the ‘initial’ streamlines density. The goal of SIFT2 is to determine a vector of weighting coefficients F, such that when the contribution of each streamline to the model is weighted according to its value in this vector, the streamlines densities match the FOD lobe integrals throughout the image. Note here that, in contrast to the original SIFT method, in SIFT2 all streamlines have a non-zero contribution to the model. The cost function to be minimised is expressed as follows:

s:jsl j N 0

196

l¼1

L h X

D

189 190

222 223

l¼1

186 187

221

L X ½PMl  FDl 

F

177 178

O

175 176

For each seeding mechanism, an additional reconstruction of 100 million streamlines was also produced, in order to compare SIFT2 operating on 10 million streamlines to the case where SIFT is used to reduce 100 million streamlines to 10 million (e.g. Smith et al., 2015).

R O

174

P

171 172

3

A¼

L h i 1X PMl  FDl 2 N l¼1

238 239 240 241 242 243

ð6Þ 245

(N is the total number of streamlines). For the results demonstrated in this article, we implemented two different expressions for the 246 regularisation function freg(s), defined as follows: 247 f tik ðsÞ ¼ F s 2 f aTV ðsÞ ¼

ð7Þ 249

  PM  js j X l l    Γ 2 F s ; F lmean PM l0  sl0  l:jsl j≠0 l0 :jsl0 j≠0 X

8 2 l < e F s −e F mean   > 2 l Γ F s ; F mean ¼  2 > : F s− Fl mean F lmean ¼

1

X

TD0l s:jsl j N 0

9 = ; F s N F lmean > > ; F s ≤ F lmean ;

250

ð8Þ 252 253

ð9Þ 255 256

½jsl j  F s 

ð10Þ

ftik(s) in Eq. (7) is a conventional Tikhonov regularisation, which constrains the weighting coefficients Fs to remain close to zero, such that the weighting factors remain close to unity. faTV(s) in Eq. (8) is a more advanced regularisation akin to total variation (TV) (Rudin et al., 1992), which differs from the simpler Tikhonov regularisation in two principal ways:

258

• The regularisation is asymmetrical: those streamlines with a greater weighting coefficient than the mean are heavily penalised based on the extent to which their weighting factor exceeds the mean, compared to those with lesser weighting where the regularisation is

263 264

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

259 260 Q5 261 262

265 266

4

267 268 269 270 271 272 273

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

applied to the weighting coefficient (shown in Eq. (9)) (the ‘a’ prefix in faTV(s) is indicative of this asymmetry). • It constrains streamlines traversing the same FOD lobes to have similar weights to one another (as encapsulated in Flmean, the mean streamline weighting coefficient of those streamlines traversing lobe l; see Eq. (10)), rather than constraining all weighting coefficients independently of one another.

274

283 284 285 286 287 288 289 290

x¼

L h X

PM L  ðμ  TDl − FDl Þ2

i

y¼

N X

f reg ðsÞ:

ð11Þ

s¼1

l¼1

304 305 306 307 308 309

C

E

R

R

302 303

O

300 301

C

298 299

F

N

296 297

Fig. 3 presents the results of the L-curve analysis on the synthetic phantom data. For Tikhonov regularisation (solid line), the curve approximately follows the expected L-curve trajectory: a minimal data residual with very light regularisation (albeit with a large solution norm), an L-shaped curve, then an increasing data residual as the strength of regularisation is increased, with little reduction of the solution norm. Once the regularisation is too strong, the data residual is comparable to that of the original, unmodified reconstruction (vertical dotted line). This occurs due to the nature of the synthetic phantom: nonzero streamline weighting coefficients are necessary to correct the errors in reconstructed streamlines density; but once the regularisation strength is too powerful, the SIFT2 algorithm is hindered from correcting this error. Conversely, the aTV regularisation (dashed line) is still capable of correcting the gross reconstruction density errors even with very strong regularisation, with only a corresponding mild detriment to the data residuals: individual streamlines are constrained to have a similar weight

U

295

to other streamlines within the same bundle, so the gross reconstruction density within each independent bundle can still be corrected. Although this mode of regularisation does not obey the typical L-shaped curve, the lack of horizontal plateau in the ‘L’ shape (which corresponds to an increase in data residuals with no decrease in solution complexity) is actually desirable behaviour. Indeed, the results indicate that in this case, there is no maximal value of λ above which the solution breaks down; rather, the strength of regularisation can be tuned to provide a compromise between variations in streamline weights and residual data errors. Fig. 4 shows some examples of reconstructions provided by SIFT2 generated in order to produce the L-curves in Fig. 3. Each streamline is coloured according to its weighting coefficient; the un-affected reconstruction is shown at top-left with an equal weighting of 1.0 for all streamlines. The SIFT solution is shown bottom-left with all retained streamlines coloured according to their effective weight given the reduction in total streamlines density, and removed streamlines are hidden.

T

Results 293 294

Fig. 3. L-curve analysis for the synthetic phantom data. Solid line: SIFT2 with Tikhonov regularisation; labels indicate value of λ. Dashed line: SIFT2 with aTV regularisation; labels indicate value of λ. Dotted line: data cost of initial streamlines reconstruction. Axes are scaled equally.

E

292

O

281 282

R O

279 280

P

277 278

The profile of the asymmetrical regularisation function Γ2(Fs, Flmean) in Eq. (7) is presented in Fig. 2, as a function of both weighting coefficient Fs, and weighting factor eFs. Global optimisation of the cost function f in Eq. (5) is prohibitively expensive; not only are there large numbers of both FOD lobes (L) and streamlines (N), but the latter greatly exceeds the former, resulting in a highly under-determined system. Our proposed optimisation method is described in detail in Appendix 2. To test the effect of regularisation on the solution, we performed an L-curve analysis by running the optimisation algorithm across a range of regularisation coefficients λ (range 10−12–103), for each of the two regularisation functions, on both synthetic phantom and in vivo data. Each point of the L-curve represents a single execution of the SIFT2 algorithm with a particular regularisation: the x- and y-axes positions are determined by the data residuals and regularisation residuals respectively, as follows:

D

275 276

Fig. 2. The function Γ2(Fs, Flmean) used for asymmetrical streamline weight regularisation, plotted as a function of (left) streamline weighting coefficient Fs, and (right) streamline weighting factor eFs. For both plots a value of Flmean = 0.0 was used.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327

5

P

R O

O

F

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363

E

T

C

• For Tikhonov regularisation, increasing the regularisation coefficient above some critical value results in a sharp increase in the data residuals, as SIFT2 is hindered from performing the necessary modifications to correct the gross errors in reconstruction density. • For aTV regularisation, the curve follows a somewhat diagonal trajectory; this suggests that the regularisation coefficient λ may be tuned for specific applications, providing a balance between the fit to the image data and the complexity of the solution. Although it may appear

364 365 366 367 368 369 370

E

335 336

R

333 334

R

331 332

that this mode of regularisation is inferior to Tikhonov regularisation (due to the elevated regularisation strength for equivalent data residual), we note that by its very design the aTV regularisation equation penalises large streamline weights much more severely than does Tikhonov regularisation, so a direct comparison between these curves is limited.

N C O

329 330

With too little regularisation (λ = 10−6), streamlines traversing the same pathways may take drastically different weighting coefficients; this is undesirable as a large fraction of the streamlines density in a voxel may be encapsulated within a small number of streamlines, reducing the ‘effective’ density of the reconstruction. With increased regularisation (λ = 10−3), the prevalence of streamlines with ‘extreme’ weighting values is reduced; streamlines within each bundle tend to have more similar weighting coefficient values, but some variability in weights is still used to correct reconstruction density errors due to random streamline seeding and propagation. The properties of those reconstructions with heavy regularisation (λ = 1) depend on the mode of regularisation employed: with Tikhonov regularisation, streamlines are prevented from reaching the necessary weighting values to correct the reconstruction density errors, resulting in a large data residual and small solution norm; whereas with aTV regularisation, the streamlines can attain the necessary values, but are restricted from taking values that differ significantly from the other streamlines within the same bundle. In the solution provided by the original SIFT method (bottom left of Fig. 4), the streamlines density is approximately equal in the two bundles, though the total number of streamlines in the image is clearly reduced compared to the SIFT2 solutions. Because SIFT has reduced freedom in terms of the mechanisms by which the errors in reconstruction density are corrected, it cannot achieve the same data fit quality as can SIFT2. Fig. 5 presents the results of the L-curve analysis as performed on in vivo data. The curves are similar to those observed in Fig. 3 for the synthetic phantom data:

U

328

D

Fig. 4. Streamlines reconstructions for the synthetic phantom data. Each streamline is coloured according to its weighting coefficient (see colour bar at bottom of figure); exception is bottom left (SIFT reconstruction) where streamlines removed by the filtering algorithm are shown in black. Parameters for each plot include the regularisation type (shown as a subscript for the parameter λ) and coefficient; residual data error; and regularisation residual.

Fig. 5. L-curve analysis for in vivo data. Line style specifies algorithm and/or regularisation type as in Fig. 3 (solid line: SIFT2 with Tikhonov regularisation; dashed line: SIFT2 with aTV regularisation; dotted line: data cost of initial reconstruction); colour specifies seeding type (red: WM seeding; green: GMWMI seeding; blue: dynamic seeding). For SIFT2 curves, dots appear where the regularisation coefficients are a power of 10, with larger dots at values: 1e−12, 1e−9, 1e−6, 1e−3, 1e0, 1e3. Note that the x-axis is stretched by a factor of 5 compared to the y-axis to highlight features of the curves.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

394

F

O

R O

392 393

D

390 391

E

388 389

T

386 387

C

384 385

E

382 383

R

380 381

R

378 379

O

377

C

375 376

Fig. 7 further highlights the differences between these reconstructions, by examining the underlying reconstructions more closely. Instead of showing the streamlines density at the native DWI acquisition resolution, this figure presents track density images (TDIs; Calamante et al., 2010) generated at super-resolution (500 μm isotropic) to highlight sub-voxel details (Calamante et al., 2011). Note that in the case of SIFT2, it is not a streamline count that determines the intensity of each voxel, but the sum of streamline weighting factors; that is, the contribution of each streamline to the TDI is weighted by the crosssectional area multiplier determined for that streamline by SIFT2. Unlike Fig. 6, here each image is windowed according to the maximum intensity in the slice. Firstly, the decrease in total streamlines density following SIFT (bottom left) is clearly evident in the appearance of the TDI at high resolution; although the fibre densities are accurate at the acquired voxel resolution, the number of streamlines remaining is inadequate to fully sample the volume of each voxel, with many grid elements in the white matter containing zero or very few streamlines (see colour bar scale). Secondly, with low levels of regularisation, the TDIs from SIFT2 appear quite dark: this is because a small number of streamlines within the slice tend to be assigned extremely large weighting factors (evident as bright spots or streaks on the images), which subsequently affects the intensity windowing of the image. This typically occurs due to noisy

N

373 374

Furthermore, this plot demonstrates that GMWMI (green lines) seeding always outperforms the more naïve WM-based seeding (red lines); but the proposed dynamic seeding (blue lines) always has equal or superior performance to GMWMI seeding. Fig. 6 presents a small selection of the in vivo reconstructions generated in order to produce the L-curves in Fig. 5. For brevity, dynamic seeding was used for all reconstructions shown. The background image is the sum of TDl for all FOD lobes within each voxel, while each lobe is coloured according to its direction and scaled according to the underlying streamlines density. The ‘target reconstruction’ (segmented FOD lobes within each voxel; lengths are FOD lobe integrals) is shown at the top of Fig. 6 for reference. It is immediately clear that the original tractogram (middle left) does not possess a reconstruction density comparable to that of the target fibre densities as estimated from the FODs (top image). Both the original SIFT algorithm (bottom left) and SIFT2 with low to moderate regularisation (middle two columns) provide a reconstruction density similar in appearance to these fibre densities. However, if the regularisation is too strong (right-most column), the algorithm is restricted from performing the necessary perturbations in streamline weights necessary to correct the reconstruction density errors; the images with heavy regularisation therefore share some features with the original reconstruction (though the effect is less pronounced with aTV regularisation than with Tikhonov regularisation).

U

371 372

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

P

6

Fig. 6. Example reconstructions for in vivo data. Background image = total streamlines density in voxel; foreground = fibre populations within each voxel, coloured according to orientation, length scaled by underlying streamlines density. Top image (‘Target fibre density image’): sum of FOD lobe integrals in each voxel, scaled by the WM fraction (an ideal streamlines reconstruction should therefore match this image exactly); middle left: original streamlines reconstruction; bottom left: filtered reconstruction provided by SIFT; all other images: SIFT2 with varying regularisation mode and coefficient. Image intensities are scaled identically between images. Parameters for each plot include: either the method used, or regularisation type (shown as a subscript for the parameter λ) and coefficient for SIFT2; residual data error; and regularisation residual.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418

7

P

R O

O

F

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454

T

C

E

426 427

R

424 425

R

422 423

FOD lobes that are very poorly tracked, but do have at least one streamline assigned to them; where this occurs, it necessitates assigning a very large weight to those streamlines in order to correct the error in reconstruction density, with large weights for those streamline having concomitant effects along the rest of their lengths. The effect is still prominent even at moderate levels of Tikhonov regularisation; this is in part due to the exponential basis used (a streamline weighting coefficient does not need to be particularly large in order to acquire a very large weighting factor). The asymmetrical nature of the proposed total variation-like regularisation is designed to address this issue, by penalising streamlines with weights larger than unity more severely than those with weights smaller than unity. With moderate aTV regularisation (bottom row, second from right), the contrast is relatively flat (in accordance with the fibre density estimates from the diffusion model; see top of Fig. 6), but the overall reconstruction density is much larger than that with SIFT (see associated colour bar scale); note that this is not simply a multiplicative factor applied to the streamlines density, but reflects the superior sub-voxel volume coverage achieved by SIFT2 by retaining all streamlines in the reconstruction. Close examination of this image also highlights some white matter pathways that are poorly or inconsistently reconstructed by the tracking algorithm (particularly in the superficial white matter), by appropriately boosting the density contributed by the few streamlines that were successful in reconstructing such pathways. Similarly to the results shown in Fig. 6, when the regularisation is too strong (rightmost column), the spatial distribution of streamlines density is not as homogeneous as the target fibre densities, and shares some spatial features with the original reconstruction (top left), as the algorithm is hindered from applying the necessary modifications to the streamline weights. These images also have a smoother appearance, as the fibre density is distributed more evenly between all streamlines rather than selectively (which is comparable to keeping some streamlines and removing others), and many non-dominant pathways are not appropriately boosted (which can otherwise increase the visual complexity of the image).

N C O

420 421

U

419

E

D

Fig. 7. Super-resolution track density images of the reconstructions shown in Fig. 6. TDIs calculated at 500 μm resolution, aligned with diffusion image voxel grid. For SIFT2 reconstructions, streamlines contribute to the TDI based on their assigned weighting factor eFs. Each TDI is independently windowed according to the maximal intensity in the displayed slice; all greyscale maps are linear.

These observations are further supported by the streamline weighting coefficient histograms presented in Fig. 8. With inadequate regularisation (second column from left), individual streamlines may take either a very small or a very large value. As the regularisation coefficient is increased, the histograms narrow; though at λ = 10−3 it is evident that aTV regularisation has a shorter ‘tail’ at large values than does Tikhonov regularisation, due to its asymmetrical nature; this is the mechanism that prevents individual streamlines from taking excessively large weights, and therefore dominating the density maps as shown in Fig. 7. Also, for low levels of regularisation the histograms contain two distinct peaks (this is akin to those streamlines removed and retained in SIFT), whereas with high regularisation a single peak is observed; the two peaks merge into one at approximately λ = 10− 1 (data not shown). An equivalent histogram for the original SIFT method is also shown (bottom left of Fig. 8). Streamlines removed from the reconstruction are shown with a weighting coefficient of −∞ (weighting factor of zero), and retained streamlines are shown with a weight equivalent to the increase in μ during filtering (in this case, around 41, or e3.7). This further demonstrates why the filtered tractogram provided by SIFT appears less ‘dense’ than the weighted tractogram provided by SIFT2. Fig. 9 further highlights the key benefit that SIFT2 provides over a reconstruction filtered using SIFT. The left-most image shows an axial image of the posterior half of the left hemisphere, at the level of the forceps major of the corpus callosum; the target fibre densities (as provided by spherical deconvolution) are overlaid. In the initial streamlines reconstruction (second column), all streamlines are interpreted as possessing unity weight. The fibre densities estimated from this reconstruction (top row) are clearly different to those estimated from the diffusion model. After applying the SIFT algorithm (third column), the reconstructed fibre densities are comparable to those estimated from the diffusion data, but the remaining streamlines reconstruction itself is relatively sparse; individual streamlines are clearly visible; and due to the decrease in the proportionality coefficient during filtering, each remaining streamline effectively possesses a very large weighting factor (around 41, or e3.7). Using the SIFT2 algorithm (rightmost column), the

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

R O

O

F

8

D

E

T

C

E

502 503

R

500 501

R

498 499

O

496 497

applied to the combined reconstruction (top right), the dense tracking from this region of interest is lost, as the gross over-representation of the pathways emanating from this region can only be corrected through gross removal of streamlines (filtering was performed to convergence in this case; Smith et al., 2013). With SIFT2 (bottom right), all 20 million streamlines are retained, but those created through targeted seeding tend to be assigned considerably smaller weights than those created from whole-brain seeding; as a consequence, the very dense sampling of the plausible pathways emanating from the structure of interest is retained (evident as the apparent high SNR appearance of the relevant TDI in the corresponding pathways), but the overall streamlines density contrast is still consistent with the fibre density estimates obtained from the diffusion model.

C

494 495

reconstructed fibre densities are again a good match for those obtained through spherical deconvolution; but the entire initial tractogram is retained by the optimization process, such that the reconstruction provided has both quantitative properties and highly dense sampling of the underlying fibre orientation field. Fig. 10 demonstrates how SIFT2 may be used to combine the quantitative properties that can be assigned to a whole-brain reconstruction, with targeted tractography to very densely sample the structural connections emanating from a region or regions of interest, whilst maintaining the beneficial attributes of both approaches. The initial reconstruction (centre) contains a combination of streamlines seeded specifically from the left precentral gyrus (top left), and those generated from whole-brain seeding (bottom left). If the original SIFT method is

N

492 493

U

491

P

Fig. 8. Streamline weighting coefficient histograms of the reconstructions shown in Figs. 6 and 7. x-axis = streamline weighting coefficient; y-axis = streamline count. Note that the maximum histogram bin density (y-axis range) varies between plots, while the x-axis is fixed between plots. In the original unmodified reconstruction (top left), all streamlines are interpreted to have identical unity weighting, corresponding to a zero weighting coefficient. For the SIFT algorithm (bottom left), streamlines removed by the algorithm are interpreted to have a weighting factor of 0, corresponding to a weighting coefficient of −∞; for those streamlines retained by the algorithm, the effective weighting factor was determined using the ratio between the proportionality coefficient before and after SIFT (41 in this case), then converted to an effective weighting coefficient (loge(41) = 3.7).

Fig. 9. Comparison between reconstructions provided by SIFT and SIFT2. Top row: fibre densities as estimated from either the raw image data (leftmost), or the streamlines reconstruction. Bottom row: streamlines reconstructions; each streamline is coloured according to its weighting coefficient (this is zero for all streamlines in the initial reconstruction, and loge(41) = 3.7 for all streamlines retained in the SIFT reconstruction in this particular instance). From leftmost to rightmost column: Fibre densities as estimated by spherical deconvolution; initial streamlines reconstruction; reconstruction filtered using SIFT; reconstruction with streamline weights determined using SIFT2.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

504 505 506 507 508 509 510 511 512 513 514 515 516

9

D

P

R O

O

F

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

C

T

E

Fig. 10. Comparison between performance of SIFT and SIFT2 algorithms on a dataset generated using a combination of whole-brain and targeted streamline seeding. Super-resolution Track Density Images (250 μm resolution); coronal slice aligned approximately with left precentral gyrus; all images windowed independently according to maximal image intensity within slice. Left column: Original streamlines reconstructions; top = targeted seeding (seed = left precentral gyrus), bottom = whole-brain seeding. Middle: Combined tractogram. Right: Processed reconstructions; top = SIFT, bottom = SIFT2.

Discussion

518

526

We have demonstrated that, as was the case for SIFT, the SIFT2 algorithm is capable of modifying a streamlines reconstruction such that the streamlines volume estimates are an accurate representation of the underlying biological fibre densities as estimated from the raw diffusion data. This is a critical step in performing quantitative evaluations of white matter connectivity using diffusion MRI, to ensure that any observations made from these data are reflective of the underlying biological connectivity and not resultant from the inadequacies of the reconstruction modality.

527

SIFT2 vs. SIFT

528

The key benefit of the SIFT2 algorithm over the previously-published SIFT method is that it retains the entire generated tractogram for further processing; that is, it does not rely on the removal of streamlines to achieve a good fit with the diffusion model. In addition to the processing efficiency benefits that this provides (as it is not necessary to generate a number of streamlines greatly in excess of the actual desired number), it means that the reconstruction to be used in subsequent processing can still contain a very large number of streamlines, which should be beneficial for the stability of quantitative analyses. The quantitative properties of a whole-brain reconstruction can even be combined with targeted tracking to provide very densely sampled, yet still quantitative, estimates of structural connectivity between user-defined regions (e.g. Fig. 10). Furthermore, in practice the application of SIFT requires a compromise between execution time, extent of filtering (which directly influences the quality of the resulting fit with the fibre densities), and density of the residual reconstruction; SIFT2 instead requires setting

524 525

529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544

R

R

522 523

N C O

521

U

519 520

E

517

regularization parameters, which have a much more direct and interpretable causal influence on the resultant reconstruction. Despite these benefits, we believe that the original SIFT method may still be applicable in a number of contexts. For instance, the use of SIFT2 requires that any subsequent processing be compatible with the definition of a cross-sectional multiplier for each streamline, which may not be possible in all applications (though in our experience it is a relatively trivial modification to any method). Another example is in visualization: highly dense streamlines reconstructions are difficult to display, even without the confounding factor of variable cross-sectional areas; SIFT provides a means to display streamlines reconstructions directly, with removal of the unwanted density fluctuations present in unfiltered streamlines reconstructions. It may also be beneficial to use a combination of these two methods when preparing streamlines data for quantitative analysis. For instance, consider an erroneous streamline that traverses a loop within the brain white matter: SIFT would almost certainly remove such a streamline, as it would likely be repeatedly contributing to an over-defined pathway, whereas SIFT2 would simply reduce the weight of such a streamline. By using these two methods in conjunction, partial filtering by application of SIFT could remove such erroneous trajectories and greatly reduce the over-definition of those pathways most easily tracked, and SIFT2 could then correct the residual reconstruction density errors, making use of a greater fraction of the initial reconstruction than would occur if using SIFT alone. The efficacy of this combined approach is however beyond the scope of this study. We note that, as is the case with SIFT, the SIFT2 method relies on the premise that the estimates of fibre density derived from the FOD are biologically accurate. This interpretation is based on a fundamental imaging property exploited through spherical deconvolution: that the diffusion signal measured from a fibre population is proportional to the MR-visible volume of that population; spherical deconvolution

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576

602 603

In addition to the SIFT2 algorithm, we have also presented the concept of dynamic seeding, where a tractogram is compared to the image data as streamlines are generated, and this information used to determine appropriate streamline seed points as the reconstruction process continues. This method therefore introduces some degree of global information into the otherwise locally-greedy streamlines tractography process. We have demonstrated that this seeding mechanism outperforms homogeneous seeding either throughout the WM or from the GM–WM interface (Fig. 5); however it is not alone adequate to fully correct the reconstruction biases that arise in whole-brain streamlines tractography. It is therefore ideally suited to be used in conjunction with methods such as SIFT and SIFT2, by reducing the degree to which the initial reconstruction must be perturbed in order to provide an adequate fit to the diffusion data.

604

Regularisation

605

In Eqs. (5)–(10), two independent modes of regularisation were defined: a Tikhonov regularisation that constrains the streamline coefficients to be close to zero (i.e. the weights should be close to unity), and a more advanced regularisation akin to total variation, that constrains the streamline weights to be similar to the weights of other streamlines traversing the same FOD lobes. In practise these two modes of regularisation perform quite differently, particularly when the regularisation is comparatively strong. This is due to the inherent differences between the two as described in the Methods section:

598 599 600 601

606 607 608 609 610 611 612 613 614

O

C

630 631 632 633 634 635 636 637 638 639 640

649

One difficulty that arises due to the operation of the SIFT2 method is the treatment of those FOD lobes that contain very little streamlines density compared to the estimated fibre density. This situation can arise due to noise in the diffusion acquisition, inadequacies of the diffusion model, the ill-posed nature of the tractography process, or the effect of isotropic signal contamination at the GM–WM interface on spherical deconvolution (Roine et al., 2014). In the original SIFT method, any streamline traversing such an FOD lobe would simply be slightly more likely to be retained (as its removal would worsen the model fit in that particular FOD lobe). In SIFT2, what occurs instead is that the optimization process tries to assign a very large weighting to those streamlines in order to match the streamlines density with the fibre density. Concomitantly, streamlines that share a common pathway with such a streamline elsewhere along its length may be assigned a lesser weight, as the former erroneous streamlines ‘claim’ a large fraction of the fibre density in all of the FOD lobes they traverse. This issue is largely addressed by the proposed aTV regularization, which heavily penalizes the assignment of very large weights to individual streamlines. This is still however an imperfect solution: the reconstructed streamline density for an FOD lobe traversed by a solitary

Table 1 Execution times, peak memory usage, and other statistics, for a selection of reconstructions provided by SIFT and SIFT2. All reconstructions generated using an Intel Core i7-4790 processor with 32 GB of RAM.

U

N

629

Poorly-reconstructed FOD lobes

t1:1 t1:2 t1:3

619 620

627 628

642 643

621 622

617 618

625 626

Based on our L-curve analyses, there is not a single ideal value for the regularization parameter λ. Instead, the optimal value may depend on the intended application: a smaller value may provide a superior fit to the image data, at the expense of potentially allowing individual streamline to obtain very large weights; whereas a larger value may provide a more visually pleasing result, at the expense of a slightly inferior fit to the image data.

• By definition, streamlines in different pathways throughout the image must acquire different weights in order to correct the reconstruction density errors in question. When the aTV regularisation is strong, it reduces the variations in streamline weights within these pathways, but still allows gross errors in reconstruction density between different pathways to be corrected. With Tikhonov regularisation, on the other hand, once a certain threshold in regularisation strength is exceeded, the streamlines are no longer able to obtain the necessary

615 616

623 624

641

T

596 597

C

594 595

E

592 593

R

586 587

R

584 585

F

590 591

583

O

Dynamic seeding

581 582

R O

589

579 580

weighting coefficients needed to correct these density errors; this is most evident in Figs. 3 and 4, where increasing the Tikhonov regularisation strength above a certain level prevents the streamlines from obtaining the weights necessary to correct the errors in streamlines reconstruction density, even in a very simple synthetic phantom. • The asymmetrical nature of the aTV regularisation penalises streamlines with very large weights more heavily that streamlines with very low weights. This helps in preventing individual streamlines from taking very large weighting factors, which is detrimental to the ‘effective’ density of the reconstruction (a reconstruction with very few streamlines with large weights, and many streamlines with small weights, is comparable to a reconstruction that simply has very few streamlines). The effect is most pronounced in Fig. 7, where the use of Tikhonov regularisation results in bright ‘streaks’ in the TDIs corresponding to individual streamlines with very large weights. The asymmetrical nature of the proposed aTV regularisation ensures that the fibre density represented by each FOD lobe is ‘shared’ relatively equally between those streamlines traversing it.

E

588

simply transforms that data such that the amplitude along a particular direction is proportional to the measured signal that is oriented in that direction. Confounding factors, such as incomplete correction of bias field or breakdown of the canonical response function assumption, will affect the accuracy of this measure; we note however that both SIFT and SIFT2 are in fact independent of the particular mechanism used to derive the fibre volumes in the target reconstruction, errors in this quantification are likely significantly less than the errors in density of streamlines reconstruction without such methods, and that alternative approaches (e.g. predicting the raw DWI signal from the tractogram; see Comparison With Other Methods section) will invariably suffer from the same limitations.

P

577 578

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

D

10

t1:4 t1:5

Seeding mechanism

Streamlines generated

Post-processing algorithm

Streamlines remaining

Cost function (data only)

Peak memory usage

Tracking time

Post-processing time

Total time

t1:6 t1:7 t1:8 t1:9 t1:10 t1:11 t1:12 t1:13 t1:14

WM GMWMI Dynamic WM GMWMI Dynamic WM GMWMI Dynamic

10 million 10 million 10 million 100 million 100 million 100 million 10 million 10 million 10 million

SIFT SIFT SIFT SIFT SIFT SIFT SIFT2 SIFT2 SIFT2

456,645 575,755 584,469 10,000,000 10,000,000 10,000,000 10,000,000 10,000,000 10,000,000

702.9 667.4 522.5 909.0 780.2 612.2 612.7 505.6 484.4

2.2 GB 1.7 GB 2.3 GB 26.2 GB 17.8 GB 21.5 GB 2.8 GB 1.9 GB 2.2 GB

2:17:44 1:09:31 1:31:49 22:57:50 11:34:19 15:11:42 2:17:58 1:09:27 1:31:40

0:31:38 0:36:30 0:44:50 2:51:02 2:35:54 4:37:26 0:53:08 0:34:55 1:07:49

2:49:22 1:46:01 2:16:39 25:48:52 14:14:13 19:49:08 3:11:06 1:44:22 2:39:29

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

644 645 646 647 648

650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

Computational expense

696 697

Table 1 provides execution statistics for a selection of experiments to summarize the relative performance between SIFT and SIFT2. Here we make two key observations for the real-world application of these methods:

698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713

C

692 693

• If a fixed number of initial streamlines is generated (as was the case for the results presented in this article), then the total necessary execution time and memory is comparable between SIFT and SIFT2; but SIFT2 provides a vastly more dense final reconstruction, and a marginally superior model fit. • If a fixed number of resulting streamlines is desired, then SIFT requires generation of an initial excess number of streamlines (a factor of 10 was used in this particular case), such that the desired count is achieved after streamline removal. This requires a much longer processing time (to generate the necessary excess of streamlines) and a greater amount of memory (as the FOD lobe visitations for all generated streamlines must be stored), and still provides an inferior model fit to SIFT2 (suggesting that the extent of filtering may not be adequate for quantitative analysis).

E

690 691

R

688 689

R

686 687

N C O

684 685

714 716 717 718 719 720 721 722 723 724 725 726 727 728

Independently of these computational advantages, we believe that SIFT2 should experience a greater uptake in the research field due to the difference in philosophy between its operation and that of SIFT. The intrinsic inefficiency of discarding the vast majority of data generated has in our experience been a significant barrier to the adoption of SIFT amongst the user community. SIFT2, on the other hand, makes use of all streamlines data generated, and can produce relatively dense quantitative reconstructions on standard desktop computing hardware.

U

715

Conclusion Providing useful measures of white matter connectivity using diffusion MR imaging is an important step in characterising the healthy structure of the human brain, as well as its perturbation in disease. The novel SIFT2 method presented here provides a quantitative and biologically-meaningful measure of structural connectivity throughout

731 732 733 734

Li et al., 2012

736

Acknowledgments

737

We are grateful to the National Health and Medical Research Council (NHMRC) of Australia, the Australian Research Council, and the Victorian Government's Operational Infrastructure Support Program for their support.

738 739

F

695

682 683

735 Q6

O

694

Subsequent to the SIFT2 method first being proposed (Smith et al., 2014), a number of alternative methods have been published that also derive appropriate weights for individual streamlines in order to match a tractogram to the diffusion image data (Daducci et al., 2014; Lemkaddem et al., 2014; Pestilli et al., 2014). Unlike SIFT2, these operate on the diffusion image data directly rather than estimated FODs; this may prevent the issues related to noisy FOD lobes mentioned above, but requires the definition of an appropriate model for predicting the image information from streamlines data (which we note is comparable to defining a canonical response function for spherical deconvolution), and may be more computationally expensive. They also use generic linear equation solvers, which can be highly efficient but may preclude the use of advanced regularization. A comparison of these techniques is beyond the scope of this study.

Uncited reference

740 741

Appendix A1. Whole-brain streamlines tractography using dynamic 742 seeding 743

R O

681

677

729 730

In this section we describe a novel approach to streamlines seeding, based on the SIFT model (Smith et al., 2013). In brief, reconstructed streamlines are mapped to the traversed FOD lobes immediately after generation, and the relative differences between FOD amplitudes and streamlines densities are used to dynamically alter the probability of seeding from each FOD lobe. Fig. A1 demonstrates the software implementation of this framework. The implementation is modular (i.e. any appropriate streamlines algorithm or target fibre density image may be used) and thread-safe. Once a streamline is generated and written to the output track file, it is then passed to an additional processing block that maps the streamline to the diffusion voxel grid, capturing both the streamline length and tangent through each voxel (Smith et al., 2013). This is then used to update both the reconstructed streamlines density attributed to each FOD lobe, and the proportionality coefficient μ within the SIFT model that performs the global scaling between FOD amplitude and streamlines density (this parameter must be updated dynamically to reflect the continual increase in streamlines density). The implemented dynamic streamline seeding mechanism uses a rejection sampling framework. Initially, all FOD lobes are assigned a uniform initial seeding probability (e.g. 10−3). Each time an FOD lobe is selected at random for potential seeding, its ‘cumulative seeding probability’ pcl is updated based on the seeding probabilities applied before and at the previous update for that lobe:

P

Comparison to similar methods

675 676

the brain, but using a mechanism that is more computationally efficient and logically direct than its predecessor SIFT algorithm. We believe it is therefore a crucial processing step in any neuroscientific analysis where the researcher aims to use streamlines tractography in a quantitative manner. This software is available as part of the MRtrix3 software package (http://github.com/jdtournier/mrtrix3).

D

680

673 674

T

678 679

streamline will still be much larger than a similarly-sized FOD lobe that is not traversed at all, introducing instability into the quantification of connectivity. Additional heuristics can be used (and indeed have been implemented) to assist in excluding such FOD lobes from the optimization process, but were omitted from this manuscript for brevity. The SIFT2 method would also benefit from higher SNR during acquisition, improved models for deriving the FOD (e.g. Jeurissen et al., 2014; Roine et al., 2015), and improvements to tractography algorithms.

E

671 672

11

pcl ¼

  Nprev  pcl þ N−N prev  pl N

744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767

ðA1Þ 769

(N: number of streamlines generated; Nprev: number of streamlines generated when lobe was previously updated; pl: seeding probability previously applied for lobe l). Then, the seeding probability pl for lobe l is updated based on the current ratio Rl between reconstructed streamlines density and estimated fibre density in lobe l as follows: Rl ¼

pl ¼

μ  TDl FDl 8 < :

min 1:0;

ðA2Þ

pcl Nt arget −N  Rl ;  Rl Nt arget −N 0:0;

Rl b 1:0 Rl ≥ 1:0

9 = ;

770 771 772 773 774

776 777

ðA3Þ 779

(μ: SIFT proportionality coefficient; TDl: track density in lobe l; FDl: fibre density of lobe l as estimated from the FOD; Ntarget: target stream- 780 line count). 781

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

 Nt arget ¼ min 2  N;

Nuser



ðA4Þ

792

797

(Nuser: user-requested total number of streamlines). Once an FOD lobe has been selected for streamline seeding, a random position within that voxel is determined; and if ACT is used, the seed position is verified against the segmented anatomical image, with revision of the seed point to the grey matter–white matter interface if necessary (Smith et al., 2012).

798

2. SIFT2 optimisation algorithm

799

The goal of the SIFT2 algorithm is to find a vector of streamline weighting coefficients F that minimises the sum of squared differences between the fibre densities as estimated from the diffusion model, and the corresponding reconstructed streamlines densities. All relevant symbols and expressions are presented in the Methods–SIFT2 section; the primary cost function is reproduced here for clarity:

793 794 795 796

800 801 802 803 804

f ¼

L h X

N h i i X PM l  ðμ  TDl −FDl Þ2 þ A  λreg  f reg ðsÞ :

ðA5Þ

s¼1

l¼1

806

0 : l:jcl j≠0 TDl

811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828

ðA6Þ

(|cl| is the length of the candidate streamline c in FOD lobe l; TD0l is the initial streamlines density in lobe l when every streamline has a weighting coefficient of zero (weighting factor of one)). Here the streamline density in each traversed FOD lobe has been broken into four components: the first term describes the streamline density attributed to that lobe given the current weighting coefficient vector F; the second term subtracts the track density currently provided to that lobe by the candidate streamline; the third term is the density

T

C

O

R

R

E

C

809 810

9

2 # = δθ  PMl  μ  TDl −jcl j  e F c þ jcl j  e F c þΔ F c þ l  ΔF c −FDl þ A  λ  f reg ðcÞ ; δF c

N

808

Generic linear system solvers are infeasible in this context due to the high dimensionality of the problem; the number of discrete FOD lobes L is ~250,000 at a coarse imaging resolution of 2.5 mm, and the number of

arg minΔ F c 8 " < X jc j l

U

807

F

789 790

O

788

R O

786 787

streamlines N can range from 106 to 109. Independent optimisation of each streamline weighting coefficient is not appropriate; large differences between streamlines densities and FOD lobe integrals would be corrected by substantially modifying the weighting factor of individual streamlines, but the cumulative effect of many streamlines correcting the same errors would be significant overshoot of the target streamlines densities. A naïve gradient descent also performs poorly: some regions of the image converge very quickly, such that a very small step size must be used to prevent overshoot and/or instability in these regions; not only is the determination of an appropriate step size ill-posed for such a high-dimensional problem, its very small magnitude makes the convergence rate of most of the reconstruction prohibitively slow. Our proposed approach borrows aspects from gradient descent and the expectation–maximisation algorithm (Dempster et al., 1977). It takes into account the potential for interactions between the weighting coefficient of different streamlines, whilst maintaining linear performance characteristics as a function of both N and L. For an individual streamline (henceforth referred to as the candidate streamline ‘c’), we seek to find a change ΔFc in the weighting coefficient Fc that satisfies the following expression:

P

784 785

This update is derived in such a way that Rl → 1.0 as N → Ntarget for all l. The seeding algorithm therefore attempts to correct for the under- or over-reconstruction of all FOD lobes throughout the image, by manipulating the probability of seeding from each lobe as new information regarding the density of streamlines reconstruction arrives, with the intent of providing a tractogram where the reconstructed streamlines density matches the estimated fibre density perfectly. We also use the following modification to provide more rapid modulation of seeding probabilities at the earlier stages of reconstruction:

D

782 783

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

E

12

Fig. A1. Flow chart for ‘dynamic’ streamline seeding framework. Rectangles: processes, parallelograms: data. Note that the framework forms a closed loop: streamlines are generated based on selected seed points; these streamlines are used to update a difference image that compares reconstructed streamlines densities to FOD lobe integrals; and this image then influences the selection of subsequent streamline seed points.

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

830 831 832 833 834 835 836 837

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866

X δTDl ΔF s δθl ¼  : δF c s:js j≠0;s≠c δF s ΔF c

T

ðA7Þ

l

868

875

C

E

R

873 874

  Δ F s  e F c   Δ F ≈ e F s : c

878 879 880 881 882 883 884 885

ðA8Þ

That is, we expect that the relative changes in streamline coefficients should be proportional to the inverse ratio of the streamline weighting factors; this is because streamlines that already possess a large weight will have a greater influence on the overall streamlines density for a given coefficient change, and are therefore likely to experience smaller perturbations. Note that any relative sign between ΔFs and ΔFc is ignored in this derivation, thus providing an over-damped solution. This reduces Eq. (A7) to the following:

U

877

R

871 872

The first term in this expression is a trivial partial derivative that can be determined from Eq. (3). The second fraction requires knowledge of the relative changes in the weighting coefficients of all streamlines for the current iteration; however this is precisely the information to be determined from Eq. (A6), and is therefore not known. Although a number of advanced heuristics were trialled, our current implementation (which proved to be superior in terms of both stability and speed of convergence) uses the following simple expression:

N C O

869 870

  δθl ¼ e F c  TD0l −jcl j δF c

ðA9Þ

887 888 889 890 891 892 893

The optimisation algorithm therefore proceeds in the following manner: 1. Update the streamline weighting coefficient vector F using line searches (Eq. (A6)). 2. Re-calculate the FOD lobe streamline densities TDl (Eq. (3)) and the mean streamline weighting coefficient for each FOD lobe Fmean l (Eq. (10)).

899

Andersson, J.L., Skare, S., Ashburner, J., 2003. How to correct susceptibility distortions in spinecho echo-planar images: application to diffusion tensor imaging. NeuroImage 20, 870–888. Calamante, F., Tournier, J.-D., Jackson, G.D., Connelly, A., 2010. Track-density imaging (TDI): super-resolution white matter imaging using whole-brain track-density mapping. NeuroImage 53, 1233–1243. Calamante, F., Tournier, J.-D., Heidemann, R.M., Anwander, A., Jackson, G.D., Connelly, A., 2011. Track density imaging (TDI): validation of super-resolution property. NeuroImage 56, 1259–1266. Cercignani, M., 2010. Strategies for Patient-Control Comparison of Diffusion MR data. Oxford University Press. Daducci, A., Dal Palú, A., Lemkaddem, A., Thiran, J.-P., 2014. COMMIT: convex optimization modeling for micro-structure informed tractography. IEEE TMI http://dx.doi.org/10. 1109/TMI.2014.2352414. Dale, A.M., Fischl, B., Sereno, M.I., 1999. Cortical surface-based analysis: I. Segmentation and surface reconstruction. NeuroImage 9, 179–194. Dempster, A., Laird, N., Rubin, D., 1977. Maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. 39, 1–38. Desikan, R.S., Segonne, F., Fischl, B., Quinn, B.T., Dickerson, B.C., Blacker, D., Buckner, R.L., Dale, A.M., Maguire, R.P., Hyman, B.T., Albert, M.S., Killiany, R.J., 2006. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage 31, 968–980. Frackowiak, R., Friston, K., Frith, C., Dolan, R., Price, C., Zeki, S., Ashburner, J., Penny, W., 2003. Human Brain Function. Academic Press. Jeurissen, B., Tournier, J.-D., Dhollander, T., Connelly, A., Sijbers, J., 2014. Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion MRI data. NeuroImage 103, 411–426. Jones, D.K., Knosche, T.R., Turner, R., 2013. White matter integrity, fiber count, and other fallacies: the do’s and don’ts of diffusion MRI. NeuroImage 73, 239–254. Lemkaddem, A., Skiöldebrand, D., Dal Palú, A., Thiran, J.-P., Daducci, A., 2014. Global tractography with embedded anatomical priors for quantitative connectivity analysis. Front. Neurol. 5, 232. Li, L., Rilling, J.K., Preuss, T.M., Glasser, M.F., Hu, X., 2012. The effects of connection reconstruction method on the interregional connectivity of brain networks via diffusion tractography. Hum. Brain Mapp. 33, 1894–1913. Mugler, J.P., Brookeman, J.R., 1990. Three-dimensional magnetization-prepared rapid gradient-echo imaging (3D MPRAGE). Magn. Reson. Med. 15, 152–157. Patenaude, B., Smith, S.M., Kennedy, D.N., Jenkinson, M., 2011. A Bayesian model of shape and appearance for subcortical brain segmentation. NeuroImage 56, 907–922. Pestilli, F., Yeatman, J.D., Rokem, A., Kay, K.N., Wandell, B.A., 2014. Evaluation and statistical inference for human connectomes. Nat. Methods 11, 1058–1063. Reese, T., Heid, O., Weisskoff, R., Wedeen, V., 2003. Reduction of eddy-current-induced distortion in diffusion MRI using a twice-refocused spin echo. Magn. Reson. Med. 49, 177–182. Roine, T., Jeurissen, B., Perrone, D., Aelterman, J., Leemans, A., Philips, W., Sijbers, J., 2014. Isotropic non-white matter partial volume effects in constrained spherical deconvolution. Front. Neuroinformatics 8. Roine, T., Jeurissen, B., Perrone, D., Aelterman, J., Philips, W., Leemans, A., Sijbers, J., 2015. nformed constrained spherical deconvolution (iCSD). Med. Image Anal. http://dx.doi. org/10.1016/j.media.2015.01.001. Smith, S.M., 2002. Fast robust automated brain extraction. Hum. Brain Mapp. 17, 143–155. Smith, S.M., Jenkinson, M., Woolrich, M.W., Beckmann, C.F., Behrens, T.E., Johansen-Berg, H., Bannister, P.R., De Luca, M., Drobnjak, I., Flitney, D.E., Niazy, R.K., Saunders, J., Vickers, J., Zhang, Y., De Stefano, N., Brady, J.M., Matthews, P.M., 2004. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage 23, S208–S219. Smith, R.E., Tournier, J.-D., Calamante, F., Connelly, A., 2012. Anatomically-constrained tractography: improved diffusion MRI streamlines tractography through effective use of anatomical information. NeuroImage 62, 1924–1938. Smith, R.E., Tournier, J.-D., Calamante, F., Connelly, A., 2013. SIFT: spherical-deconvolution informed filtering of tractograms. NeuroImage 67, 298–312. Smith, R.E., Tournier, J.-D., Calamante, F., Connelly, A., 2014. SIFT2: enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography. Proc. ISMRM 22, 0278. Smith, R.E., Tournier, J.-D., Calamante, F., Connelly, A., 2015. The effects of SIFT on the reproducibility and biological accuracy of the structural connectome. NeuroImage 104, 253–265. Sporns, O., Tononi, G., Kötter, R., 2005. The human connectome: a structural description of the human brain. PLoS Comput. Biol. 1 (4), e42. Tax, C.M., Jeurissen, B., Vos, S.B., Viergever, M.A., Leemans, A., 2014. Recursive calibration of the fiber response function for spherical deconvolution of diffusion MRI data. NeuroImage 86, 67–80. Tournier, J.-D., Calamante, F., Connelly, A., 2007. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution. NeuroImage 35, 1459–1472.

900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974

F

847 848

References

O

845 846

R O

844

894 895

P

842 843

3. Re-calculate the value of the cost function f (Eq. (A5)) based on the updates to F and Fmean . l 4. Iterate until the decrease in the cost function between iterations falls below some threshold, or some maximum number of iterations is reached.

D

840 841

provided by the candidate streamline to that lobe under perturbation of its weighting coefficient; and the fourth is a ‘correlation term’ that will be described later. The contribution from each lobe traversed by the candidate streamline to this expression is weighted by the fraction of the streamline density in that lobe attributed to the candidate streamline; this ensures appropriate scaling between data and regularisation terms, and gives appropriate emphasis to different lobes along the length of the candidate streamline (depending on both the total number of streamlines traversing each lobe, and the length of the streamline within each voxel). The final term in the expression describes the cost of regularisation under perturbation of the candidate streamline weighting coefficient Fc (see Eqs. (7)–(8)). Given an appropriate expression for δθl/δFc, this equation can be minimised independently for each individual streamline using a line search algorithm (we found the iterative Newton method to give the best combination of speed and accuracy). ΔFc may be limited to a finite range (e.g. [− 1,+ 1]) to mitigate unwanted large perturbations between iterations; streamlines may also be rejected entirely if their weighting coefficient falls below some threshold, or used to exclude noisy FOD lobes from the optimisation process if their coefficient exceeds some threshold. The challenge then in the solution of Eq. (A6) is the definition of the correlation term δθl/δFc. This term indicates a change in the streamline density TDl that is related to a change in the candidate streamline weighting coefficient Fc, but is not the actual change in streamline density in FOD lobe l from the candidate streamline c. This term attempts to encapsulate the projected change in the streamline density of FOD lobe l due to changes in the weighting coefficients of all other streamlines traversing that lobe. This term is defined as follows:

E

838 839

13

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

896 897 898

14

975 976 977 978 979 980 981

R.E. Smith et al. / NeuroImage xxx (2015) xxx–xxx

Tournier, J.-D., Calamante, F., Connelly, A., 2010. Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions. Proc ISMRM 18, 1670. Tournier, J.-D., Mori, S., Leemans, A., 2011. Diffusion tensor imaging and beyond. Magn. Reson. Med. 65, 1532–1556. Tournier, J.-D., Calamante, F., Connelly, A., 2012. Mrtrix: diffusion tractography in crossing fiber regions. Int. J. Imaging Syst. Technol. 22, 53–66.

Tustison, N., Avants, B., Cook, P., Zheng, Y., Egan, A., Yushkevich, P., Gee, J., 2010. N4ITK: improved N3 bias correction. IEEE TMI 29, 1310–1320. Zhang, Y., Brady, M., Smith, S., 2001. Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE TMI 20, 45–57.

U

N

C

O

R

R

E

C

T

E

D

P

R O

O

F

987

Please cite this article as: Smith, R.E., et al., SIFT2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography, NeuroImage (2015), http://dx.doi.org/10.1016/j.neuroimage.2015.06.092

982 983 984 985 986