Sensory‐motor network topology in multiple sclerosis: Structural connectivity analysis accounting for intrinsic density discrepancy

Abstract Graph theory and network modelling have been previously applied to characterize motor network structural topology in multiple sclerosis (MS). However, between‐group differences disclosed by graph analysis might be primarily driven by discrepancy in density, which is likely to be reduced in pathologic conditions as a consequence of macroscopic damage and fibre loss that may result in less streamlines properly traced. In this work, we employed the convex optimization modelling for microstructure informed tractography (COMMIT) framework, which, given a tractogram, estimates the actual contribution (or weight) of each streamline in order to optimally explain the diffusion magnetic resonance imaging signal, filtering out those that are implausible or not necessary. Then, we analysed the topology of this ‘COMMIT‐weighted sensory‐motor network’ in MS accounting for network density. By comparing with standard connectivity analysis, we also tested if abnormalities in network topology are still identifiable when focusing on more ‘quantitative’ network properties. We found that topology differences identified with standard tractography in MS seem to be mainly driven by density, which, in turn, is strongly influenced by the presence of lesions. We were able to identify a significant difference in density but also in network global and local properties when accounting for density discrepancy. Therefore, we believe that COMMIT may help characterize the structural organization in pathological conditions, allowing a fair comparison of connectomes which considers discrepancies in network density. Moreover, discrepancy‐corrected network properties are clinically meaningful and may help guide prognosis assessment and treatment choice.


| INTRODUCTION
Graph theory and network modelling have been applied to characterize structural motor network topology in multiple sclerosis (MS), demonstrating a reduced motor network efficiency through the quantification of structural damage in white matter (WM) bundles connecting pairs of cortical and subcortical grey matter (GM) regions (Pardini et al., 2015). More broadly, graph analysis of the structural connectome (Sporns, Tononi, & Kötter, 2005) (i.e., the set of whitematter pathways between pairs of GM regions) has been successfully used to discriminate MS patients from healthy controls (HCs) and to classify MS clinical phenotypes (Kocevar et al., 2016;Li et al., 2013;Llufriu et al., 2017;Nigro et al., 2015). However, such between-group differences may be primarily driven by discrepancy in network density (van Wijk, Stam, & Daffertshofer, 2010), which is likely to be reduced in pathologic conditions as a consequence of macroscopic damage and fibres loss. Thus, resulting in a less accurate tracking of streamlines (Ozturk et al., 2010). In the framework of graph analysis, methods such as the minimum spanning tree have been applied to account for differences in density, by reducing networks to a backbone structure insensitive to alterations in connection strength or linked density (Tewarie, van Dellen, Hillebrand, & Stam, 2015). An alternative and indirect way to deal with group differences in density is to extract connectivity metrics from an atlas of bundles built from healthy subjects keeping network density constant (Pagani et al., 2019). Tracing fibres in HC offers the additional advantage to avoid inaccuracy in tract reconstruction related to the presence of WM lesions. Therefore, in MS studies, tractography is often performed in the control group (or a subset of it), and the reconstructed tracts are subsequently registered to patients' data to derive the metrics of interest (Pagani et al., 2019;Pardini et al., 2015;Steenwijk et al., 2015). Although the underlying idea is the same, its implementation is slightly different in each of these works. Pagani et al. (2019) first coregistered the diffusion tensor images of HCs to the standard Montreal Neurological Institute (MNI) space, then they used the average of those data to perform tractography saving only the tracts connecting pairs of cortical areas with more than five streamlines as voxel maps.
Finally, they registered all the remaining subjects to MNI space and they used the common tractogram to compute the individual connectomes. Pardini et al. (2015) instead performed tractography in each individual healthy subject's space and then registered the recovered track density images to the MNI space to create populationaveraged maps for each tracts of interest. They then coregistered these maps to each subject involved in the study to compute the connectomes. Finally Steenwijk et al. (Steenwijk et al., 2015) implemented a similar method of Pardini et al. (2015), but they computed for each subject and tract separately the average of weighted lesion volume and weighted average of fractional anisotropy (FA) in normal appearing WM. When tractography is conducted directly in MS patients, an FA threshold is set during fibre reconstruction and a minimum number of fibres are selected to define single bundles in order to reduce the risk of false-positive connections (Nigro et al., 2015;Shu et al., 2011). The shortcoming of this approach is the drastic reduction in reconstructed fibres, especially in those bundles that are rich in crossing fibres (Sinke et al., 2018). More recently, a sphericaldeconvolution-informed filtering of tractograms (SIFT; R. E. Smith, Tournier, Calamante, & Connelly, 2013) has been employed to reduce reconstruction bias and improve biological plausibility (Koubiyr et al., 2019), but the accuracy of SIFT application to pathological brains is still under debate (Zalesky, Sarwar, & Ramamohanarao, 2020).
Furthermore, the characterization of the structural connectome in MS has to take into account the impact of WM lesions on connectivity which is usually assessed through correlation analysis between graph metrics and lesion loads (He et al., 2009;Romascano et al., 2015). A more specific disconnection analysis can also be conducted, quantifying dedicated graph measures that estimate the indirect, compensatory connections between two regions developed after the transection of the direct connection between them (Li et al., 2013). More recently, the impact of macroscopic lesions on structural connectivity was modelled by assuming transection of all fibres passing through WM lesions (Pagani et al., 2019).

Finally, the quantification of the connection strength in structural
connectomes is an open issue. Typically, the connection strength between each pair of grey-matter regions is 'quantified' by counting the number of streamlines connecting them, that is, streamline count, but this approach is not quantitative (Jones, Knösche, & Turner, 2013). Microstructure-informed tractography (Daducci, Dal Palu, Descoteaux, & Thiran, 2016) was recently proposed as a means to improve the estimation of structural connectivity by combining tractography with local microstructural features of the tissue and fitting the actual contributions of the streamlines to the measured diffusion magnetic resonance imaging (MRI) data. These contributions do not allow to estimate the microscopical fibre count, but this approach has the potential to provide a more 'physically quantitative' assessment of the connectivity than the simple streamline count. In fact, as the contributions of the streamlines (or weights) are estimated such that they explain the diffusion MRI data, and the connectivity is 'physically quantified' based on these weights. This possibility to extract more 'quantitative' metrics from the reconstructed connectomes may allow for a fair comparison of network properties despite density discrepancies. However, to the best of our knowledge, this approach has never been proposed in clinical studies.
In this proof of concept study, we investigated the topology of the 'physically quantitative' sensory-motor network (SMN) (i.e., the network whose weights are estimated through microstructure informed tractography) in MS using the convex optimization modelling for microstructure-informed tractography (COMMIT) (Daducci, Dal Palu, Lemkaddem, & Thiran, 2013

| Diffusion MRI processing
Diffusion MR images were corrected for motion and eddy currents (Andersson & Sotiropoulos, 2016) using FSL. To perform whole brain anatomically constrained tractography (R. E. Smith, Tournier, Calamante, & Connelly, 2012), we first coregistered the T1 and diffusion images using FMRIB's linear image registration tool (FLIRT) (Jenkinson, Bannister, Brady, & Smith, 2002) of FSL with boundarybased cost function (Greve & Fischl, 2009). Then we computed the fibre orientation distribution functions using the multishell multitissue-constrained spherical deconvolution approach (Jeurissen, Tournier, Dhollander, Connelly, & Sijbers, 2014;Tournier, Calamante, & Connelly, 2007) and generated 1 million streamlines using the iFOD2 (Tournier, Calamante, & Connelly, 2010) tractography algorithm implemented in MRtrix (http://www.mrtrix.org). In light of the discussion in Zalesky et al (2020), we processed the resulting tractograms using the COMMIT (Daducci et al., 2013(Daducci et al., , 2015 with stick and zeppelin ball model (Alexander et al., 2010). COMMIT is a powerful framework that allows to decompose a signal in contributions coming from different compartments. The main assumption of the framework is that the contribution of a streamline is constant along its path, while the remaining components can be different in each voxel. In this case, we applied COMMIT to diffusion MR signal and we decomposed the signal in intra-axonal, extra-axonal and isotropic contributions according to the stick and zeppelin ball model (Alexander et al., 2010). Indeed, with this model, we imposed that the intraaxonal diffusion signal was constant along each tract and (when needed) we indirectly accounted for the presence of free water due to a lesion with the zeppelin and ball compartments.
Finally, for each subject, both the raw (i.e., obtained using the number of streamlines as entries) and the COMMIT-weighted connectomes (i.e., obtained using COMMIT weights as entries) were built using the motor network parcellation described above ( Figure 2). As entries (a ij ) of COMMIT-derived matrices, we used the weighted average intra-axonal signal contribution of each bundle: where i and j are the indices of ROIs connected by the bundle, N ij is bundle's number of streamlines, x k ij is the weight of the streamline k obtained by COMMIT and l k its length. In this way, each entry contained the total signal fraction associated to the bundle, which was given by the weighted average of the streamline contribution (obtained with COMMIT) multiplied by its length and divided by the average length of the bundle.
In light of the recent results showed in Buchanan et al. (2020), in Data S1, we also report additional results obtained by thresholding the number of streamlines in the raw connectomes according to two widely used techniques: proportional and consistency thresholding.
For further details, we recommend readers to refer Data S1.

| Graph analysis
As it was done in previous works (Pagani et al., 2019;Pardini et al., 2015;Steenwijk et al., 2015), for each subject we computed six global network measures from the obtained connectomes using the brain connectivity toolbox (Rubinov & Sporns, 2010): modularity and density (corresponding to the fraction of present connections to possible connections). For each node of the subjects' connectome we also computed local efficiency and nodal strength to investigate which node of the SMN was more affected by the disease.

| Statistical analysis
All analyses were performed using Statistical Package for Social Science (SPSS V.25.0).
Between-group comparisons were performed via analysis of covariance analysis, entering age and gender as covariates. In order to assess differences in density estimation related to the application of COMMIT, we performed between-group comparisons both on results from raw connectomes and COMMIT-weighted connectomes and repeated the analysis entering density as additional covariate.
The relationship between network global properties, T2 lesion load and GM atrophy were tested via partial correlation accounting for age and gender.
F I G U R E 1 Motor network hubs used in our analysis in a representative healthy subject. The primary sensory-motor cortex (S-M1) is shown in red; the secondary motor cortex (M2) in green; the secondary sensory cortex (S2) in light blue; the posterior associative sensory cortex (AS Sens C) in yellow; the prefrontal cortex (PFC) in blue; the deep grey matter (Deep GM) in pink (for the right hemisphere) and orange (for the left hemisphere) and the cerebellum in purple The relationship between network properties and clinical disability was tested with stepwise regression models, entering age and gender in the first block and network global/local properties in the second block.
Results were considered significant for p < .05 (Bonferroni   Table 1. After Bonferroni correction for multiple comparisons, F I G U R E 2 Matrix representation of the connectomes obtained with the two different methods: counting the number of streamlines connecting two pairs of grey matter regions (top); or assigning the quantitative measures obtained with COMMIT (bottom). For both method we report the average connectomes obtained for the two groups of subjects: healthy controls (left) and PMS patients (right). In both cases (raw and COMMIT), the pattern of connections is similar, but while in the upper case the information contained in the connectomes is nonquantitative, in the bottom ones it represents the intra-axonal signal fraction associated to each connection. We also observe that some interhemispheric connections present in the raw connectomes disappear after the application of COMMIT. COMMIT, convex optimization modelling for microstructure informed tractography; PMS, progressive multiple sclerosis modularity, global efficiency and mean strength were significantly different between the two groups of subjects when accounting for age and sex. When controlling also for density only the difference in modularity was still present. Of note, no significant differences in density were identified between the two groups.
Mean values and SDs of the local network metrics are reported in Table 2 (strength) and Table 3 (efficiency). After Bonferroni correction for multiple comparisons, significant differences were identified in five nodes in terms of strength and in nine nodes in terms of efficiency between the two groups of subjects when accounting for age and sex. When controlling also for density, significant differences were still identified in two nodes in terms of strength and in six nodes in terms of efficiency.

| COMMIT-weighted connectomes
Mean values and SDs of the global network metrics are reported in Table 4. After Bonferroni correction for multiple comparisons, all the explored metrics, except the clustering coefficient, were significantly different between the two groups of subjects when controlling for age and sex. When controlling also for density the difference in assortativity disappeared.
Mean values and SDs of the local network metrics are reported in Table 5 (strength) and Table 6 (efficiency). After Bonferroni correction for multiple comparisons, significant differences were identified in six nodes in terms of strength and in seven nodes in terms of efficiency between the two groups of subjects when accounting for age and sex. When controlling also for density significant differences were still identified in the same nodes in terms of strength and in seven nodes in terms of efficiency.

| Relationship between raw connectome global properties, WM lesions and GM atrophy
Accounting for age and gender, significant correlations were identified between T2 lesion volume and global efficiency (r = −.655, p < .0001),

| Clinical impact of COMMIT-weighted connectome abnormalities
The models including demographic variables and network global properties accounted for 27% to 35% of variance in 9HPT scores (for

| DISCUSSION
Notwithstanding all previous efforts in investigating structural connectivity and disconnection in MS, in this study we propose a methodological approach-COMMIT-that accounts for the presence of lesions and fibres loss and provides a means to directly compare connectomes with different density. Note: All values are expressed as mean ± SD; ANCOVA age and gender corrected (p a ), ANCOVA age, gender and density corrected (p b ). Statistically significant p values after Bonferroni correction are highlighted in bold.
Thanks to its capability of decomposing the intrinsic signal contribution of each streamline in the tractogram, COMMIT may represent an effective method to cope with density discrepancies between healthy subjects and patients. The main idea behind this method is to assume that one (or more) microstructure feature does not vary along the length of a tract and therefore it is possible to effectively estimate its value for the entire tract (rather than only voxel-wise). This estimation is done simultaneously for all the streamlines by fitting them to a map related to the selected microstructure feature. If only diffusionweighted magnetic resonance imaging data are available, it is reasonable to assume that the intra-axonal diffusion signal is constant along the tract and COMMIT uses any predefined microstructural model to estimate it. Similarly to what was recently found in Lipp et al. (2019), using the recently introduced multishell multitissue spherical deconvolution (Jeurissen et al., 2014) and the probabilistic algorithm (Tournier et al., 2010) to generate streamlines, we were able to propagate the tracking also inside MS lesions to build the input tractograms.
We then applied COMMIT to decide if a streamline passing through a lesion is essential to explain the signal or not and consequently keeps or discards it to construct the final tractogram. In the present work, we employed as microstructural model the stick and zeppelin ball model (Alexander et al., 2010) which indirectly accounts for the presence of free water due to a lesion with the zeppelin and ball compartments. Finally, to construct the COMMIT-weighted connectomes, we chose not to use the traditional number of streamlines connecting two cortical regions of interest (streamlines count), which was shown not to be quantitative (Jones et al., 2013). Conversely, we considered the more informative total signal fraction associated to the bundle, which is given by the weighted average of the streamline contribution (obtained with COMMIT) multiplied by its length and divided by the average length of the bundle. This approach offers two main advantages. First, by forcing fibre tracking within lesions and subsequently filtering them according to the signal preservation along the streamline, COMMIT retains in the tractogram only fibres whose microstructure is not irredeemably damaged by lesions or subtle inflammatory/ neurodegenerative processes ongoing in the normal appearing WM (Lassmann, 2018). Thus, producing a weighted network composed by 'healthy' and partly damaged fibres whose signal is not irreversibly compromised and can be fitted with a stick. As a consequence of COMMIT's filtering, in the COMMIT-weighted connectomes, we observed a reduction in density in comparison with the raw connectomes both in patients and controls (Figure 3). A number of implausible connections, related to tractography intrinsic limitations, as well as the fact that our control group presumably presented agerelated subtle WM abnormalities, were removed in HCs. As expected though, the number of implausible connections removed in patients was even higher, which explains why differences in terms of density between patients and controls became apparent only after COMMIT application. Second, by giving the possibility to compare more 'quantitative' metrics rather than measures derived from the nonquantitative streamline number (Jones et al., 2013), COMMIT offers the possibility to assess differences in network properties beyond changes driven by density discrepancy. This is supported by the results of our between-group comparison, which shows that, while topology differences identified with standard tractography were mainly driven by density, differences in global and local properties derived from the COMMITweighted connectomes were insensitive to density correction . Finally, it is worth highlighting that although COMMIT estimates the actual weight of the edges in the network by fitting the corresponding streamlines to the white-matter signal, normalization may still be required to account for ROI size differences in the chosen parcellation (Sotiropoulos & Zalesky, 2019). In fact, larger ROIs may be connected with more streamlines simply because of their size.
Note, however, that this applies to raw and COMMIT-weighted connectomes alike, and hence it does not bias our results. Future studies will investigate the possibility to use COMMIT to account also for this aspect.
Differences in connectome global properties estimated after COMMIT application suggest that also the COMMIT-weighted connectome presents the topology abnormalities previously described in MS (Kocevar et al., 2016;Li et al., 2013;Llufriu et al., 2017;Nigro et al., 2015;Pardini et al., 2015). Indeed, the COMMIT-weighted SMN was less efficient, more dispersed and weaker in MS than in HC, supporting the notion that also seemingly intact connections are not sufficient to preserve brain structure. As COMMIT retains also con- COMMIT-weighted SMN global properties showed strong to moderate associations with WM lesion load and atrophy, confirming that brain topological organization is related to the accrual of macrostructural damage (Pagani et al., 2019), with lesion load playing a predominant role in PMS (Steenwijk et al., 2015).
Boxplots showing the differences in global network measures between HCs (white) and PMS patients (grey) for both raw and COMMIT tractograms. We observe that after the application of COMMIT the differences between HC and PMS patients are more pronounced. Also, the presence of outliers is often mitigated when COMMIT is applied. COMMIT, convex optimization modelling for microstructure informed tractography; HCs, healthy controls; PMS, progressive multiple sclerosis F I G U R E 4 Barplot showing the local efficiency of all the hubs of the motor network for both raw and COMMIT connectomes. The statistically significant differences between HCs in white and PMS patients in grey and accounting for discrepancies in age, sex and density are marked with an asterisk. COMMIT, convex optimization modelling for microstructure informed tractography; HCs, healthy controls; PMS, progressive multiple sclerosis for manual dexterity performance, highlighting the importance of motor planning for the execution of fine motor movements, while efficiency of associative sensory cortex was significantly correlated with the ambulation performance. Interestingly, it seems that the efficiency of integrative rather than primary areas is particularly relevant for clinical function preservation within the weighted connectome, highlighting the compensatory role of these regions in advanced disease stages.

| CONCLUSIONS
Topology differences identified with standard tractography in MS seem to be mainly driven by density, which, in turn, is strongly influenced by the presence of lesions, suggesting caution when interpreting between-group differences in connectome properties. Moving from a qualitative towards a more 'quantitative' appraisal of the brain structural connectome, COMMIT application allowed the identification of a significant difference in density between patients and HC and the exploration of network topology in the COMMIT-weighted connectome. Differences observed in network global and local properties suggest that preserved connections undergo a topological reorganization in MS. Within such reorganization of the brain connectome, decreased local efficiency in key areas of the SMN represent the most relevant correlates of motor disability. Based on these results, we believe that COMMIT may help characterize the topological organization of structural networks in pathological conditions, allowing a fair comparison of connectomes which takes into account discrepancies in network density. More importantly, our study shows that discrepancy-corrected network properties are clinically meaningful and, therefore, may help guide prognosis assessment and treatment choice.

ACKNOWLEDGMENTS
This study was supported in part by grants from NMSS (RG 5120A3/1) and Teva Neuroscience (CNS-2014-221). This work was supported by the Rita Levi Montalcini Programme of the Italian Ministry of Education, University and Research (MIUR).

DATA AVAILABILITY STATEMENT
The brain MRI images used were obtained from the Mount Sinai Hospital of New York and could available from the corresponding author upon reasonable request.
F I G U R E 5 Barplot showing the strength of all the nodes of the motor network for both raw and COMMIT tractograms. The statistically significant differences between HCs in white and PMS patients in grey and accounting for discrepancies in age, sex and density are marked with an asterisk. With the application of COMMIT differences in the left associative sensory cortex, sensory-motor and right deep grey matter strength appears. COMMIT, convex optimization modelling for microstructure informed tractography; HCs, healthy controls; PMS, progressive multiple sclerosis