Mapping the subcortical connectome using in vivo diffusion MRI: feasibility and reliability

Tractography combined with regions of interest (ROIs) has been used to non-invasively study the structural connectivity of the cortex as well as to assess the reliability of these connections. However, the subcortical connectome (subcortex to subcortex) has not been comprehensively examined, in part due to the difficulty of performing tractography in this complex and compact region. In this study, we performed an in vivo investigation using tractography to assess the feasibility and reliability of mapping known connections between structures of the subcortex using the test-retest dataset from the Human Connectome Project (HCP). We further validated our observations using a separate unrelated subjects dataset from the HCP. Quantitative assessment was performed by computing tract densities and spatial overlap of identified connections between subcortical ROIs. Further, known connections between structures of the basal ganglia and thalamus were identified and visually inspected, comparing tractography reconstructed trajectories with descriptions from tract-tracing studies. Our observations demonstrate both the feasibility and reliability of using a data-driven tractography-based approach to map the subcortical connectome in vivo.


Introduction
18 A brain network is comprised of bundles of axons, which form the structural pathways 19 (also referred to as tracts or connections), that allow transfer of information between the different 20 regions (Sotiropoulos and Zalesky, 2019) and facilitate the performance of complex functions 21 (Klingberg et al., 1999;Mesulam, 1998). Axons can be computationally reconstructed (represented 22 as a streamline) using diffusion magnetic resonance imaging (dMRI), a non-invasive technique 23 sensitive to the direction of water motion (Bammer, 2003;Conturo et al., 1999). As axons are 24 bundled together, water molecules will preferentially diffuse parallel to the axonal trajectory, 25 which can then be detected using dMRI to enable an in vivo estimation of tract trajectories. This 26 process, known as tractography, first estimates the diffusion orientations within all imaging voxels 27 before traversing from a starting seed location until termination criteria are met (e.g. quantitative 28 value drops below defined thresholds) (Sotiropoulos and Zalesky, 2019). Additionally, regions of 29 interest (ROI) can be used to define inclusion and exclusion criteria to constrain tract trajectories 30 and facilitate identification of connections between terminal regions. 31 Mapping the human connectome is an important, non-trivial task that contributes to 32 disentangling the network organization of the brain and increased understanding of changes in 33 healthy aging or due to disease (Sporns, 2011;Toga et al., 2012). To date, much of the work 34 studying structural connectivity using dMRI has focused on the cortico-cortical (between regions 35 of the cortex) and cortico-subcortical (between cortex and subcortex) tracts, resulting in the 36 development of a number of structural connectivity atlases. Such connectivity can be described as 37 the cortical connectome. Examples of such atlases include the Johns Hopkins University white 38 matter atlas, which identified a number of cortico-cortical white matter tracts (Hua et al., 2008; 39 Mori et al., 2005;Wakana et al., 2007), and the Oxford thalamic connectivity atlas, which aimed 40 important to motor control (Gallay et al., 2008;Sommer, 2003), as well as cognition and emotion 48 (Hollander et al., 2015). Accordingly, connections between the subcortical structures are integral 49 and have been studied extensively in non-human primate (NHP) studies of the motor network 50 to map its connections. One such example involved the injection of an anatomical tracer at the 55 ventral pallidum, which determined projections to the subthalamic nucleus (STN), as well as the 56 hypothalamus and brainstem (Haber et al., 1993), 57 Studies that attempt to more comprehensively identify the subcortical connections non-58 invasively via tractography, that is to map the subcortical connectome, have been limited. The 59 scarcity of subcortical connectome studies is in part due to the difficulty of tracking the 60 connections in a compact region where the underlying diffusion signal is complicated by multiple 61 diffusion orientations arising from numerous intersecting connections and structures with low 62 anisotropy. One previous study demonstrated the ability to map connections between the basal 63 8 unrelated datasets. The use of a group average response function minimizes biases in FOD maps 149 (Raffelt et al., 2012), improving the comparability of tractography within datasets with observed 150 differences attributed to the underlying diffusion data of an individual. Prior to performing 151 tractography, multi-tissue informed log-domain transformed normalization was performed 152 (Raffelt et al., 2017) on the FOD maps. 153 As the primary diffusion orientation is also reflected in FODs, major white matter 154 connections (e.g. the corticospinal tract (CST)) passing through the subcortical region will hinder 155 the ability to identify subcortical trajectories. To traverse trajectories along non-primary diffusion 156 orientations, the iFOD2 probabilistic algorithm (Tournier et al., 2010) with a step-size of 0.35mm 157 and maximum angle of 45° between successive steps was used. Random seeding was performed 158 throughout the brain until 20 million streamlines, constrained to the subcortical region with the 159 previously created exclusion mask, were selected. The chosen parameters are comparable to what 160 is typically used in iFOD2 algorithms to perform whole-brain tractography with a noted decrease 161 in step-size (from 0.5 × to 0.25 × ) to sample more frequently along a 162 streamline's trajectory. 163 Following tractogram creation, each streamline was assigned a weighting to reflect its 164 contribution to the underlying diffusion signal using the updated spherical-deconvolution informed 165 filtering of tractograms (SIFT2) technique enabling the assessment of tract densities (Smith et al., 166 2015). Using MRtrix3, structural connectivity was established by identifying the nearest 167 subcortical label within a 1.5mm radius at each terminal end of a given streamline. Due to the low 168 anisotropy within gray matter, streamlines whose trajectories intersect other subcortical labels 169 prior to reaching the terminal structures were discarded (see Discussion). Furthermore, the CST, 170 which represents a dominant tract passing in proximity to many subcortical connections, was 9 separately identified in order to visually assess its influence on derived tracts. Identification of the 172 CST was performed using the brainstem and segmentations of both pre-and post-central gyri 173 identified by FreeSurfer as inclusion regions of interest. Generation of the CST was performed 174 until 500 streamlines were identified in each hemisphere. Similar to the connectivity of the 175 subcortical connectome, streamlines had to terminate within a 1.5mm radius of these 176 segmentations to be considered a part of the CST. 177

Assessment of reliability and accuracy
Additionally, intraclass correlation (ICC) was computed for the tract densities between the two 203 datasets as a metric of consistency using a two-way, mixed effects model (McGraw and Wong, 204 1996). In this model, the "raters" (column factor) were the tract densities and the "targets" (row 205 factor) were the test and retest session connectivity. A paired t-test was also conducted between 206 average TDs of the test and retest sessions. To compare average connectivity of the basal ganglia 207 with average connectivity between the basal ganglia and thalamus, an one-way ANOVA was 208 performed. The impact of ROI volume on TD was assessed using ordinary least squares multiple 209 regression, treating the average ROI volume across subjects as an independent variable and TD as 210 the dependent variable. Further, Spearman's correlation was performed between average TD and 211 the absolute percent change between test and retest sessions. 212 213   A TD map was first created for each tract identified in the test and retest sessions by  214 identifying streamlines passing through each voxel. The sum of streamline weights were assigned 215 to corresponding voxels. Following assignment of streamline weights, the fraction of the tract (a 216 value between 0 and 1) passing through a voxel is determined from the TD map and used to 217 compute the overlap between tracts from the weighted Dice similarity coefficient (wDSC; In addition to comparing the tract overlap, a Spearman's correlation was also computed 225 between the average TD across the two datasets. Similar to TD, a one-way ANOVA was performed 226 to compare the wDSC for connectivity of the basal ganglia with connectivity between the basal 227 ganglia and thalamus. 228

229
The previously computed connectivity matrices were thresholded in an effort to retain the 230 most reliable tracts, while discarding potential false positive connectivity. As noted in previous 231 studies, defining a threshold is a non-trivial task (Shadi et al., 2016;Zhang et al., 2018). If the 232 chosen threshold is too low, tracts that are not reliably identified may remain, including those that 233 do not exist in reality, but if it is too high, legitimate connections may be discarded. Common 234 approaches include choosing an arbitrary threshold such that the majority of the subjects to be 235 analyzed retain the same connections (Li et al., 2009) or by sweeping through a range of thresholds 236 (Li et al., 2012). More recently, a test-retest metric was proposed, wherein reliability was evaluated 237 across a range of thresholds and a final threshold was selected where the change was at a minimum 238 (Zhang et al., 2018). Here, we selected our threshold by following a similar test-retest reliability 239 procedure, using the tract overlap (wDSC) as the reliability measure. To select our threshold, we 240 first identified the range where the rate of increase in tract overlap was at a minimum. We then 241 selected the first threshold where the rate of increase in tract overlap is 0 and determined the 242 corresponding average TD threshold. Supplementary Fig. 1

245
Processing and analysis of the HCP unrelated dataset followed the same workflow as 246 before with the test-retest dataset. As before, known connectivity of the subcortical connectome 247 was both visually and quantitatively assessed. TD matrices were computed for each subject as 248 before, and further separated by hemispheric connectivity to compare with previous findings. With 249 only a single acquisition session in the unrelated dataset, an average TD matrix was computed 250 across subjects, and a Pearson's correlation was performed against the average TD matrices of the 251 test and retest sessions to evaluate the similarity of the subcortical connectome. As with the test-252 retest dataset, the relationship between TD and the size of the subcortical structures was also 253 evaluated.

262
Motor network connectivity using the described tractography methods could successfully 263 recapitulate known connections as previously described in the literature (Fig. 2). Similarly, known 264 connections of both the associative and limbic network connectivity were also successfully 265 captured (Supplementary Fig. 2A and Supplementary Fig. 2B respectively). A wDSC of 0.58 was 266 selected as the final overlap threshold, which corresponded to a TD threshold of 6.5 AFD. Of the 267 known connectivity comprising the motor network, 78% (14 out of 18) of the identified 268 connections met the threshold. In the associative and limbic networks, 100% and 79% (11 out of 269 14) of the observed connectivity met the TD threshold respectively. Connectivity failing to meet 270 this threshold was commonly found between a thalamic nucleus (which was often small) and 271 another subcortical structure (see Supplementary Table 1  Identified connections were also visually inspected, examining the connected structures 278 and their trajectories. In observations of tract density, it was previously noted that basal ganglia 279 connections (e.g. non-thalamic ROI to non-thalamic ROI) were denser, while connections between 280 the basal ganglia and thalamus (e.g. non-thalamic ROI to thalamic ROI) were sparser. Visual 281 inspection of the known trajectories, reflected the previous observation of denser connectivity 282 between basal ganglia structures, which are also shorter and more direct. Conversely, connections 283 between the basal ganglia and thalamus were sparser with longer and more curved trajectories. 284 These longer trajectories increased the potential for intersecting GM structures between the basal 285 ganglia and thalamus as was the case for connections between the ventrolateral anterior nucleus of 286 the thalamus (VLa) and GPi (Fig. 3), as well as between STN and globus pallidus externa (GPe) / 287 GPi (Fig. 4). It was observed that certain thalamic nuclei were more difficult to reach, as 288 trajectories would have to pass through other surrounding thalamic nuclei. Some spurious 289 streamlines were also noted (e.g. streamlines that looped in the brainstem). Full descriptions of 290 known subcortical connections can be found in Supplementary Table 1. 291 15 the same ones that demonstrated a low TD. For connections between basal ganglia structures, a 299 low wDSC was only found between the caudate and amygdala, where the tract was sparse and 300 trajectories would have had to pass through other GM structures (e.g. putamen, GPe, GPi). 301

Reliability of known subcortical connections
Additionally, lower overlap was observed in the connections between the basal ganglia and 302 thalamus, in particular connections to the mediodorsal nuclei of the thalamus (MD), which was 303 more difficult to reach and in which trajectories also had to potentially traverse other nuclei of the 304 thalamus. Despite the overlap observed in a few connections, good overall reliability was 305 demonstrated for connectivity of each network, with similar measurements for each hemisphere. 306

307
Connectivity matrices were created for both test and retest sessions between subcortical 308 structures for all subjects. Visual assessments were first performed, followed by quantitative 309 evaluation of all intra-hemispheric subcortical connections. Reliability of the subcortical 310 connectome was also evaluated and noted to be similar to what was previously assessed for known 311 connections. 312

313
Connectivity matrices for test and retest sessions were created from the computed TD 314 between subcortical structures for all subjects. A visual assessment of the computed matrices was 315 first performed, followed by a quantitative evaluation of the computed TDs. Matrices were 316 observed to be similar across subjects and test-retest sessions ( Supplementary Fig. 3). Connections 317 between basal ganglia structures were often denser (TD = 2.33 log(AFD) and 2.24 log(AFD)) in 318 left and right hemispheres respectively, averaged across test and retest sessions) than connections 319 16 between the basal ganglia and thalamus (TD = 1.34 log(AFD) and 1.26 log(AFD)) in left and right 320 hemispheres respectively, averaged across test and retest sessions). The difference between the 321 two groups was also corroborated with a one-way ANOVA for both the left (F = 19.19, p < 0.05) 322 and right (F = 3.75, p < 0.05) hemispheres. Fig. 5A demonstrates the average tract densities across 323 the test-retest session. 324 The influence that the volume of terminal subcortical structures had on TD was also 325 assessed. By plotting the average TD against the volume of the two terminal subcortical structures 326 ( Fig. 5B), the average tract density was observed to increase as the volume of one of the two 327 structures increased. Performing an ordinary least squares regression, we identified a positive 328 linear relationship between the TD and the size of the two subcortical structures (r 2 = 0.210, p < 329 0.05). 330

331
Using the previously computed connectivity matrices, an average TD matrix across 332 subjects was created for the test and retest sessions independently (Fig. 6A). Individual subject 333 matrices, as well as average session matrices were visually inspected, and minimal differences 334 were observed between test and retest sessions. A linear relationship was identified between the 335 TDs of the test and retest sessions (Supplementary Fig. 3A; ρ = 0.997, p < 0.05), with greater 336 variability observed between sessions when the TD was low. Connectivity was further divided by 337 hemisphere (i.e. intra-hemispheric left vs right) and a box plot was created to visually compare 338 test and retest tract densities (Fig. 6B). No differences in hemispheric connectivity were observed 339 between test and retest sessions, which was further corroborated after performing a paired t-test 340 between average test and retest TDs (t=1.52, p = 0.264 and t=1.06, p = 0.293 for left and right 341 hemispheres respectively). 342 To quantify the consistency of TD between test and retest sessions, we computed the 343 percent change of corresponding subcortical connections between sessions finding on average a 344 percent change of 36% and 32% for the left and right intra-hemispheric TD respectively (Fig. 6C). 345 As previously noted, in test-retest pairs where TD was low, greater variability was observed. 346 Correspondingly, a greater absolute percent change was more likely to be associated with a sparser connectivity respectively between test and retest TDs (Fig. 6D). 352 Voxel-wise spatial overlap of tracts (calculated via wDSC), was also computed between 353 test-retest pairs as another measure of reliability. We observed an average wDSC of 0.46 and 0.45 354 for the left and right intra-hemispheric connectivity respectively. We also plotted and performed 355 a Spearman's correlation between wDSC and average TD across both sessions where wDSC is 356 expected to increase with TD before plateauing. As expected, we identified a sigmoid relationship 357 (Supplementary Fig. 3C; ρ = 0.950, p < 0.05) between the two (Fig. 3E), with good overlap 358 achieved at a TD around 2.0 log(AFD) and reaching maximum overlap at approximately around 359 2.5 log(AFD). The overlap remained low while the TD was less than 0 log(AFD) before slowly 360 increasing until the overlap began to peak at a log-transformed TD of around 2 log(AFD). As 361 wDSC was highly correlated with TD, we further separated and evaluated the wDSC to the two 362 previously identified groups. As with observations from known subcortico-subcortical 363 connections, we noted better overlap in basal ganglia connectivity (average wDSCs = 0.75 and 364 0.72 for left and right hemispheres respectively) than in connectivity between the basal ganglia 365 18 and thalamus (average wDSCs = 0.46 and 0.44 in left and right hemispheres respectively). The 366 difference observed between the two groups was also supported with a one-way ANOVA for both 367 the left (F = 33.53, p < 0.05) and right (F = 41.64, p < 0.05) hemispheres. 368 thalamus, as well as between the STN and both GPe and GPi ( Fig. 3 and Fig. 4), which have both 401 been well studied. Full details regarding our observations are found in Supplementary Table 1.  402 First, focusing on the connection between GPi and VLa/VLp, one plausible trajectory 403 observed was the ansa lenticularis (Fig. 3). Similar origins of the ansa lenticularis observed in this 404 study from tractography reconstruction have been previously described from tract-tracing studies 405 in NHPs, with similar projections from the GPi to the VLa and VLp (Gallay et al., 2008). The 406 connection has also been described by (Parent and Carpenter, 1996), noting a trajectory that "forms 407 a well-defined bundle on the ventral surface of the pallidum…" curving around the posterior limb 408 of the internal capsule before continuing posteriorly. Along this trajectory, the ansa lenticularis is 409 20 known to converge with the lenticular fasciculus in the fields of Forel to form the thalamic 410 fasciculus, which continues to VLa and VLp. While we were able to note the termination in the 411 VLa and VLp in our observations, we were unable to delineate the transition from ansa lenticularis 412 to thalamic fasciculus. Further, sparse connections were observed to cross the region of the internal 413 capsule to connect the GPi with VLa and VLp, which may be part of the lenticular fasciculus. 414

Observations in HCP unrelated dataset
Another notable connection observed was between the STN and globus pallidus, including 415 both the internus (GPi) and externus (GPe) segments (Fig. 4) (Cousineau et al., 2017). The same study also examined reliability of 493 cortico-subcortical connections using an ROI defining a general cortical region (e.g. sensorimotor 494 cortex, associative cortex, limbic cortex) to an ROI defining a subcortical structure (e.g. caudate, 495 putamen, thalamus), where poor reliability of cortico-subcortical connectivity was noted. The poor 496 reliability was attributed to a combination of the quality of atlas used to define ROIs, partial 497 volume effects, motion and the low resolution of the data, as well as the difficulty of performing 498 tractography to the subcortical brain regions where structures are generally smaller (Cousineau et 499 al., 2017). 500 24 Overall, we demonstrated that the methods employed in the present study can produce 501 subcortical connectomes with comparable reliability to those that have been used to study cortico-502 cortical and cortico-subcortical connectomes. Although we observed worse reliability in 503 connectivity between the basal ganglia and thalamus, we believe this is primarily attributed to the 504 sparse connections resulting from imposed GM constraints. Further some nuclei of the thalamus 505 were more difficult to reach, with other nuclei present along the expected trajectory, while others 506 were smaller in size. Some changes to the tractography algorithm (e.g. angle, maximum streamline 507 length, etc.) or further optimization of constraints may be required to improve the overall 508 connectivity with the thalamus. In subsequent work, the described framework can be leveraged to segmentations, ROIs may be defined by the underlying nuclei, whereas in structural connectivity-533 based segmentations, ROIs may be related to regional connectivity. With many segmentation 534 schemes are readily available and there are multiple considerations to contemplate, it is important 535 to note that choice of segmentation is often dependent on the aims of the specific study (Arslan et 536 al., 2018). 537 Segmentation accuracy is also important for capturing the true underlying subcortical 538 structure, with size influencing the reconstructed connectivity as has been previously noted 539 (Sotiropoulos and Zalesky, 2019) and also observed in this study. An ROI larger than the structure, 540 especially in the small subcortical region, may overlap with other structures or extend into the WM 541 or CSF. On the other hand, an ROI smaller than the structure may exclude connectivity that does 542 not reach the boundaries, although some of this is alleviated by using the radial search strategy 543 employed. Due to the relationship between ROI size and tract density (see Fig. 5B), the ability to 544 identify connections in small structures is challenging and some expected connections may remain 545 unidentified. Nonetheless, the segmentations present incorporated data from histology and largely 546 reflect the underlying anatomy. 547 We chose to include a GM exclusion criteria in our tractography algorithm, removing 548 connectivity passing through other subcortical structures along its trajectory. This choice was made ROIs, the ability for tracts to pass through GM will be challenging due to the reduced anisotropy 563 in GM (e.g. in the pallidum). 564

27
In a similar manner, the inclusion of WM priors as wayward ROIs may improve anatomical 565 accuracy. With a data-driven approach to identifying the subcortical connectome, we had observed 566 the presence of major tracts (e.g. CST), spurious streamlines, and in some instances, multiple 567 trajectories between two subcortical structures. By using a WM prior, trajectories from major tracts 568 and spurious streamlines could be filtered, while individual trajectories can be isolated. In a 569 previous study of tractography reproducibility, a suggestion was made to include the use of 570 anatomical priors as guidance to improve identification of connectivity (Maier-Hein et al., 2017). identified connections in vivo is a known challenge, given the limited ability to compare to ground 585 truth trajectories, which have been conventionally identified using tract-tracing in experimental 586 28 animals. While the most accurate comparisons would be performed between tract-tracing and 587 tractography on the same brain, this is not feasible in humans. Fortunately, connections between 588 regions are highly similar across different primates (Grisot et al., 2021). In the current study, we 589 limited our investigation to known connections between ROIs of the basal ganglia and the 590 thalamus in order to compare our observations with previously described trajectories. As a result, 591 we did not explore the complete network circuitry to other regions of the brain (e.g. brainstem, 592 cerebellum, cortex, etc). Some of these unexplored regions contain important nodes, such as 593 connectivity with the hypothalamus (Haber et al., 1993) and pedunculopontine nucleus of the 594 brainstem (DeLong and Wichmann, 2007). Future work should explore the network circuitry more 595 comprehensively, which should be increasingly feasible with increasing availability of brain 596 atlases. 597 598

599
In this study, we demonstrated that identifying the subcortical connectome using a data-600 driven probabilistic approach with in vivo tractography was both feasible and reliable, with a 601 particular focus on the assessment of known connections that have been previously described. The data used in this study are available as part of the publicly available Human 616 Connectome Project S1200 release (https://humanconnectome.org/study/hcp-young-adult). 617 Analysis was performed with JupyterLab using Python (version 3.8). Code and results are 618 available upon request and will be made openly available upon acceptance for publication. 619