Towards deep learning for connectome mapping: A block decomposition framework

We propose a new framework to map structural connectomes using deep learning and diffusion MRI. We show that our framework not only enables connectome mapping with a convolutional neural network (CNN), but can also be straightforwardly incorporated into conventional connectome mapping pipelines to enhance accuracy. Our framework involves decomposing the entire brain volume into overlapping blocks. Blocks are sufficiently small to ensure that a CNN can be efficiently trained to predict each block's internal connectivity architecture. We develop a block stitching algorithm to rebuild the full brain volume from these blocks and thereby map end-to-end connectivity matrices. To evaluate our block decomposition and stitching (BDS) framework independent of CNN performance, we first map each block's internal connectivity using conventional streamline tractography. Performance is evaluated using simulated diffusion MRI data generated from numerical connectome phantoms with known ground truth connectivity. Due to the redundancy achieved by allowing blocks to overlap, we find that our block decomposition and stitching steps per se can enhance the accuracy of probabilistic and deterministic tractography algorithms by up to 20-30%. Moreover, we demonstrate that our framework can improve the strength of structure-function coupling between in vivo diffusion and functional MRI data. We find that structural brain networks mapped with deep learning correlate more strongly with functional brain networks (r = 0.45) than those mapped with conventional tractography (r = 0.36). In conclusion, our BDS framework not only enables connectome mapping with deep learning, but its two constituent steps can be straightforwardly incorporated as part of conventional connectome mapping pipelines to enhance accuracy.


Introduction
Mapping the human connectome is a major goal in neuroscience (Bullmore and Bassett, 2011, Behrens and Sporns, 2012, Fornito et al., 2013, Sporns et al., 2005).Diffusion-weighted magnetic resonance imaging (dMRI) is currently the only non-invasive imaging modality that enables in vivo structural connectome mapping.Individual fiber bundles can be reconstructed from dMRI data using a fiber tracking process called tractography (Mori et al., 1999, Catani et al., 2002).The result of this process is millions of streamlines that trace out the trajectories of white matter fiber bundles (Côté et al., 2013a, Neher et al., 2015, Maier-Hein et al., 2017).The number of streamlines interconnecting each pair of regions comprising a predefined cortical parcellation (Zalesky et al., 2010) is typically enumerated to provide a pairwise measure of interregional connectivity and arranged to form a connectivity matrix (Sotiropoulos and Zalesky, 2017).Network analyses of these connectivity matrices have yielded new insights into brain disorders (Bassett and Bullmore, 2009, Crossley et al., 2014, Fornito et al., 2015), cognition (van den Heuvel and Sporns, 2013b) and neurodevelopmental processes (Khundrakpam et al., 2012, Cropley et al., 2016).
The accuracy of connectome reconstruction can be degraded by spurious connections between pairs of regions that are not genuinely connected by a fiber bundle (false-positive connections) as well as missing connections (false-negative connections) (Reveley et al., 2015, Thomas et al., 2014a, Sotiropoulos and Zalesky, 2017, Calabrese et al., 2015, Côté et al., 2013a).The impact of falsepositive connections is particularly detrimental to the topological analysis of binary connectomes (Zalesky et al., 2016).While the use of weighted connectomes can overcome some of these challenges, numerous applications in network neuroscience require knowledge of whether pairs of regions are connected or not, independent of connectivity strength.This simplification can speedup computation and aid the interpretation of some graph-theoretic measures.For example, generating null networks for binary connectomes can be efficiently and straightforwardly achieved with the Maslov-Sneppen rewiring algorithm, whereas extending this algorithm to weighted networks is computationally expensive and mandates several assumptions (Rubinov and Sporns, 2011).
Evaluation of connectome mapping methodologies using numerical connectome and fiber phantoms (Neher et al., 2015, Côté et al., 2013b, Girard et al., 2014, Fillard et al., 2011) have identified an inherent trade-off between false-positive and false-negative connections.In particular, connectome mapping methods based on probabilistic tractography are typically hampered by an abundance of false-positive connections, leading to reduced specificity in connectivity localization due to streamline dispersion, while deterministic tractography can yield connectomes with missing connections (Zalesky et al., 2016).Alternative methods have been successfully used to identify missing connections (Markov et al., 2013).While conclusions drawn from studies utilizing numerical connectome phantoms are not necessarily entirely representative of performance on in vivo data, mounting evidence derived from these studies indicates that mapping connectomes with both high specificity and high sensitivity is very challenging with current tractography algorithms (Maier-Hein et al., 2017, Thomas et al., 2014b).Multi-fiber deterministic tractography algorithms have been shown to yield the most accurate reconstruction of binary connectomes, but with an F-measure of approximately 0.4 (Sarwar et al., 2018), the scope for improvement in reconstruction accuracy remains substantial.
Although tractography has inherit limitations and challenges, it remains the only non-invasive method to map structural connectivity using in vivo data (Schilling et al., 2018).Tractography has yielded fundamental insights into white matter connectivity architecture and has proven crucial to the emergence of the field of connectomics (Sporns et al., 2005, Hagmann et al., 2008).The overarching aim of the current study is to provide a framework to improve the accuracy of connectome mapping with tractography.
Machine and deep learning techniques can potentially improve the accuracy of fiber tracking and connectome reconstruction.These techniques have recently demonstrated utility in learning the mapping between dMRI signals and local fiber orientations.This has given rise to the concept of machine (and deep) learning tractography, whereby a neural network is trained to successively predict the direction in which to propagate streamlines based on sequences of dMRI signals (Neher et al., 2017).A key advantage of deep learning approaches is that they do not mandate a specific diffusion model or parameterization of the diffusion signal and may thus provide improved robustness to variation in data acquisition schemes and noise.Numerous network architectures have been considered in relation to this new fiber tracking paradigm, including recurrent (Poulin et al., 2017, Benou andRiklin-Raviv, 2018), convolutional (Koppers et al., 2017) and multi-layer perceptron (Wegmayr et al., 2018) architectures.
In the present study, we aim to train a convolutional neural network (CNN) to learn the mapping between dMRI signals and whole-brain connectivity.Connectomes are typically mapped based on the locations of streamline endpoints reconstructed with tractography.If the endpoints of sufficiently many streamlines reside within two separate regions, those regions are deemed to be structurally connected, irrespective of the geometrical properties of the streamlines (Maier-Hein et al., 2017, Côté et al., 2013a).Streamline geometry is therefore not explicitly modelled when mapping connectomes, particularly connectivity matrices.For this reason, we seek to investigate whether deep learning can be applied to learn the relation between dMRI signals and connectivity matrices.
Learning the relation between dMRI signals and whole-brain connectivity is a challenge that is of substantially higher dimensionality than current applications of machine and deep learning in dMRI.
The input space is defined by a series of three-dimensional images acquired across possibly hundreds of gradient directions and strengths, each comprising thousands of relevant white-matter voxels (pixels); while the output space spans thousands of unique elements comprising a connectivity matrix.
Given the high dimensionality of the input and output spaces, we first aim to develop a tractable decomposition framework in which the dMRI data are decoupled into thousands of blocks, where each block delineates a contiguous sub-image comprising × × voxels.In this way, each row/column of a block's local connectivity matrix corresponds to one of the − − 2 whitematter voxels residing at its boundary.For example, a whole-brain dMRI dataset with dimensions of 150 × 150 × 150 × 120 requires 100 million inputs (voxels representing brain volume), which is computationally infeasible from the perspective of a CNN.Block decomposition reduces the input space to 216 when using blocks of width 6 voxels enabling efficient training and prediction.Block size is determined to ensure that the chosen deep learning architecture can be tractably trained to learn the mapping between the dMRI sub-images and intra-block connectivity.The local connectivity matrices estimated for each block are not of direct interest, but rather serve as buildings blocks to estimate connectivity.In particular, we propose a simple block stitching algorithm to stitch together the local connectivity matrices and thereby yield a whole-brain connectivity matrix.The motivation for using blocks is to minimize the error introduced in voxel-level estimations.
In this study, we train a CNN to predict an arbitrary block's local (internal) connectivity matrix based on the block's minimally pre-processed dMRI signals.Furthermore, to evaluate the performance of our block stitching algorithm in isolation of CNN performance, we directly reconstruct each block's local connectivity matrix with streamline tractography.In this case, we find that our block stitching algorithm reconstructs whole-brain connectivity matrices that are more accurate than those yielded with conventional connectome mapping pipelines based on deterministic and probabilistic tractography.This is due to the added robustness to noise achieved by allowing blocks to overlap with each other.The first contribution of this paper is thus the proposal of a new connectome mapping methodology based on block decomposition and stitching.The second contribution is to show how this block decomposition framework could facilitate tractable connectome reconstruction with deep learning.Figure 1 provides a schematic overview comparing our proposed framework to classic tractography-based connectome mapping.

Materials and Methods
This section is organized as follows: We first describe the diffusion MRI (dMRI) and resting-state functional MRI (rsfMRI) data that was used in this study.Note that to evaluate the accuracy of the structural connectivity matrices mapped in this study, we considered the extent to which they correlated with functional connectivity matrices inferred from rsfMRI data (Honey et al., 2007), where stronger structure-function coupling was taken as an indicator of better performance.Next, we consider the estimation of each block's local (internal) connectivity matrix.This was achieved with: i) classic streamline tractrography; and, ii) by training a CNN to learn the mapping from dMRI subimages and their internal connectivity.At this point, we describe the details of the CNN architecture utilized in this study.We then describe in detail the probabilistic and deterministic variants of our block decomposition and stitching (BDS) framework for connectome mapping and define the concepts of block and stride.Finally, we describe how the accuracy of the whole-brain structural connectivity matrices reconstructed with BDS was evaluated, including: i) the use of simulated numerical connectome phantoms with known ground truth connectivity matrices; and, ii) the evaluation of structure-function coupling based on functional connectivity matrices inferred from rsfMRI data.Connectivity matrix reproduced from (van den Heuvel and Sporns, 2013a)
The MRI acquisition protocols and minimal pre-processing pipelines are described in detail elsewhere (Glasser et al., 2013, Sotiropoulos et al., 2013, Smith et al., 2013).In brief, for dMRI, 90 gradient directions were acquired at each of three b-values (1000, 2000 and 3000 s/mm 2 ) using spin-echo planar imaging.Minimal pre-processing included correction of eddy-current, head motion and gradient-nonlinearity distortions, followed by transformation of the corrected dMRI volumes to native structural space and rotation of gradients vectors.For rsfMRI, multiband gradient-echo planar imaging was acquired for a period of 14.4 minutes (TR = 720ms, multiband factor = 8), with eyes open and relaxed fixation on a projected bright cross-hair on a dark background.Four runs were acquired, two runs (right-to-left and left-to-right phase encodings) in one session and two in another session.Minimal pre-processing included data cleaning with ICA+FIX (Salimi-Khorshidi et al., 2014), correction of distortion in the phase-encoding direction with TOPUP (Andersson et al., 2003) and nonlinear transformation to Montreal Neurological Institute (MNI) standard space.The four runs were temporally demeaned and then concatenated, effectively yielding a total of 4800 volumes for each individual.

Mapping Block Connectivity
The first step of our framework involves uniformly decomposing the dMRI images into thousands of overlapping blocks, where each block delineates a contiguous sub-image comprising × × voxels.In the current study, blocks were uniformly positioned throughout the entire brain volume.
Further details pertaining to block delineation and positioning are provided in Section 2.3.1.We aimed to map each block's local connectivity architecture, using the local dMRI data that it encapsulated.More specifically, we aimed to represent each block's local connectivity with a × symmetric connectivity matrix, where = − − 2 denotes the number of voxels residing on the periphery (faces) of each block.Connectivity matrices for each block were mapped using: i) classic streamline tractography; or, ii) CNN.

Tractography-based mapping of block connectivity
In the absence of a ground truth for the in vivo data, conventional streamline tractography was used to map intra-block connectivity (i.e.map a connectivity matrix for each block).More specifically, we estimated multiple fiber orientations for each white matter voxel using constrained spherical deconvolution (CSD) (Tournier et al., 2007) and then used deterministic fiber tracking to propagate streamlines within the volume encapsulated by each block.We opted for CSD-based deterministic tractography on the basis of recent recommendations (Sarwar et al., 2018).A total of 1000 streamlines were initiated from randomly chosen coordinates within each block (refer to Section 2.4.4 for implementation details of tractography).
For each streamline intersecting the periphery of a given block, we determined whether the streamline intersected at least two of the peripheral voxels positioned at the block's periphery.Peripheral voxels are voxels that are intersected when streamlines enter or exit a block's internal volume.
Streamlines satisfying this criterion were used to map the block's connectivity matrix.In particular, element , of the block's connectivity matrix enumerated the total number of streamlines intersecting peripheral voxels and .For streamlines that intersected more than two peripheral voxels comprising a given block, we only considered the two peripheral voxels that were intersected when the streamline first entered and last exited the block's internal volume.It should be noted that connections between adjacent voxels were discarded.Streamlines that did not traverse the interior of the block were also discarded.Finally, each block's connectivity matrix was binarized, such that element , was set to 1 if one or more streamlines were found between peripheral voxels and , otherwise this element was set to 0.
Figure 2(a) shows a schematic of the block connectivity mapping process based on streamline tractography.For illustrative purposes, we only show two spatial dimensions, and thus blocks effectively become squares in this scheme.

CNN-based mapping of block connectivity
A CNN was trained to learn the mapping between the minimally pre-processed dMRI data encapsulated by a block and the block's connectivity matrix.The aim of the CNN was therefore to learn the connectivity mapping process within a circumscribed spatial locale, from dMRI data to connectivity matrix.In the current study, each training sample comprised the dMRI data encapsulated by a block (input) and the block's local connectivity matrix (output).Trainings samples were drawn from randomly selected blocks in 10 dMRI datasets.The local connectivity matrix for each sampled block was mapped using streamline tractography, as described above (Section 2.2.1).In this way, the connectivity mapping derived with streamline tractography served as the ground truth for which the CNN was trained to predict.The CNN used the dMRI data to directly determine each block's connectivity matrix, eliminating the need for modelling of fiber orientations, tractography and assignment of streamlines to pairs of nodes.The details of the CNN architecture used in the study can be found in Appendix.

Blocks and strides
The first step in our block decomposition and stitching (BDS) algorithm involves decomposing the dMRI images into blocks of dimension × × and estimating intra-block connectivity independently for each block.The concept of stride is used to determine the extent of overlap between neighboring blocks.Denoted with , stride defines the extent of separation between two adjacent blocks, i.e. two blocks are separated by voxels along a given plane.For a stride of one, each brain voxel is associated with a unique block; for a stride of two, only every second voxel along any given plane is associated with a unique block, and so on.Therefore, shorter strides and/or larger block dimensions result in increased overlap between neighboring blocks, which in turn provides a degree of redundancy to counteract the effects of noise.Nevertheless, while short stride lengths and large block dimensions can improve robustness to noise, they increase the computational burden.Therefore, stride length and block dimension are important parameters that enable a trade-off between accuracy and running time.(blocks appear as squares in two dimensions), where the blue square indicates the boundaries of one such representative block.A 20 × 20 connectivity matrix was then mapped for this representative block, where each row/column indexes one of the block's peripheral voxels, labelled from 1 to 20.The representative block encapsulates a portion of a single fiber bundle that comprises streamlines intersecting peripheral voxels 3, 13 and 14.Accordingly, elements (3,13), (3,14) and (13,14) of the block's connectivity matrix (rightmost) were set to one.(b) Schematic of the block decomposition process.The 14 × 14 image (left) is decomposed into 9 blocks (right) of size 6 × 6 using a stride of 4. This stride length results in an overlap of 12 voxels between consecutive blocks.Blocks are indexed using row ( ) and column ( ) subscripts.Colored boxes represent the mapping of voxels to blocks and their locations.(c) Schematic of the block stitching process used to estimate connectivity.The block with index is chosen as the seed block (red) and the connection between peripheral voxels 4 and 14 (yellow) within this block is selected as the seed connection, where denotes the connection's index.The unit vector, denoted by , determined for this seed connection points in the direction of the next-selected block (blue), denoted .For Deterministic-BDS, unit vectors for all connections comprising the candidate block (blue) are estimated, and the connection with unit vector that forms the most acute angle with the preceding block's unit vector is selected (green line) (d) Schematic of iterative block stitching resulting in a block chain.In this simplified example, each block comprises a single connection (i.e.= 1).Block is determined based on the unit vector within block .A block chain is an ordered set of blocks ( , , , … ) that defines and connection between gray matter regions.
Once the dMRI image is decomposed into blocks based on a desired stride length, the next step is to map an intra-block connectivity matrix for each block.Intra-block connectivity mapping is performed independently for each block, as described in Section 2.2.We use ! to denote the connectivity matrix mapped for the block with index .Finally, blocks are stitched together to form chains of blocks, or block chains (Figure 2(c,d)), which characterize end-to-end connectivity between gray matter regions and can be used to map a connectivity matrix.Block chains are initiated from a seed block and iteratively propagated block-by-block in much that same way as a conventional tractography streamline is propagated.In the following sections, we describe the selection of seed blocks (Section 2.3.2) and propagation of block chains (Section 2.3.3).

Seeding methodology
In the current study, block chains were successively seeded (initiated) from either: i) all possible blocks; or, ii) only blocks that resided within gray matter regions comprising the cortical parcellation.
These two seeding methodologies are reminiscent of the distinction between whole-brain and gray matter seeding strategies used for conventional streamline tractography (Li et al., 2012).Once a given seed block is selected, we initiate a block chain for each connection comprising the block.Therefore, each seed block can potentially be the origin of multiple block chains.Given that intra-block connectivity is inherently undirected, for each connection, a block chain is tracked in both diametrically opposing directions associated with each connection.

Propagating block chains
Block chains are propagated in much the same way as streamlines are propagated under conventional tractography.While a conventional streamline is parameterized by an ordered set of spatial coordinates, a block chain is an ordered set of blocks ( , , , … ; Figure 2(d)) that can potentially define a connection between gray matter regions.For a given seed block, denoted , we propagate a separate block chain for each connection comprising the block's connectivity matrix ! .We use = 1, … , to index the connections within a given block (i.e. the non-zero elements of the connectivity matrix).For each such connection in , we determine a unit vector probabilistic choice that shares the same rationale as probabilistic tractography.

Deterministic criterion:
For each new block , we select the connection from the set of all candidate connections # $ % ,…,& that minimizes the angle formed with the connection in the current block .This criterion ensures that a block chain does not form sharp angles that might be considered biologically unrealistic.Given that connections are inherently undirected, the unit vector assigned to connection in block is oriented such that it is consistent with the direction of .The block chain is terminated at block if: i) the connectivity matrix associated with block comprises no connections; or, ii) the angle formed between exceeds a prescribed angle threshold, denoted '.We refer to block chain propagation under this deterministic criterion as deterministic block decomposition and stitching (Det BDS).

Probabilistic criterion:
For each new block , we first determine the angle formed between each candidate connection # $ % ,…,& and the current direction .We then randomly sample a connection from the subset of candidate connections that satisfy the prescribed angle threshold '.The block chain is terminated at block if no candidate connections satisfy the angle threshold.We refer to block chain propagation under this probabilistic criterion as probabilistic block decomposition and stitching (Prob BDS).
The CNN is trained to predict each block's local connectivity matrix (Section 2.2.2), which is then used to generate block chains using Det BDS or Prob BDS, referred as Det CNN-BDS and Prob CNN-BDS, respectively.CNN predicts the probability of every connection in a block, which can be utilized to generate block chains.Similar to Det BDS, from the set of all candidate connections # $ % ,…,& that satisfy the prescribed angle threshold ', the connection with the highest probability predicted by CNN is selected.We refer to block chain propagation under this deterministic criterion as Det-2 CNN-BDS.The nomenclature used to label each variation of our BDS framework is listed in the Appendix.The block chain generation process is described using pseudocode in Supplementary Table S1.

Mapping Connectomes
Once a sufficient number of block chains have been generated, they can be used to map a whole-brain connectivity matrix for a given cortical parcellation.For each block chain with starting and ending blocks that reside within the confines of the cortical parcellation (gray matter), we first identify the peripheral voxels in each of these two blocks that are associated with the intra-block connection used by the block chain.If one peripheral voxel resides in cortical region and the other in region , we increment the block chain count at element ( , ) of the connectivity matrix.In this way, the wholebrain connectivity matrix quantifies the number of block chains that are found to interconnect each pair of regions.Block chains with starting and/or ending blocks that terminate prematurely before reaching the cortical parcellation are discarded.

Evaluation Methodology
We applied our BDS framework to map connectivity matrices for: i) simulated connectome phantoms (Section 2.4.2); and ii) in vivo dMRI data acquired in healthy adults (Section 2.4.3).The connectivity matrices mapped for the simulated phantoms were benchmarked against the known ground truth and compared to connectivity matrices mapped with conventional streamline tractography.For the in vivo data, we tested the extent of correlation between functional connectivity inferred from rsfMRI and the connectivity matrices mapped with BDS as well as conventional tractography (structural connectivity).This correlation was computed across all pairs of regions.Structural connectivity matrices that showed a stronger correlation with functional connectivity were deemed more accurate than those yielding weaker structure-function coupling.

Evaluation using numerical connectome phantoms
Using an established framework (Sarwar et al., 2018), we simulated a range of two-dimensional connectome phantoms that recapitulated key properties of the human connectome.In brief, each phantom comprised multiple tubular fiber bundles of varying diameters and lengths.The fiber bundles were positioned within the confines of circle, such that both ends of each fiber terminated at distinct locations on the circle's circumference.The circle's circumference was uniformly segmented into 20 areas (nodes), enabling each fiber bundle to be assigned to a pair of areas.In this way, a ground truth connectivity matrix could be mapped for each phantom, where rows/columns of the ground truth connectivity matrix denote areas on the circle's circumference.The procedure used to simulate fiber bundle trajectories is described in detail elsewhere (Sarwar et al., 2018).
We simulated dMRI signals for each voxel within the circle defining our phantom's spatial extent.
The dMRI signal for each voxel comprised a linear combination of isotropic ( ( ) and anisotropic ( ) ) signal components, respectively modeled using the ball (Behrens et al., 2003) and zeppelin (Alexander, 2008, Panagiotaki et al., 2012) models, with volume fraction f and 1-f respectively.In particular, the dMRI signal for the ith voxel for diffusion-weighted gradient * and b-value b was given by, .2 1 ) ; ( = -2(3 ; ) = ∑ 5 -2(6 7 8 9 : 6 ; <= ;∈? @ (1) where d is the diffusivity of the isotropic compartment, A + is the set of all fibers inside the ith voxel, B ; C is the prolate tensor defining the anisotropic diffusion at point = of path D ∈ AE, ./ is the echo time and .2 is the relaxation time.This follows the same methodology that we have described in detail elsewhere (Sarwar et al., 2018).The only difference is that in the current study, we additionally modelled the effect of echo and relaxation times (Neher et al., 2014, Ferizi et al., 2015).
Table 1 lists parameter settings of the dMRI signal model and image dimensions of the connectome phantom.We used a volume fraction of , = 0.2 (Neher et al., 2014).We evaluated a range of alternative volume fraction values elsewhere (Sarwar et al., 2018) and found that the relative performance of tractography algorithms was not substantially affected by this parameter choice.

Phantom-based performance evaluation
BDS and conventional streamline tractography were used to map connectivity matrices based on the dMRI data simulated for a range of connectome phantoms.Under BDS, element ( , ) of each connectivity matrix stored the number of block chains that interconnected nodes and , while this element stored the streamline count in the case of conventional tractography.Connectivity matrices were first binarized with respect to a block/streamline count threshold, such that any matrix element exceeding the count threshold was set to one, while all other elements were set to zero.The count threshold was incrementally increased from zero to the largest block/streamline count in the connectivity matrix under consideration.This yielded a range of binary connectivity matrices that varied in connection density.Each binary connectivity matrix was then benchmarked against the corresponding ground truth connectivity matrix (Sarwar et al., 2018).In particular, we evaluated the true positive (TP), false positive (FP) and false negative (FN) rates as well as the F-measure as a function of the block/streamline count threshold.The F-measure was computed as 2TP / (2TP + FP + FN).Due to symmetry, these measures were evaluated using only the upper triangular elements of the connectivity matrices.(Hagmann et al., 2008, Honey et al., 2009, Koch et al., 2002, Seguin et al., 2018), and thus the extent to which a structural connectivity matrix correlates with functional connectivity can provide a proxy measure of connectome mapping performance.Using this criterion, we compared connectivity matrices mapped with BDS and conventional tractography.It is important to note that the extent of structure-function coupling can also be impacted by the pipeline used to map functional connectivity.In this study, we used Spearman's rank correlation coefficient to assess the extent of correlation between a pair of structural and functional connectivity matrices.Correlation was computed across all upper triangular elements of the connectivity matrices.
Functional connectivity matrices were mapped using an established pipeline that is described elsewhere (Seguin et al., 2018).In brief, the Pearson correlation coefficient was used to compute functional connectivity between all pairs of regions comprising the same cortical parcellation used to map the structural connectivity matrices.This involved forming regionally-averaged time series for each region.These time series spanned approximately 30 minutes, as a consequence of temporal concatenation across four runs (Section 2.1).A structural connectivity matrix was mapped for each of ten individuals.The ten matrices were then averaged and log-transformed to obtain a group-level structural connectivity matrix.Density-based thresholding was applied to the matrices and the correlation between log-transformed structural and functional connectivity was calculated for the whole brain and right hemisphere.
Connectivity matrices were thresholded and binarized using the same process described above.
A spherical harmonic order of 8 was used for CSD, and 2 million streamlines were generated for each in vivo dataset.Default parameters were used for step size (0.5 × H=-IE1for ifod2 and 0.1 × H=-IE1for rest), angle threshold (9 K × IL-D IE1-/ H=-IE1-) and FOD/FA threshold (0.1).Streamlines were uniformly seeded from the white matter mask, where white matter boundaries were dilated by 1 voxel to fill potential gaps between gray and white matter.This ensured that streamlines terminate within gray matter.

Results
We developed a block decomposition and stitching (BDS) framework to: i) enable scalable connectome mapping with deep learning; and, ii) improve the accuracy of conventional connectome mapping pipelines based on streamline tractography.Here we demonstrate and evaluate the performance of our proposed framework in both these applications.First, we evaluate BDS in the case where intra-block connectivity is mapped using conventional streamline tractography.This enabled us to evaluate the performance of the decomposition and stitching steps of our framework, in isolation of CNN performance.We evaluated the impact of key parameter choices inherent to BDS, including stride length, seeding methodology and angle threshold.Furthermore, we compared the connectivity matrices mapped with BDS to those generated with deterministic and probabilistic tractography.Performance was evaluated by benchmarking the reconstructed connectivity matrices to known ground truths (phantom data) and assessing the extent to which reconstructed structural connectivity correlated with functional connectivity (in vivo data).Second, we demonstrate how BDS can provide a scalable framework to enable connectome mapping with a CNN.In particular, we train a CNN to learn the mapping between dMRI signals and intra-block connectivity and evaluate the accuracy with which this mapping can be learned.Block size was determined to ensure that the chosen deep learning architecture can be tractably trained to learn the mapping between the dMRI sub-images and intra-block connectivity.

Evaluation of stride length, seeding methodology and angle threshold
We generated realistic connectome phantoms that recapitulated key properties of the human connectome (Sarwar et al., 2018).In particular, 100 two-dimensional phantoms were generated, each comprising 20 nodes positioned on the circumference of a circle as well as 10 three-dimensional phantoms, each comprising 40 nodes positioned on the surface of a sphere.Fiber bundle geometry and trajectories differed between each phantom.For each phantom, dMRI signals were simulated for each voxel (Section 2.4.1).The simulated dMRI images were decomposed into blocks of dimensions 6 × 6 (two-dimensional) and 4 × 4 × 4 (three-dimensional).Intra-block connectivity was mapped with tractography (Section 2.2.1).Blocks were then stitched together to yield block chains, which in turn were used to map connectivity matrices that were benchmarked to the ground truth connectivity matrices for each phantom.Blocks chains comprising fewer than three blocks were discarded.
We aimed to evaluate the impact of stride length, seeding methodology and angle threshold on the accuracy of connectome mapping with BDS.In this section, intra-block connectivity was mapped using Det CSD and block chains were propagated using Det BDS.The default parameters were an angle threshold of 20°, a stride of =1 and ROI-based seeding.In the following sections (Section 3.1.1-3.1.3),we investigate BDS performance to variation in these default parameter values (Supplementary Figures S2-S3).We report performance as ensemble averages across the set of all connectome phantoms.

Varying stride
Figure 3 shows the effect of variation in stride length on the reconstruction of block chains for a twodimensional connectome phantom.As stride length is increased, the extent of spatial overlap between neighboring blocks is reduced, which in turn reduces spatial redundancy and robustness to noise effects.It can be seen that detailed features of fiber bundle geometry and curvature are no longer recovered with increasing stride length.Figure 4(a) shows the accuracy of connectome reconstruction with Det BDS for stride lengths that vary from 1 to 4 (see Supplementary Figure S1(a) for Prob BDS).
As expected, accuracy tended to decrease with increasing stride length and the maximum F-measure was achieved with the shortest stride lengths of 1 and 2. Figure 4(a) also suggests that increasing stride length tended to compromise the sensitivity (true positives) of connectome reconstruction more than specificity.We conclude that a stride length of 1 is optimal with respect to the accuracy of connectome reconstruction, although longer strides enable accuracy to be traded for improvements in computational tractability.

Variation in seeding methodology
With conventional streamline tractography, the choice of voxels from which streamlines are seeded is crucial and can substantially influence connectome reconstruction (Catani et al., 2002, Girard et al., 2014, Mori and Van Zijl, 2002).We aimed to determine the extent to which the seeding of block chains under BDS impacts the accuracy of connectome reconstruction.As described in Section 2.3.2,we considered two distinct block chain seeding methodologies: i) uniformly seeding from all blocks (whole-image seeding); and, ii) seeding only from blocks comprising the cortical parcellation (ROIbased seeding).These two methodologies are broadly equivalent to seeding strategies that are currently used to seed streamlines under conventional tractography.These experiments were repeated for three-dimensional connectome phantoms that provided a more realistic representation of the brain's spatial embedding compared to the two-dimensional phantoms investigated above.We found that stride length, seeding methodology and angle threshold demonstrated the same behavior for the three-dimensional phantoms (Supplementary Figure S3).
Angle thresholds of 15 o -30 o yielded comparable accuracies for both two-and three-dimensional phantoms.Therefore, for all subsequent experiments, we use a fixed stride length of 1, ROI-based seeding and a fixed angle threshold of 20 o , unless otherwise noted.It should be noted that the angle threshold and lower stride aids in mitigating inconsistent intra-block connections and block-chains (Figure S4).

Evaluation of connectome reconstruction accuracy
Having found the optimal parameter settings for block chain propagation based on our BDS framework, we next sought to compare the accuracy of connectome reconstruction under BDS relative to conventional streamline tractography methods.Comparisons were performed using the same two-and three-dimensional connectome phantoms with known ground truth connectivity matrices that were used above.The accuracy of connectome reconstruction was quantified using ROC curves as well as the F-measure evaluated as a function of the block chain (or streamline) threshold.
Following an established evaluation methodology (Sarwar et al., 2018), we compared four streamline tractography algorithms (Det Tensor, Prob Tensor, Det CSD, Prob CSD; see Section 2.4.4) relative to deterministic and probabilistic variants of our BSD framework (Det BSD, Prob CSD).Importantly, for BSD, intra-block connectivity for each block was mapped using streamline tractography, and thus we aimed to quantify the improvement in accuracy that can be achieved simply by adding a block decomposition step to conventional connectome mapping pipelines.For the conventional tractography methods, we generated 2 million streamlines per phantom by seeding streamlines throughout the entire image (Section 2.4.4).In particular, for the three-dimensional phantom (Figure 5(b)), BDS yielded a 23% improvement in the F-measure compared to the best performing conventional tractography method (Det CSD).More specifically, Det BDS yielded an F-measure of 0.52 without any thresholding, whereas Det CSD yielded 0.4 (leftmost points on the F-measure curve).For the two-dimensional phantom (Figure 5(a)), BDS yielded a 30% improvement in the F-measure compared to the best performing conventional method.This suggests that the redundancy introduced by decomposing dMRI images into overlapping blocks can substantially improve the accuracy of connectome reconstruction, albeit at the expense of an increased computational burden.
We found that thresholding improved the accuracy of some of the BDS and conventional tractography methods, particularly the probabilistic variants.This was most likely owing to the abundance of false-positive connections mapped with probabilistic tractography and the probabilistic stitching method (Prob BDS), which can be alleviated to a certain extent by eliminating the weakest connections with thresholding.However, no amount of thresholding yielded an F-measure that exceeded the case of Det BDS without any thresholding.In other words, while a moderate level of thresholding (<20%) can improve the performance of probabilistic methods, this improvement is not sufficient to rival the performance of deterministic methods without any thresholding.Of the four conventional streamline tractography methods evaluated, we found that Det CSD yielded the highest F-measure for both the two-and three-dimensional connectome phantoms, whereas Det BDS was the best performing BDS variant (Figure 5).
As expected, tractography methods steered by the diffusion tensor yielded the least accurate connectome reconstructions (Figure 5).This is due to the abundance of crossing fiber bundles instantiated within the connectome phantoms, which cannot be adequately resolved with the diffusion tensor model.Indeed, we ensured that the proportion of voxels comprising crossing fibers in the connectome phantoms was matched to the same proportion in genuine dMRI data acquired in humans (Sarwar et al., 2018).This highlights the utility of fitting a multi-fiber model to each voxel that explicitly accounts for distinct fiber populations.We also found that deterministic tractography and the deterministic variant of BDS outperformed their probabilistic counterparts for both the two-and three-dimensional phantoms.This may be due to a reduction in specificity owing the spatial dispersion of streamlines that is associated with probabilistic methods.

In vivo dMRI data
Based on the evaluation performed on two-and three-dimensional connectome phantoms, we found that our proposed BDS framework outperformed conventional probabilistic and deterministic connectome mapping pipelines in terms of the accuracy of connectome reconstruction.However, performance and conclusions derived from simulated phantoms do not necessarily generalize to in vivo dMRI data.In this respect, while the connectome phantoms used in this study preserved the fiber complexity and topological properties of the human connectome, they were based on relatively simple fiber geometries as well as a somewhat unrealistic spherical spatial embedding of the cortex.
To address these limitations, we next applied our BDS framework to in vivo dMRI data provided by the Human Connectome Project (HCP; Section 2.1).In the absence of ground truth connectivity maps for the in vivo data, we evaluated the accuracy of connectome reconstruction based on an independent modality; namely, functional connectivity inferred from resting-state fMRI (rsfMRI) data.In particular, we assumed that the strength of correlation between structural (dMRI) and functional (rsfMRI) connectivity across the set of all pairs of regions (connections) served as a proxy measure of connectome reconstruction accuracy.It is well established that structural and functional connectivity are coupled (Honey et al., 2009), and thus any reduction in structure-function coupling can potentially indicate inferior mapping of structural connectivity, or the effects of pathology (Cocchi et al., 2014, Zhang et al., 2011).The latter was not the case here given that we only analyzed healthy adults.A functional connectivity matrix was mapped for each of ten individuals based on a cortical parcellation comprising 84 regions (Section 2.4.3).Structural connectivity matrices were mapped for the same individuals using: i) conventional streamline tractography (Det CSD; Section 2.2.1); and, ii) the deterministic and probabilistic variants of our BDS framework (Det & Prob BDS).For BDS, we uniformly seeded block chains using ROI-based seeding and a stride length of 1 with an angle threshold of 20 K .
Figure 6 provides a qualitative comparison of whole-brain tractograms as well as two fiber bundles reconstructed for a representative individual with conventional streamline tractography (upper panel) and our BDS framework (lower panel).While the primary purpose of BDS is connectivity matrix estimation, it can be seen that the reconstructed fiber bundle trajectories, as represented by block chains, provide a faithful approximation of fiber bundle geometry and accord reasonably well with the streamlines derived from conventional tractography.However, due to the discretization of intra-block connectivity, the spatial smoothness of each block chain is reduced compared to streamlines, resulting in a somewhat jagged path from one block to another.
Figure 7(a-c) shows group-averaged structural connectivity matrices estimated with conventional streamline tractography (Det CSD) as well as the deterministic and probabilistic variants of our BDS framework.It can be seen that BDS yields connectivity matrices that are typically greater in connection density than conventional tractography, although use of a probabilistic method would have likely reduced this discrepancy.All connectivity matrices show relatively strong off-diagonal elements, corresponding to prominent structural connectivity between contralateral homologues.
Spearman's rank-order correlation coefficient was computed to assess the extent of structurefunctional coupling for each structural connectivity matrix.The correlation coefficient was computed across all pairs of regions, since an absence of a direct structural connection does not imply that the pair of regions is not functionally connected (Honey et al., 2009).Subcortical regions were excluded due to the inherent difficulty in parcellating the subcortex, reduced subcortical signal-tonoise ratio and the inherent challenges of mapping subcortical white matter connections.Furthermore, given that the range of structural connectivity values spanned several orders of magnitude, these values were logarithm-transformed.
We evaluated the correlation coefficient assessing the strength of structure-function coupling as a function of the streamline and block chain threshold.Increasing this threshold yielded sparser structural connectomes.Thresholding was performed to eliminate connections comprising the fewest numbers of streamlines, under the assumption that these connections are spurious (van den Heuvel et al., 2017, Roberts et al., 2017).The functional connectivity matrices were not subjected to any thresholding.Figure 8(a) shows the correlation coefficient as a function of connection density for Det CSD (blue), Det BDS (red) and Prob BDS (orange) for whole brain, while Figure 8(b) shows the same analysis restricted to the right hemisphere.It can be seen that conventional tractography yields higher structure-function coupling than BDS for relatively low connection densities (<5-10%).However, for higher connection densities both variants of BDS yield significantly higher coupling compared to conventional tractography (p<0.001).This difference is particularly evident for the analysis restricted to the right hemisphere, where for a connection density of 55%, BDS yields a structure-function coupling of r = 0.62 (Det BDS) and r = 0.6 (Prob BDS), whereas conventional tractography yields r = 0.45.Connection densities exceeding 45-55% could not be assessed for Det CSD due to the inherent sparsity of the connectivity matrices.
We defined connectivity strength as the number of streamlines/block-chains per connection.The connection density threshold discards connections between pair of nodes, whereas the streamline/block-chain threshold only discards connections with strength that is less than the threshold selected.Therefore, the analysis was repeated for correlation coefficient as a function of streamline/block chain threshold (Supplementary Figure S5).The results were consistent under varying streamline/block-chain thresholds i.e. structure-function coupling for BDS was superior to Det CSD.
Therefore, BDS yielded structural connectomes that were more strongly coupled with functional connectivity for the majority of connection densities investigated, compared to connectomes mapped with conventional methods.It is important to remark that while the extent of structure-function coupling is not a direct measure of the accuracy of connectome reconstruction, the coupling strength across these two modalities provides a useful surrogate marker in the absence of a ground truth connectivity matrix.

Connectome mapping with deep learning
For all the results presented thus far (i.e.Det BDS and Prob BDS), intra-block connectivity was mapped directly using conventional streamline tractography.This enabled us to establish the utility of the block decomposition step per se as an augmentation to enhance the accuracy of conventional streamline-based connectome mapping pipelines.In this section, we demonstrate the utility of the BDS framework in facilitating tractable connectome mapping with deep learning.In particular, we now train a convolutional neural network (CNN) to learn the mapping from the dMRI images to intrablock connectivity.In this section, we introduce the nomenclature Det CNN-BDS and Prob CNN-BDS to denote the deterministic and probabilistic variants of our BDS framework, respectively, in the case where intra-block connectivity is predicted based on a CNN.

Phantom connectome mapping
To train the CNN to learn the mapping from dMRI data to intra-block connectivity, we simulated dMRI data for ten connectome phantoms and decomposed each image into blocks based on a stride length of 1. Intra-block connectivity was extracted using conventional tractography (Section 2.2.1).
Increasing the number of training data with inclusion of additional connectome phantoms had no significant effect on performance and only resulted in computational overhead in training the CNN.
We trained a separate CNN for the two-and three-dimensional phantoms, using block dimensions of 6 × 6 and 4 × 4 × 4, respectively.To improve the learning process, we found that excluding the feed-forward layer (layer 7-Table A1) for the three-dimensional CNN was beneficial.It should be noted that although the underlying structural connectivity was pre-defined for the phantoms, the training dataset generation was independent of the underlying connectivity.
After training the CNN, we repeated the evaluations described in Sections 3.1-3.2based on the same connectome phantoms, but now with intra-block connectivity predicted using the CNN.We found that the optimal parameter settings for BDS with CNN-based prediction of block connectivity were largely consistent with those reported in Sections 3.1.1-3.1.3.More specifically, i) both seeding methodologies (ROI-based and whole image) yielded similar results; ii) the impact of variations in stride length and angle threshold were highly comparable; and, iii) connectomes reconstructed with BDS (both CNN and streamline-based) were generally more accurate compared to those mapped with conventional streamline tractography, as indexed by the F-measure (Figure 9).Details can be found in Supplementary Figures S6-S11.

Human connectome mapping
Next, we applied the CNN trained using the three-dimensional connectome phantoms to predict intrablock connectivity in the in vivo dMRI data, with the goal of using a CNN to map the human connectome.However, we found that the CNN trained on the phantom data generally yielded inaccurate predictions when applied to the in vivo dMRI data.This was most likely due to complex fiber bundle geometries in the human brain (e.g.fanning and kissing fiber configurations) that were not adequately represented in our relatively simple connectome phantoms.The absence of these complexities in the training data may have precluded the CNN from learning the full extent of the mapping from dMRI data to intra-block connectivity, resulting in the poor performance observed when applied to white matter regions encompassing particularly complex fiber architectures.This motivates the need to develop more realistic connectomes phantoms in the future, possibly scaling up and improving the level of detail represented in existing phantoms (Fillard et al., 2011, Maier-Hein et al., 2017, Girard et al., 2014).
In the present study, to provide a proof-of-principle demonstration of how our proposed BDS framework facilitates mapping of the human connectome with deep learning, we used in vivo dMRI data to train a CNN.In particular, the training data comprised sub-images (blocks) from the in vivo dMRI data (input) and corresponding intra-block connectivity matrices (output) that were mapped with conventional streamline tractography.The training data derived from the in vivo dMRI data comprised one million blocks that were uniformly sampled from the white matter volume of ten healthy adults, using stride lengths ranging between 2 and 5 (Section 2.1).The input space for each training example was the dMRI signal for 60 diffusion-weighted gradients for each voxel comprising a particular dMRI sub-image (block), while the output space comprised the intra-block connectivity matrix, as mapped with conventional streamline tractography (Det CSD; Section 2.2.1).We limited the input space to 60 diffusion-weighted gradients to reduce the computational overhead.Once the CNN was trained, we evaluated its performance using dMRI data acquired in an independent group of 10 healthy adults.Table 2 shows the performance of the CNN averaged across 100,000 unseen blocks sampled from these individuals.To serve as a benchmark, Table 2 also shows the performance of the CNN trained using the dMRI data simulated from the connectome phantoms.As alluded to above, the CNN trained on simulated data performed poorly when applied to in vivo data, however, reasonable prediction accuracies were achieved for the in vivo data when training was performed on the in vivo data.
Table 2. Accuracy of intra-block connectivity matrix prediction based on convolutional neural networks (CNNs).CNNs were trained using dMRI data simulated from two-and three-dimensional connectome phantoms as well as in vivo dMRI data.Training examples comprised dMRI data extracted from a particular block/sub-image (input space), paired with the block's connectivity matrix (output space).The ground truth intra-block connectivity matrix was known for the connectome phantoms, while conventional streamline tractography was used to map intra-block connectivity for the in vivo dMRI data.Each CNN was evaluated using unseen data and the performance measures shown represent averages across 100,000 test blocks.

Discussion
In this study, we developed and evaluated a new framework to facilitate the application of deep and machine learning techniques to the task of connectome mapping.We demonstrated that our framework not only enables tractable and scalable deep learning in connectomics, but can also be used as a simple augmentation to enhance the accuracy of conventional connectome mapping pipelines that use streamline tractography.Based on evaluations performed with two-and threedimensional phantoms, we demonstrated that incorporating our framework into typical connectome mapping pipelines can improve the accuracy of reconstructed connectivity matrices by up to 20-30%, albeit with some additional computational burden.This improvement is owing to the robustness to noise and partial volume effects that is achieved by decomposing the dMRI data into overlapping blocks (sub-images), thereby introducing redundancy in the dMRI data among neighboring blocks.
Training a deep learning architecture to learn the mapping from dMRI data to a whole-brain connectivity matrix representation of the human connectome is currently intractable due to the high dimensionality of the data produced by dMRI.To establish tractability, the principle of our BDS framework is to spatially decompose the dMRI data into blocks.We demonstrated that the mapping from the dMRI data encapsulated by a block, to the block's connectivity architecture, can then be tractably learned under this decomposition approximation.In this study, we trained a CNN to learn this mapping, where the training data was derived from conventional streamline tractography.While this enabled a proof-of-principle demonstration of how BDS can facilitate mapping of the human connectome with deep learning techniques, future work should focus on developing training data that accurately recapitulates the complex geometry of white matter fibers.This would avoid the need for the CNN to learn intra-block connectivity from conventional streamline tractography.
Several recent studies have developed deep learning techniques to improve the performance of tractography.For example, the utility of deep learning has been demonstrated in learning the mapping between dMRI signals and local fiber orientations (Koppers et al., 2017), predicting the fiber orientation distributions (Benou and Riklin-Raviv, 2018) and directions (unit vectors) (Poulin et al., 2017, Wegmayr et al., 2018) from dMRI data.These tractography methods involve training a neural network to successively predict the direction in which to propagate streamlines based on sequences of dMRI signals.These methods estimate fiber orientations within each voxel, whereas in our BDS framework, the orientation used to propagate is collectively determined by a block.
Initially, we used connectome phantoms with known ground-truth connectivity to evaluate the performance of BDS.We used simulated dMRI data derived from two-and three-dimensional phantoms to determine the impact of variations in the stride length, seeding methodology and angle threshold parameters inherent to BDS.Variation in stride length, from 1 to 4, was used to study the importance of redundant voxels among consecutive blocks for generation of streamlines.The Fmeasure was primarily considered for evaluation of the numerical connectomes because it was not dominated by the high number of true negatives.The F-measure was maximized for a stride length of 1 because this stride resulted in the maximum number of shared voxels among consecutive blocks, which can aid in reconstructing complex fiber geometries and crossing fibers.Increasing the stride length resulted in reduction of the extent of overlap and redundancy between neighboring blocks, which yielded a greater number of false-negative connections and a reduction in F-measure.
We evaluated two seeding methodologies, seeding from the blocks that coincided with the underlying cortical parcellation (ROI-based) and seeding from all possible blocks (whole-image seeding).We found that both seeding methodologies yielded comparable accuracy (Figure 4(b)).However, wholeimage seeding yielded connectomes with a greater number of true-positive and false-positive connections, compared to ROI-based seeding, and thus the F-measure for both methods was comparable.ROI-based seeding is preferred over whole block-image seeding, as it is computationally less expensive.
We finally evaluated the impact of variation in angle thresholds.We found that smaller values of the angle threshold θ = 10 K ) resulted in failure to reconstruct many connections (false negatives), whereas a high threshold (θ > 30 K ) yielded a lower F-measure (Figure 4(c) and Supplementary Figure S2-S3).Taking both two-and three-dimensional phantoms into consideration, we recommended an angle threshold between 15 K -30 K .
BDS was found to improve the accuracy of connectomes mapped with deterministic and probabilistic streamline tractography.In particular, both the deterministic and probabilistic variants of BDS (Det and Prob BDS) outperformed Det CSD, which was previously shown to yield the most accurate reconstructions relative to a variety of other tractography methods (Sarwar et al., 2018).We found that many streamlines terminated before reaching a node, which led to low TPR and FPR (< 60%specifically for Det CSD).This behavior was observed in previous studies as well (Côté et al., 2013a, Sarwar et al., 2018, Girard et al., 2014).While increasing the angle threshold led to an increase in the total number of connections detected with Det and Prob CSD (Supplementary Figure S12), according to the F-measure, the new BDS framework consistently outperformed conventional tractography.
It is important to remark that conclusions derived from phantoms do not necessarily generalize to in vivo data.To address this limitation, we also assessed performance based on the extent of coupling between structural and functional connectivity, using in vivo functional and diffusion MRI data.For the majority of connection densities investigated, the strength of structure-function coupling was greater for structural connectomes mapped with BDS, compared to Det CSD (Figure 8).Interestingly, we found that Det CSD was associated with higher structure-function coupling than BDS variants for low connection densities (< 15%).This suggests that conventional methods reconstruct the strongest and most reliable connections with higher accuracy than BDS, but the reconstruction of weaker connections (i.e.connections that only appear at high connection densities) is better performed with BDS.In other words, BDS is particularly suited to the reconstruction of weak connections that are most susceptible to the effects of noise, and thus benefit the most from the redundancy achieved through the block decomposition process.
With conventional tractography, spurious deviations in streamline trajectories are due to inaccurately estimated fiber orientations within one or more voxels traversed by the streamline, leading to the reconstruction of false-positive and/or false-negative connections.Under our BDS framework, the orientation used to propagate and steer a block chain is collectively determined by a block ( voxels), unlike conventional tractography, where a single voxel determines the next propagation step.
Using a block of many voxels rather than a single voxel introduces a degree of redundancy in block chain propagation, which endows the propagation process with robustness to noise and partial volume effects.This robustness comes at the cost of discretization in the block chain trajectories, resulting in somewhat jagged trajectories and reduced spatial specificity.However, in connectomics, the accuracy of connectivity matrix estimation is typically a more important consideration than smooth reconstruction of fiber bundle geometry.By forgoing the spatial resolution of fiber trajectory localization caused by this discretization error, we can yield connectivity matrices that are substantially more accurate than established connectome mapping methods.
Several limitations of our proposed BDS framework require consideration, together with general challenges of mapping connectomes with deep learning.First, BDS is substantially more computationally expensive than conventional connectome mapping pipelines and yields inferior reconstruction of fiber bundle trajectories, due to discretization error.Some applications of tractography mandate precise localization of fiber bundle geometry, length and curvature, which are properties that can be straightforwardly estimated using conventional tractography algorithms (Poupon et al., 2008, Fillard et al., 2011, Neher et al., 2015, Tournier et al., 2008, Neher et al., 2014, Côté et al., 2013a).In contrast, BDS is not well suited to the estimation of these attributes.Second, blocks can comprise a heterogeneous composition of voxels spanning multiple tissue classes, including gray and white matter as well as cerebrospinal fluid (CSF).While this can potentially exacerbate partial volume effects, in principle, a CNN should be able to learn to distinguish the connectivity architecture of different tissue classes within a given block, including the absence of connectivity in CSF.Indeed, empty intra-block connectivity matrices were predicted for blocks positioned within CSF and connection density was typically reduced for blocks comprising a mix of white and gray matter, compared to blocks positioned entirely within white matter.Therefore, we can infer that dense intra-block connectivity matrices indicate white matter, whereas sparse or empty matrices indicate the presence of CSF and/or gray matter.Third, in the absence of a ground truth connectivity matrix for the in vivo dMRI data, we evaluated performance based on the extent of structure-function coupling.It is important to remark that structure-function coupling is not an explicit measure of connectivity matrix accuracy and may be influenced by other methodological factors.
Fourth, our findings point to the pressing need for the development of realistic connectome phantoms that can be used to generate training data with known ground truth connectivity.We found that training data comprising simulated dMRI data derived from existing two-and three-dimensional phantoms yielded CNNs that typically performed poorly when applied to in vivo data.Therefore, in this study, we used training data derived from in vivo dMRI data, where the ground truth intra-block connectivity was approximated with conventional streamline tractography.While this was sufficient to demonstrate the feasibility of connectome mapping with deep learning and our BDS framework, it remained limited by the shortcomings of streamline tractography, since the CNN most likely learned these shortcomings.Fifth, in this study, we only investigated performance using high-quality HCP data.The performance that we report here might be unattainable for low-quality data (low SNR), fewer diffusion-weighted gradients and/or fewer b-values.Finally, the proposed BDS algorithm reconstructs block-chains in a locally greedy manner using step-by-step block-chain propagation, analogous to conventional local tractography algorithms.To overcome the limitations inherent to this locally greedy paradigm (Zalesky, 2008), future work should focus on reconstructing block chains in a globally optimized manner based on global fits of the reconstructed block-chain trajectory to the dMRI signal.Global tractography methods (Mangin et al., 2013, Jbabdi et al., 2007, Sherbondy et al., 2009, Pestilli et al., 2014) have proven useful in alleviating the limitations inherent to locally greedy algorithms and we envisage that this would also be the case for BDS.
In conclusion, we developed and evaluated a new framework to enable tractable connectome mapping with deep learning.The task of connectome mapping is challenging and limited by numerous shortcomings of current diffusion models and tractography algorithms.The advantage of deep learning is that the mapping from an individual's dMRI data to their structural connectivity matrix can be learned, in the absence of any explicit diffusion model.We also demonstrated that our framework not only facilitates deep learning in connectome mapping applications, but can also be used to improve the accuracy of existing tractography-based connectome mapping pipelines.We consider this to be an immediate benefit of our framework, which can be straightforwardly incorporated into existing connectome mapping pipelines.In particular, we showed that block decomposition can improve the accuracy of connectivity matrix reconstruction by up to 30%, relative to the best performing streamline tractography algorithm.While performance gains demonstrated with connectome phantoms do not necessarily hold for in vivo data, we showed that our framework can improve the strength of structure-function coupling in the human connectome, which is a proxy measure of performance.
Code that implements our BDS framework is freely available at https://github.com/sarwart/BDS/

Figure 1 .
Figure 1.Block decomposition and stitching (BDS) framework for tractable connectome mapping with deep learning (lower arc), compared with conventional connectome mapping based on whole-brain tractography (upper arc).Deep learning is used to learn the mapping between dMRI data and whole-brain connectivity, in the absence of any intermediate steps involving streamline propagation.To ensure tractability, the dMRI images are first decomposed into spatially overlapping sub-images, referred to as blocks.A convolutional neural network (CNN) is trained to learn the mapping between these subimages and their local (internal) connectivity.A block stitching algorithm is then applied to stitch neighboring blocks together and yield a whole-brain connectivity matrix.Each block's local connectivity matrix can be determined by performing tractography locally within the volume encapsulated by the block or by training a CNN to learn this mapping.Connectivity matrix reproduced from (van den Heuvel and Sporns, 2013a) Figure 2(b) shows an example of block decomposition in the case of = 6 for= 4 .Images corresponding to different diffusion-weighted gradients and strengths are decomposed using the same stride and block dimension.Henceforth, we use to denote the index of the th block and to denote stride length.

Figure 2 .
Figure2.Intra-block connectivity mapping and block decomposition and stitching (BDS) demonstrated on a twodimensional dMRI phantom.(a) Schematic of intra-block connectivity mapping.Streamlines were mapped for a twodimensional dMRI numerical phantom using tractography.The image was then decomposed into overlapping 6 × 6 blocks (blocks appear as squares in two dimensions), where the blue square indicates the boundaries of one such representative block.A 20 × 20 connectivity matrix was then mapped for this representative block, where each row/column indexes one of the block's peripheral voxels, labelled from 1 to 20.The representative block encapsulates a portion of a single fiber bundle that comprises streamlines intersecting peripheral voxels 3, 13 and 14.Accordingly, elements (3,13), (3,14) and (13,14) of the block's connectivity matrix (rightmost) were set to one.(b) Schematic of the block decomposition process.The 14 × 14 image (left) is decomposed into 9 blocks (right) of size 6 × 6 using a stride of 4. This stride length results in an overlap of 12 voxels between consecutive blocks.Blocks are indexed using row ( ) and column ( ) subscripts.Colored boxes represent the mapping of voxels to blocks and their locations.(c) Schematic of the block stitching process used to estimate connectivity.The block with index is chosen as the seed block (red) and the connection between peripheral voxels 4 and 14 (yellow) within this block is selected as the seed connection, where denotes the connection's index.The unit vector, denoted by , determined for this seed connection points in the direction of the next-selected block (blue), denoted .For Deterministic-BDS, unit vectors for all connections comprising the candidate block (blue) are estimated, and the connection with unit vector that forms the most acute angle with the preceding block's unit vector is selected (green line) (d) Schematic of iterative block stitching resulting in a block chain.In this simplified example, each block comprises a single connection (i.e.= 1).Block is determined based on the unit vector within block .A block chain is an ordered set of blocks ( , , , … ) that defines and connection between gray matter regions.
Figure 4(b) shows ROC curves and the F-measure as a function of the block chain threshold for both seeding methodologies.While whole-image seeding yielded greater numbers of true and false positive connections, compared to ROI-based seeding, consideration of the F-measure indicates a negligible difference between the two methodologies.Therefore, given that both methodologies performed equally well with respect to the accuracy of connectome reconstruction, we recommend using ROI-based seeding because it is less computationally expensive than whole-image seeding.The results for Prob BDS are presented in Supplementary Figure S1(b), which provides further evidence indicating that ROI-based seeding outperforms whole-image seeding.3.1.3.Variation in angular threshold The angle threshold θ was varied from 10 o to 50 o in increments of 5 o .Figure 4(c) shows ROC curves and the F-measure as a function of the block chain threshold for each of these angle thresholds.The most stringent angle threshold (10 o ) yielded poor sensitivity and a relatively low F-measure, compared to more lenient thresholds.Of all the increments tested, we found that an angle threshold of 20 o yielded the highest F-measure.Thresholds exceeding 45 o led to a substantial increase in falsepositive connections, thereby compromising the specificity of connectome reconstruction.Based on the F-measure, in the absence of any thresholding, we found that Det BDS outperformed Prob BDS because the additional true-positive connections identified with Prob BDS were overshadowed by substantial increases in false-positive connections (Figure 5 and Supplementary Figure S2(b)).The variation in angle threshold θ had a consistent effect on both seeding methodologies for both algorithms (Det and Prob BDS; Supplementary Figure S2(a) and S2(c)).

Figure 3 .
Figure 3. Representative block chains reconstructed for a two-dimensional connectome phantom using stride lengths of = 3 (a) = 4 (b) = 5 (c).The non-diffusion weighted image for the connectome phantom is shown at the bottom-left of Panel (a).Blocks are indicated with squares boxes and the colored lines within each block denote intra-block connections inferred from conventional streamline tractography.Block dimensions are 6 × 6, and thus each block comprises 36 voxels.Although blocks spatially overlap in practice, the extent to which is determined by stride length, for visualization purposes, block are shown consecutively without any overlap.Shorter stride lengths increase the redundancy between neighboring blocks.This redundancy improves robustness to noise effects and can be seen to improve the resolution of fiber bundle reconstruction.

Figure 4 .
Figure 4. Evaluation of stride length (a), seeding methodology (b) and angle threshold (c) on the accuracy of connectivity matrix reconstruction with the block decomposition and stitching (BDS) framework.Connectivity matrices were mapped for two-dimensional connectome phantoms with known ground truth connectivity.Intra-block connectivity was mapped using conventional streamline tractography (Det CSD) and blocks were stitched using the deterministic variant of BDS (Det BDS).Accuracy was quantified using receiver operating characteristic (ROC) curves (top row) and the F-measure (bottom row), where curves are parameterized by the block chain threshold.For each ROC curve, the rightmost point on each curve corresponds to a threshold of zero, while the most severe threshold corresponds to the leftmost point.Unless indicated otherwise, a stride length of = 1, ROI-based seeding methodology and an angle threshold of ' = 20 K was used.Inset (bottom row) characterizes performance for block-chain thresholds O0.01In the absence of any thresholding, we found that Det BSD yielded the most accurate connectome

Figure 5 .
Figure 5. Accuracy of connectome reconstruction compared between the block decomposition and stitching (BDS) framework and conventional streamline tractography methods.Comparisons were performed using two-(a) and threedimensional (b) connectome phantoms with known ground truth connectivity matrices.The accuracy of connectome reconstruction was assessed using receiver operating characteristic (ROC) curves (left) as well as the F-measure evaluated as a function of the streamline/block-chain threshold (right).In the absence of any thresholding, the deterministic variant of BDS (Det BDS) yielded the most accurate connectome reconstruction, as determined by the F-measure.BDS methods: Deterministic BDS (Det BDS, magenta), Probabilistic BDS (Prob BDS, yellow).Streamline methods: Diffusion tensorbased deterministic (Det Tensor, black) and probabilistic (Prob Tensor, green) tractography, Constrained spherical deconvolution-based deterministic (Det CSD, blue) and probabilistic (Prob CSD, red) tractography.Inset (right) characterizes performance for streamline/block-chain thresholds O0.01

Figure 6 .
Figure 6.Qualitative comparison of whole-brain tractograms (a) as well as two selected fiber bundles (b,c) reconstructed for a representative individual with conventional streamline tractography (upper panel) and the block decomposition and stitching (BDS) framework (lower panel).Streamlines were generated using deterministic tractography and constrained spherical deconvolution (Det CSD), whereas block chains were generated using the deterministic variant of BDS.The fiber bundles shown correspond to the cortico-spinal tract (a) and a portion of commissural fibers traversing the corpus callosum (b).Many block-chains comprised of the same number and location of blocks, which are represented here by tubes with a larger diameter.Images were generated with DSI-studio (http://dsi-studio.labsolver.org/).

Figure 8 .
Figure 8. Assessment of structure-function coupling for structural connectomes mapped with the block decomposition and stitching (BDS) framework (a,b) and BDS with intra-block connectivity mapped with a convolutional neural network (CNN-BDS; c-d).Spearman's rank-order correlation coefficient was computed to assess the extent of structure-functional coupling for each structural connectivity matrix.Functional connectivity matrices were mapped using resting-state functional MRI data.The correlation coefficient was computed across all pairs of regions, excluding subcortical regions (a,c) and only regions comprising the right hemisphere (b,d).The correlation coefficient is shown as a function of the connection density of the structural connectivity matrix, as adjusted with a streamline or block chain threshold.See Appendix for a detailed description of each BDS variant listed in the legends.

Figure 9 .
Figure 9. F-measure quantifying the accuracy of connectome reconstruction compared between the block decomposition and stitching (BDS) framework and conventional streamline tractography methods (a) deterministic methods (b) probabilistic methods.All results pertain to two-dimensional connectome phantoms.

Figure 7
Figure 7(d-f) shows group-averaged structural connectivity matrices reconstructed with the CNN, as ) shows the propagation of a block chain according to this rule in which the new block is determined by the direction of the connection used in the current block.When arriving at a new block , we need to select a connection (i.e.new direction) from the set of candidate connections in the new block # $ % ,…,& .We propose two alternative criteria to select a connection within each new block.The first criterion is based on deterministic selection of a connection, while the second is a , calculated as (p i -p j ) / | (p i -p j ) |, where p i and p j is the spatial location of node i and j respectively (Figure2(c)).For each connection in !, an attempt is made to propagate a block

Table 1 .
Single-shell dMRI Acquisition and Signal Parameters