Differential Tractography as a Track-Based Biomarker for Neuronal Injury

Diffusion MRI tractography has been used to map the axonal structure of human brain, but its ability to detect neuronal injury is yet to be explored. Here we report differential tractography, a new type of tractography that utilizes repeat MRI scans and a novel tracking strategy to map the exact segment of fiber pathways with a neuronal injury. We examined differential tractography on multiple sclerosis, Huntington disease, amyotrophic lateral sclerosis, and epileptic patients. The results showed that the affected pathways shown by differential tractography matched well with the unique clinical symptoms of the patients, and the false discovery rate of the findings could be estimated using a sham setting to provide a reliability measurement. This novel approach enables a quantitative and objective method to monitor neuronal injury in individuals, allowing for diagnostic and prognostic evaluation of brain diseases.

Greater specificity and sensitivity could be achieved by aggregating voxel-wise anisotropy changes 47 diffeomorphic reconstruction (QSDR) (Yeh and Tseng, 2011), a method that generalized GQI to accept 162 spatial transformation in the reconstruction. QSDR allowed us to simultaneously reconstruct and 163 transform SDF from the follow-up scan to the space of the baseline scan using the following formula: 164 Ψ ( , ) = ∑ ( ( ), ) 6 ( )〈 ( ), ( ) 〉 (4) 165 where ( ) transforms spatial coordinate r from the space of the baseline scan to that of follow-up scan. 166 ( ( ), ) is the diffusion-weighted signals at coordinate ( ). ( ) is the Jacobian matrix at the 167 same coordinate that rotates the unit vector . There other variables follow the same notations in Eq. 168 (3). 169 Since the scans were from the same subject, we assumed there was only "rigid body" transformation (i.e. 170 only rotation or translocation) between the scans, and the transformation was a simple matrix-vector 171 multiplication. Please note that this assumption could be violated if there was a large tissue distortion 172 due to edema or tissue removal, and a nonlinear spatial registration should be used in QSDR to handle 173 this problem. The rigid body transformation matrix was obtained by linear registering the b0 images (or 174 the sum of all diffusion-weighted images). We used the negative of the correlation coefficient between 175 the images from baseline and follow-up scans as a cost function to calculate the transformation matrix. 176 The cost function was minimized using a gradient descent method. The rotation matrix of the rigid body 177 transformation was used as the Jacobian for Eq. (4). 178 Please note that the SDFs calculated from Eq. (3) and (4) have "arbitrary units". Therefore, the Z1 179 constant in Eq. (4) had to be scaled to match the same unit of Z0 in Eq.
Equation (6) and (7) provide the anisotropic component of SDF (termed anisotropy hereafter) to 186 minimize the effect of free water diffusion (Yeh et al., 2013). It is noteworthy that this anisotropy 187 measurement has a different physical meaning from the fractional anisotropy (FA) calculated in DTI. FA 188 is a ratio between 0 and 1 calculated from diffusivities and has no unit. The anisotropic SDF has the 189 same physical unit of the SDF, which is the spin density of diffusing water. 190 Tracking differences in the SDF

191
To track differences along the existing fiber pathways, we first determined the local fiber orientations 192 using the peaks on the sum of Ψ ( , ) and Ψ ( , ), and then the anisotropy estimated from the 193 summed SDF was used to filter out noisy fibers and to define the termination of the white matter tracks, 194 as what had been done in conventional tractography (Yeh et al., 2013). The percentage difference in the 195 anisotropy between baseline and follow-up scans was then calculated (Fig. 1c): 196 The percentage changes in the anisotropy, Ψ ( , ) , can have positive values (blue SDFs in Fig. 1c), 198 which indicates an increase in the density of anisotropic diffusion, or negative values (red SDFs in Fig.  199 1c), which indicates a decrease in the density of anisotropic diffusion. 200 An additional tracking-the-differences criterion was added to the fiber tracking algorithm to track the 201 exact segment with a decrease or an increase in the anisotropy greater than a change threshold, 202 Specifically, to track pathways with an increase of anisotropy, the additional criterion checked whether 203 the increase of anisotropy was greater than a predefined value of percentage change (e.g. 20%) and 204 continued tracking as long as the criterion was satisfied: 205 where is the local fiber directions used in the fiber tracking algorithm. Similarly, to track pathways with 207 decreased anisotropy, the criteria continued tracking if the decrease of anisotropy was greater than a 208 predefined value of percentage change (e.g. -20%): 209 This allowed us to track two different sets of pathways, one for increased anisotropy, one for decreased 211 anisotropy. The other existing criteria in conventional tractography (e.g. seeding strategy, propagation 212 interval, angular threshold, length constraint…etc.) remained in effect as what has been used in the 213 generalized deterministic fiber tracking algorithm (Yeh et al., 2013). It is noteworthy that the angular and 214 anisotropy thresholds in the original tracking algorithm were still used in differential tractography to 215 eliminate noisy fiber and to ensure a correct white matter coverage. The criteria (9) or (10) (termed 216 "positive change threshold" and "negative change threshold") served as additional constraints to limit the 217 findings to the exact segment of pathways with a substantial change in the anisotropy value. 218 219 In this study, the differential tractogram was obtained by placing a total of 5,000,000 seeding points in 220 the white matter. The angular threshold was randomly selected between 15 to 90 degrees. The step size 221 was 1 mm, and the anisotropy threshold was automatically determined by DSI Studio. Two iterations of 222 topology-informed pruning (Yeh et al., 2018) were applied to the tractography to remove noisy findings. 223 The above-mentioned setting was regularly used in conventional tractography. We tested differential 224 tractography with different values of change threshold (5%, 10%, 15% …. 50%) and length threshold (5 225 mm, 10 mm, 15 mm, … 50 mm). Track with lengths shorter than the length threshold was discarded, and 226 the results of different length threshold and change threshold were compared to access its effect on the 227 sensitivity and specificity of differential tractography. 228 Estimating false discovery rate using a sham setting 229 Figure 2 illustrates the experimental design that allows for estimating the false discovery rate (FDR) of 230 differential tractography in individual scans. As shown in Fig. 2a, the baseline scan is compared with the 231 follow-up scan (upper row) using differential tractography to reveal tracks with decreased anisotropy (Fig.  232 2b, upper row). If a total of 5,000,000 tracking iterations are conducted, we can view each of them as an 233 independent hypothesis testing since each tracking iteration is done independently. The null hypothesis 234 is "there exists no track with a decrease in the anisotropy". This null hypothesis will be rejected if the 235 track length is longer than a predefined length threshold (e.g. 40 mm). Each rejected hypothesis is thus 236 regarded as a "positive finding", but here the finding can be either a true positive or false positive  Fig. 2b shows that the sham scan generates a total of 177 tracks with a length longer than 40 mm, meaning that the estimated number of false-positive findings is 177. Using 247 this information, we can calculate the FDR, which is 177/13947=0.0127 (Fig. 2c). 248

249
In this retrospective study, we did not have an additional sham scan, and a "substitute sham" approach 250 was used. We assumed that there should be no increased track integrity during the disease course, any 251 findings showing increased anisotropy can be regarded as false-positive findings. Therefore, we can use 252 the number of tracks with increased anisotropy in the follow-up scan as a substitute for the sham scan. 253 Even if tracks with increased integrity do exist due to re-myelination, it will increase our estimated 254 number of false-positive findings. This means that the FDR estimated by this substitute sham approach 255 will be an upper bound of the actual FDR. Therefore, instead of reporting "FDR=0.0127", the FDR 256 calculated using a substitute sham should be reported as "FDR ≤ 0.0127". 257 The processing pipeline for differential tractography and the quality control procedure is implemented in 258 DSI Studio (http://dsi-studio.labsolver.org). Documentation and source code to reproduce the same 259 result is also available on the website. 260

261
Neuronal injury reflected by a decrease of anisotropy 262 Figure 3a shows the intermediate results of differential tractography applied to an MS patient with optic 263 neuritis (patient #1, demographics summarized in Table 1). The baseline scan was acquired right after 264 the onset, whereas the follow-up scan was acquired 6 months after. For each fiber orientation in a voxel, 265 differential tractography compares the anisotropy differences between two MRI scans in a common 266 subject space (Fig. 1a-1c). The fiber orientations with a decrease of anisotropy greater than 30% could 267 are plotted by red sticks in Fig. 3a. In the figure, most of the differences are distributed near the primary 268 visual pathways, whereas some spurious differences are randomly distributed throughout the entire 269 whiter matter regions, most likely due to local signal variations or registration error. 270

271
To eliminate these local spurious differences, the "tracking-the-difference" algorithm is applied to the 272 track and link all local differences together into continuous trajectories, and short fragments are 273 discarded using a length threshold of 40 mm (Fig. 3b). The rationale behind this length threshold is that 274 the local random error does not propagate along fiber pathways, whereas true findings due to neuronal 275 injury will form a continuous decrease of anisotropy along the fiber bundles. A length threshold will 276 effectively differentiate between them and help eliminate false results. 277 The left inset figure in Fig. 3b shows affected tracks in directional colors (red: left-right green: 278 anterior-posterior blue: superior-inferior), whereas the tracks in the right inset figure are color-coded by 279 the percentage decrease of anisotropy suggesting the severity of neuronal injury (yellow: 0% decrease 280 red: 70% decrease). In overall, the differential tractogram in Fig. 3b reveals a heterogeneous decrease 281 of anisotropy between 20% to 50%. All findings are in the bilateral primary visual pathways or their 282 collateral connections. The location of the finding matches well with the patient's medical history of 283 visual field loss in both left and right quadrants. The topology of affected pathways seems to present a 284 ripple effect: not only the primary visual pathway is affected, but also certain collateral connections to the 285 visual cortex has shown a decrease in the anisotropy. Although this patient was fully recovered from the 286 symptoms during the follow-up scan (brief medical history in the Supplementary materials), differential 287 tractography still captures subclinical change near the bilateral optic radiation. 288 We further compare differential tractography with conventional voxel-wise statistics. Figure 3c shows the 289 axial mapping of anisotropy differences for each voxel using the same data. The red regions are voxels 290 with anisotropy decrease greater than 30%, and the numerous fragments can be observed across the 291 entire brain regions. Those fragments could be due to local random error and thus may not be true 292 findings. This illustrates a typical limitation of a voxel-based imaging biomarker, and there is no pathway 293 information to assist diagnostic evaluation. In comparison, the differential tractogram in Fig. 3b provides 294 track-based imaging biomarkers that can easily assist diagnostic evaluation by offering trajectories 295 information. The findings can be associated with an anatomical pathway to infer the functional 296

correlation. 297
Conventional tractography versus differential tractography  Table 1). The conventional 300 tractography was generated using the baseline scan, whereas the differential tractography was 301 configured to map pathways with more than a 30% decrease in anisotropy with a length threshold of 40 302 mm. The trajectories in Fig. 4a and Fig. 4b are both colored-coded with directional colors (red: left-right 303 green: anterior-posterior blue: superior-inferior). The first row shows tractography viewing from the top, 304 whereas the second rows show from the left. The conventional tractography in Fig. 4a visualizes the 305 trajectories of the entire brain fiber pathways, and there is no gross anomaly visible from the 306 tractography that may suggest a major neurological disorder. This is expected because, at the early 307 stage of multiple sclerosis, the patient usually does not present a gross structural change that can be 308 readily identified in conventional tractography. In comparison, the different tractography in Fig. 4b  309 pinpoints the location of affected pathways in the bilateral primary visual pathway near the visual cortex. 310 The location matches well with the patient's disease presentation of optic neuritis, whereas conventional 311 tractography in Fig. 4a shows no clue to this critical information. 312 Reliability assessment

313
We further use patient #1 and a 42-year-old normal subject as the examples to show how the reliability 314 of differential tractography findings can be quantified using a sham setting. Both the follow-up scans of 315 these two subjects were acquired 6 months after the baseline scans. Figure 5a  that around 1% of the tracks shown in the differential tractogram are false results. We can use these two 345 thresholds to leverage sensitivity and specificity. For example, lower thresholds are geared toward 346 higher sensitivity to explore potential neuronal injury, whereas higher thresholds can provide a 347 confirmatory answer to the axonal damage. The optimal setting can be different based on the disease 348 condition, scan interval, and purposes (e.g. exploratory or confirmatory). 349 Differential tractography on patients with neurological diseases 350 We further apply differential tractography to patients with different neurological disorders in Fig. 6 and  351 list the FDR of these findings in Table 2 (Table 1 and Supplementary Materials). The greater severity in 361 patient #1 is reflected by a larger volume of affected pathways diffusion (patient #1: 55681.6 mm 3 , 362 patient #2: 26124 mm 3 ) and a greater decrease of anisotropy shown in the last row. This suggests that 363 differential tractography has a quantitative potential to evaluate disease severity using either the volume 364 of affected tracks or the decrease of their anisotropy. 365 The differential tractogram of HD patients #3 and #4 in Fig. 6 both shows affected pathways around the 366 stratum. This matches well with the common understanding that striatal pathways are usually involved in 367 Huntington's disease. Moreover, the differential tractogram in patient #4 has a greater involvement 368 extending to brainstem and cerebellum, suggesting a worse motor performance. This seems to match 369 the patient's medical history of a higher motor score of 64 (Table 1 and Supplementary Materials). The 370 patient also had more asymmetric dystonia, matching the asymmetry presentation of the differential 371 tractography. 372 Patient #5 in Fig. 6 is an ALS patient. It is noteworthy that this patient had mostly lower motor neuron 373 symptoms (weakness), and thus might not have positive findings in the brain. The differential tractogram 374 of this patient were obtained using a 15% negative change threshold because a threshold of 30% 375 yielded no findings. Differential tractography reveals only a minor decrease in this patient (other cases in 376 Fig. 6 have mostly greater than 30% decrease). This could be explained by the fact that the patient had 377 predominately lower motor neuron symptoms affecting predominately peripheral nerves. Therefore, the 378 findings in the central nervous system could be only subclinical. Nonetheless, as we lowered the change threshold to 15%, differential tractography showed affected pathways in the right lower corticospinal 380 pathway (blue-purple colored), superior cerebellar peduncle, and posterior corpus callosum, as shown in 381 Fig. 6. The right corticospinal pathway involvement seems to match the patient's history of left side 382 involvement, but it is noteworthy that the FDR of these findings were much higher (FDR~0.2 in Table 2 Furthermore, the last row shows that the decrease of anisotropy is mostly greater than 50%, indicating a 393 greater axonal loss due to the surgical removal of the brain tissue. 394 The last column in Fig. 6 shows the differential tractogram of a control. We applied the same settings to 395 examine how differential tractography may capture false results. The result shows a mild decrease of 396 anisotropy as presented by yellow tracks in the last row, a clue that the change may be a false positive 397 result. Furthermore, there are only 74 findings located at the prefrontal cortex, and these findings are 398 relatively insignificant compared to those of the patient population that shows thousands of findings 399 (Table 2). Moreover, the location of the findings is known to be heavily affected by the phase distortion 400 artifact, and the findings could be due to the different level of distortion between the repeat scans. 401 However, it is noteworthy that the calculated FDR will be 0 in this case because there are no findings 402 with increase anisotropy for this control subject. This result suggests that the FDR estimation still have 403 its limitation if the number of findings is sufficiently small. The interpretation of differential tractography 404 results still needs to consider the percentage decrease of anisotropy, the total number of findings, and 405 possible sources of imaging artifact. 406 In general, the findings in Fig. 6 allows us to quickly differentiates the possible locations of neuronal 407 injury and evaluate the severity. The affected pathways in MS, HD, and ALS patients show distinctly 408 different topology, allowing for differential diagnosis or prognosis evaluation. We further investigate whether high b-value signals contribute to the detection power of differential 412 tractography. In this study, we acquired 22 b-values ranging from 0 to 7,000 s/mm 2 , but most diffusion 413 MRI studies would acquire b-values between 1,000 and 3,000 s/mm 2 because higher b-values result in 414 poorer signal-to-noise (SNR) ratio. This can be observed in Fig. 7a showing the diffusion-weighted 415 images at different b-values from patient #1. It is thus of a practical consideration to limit the b-values 416 under a certain maximum value. To examine whether the high b-value acquisition is necessary for 417 differential tractography, we repeated the same differential tractography analysis but used only the 418 signals with b-values between 0 and 3,000 s/mm 2 (a total of 85 sampling directions). The reduced 419 dataset is very similar to the scheme acquired in a previous study (Wang et al., 2011). 420 Fig. 7b shows a qualitative comparison on subject #1 before and after dropping high b-value signals. 421 The location of the findings is mostly consistent; however, dropping high b-values results in substantial 422 loss of the findings (annotated by the red circle). This implies that certain stages of neuronal injury can 423 only be detected by very high b-values. The quantitative comparison is listed in Table 3. As shown in the and #3 have a dramatic FDR increase from less than 0.05 to more than 0.3, making the findings not 427 statistically reliable. It seems that differential tractography using only low b-value signals is more 428 sensitive to physiological variations, whereas differential tractography using both high and low b-value 429 signals are more reliable in detecting are more specific to neuronal injury. This could be explained by the 430 fact that the low b-value images in Fig. 7a have more signal contribution from free water diffusion in 431 ventricles and subarachnoid space, whereas high b-value images only have signals from restricted 432 diffusion in the core white matter. Excluding high b-value signals will shift the focus to free water 433 diffusion and lead to more false-positive results. The overall result supports the necessity of including 434 high b-value acquisition for differential tractography to detect neuronal injury. 435

436
Here we report a novel tractography method to reveal fiber pathways affected by a neuronal injury. We 437 found that differential tractography could serve as a track-based biomarker to provide localization of 438 neuronal injury and allow for quantifying its severity using the decrease of anisotropy and the total 439 volume of affected pathways. The estimated FDR further offered reliability information to interpret the 440 results, and the findings correlated well with clinical presentations of each individual. 441 Comparison with other methods

442
Diffusion MRI fiber tracking can be viewed as a clustering process that is directionally sensitive. By using 443 spatial relations across multiples voxels, the tracking-the-difference strategy used in differential tractography has the potential to differentiate true findings from local errors, since neuronal injury will 445 propagate along axonal fibers while local error stays locally. This tracking-the-difference strategy is 446 conceptually similar to clustering used in voxel-based morphometry (Ashburner and Friston, 2000) or 447 fMRI studies (Woo et al., 2014), which groups voxel-wise statistics into clusters to achieve greater 448 statistical power. The difference is that the clustering used in previous studies did not consider local fiber 449 orientations and would include voxels at all possible neighboring directions, whereas differential 450 tractography only allows findings along the fiber pathways. This improvement is similar to those using 451 fiber geometry to increase specificity (Raffelt et al., 2017;Zhang et al., 2018), and It is likely that this 452 structural restriction may achieve a better specificity than a regular clustering approach. 453 Optimal b-table for differential tractography

454
The diffusion MRI acquisitions played a critical role to boost the sensitivity of differential tractography. 455 Our result shows that acquisitions using only b-values lower than 3,000 s/mm 2 may have a limited 456 detection power for early neuronal injury. We also tested differential tractography on existing DTI data, 457 and the preliminary result (not reported here) also showed a substantially higher rate of false findings. 458 This is not surprising because a typical DTI acquisition only acquires only one b-value of 1,000~2,000 459 s/mm 2 at 30~60 directions, whereas in this study, we acquired 22 b-values from 0 to 7,000 s/mm 2 at 257 460 directions. It is likely that early axonal injury affects mostly restricted diffusion and can only be reliably 461 captured if a wider range of b-value is acquired with enough diffusion sampling directions. 462 Multishell-acquisition could also be used by differential tractography to detect neuronal injury, but its 463 challenges are how to ensure a homogeneous sampling density in the q-space acquisition to make the 464 acquisition "rotation invariant". In addition, most multi-shell acquisitions have an over-sampling problem 465 within each shell and under-sampling problem between shells, which can be observed by plotting the sampling points in the q-space. Consequently, the inter-shell coverage may not be sufficient enough to 467 capture a variety of diffusion patterns, while the intra-shell signals can be easily interpolated due to 468 oversampling, meaning that there is redundancy which can be further reduced to save scanning time. 469 The q-space grid sampling used in this study seems to the method of choice for differential tractography 470 because it has uniform sampling density in the q-space and covers 22 different b-values. This 471 maximized the chance to detect a variety of diffusion pattern that could be altered during the disease 472 process. Q-space imaging used to be criticized for its lengthy scanning time; however, after the 473 introduction of the multi-band sequences, the updated q-space grid sampling scheme could be acquired 474 within 12 minutes for 256 directions and 6 minutes for 128 directions, making it highly feasible for clinical 475 studies. The exact steps to reproduce these two q-space grid acquisitions are documented on the DSI 476 Studio website (http://dsi-studio.labsolver.org/). 477 Anisotropy measurement for differential tractography

478
The anisotropy used in this study is not the commonly used FA provided by DTI. FA is a voxel-wise 479 measurement, and it does not selectively quantify the anisotropy at different fiber population within the 480 same voxel. In comparison, the anisotropy in this study can quantify the anisotropy for each fiber 481 population at different directions This results in greater specificity to individual's connectivity patterns 482 (Yeh et al., 2016) and better performance in handling the partial volume effect (Yeh et al., 2013). 483

484
Another interesting result we observed in this study is that differential tractography seems to capture scans. It cannot access track integrity in a cross-sectional setting, nor does it able to detect any 503 abnormality or axonal injury prior to the baseline scans. Furthermore, differential tractography still has 504 false results if the artifact also propagates coincidentally along a fiber pathway. The parallel imaging or 505 eddy current artifact often gives rise to straight lines near the brain surface but may appear like a 506 spuriously legitimate connection. Misalignment between baseline and follow-up scans can also generate 507 a false result, and it can happen due to registration error or brain tissue shift after surgical intervention. 508 There are still other possible causes of false results. The subjects may have substantial movement in 509 the follow-up scan, but not in the study or sham scan. There may be inconsistency in image acquisition 510 between repeat scans such as changing the head coils or scanning protocol. Both scenarios can 511 produce spurious findings, and thus a series of quality control is always needed to avoid these getting 512 false results. Furthermore, the findings in differential tractography still need to be validated against 513 neuroanatomy. Spurious findings often appear near the brain surface with odd trajectories (straight 514 lines), while true findings tend to follow the trajectories of well-known neuroanatomical pathways. Prior 515 neuroanatomy knowledge may help exclude false results from true findings. 516 There are also limitations in the reliability assessment. The sham setting in this study uses only one scan, 517 and thus the calculated FDR only considered false results due to local random error (e.g. noise which 518 randomly occurs at each imaging voxel) and does not include those due to systematic errors (errors that 519 affects all imaging voxels simultaneously) such as subject movements, coil quality, or signal drift. To 520 detect these systematic errors, we introduced four quality-control routines in the Material and Methods 521 section to discard scans that had a substantial amount of systematic errors. Moreover, we did not 522 acquire an additional sham scan in this study and used the "substitute sham" approach. The substitute 523 sham approach can overestimate the number of false results due to a genuine increase of anisotropy 524 (e.g. re-myelination after recovering from a neuronal injury). The consequence is that the FDR could be 525 higher than the actual value, and the reliability of the findings could be underestimated, leading to a risk 526 of missing meaningful findings. 527 Last, the diffusion MRI protocol in this study could be further optimized. Our result showed that high 528 b-value played an important role in boosting the sensitivity of differential tractography, but the optimal 529 value range still needs further investigation. Moreover, this study did not utilize the recent multi-band 530 acquisition to reduce the scanning time, and future studies will utilize multi-band sequences to make 531 differential tractography more feasible for clinical applications (e.g. 6-minutes grid-128 and 12-minute 532 grid-256 acquisitions mentioned in the previous paragraphs). 533 Clinical applications 534 Differential tractography can be used in differential diagnosis or prognostic evaluation after treatment or 535 intervention. Neurologists can use it to differentiate the cause of a neurological disorder as patients with 536 different neurological disorders will present distinctly different spatial patterns in their affected pathways. 537 The location will provide a clue about the possible causes to resolve difficult clinical cases. This is 538 otherwise not achievable in structural MRI unless a gross lesion or atrophy is visible in the late stage of 539 neuronal injury. Another application of differential tractography is for evaluating an intervention or 540 treatment. Differential tractography can provide an objective quantitation that is directly comparable 541 across subjects and less susceptible to observer differences. It could minimize variance due to evaluator 542 differences and increase effect size in comparison with the conventional evaluation conducted by 543 patients or neurologists. This opens a gate for early treatments to restore subclinical injuries before 544 those injuries accumulate to become a major functional deficit. 545