Disordered network structure and function in dystonia: Pathological connectivity vs. adaptive responses

Primary dystonia is thought to emerge through abnormal functional relationships between basal ganglia and cerebellar motor circuits. These interactions may differ across disease subtypes and provide a novel biomarker for diagnosis and treatment. Using a network mapping algorithm based on resting-state functional MRI (rs-fMRI), a method that is readily implemented on conventional MRI scanners, we identified similar disease topographies in hereditary dystonia associated with the DYT1 or DYT6 mutations and in sporadic patients lacking these mutations. Both networks were characterized by contributions from the basal ganglia, cerebellum, thalamus, sensorimotor areas, as well as cortical association regions. Expression levels for the two networks were elevated in hereditary and sporadic dystonia, and in non-manifesting carriers of dystonia mutations. Nonetheless, the distribution of abnormal functional connections differed across groups, as did metrics of network organization and efficiency in key modules. Despite these differences, network expression correlated with dystonia motor ratings, significantly improving the accuracy of predictions based on thalamocortical tract integrity obtained with diffusion tensor MRI (DTI). Thus, in addition to providing unique information regarding the anatomy of abnormal brain circuits, rs-fMRI functional networks may provide a widely accessible method to help in the objective evaluation of new treatments for this disorder.


Introduction
Dystonia is a brain disorder characterized by sustained muscle contractions resulting in involuntary twisting movements and abnormal postures (1, 2). While frequently sporadic, primary dystonia has been identified with a variety of genetic mutations (3)(4)(5). Of these, the most common are the DYT1 (TOR1A) and DYT6 (THAP1) mutations, which are inherited as autosomal dominant traits with incomplete penetrance. Primary dystonia has often been viewed as a disorder of the basal ganglia, although considerable evidence points to widespread involvement of abnormal functional networks connecting these structures to the cerebellum, thalamus, and motor cortex (6)(7)(8)(9). Indeed, spatial covariance analysis of [ 18 F]fluorodeoxyglucose (FDG) PET scan data (10), which maps afferent synaptic activity in brain regions (11,12), has revealed a reproducible dystonia-related metabolic pattern involving the putamen, cerebellum, and supplementary motor area (SMA) in patients with inherited as well as sporadic forms of the disorder (13,14).
A limitation in the wide application of PET-based network mapping approaches has been restricted availability and patient concerns regarding radiation doses. Recently, we developed a non-invasive method to identify and validate disease networks in resting-state functional MRI (rs-fMRI) scan data in conjunction with independent component analysis (ICA) and bootstrapping resampling (15)(16)(17). This method is available on standard clinical MRI commonly used for neuroradiological exams. The ability of rs-fMRI to detect disease networks has been validated against spatial covariance analysis with FDG PET (16,17), including the ability to partition the corresponding graphs into core and periphery zones (18), and to map node-to-node connections within and between modules (19,20). Indeed, this approach has been used to examine the effects of genotype, disease progression, and treatment on functional connectivity in Vo et al. 4 the network space (19)(20)(21).
Here, we used rs-fMRI to characterize an abnormal hereditary dystonia-related pattern (H-DytRP) in clinically manifesting (MAN) DYT1 or DYT6 patients. Given that these genes are only 30-50% penetrant (3-5, 13, 22), we additionally measured H-DytRP expression in clinically non-manifesting (NM) mutation carriers, sporadic dystonia (SPOR) patients, and healthy control (HC) subjects. Likewise, we identified an analogous sporadic dystonia-related pattern (S-DytRP) in rs-fMRI scans from SPOR and used the results to evaluate similarities and differences in the topographies of the two networks and in their respective relationships to dystonia severity ratings. Lastly, we used a graph theory-based approach to determine whether the internal structure of the disease network differed for the various groups. In a recent study of Parkinson's disease PD networks, we found that pathological and adaptive connectivity responses were associated with specific changes in graph metrics within key modules (20). By applying the same strategy to dystonia, we discerned analogous connectivity patterns in patients and clinically unaffected gene carriers. This distinction provides insight into the developmental mechanism that underlies dystonia, as well as the design of new treatment strategies for this and related disorders.

Hereditary dystonia-related pattern (H-DytRP)
To identify a significant disease network associated with hereditary dystonia, we analyzed rs-fMRI scans from manifesting DYT1/DYT6 mutation carriers (MAN) and age-matched healthy control (HC1) subjects using an ICA-based algorithm (see Methods). Patients and control subjects were distinguished by expression levels for three of the 40 ICs that were identifiable in Vo et al. 5 the data (Fig. 1A), corresponding to the cerebellum (IC18), the basal ganglia and thalamus (IC9), and sensorimotor cortex and occipital association regions (IC10). The estimated coefficients (weights) on the individual components (Fig. 1B) were used to combine the individual topographies into a composite disease network, which we termed the hereditary dystonia-related pattern (H-DytRP). The salient regions that comprise this network ( Fig. 2A) are presented in Table 1. Expression levels for this network (Fig. 2B, left) were elevated in MAN patients (red bar) relative to the HC1 subjects (gray bar) used for network identification (p<0.001; permutation test, 5000 iterations). Prospectively computed network expression (Fig. 2B, right) was likewise increased in sporadic dystonia (SPOR) patients (black bar) relative to healthy control subjects (gray bar) (p=0.03; post-hoc test). H-DytRP expression was also increased in NM mutation carriers (blue bar) compared to healthy control values (p<0.001; post-hoc test).
Significant differences in network expression were not observed across the MAN, NM, and SPOR groups (p=0.52; one-way ANOVA).

Sporadic dystonia-related pattern (S-DytRP)
We used an analogous strategy to identify a sporadic dystonia-related pattern (S-DytRP) in rs-fMRI scans from SPOR patients and HC1 subjects, but a reliable frequency histogram could not be identified for IC selection. Nonetheless, the same three ICs (IC18, IC9 and IC10) were also detected by the algorithm when it was applied to the combined MAN and SPOR patient group.
Indeed, a significant S-DytRP was identified based on the same ICs as those in the H-DytRP derivation, although the weights on two of these differed for the two dystonia networks ( Fig. 1B and S3A). Specifically, the distribution of weights on IC18 (cerebellum) (left panels) was different for H-DytRP and S-DytRP (Jensen-Shannon divergence (JSD)=0.25), with a larger Vo et al. 6 coefficient (represented by the mean value of the distribution; vertical dashed lines) on the former network. The weight distributions on IC10 (SMC/OCC) (right panels) also differed (JSD=0.26), but with a larger coefficient on the latter network. By contrast, the weight distributions on IC9 (basal ganglia/thalamus) (middle panels) were similar for the two networks (JSD=0.02), with nearly identical coefficients. Thus, using the same ICs as H-DytRP but with revised coefficients, we identified a candidate S-DytRP that accurately discriminated SPOR from HC1 (p<0.005; permutation test, 5000 iterations). As expected, voxel weights on the two dystonia networks were strongly correlated (r=0.89, p<0.001; voxel-wise correlation, corrected for spatial autocorrelation), as were corresponding expression levels across the study population (r=0.91, p<0.0001; Fig. 2C). Moreover, as with H-DytRP, S-DytRP expression ( Fig. S3B) was elevated in MAN, NM, and SPOR compared to healthy control values (p<0.015, post-hoc tests).

Clinical correlates of network expression
H-DytRP expression in MAN and SPOR patients (n=20) correlated with Burke-Fahn-Marsden Dystonia Rating Scale (BFMDRS) motor ratings obtained at the time of imaging (R 2 =0.36, p=0.005; linear regression). Diffusion tensor imaging (DTI) measurements of fractional anisotropy (FA) in a prespecified subrolandic white matter (SWM) volume also correlated with motor ratings (R 2 =0.29, p=0.015) in this group. Indeed, we found that a 2-variable model based on H-DytRP expression and SWM FA in combination predicted motor ratings more accurately (R 2 =0.62, p=0.0003) and with greater generalizability than 1-variable models based on each predictor alone (AIC=109 for the 2-predictor model vs. 117 and 120 for 1-predictor models based on H-DytRP expression or SWM FA, respectively). These relationships are depicted graphically for each predictor (Fig. 3A, B), and for both together (Fig. 3C), as partial correlation Vo et al. 7 leverage plots (see Methods). Of note, similar results were also obtained using expression values for S-DytRP, rather than H-DytRP, as the network predictor (Fig. S3C).

Changes in functional connectivity
We next identified the functional connections gained in the network space for each dystonia group relative to the healthy subjects (Tables S2-S4). Given the similarity of the H-DytRP and S-DytRP topographies shown above, connectional analysis focused on the former network. To evaluate the organizational structure of the H-DytRP, we calculated the gain and loss of functional connections linking its component nodes.

Gain:
The gain in connections in the H-DytRP space (Fig. 4A, left panel) was similar for the MAN, NM, and SPOR groups (F(2,297)=0.34, p=0.71; one-way ANOVA). To delineate group differences in connectional gain within specific network subspaces, we partitioned the H-DytRP into core and periphery subgraphs using a modularity maximization algorithm ( Fig. 4B; see Methods). The core included the cerebellum, thalamus, putamen, and the left precentral and perirolandic areas, while the periphery was composed mainly of prefrontal, superior temporal, and parietal cortical regions. We found that the distribution of gained connections in the core and periphery (Fig. 4A, middle panel) differed for the MAN, NM, and SPOR groups (F(2,594)=194, p<0.0001; group × subgraph interaction effect). Indeed, group differences in connectional gain were more pronounced within the core compared to the periphery. In this subgraph (black outline), more connections were gained in MAN compared to either SPOR or NM (p<0.0001; post-hoc tests). We note that in MAN, NM, and SPOR, connectional gain was greater outside compared to inside the core (Fig. 4A, right panel). This tendency differed significantly across the groups (F(2,297)=137, p<0.0001; one-way ANOVA), with greater relative gain outside the Vo et al. 8 core in SPOR and NM compared to MAN (p<0.0001; post-hoc tests).
The individual connections gained in the MAN, NM, and SPOR groups are summarized in Tables S2A, S3A, and S4A. MAN and SPOR, the two clinically affected groups, displayed abnormal core-core connections linking the right thalamic and putamen, as well as coreperiphery connections linking the left middle frontal gyrus and precentral cortex, and the right Rolandic operculum and angular gyrus. (In the left hemisphere, abnormal connections linking the latter two regions were gained in MAN but not SPOR.) By contrast, MAN and NM, the two gene positive groups, exhibited abnormal gain in core-core connections linking the Rolandic operculum and the putamen bilaterally, and in core-periphery connections linking the left middle and inferior opercular frontal regions. Notably, gain in connections linking nodes in the cerebellar vermis and lateral hemisphere was present in all three groups. That said, gain in connections linking the thalamus to other H-DytRP core nodes, the putamen and Rolandic operculum in particular, was seen in MAN and SPOR patients, but not in NM gene carriers.
Loss: Loss of normal H-DytRP connections was also noted (Fig. 5A, left panel), which was most pronounced in NM carriers compared to MAN or SPOR patients (p<0.001; post-hoc tests). Connectional loss in the core and periphery of the network (Fig. 5A, middle panel) differed across the groups (F(2,594)=142, p<0.0001; group × subgraph interact effect). In contrast to connectional gain in the core (black outline), which was most pronounced in MAN, loss of connections in this module was greater in NM (p<0.001; post-hoc tests). Indeed, the NM group exhibited relatively greater connectional gain outside the core, but greater loss inside the core compared to the other groups (Fig. 5A, right panel).
The normal H-DytRP connections that were lost in the MAN, NM, and SPOR groups are summarized in Tables S2B, S3B, and S4B. In MAN, loss of core-core connections was limited Vo et al. 9 to those linking the left pallidum to the putamen. By contrast, loss of core-periphery connections involved those linking the left pons and cerebellar vermis, the vermis and right middle frontal gyrus, and the angular and lingual gyri, while cortico-cortical connections in the periphery were also lost in this group. In NM, by contrast, the loss of normal core-core connections was more extensive than the other groups. In particular, connections normally linking the left thalamus and putamen were not present in NM, as were those linking the vermis and right putamen, the left pallidum and thalamus, and the left motor cortex and globus pallidus. A variety of coreperiphery connections, including those linking the cerebellum with the cortical association regions, were also lost in NM, as well as cortico-cortical connections in the periphery. Finally in SPOR, connectional loss was similar to that observed in MAN, with limited involvement of core-core connections such as those linking the cerebellar vermis to the right putamen, and left motor cortical regions and the globus pallidus. More extensive loss was noted for core-periphery connections such as those linking the cerebellar vermis with the inferior frontal operculum, and the right thalamus with frontal, parietal, and occipital association regions, as well as for corticocortical connections within the periphery. In aggregate, connectional gain in the core was greatest in MAN and least in NM and SPOR. In the same subgraph, however, loss was greatest in NM but only minimal in MAN and SPOR.

Altered graph metrics in the H-DytRP core zone
To determine the effects of the connectivity changes on network function in dystonia, we quantified graph metrics in specific H-DytRP subgraphs. Given that group differences were greatest in the core, further analysis focused on this subgraph. The results for each of the metrics are summarized in Table 2. Core degree centrality (Fig. 6A), a measure of overall connectivity Vo et al. 10 in the subgraph, was greater in MAN compared to NM, SPOR, and HC1 (PCORR<0.001).
Likewise, assortativity measured in the same subgraph ( Fig. 6B) was increased in MAN compared to NM (PCORR<0.001) and HC1 (PCORR<0.05), but was reduced in NM compared to control (PCORR<0.05). Thus, while core degree centrality was similar for NM and HC1, assortativity was lower in the former group relative to the latter. Regarding the other core metrics, clustering (Fig. 6C) was reduced in MAN relative to HC1 and NM (PCORR<0.001), whereas characteristic path length (Fig. 6D) was increased in MAN relative to the other groups (PCORR<0.001). NM, by contrast, exhibited increased clustering (PCORR<0.001) without significant difference in characteristic path length in comparisons with HC1. These changes are consistent with significant reduction in small-worldness seen in the core (Fig. 6E) for MAN compared to NM and HC1 (PCORR<0.001) and the modest increase in the core (PCORR<0.05) that was present in NM relative to HC1. Graph analysis of the H-DytRP core additionally revealed the presence of a "rich club," i.e., a set of highly interconnected central nodes, within this subgraph, but only in the MAN group (Fig. 6F) and not in the NM, SPOR, or HC1 groups (PCORR<0.001). The rich-club coefficient (see Methods) was normal in the NM and SPOR groups.
Reconstructions of the H-DytRP core zone (see Methods) in MAN (Fig. 7A) revealed an abundance of connections linking high degree nodes (yellow-yellow) at the center. The gain in connections in this subgraph (red lines) was largely, but not exclusively, assortative as evidenced by those linking nodes of low degree (blue-blue) as well as high degree (yellow-yellow). In NM ( Fig. 7B), by contrast, several disassortative connections were observed in the core linking high and low degree nodes (yellow-blue); however, links between nodes of similar degree were also observed in this group. Graph visualization did not reveal prominent connectivity differences in Vo et al. 11 the core subgraph of SPOR patients (Fig. 7C); connectional gain was modest in this group as it was in NM. Subgraph reconstruction of the H-DytRP core in HC1 subjects is displayed for reference (Fig. 7D).

Discussion
In this study, we used a novel machine learning approach to identify dystonia-related functional networks in rs-fMRI data, characterized by contributions from the basal ganglia, cerebellum, thalamus and sensorimotor cortex, as well as frontal and parieto-occipital association regions, which is in agreement with previously reported network abnormalities in dystonia (23)(24)(25). The findings demonstrated striking similarity of disease networks derived from hereditary vs. sporadic dystonia patients, supporting a common nosology. Moreover, given the significant relationship between functional network expression in individual patients and dystonia motor ratings, the rs-fMRI-based disease topography may have utility as an objective imaging biomarker in clinical trials. This approach also allowed for a detailed analysis of the connections (edges) linking network regions (nodes) within and between key modules in the various groups.
In addition to identifying commonalities and differences between groups in individual connections, we found distinctive changes in network organization that were consistent with a neurodevelopmental mechanism for the disorder. The results also suggest the possibility of new treatment targets for patients with this disorder.

Gain and loss of functional connections in the H-DytRP network
As part of the study, we mapped individual functional connections that linked H-DytRP nodes in the MAN, NM, and SPOR groups but not in HC. Abnormal gain in connections linking the cerebellar vermis and hemispheres was observed in all three groups, supporting a major role for this structure in dystonia (7,9,13,26). Certain connections were gained in dystonia gene carriers irrespective of penetrance, i.e., in both MAN and NM, such as those linking the Rolandic operculum with the putamen, and the inferior frontal operculum with the middle frontal gyrus.
These abnormal connections were not detected in SPOR, suggesting that dystonic movements are mediated by other pathways in the latter group. That said, other functional connections were present in both MAN and SPOR but not in NM or HC, such as those linking the centromedian thalamus with the putamen, the Rolandic operculum with the angular gyrus, and the middle frontal gyrus with the precentral cortex. Of note, connectivity changes involving these regions have been reported previously in fMRI studies of focal dystonia (27,28). The delineation of abnormal functional connections such as these may help customize targets for therapeutic intervention in individual patients. Indeed, a recent deep brain stimulation (DBS) study found that cervical and generalized dystonia patients benefited clinically from distinct DBS stimulation pathways (25).
The study of H-DytRP connectivity in dystonia revealed parallels with recent data on the influence of genotype on the organization of Parkinson's disease networks. Specifically, in both disorders, these effects were most pronounced in core subgraphs defined independently through community detection algorithms (19). Indeed, in MAN, prominent connectional gain was evident within the H-DytRP core, whereas in the other groups, this effect was greater outside this subgraph. The shift of the abnormal connections from the core in MAN to the periphery in SPOR may be clinically relevant. The H-DytRP periphery is composed largely of frontal, parietal, and occipital association regions, and the current data suggest a preponderance of such connections in SPOR patients. In hereditary dystonia, by contrast, connectional gain was more prominent in Vo et al. 13 the core, linking cerebellar, thalamic, striatal, and motor cortical nodes within this module. As noted above, significant group differences in connectional loss were also observed in the core, but these were most pronounced in NM carriers as opposed to MAN and SPOR patients.
The meaning of the loss of normal connections in dystonia is less clear. It is tempting to associate reductions in functional connectivity as seen in the core zone of NM with microstructural changes in the integrity of thalamic outflow pathways to the motor cortex and striatum (22,29). The prominent loss of functional connections linking these core regions is therefore consistent with earlier DTI studies conducted in DYT1/DYT6 carriers which demonstrated a close relationship (13,30). Of note, prominent loss of functional connectivity was also present in MAN and SPOR, but unlike NM, these changes were predominantly observed outside the core zone. For example, functional connections in the H-DytRP periphery linked the inferior frontal operculum (BA 44) and the angular gyrus (BA 39) in healthy subjects but were not present in either of the two affected groups. Given that activity in these regions was anticorrelated in healthy individuals, it is possible that this projection normally has an inhibitory function that is lost in dystonia patients. Even so, one cannot determine from the rs-fMRI data alone whether the underlying anatomical connections are intact but functionally deactivated through interactions with other involved nodes or subnetworks, or alternatively whether the two regions are anatomically disconnected on a neurodevelopmental basis. DTI tractography can help disambiguate these possibilities. In this case, we observed a substantial reduction in the number of fiber tracts connecting these regions in MAN compared to HC subjects (31). It is possible that similar changes underlie deficits in higher order functions such as sequence learning and visual motion perception reported in dystonia patients (31)(32)(33). Further multimodal imaging studies with rs-fMRI and DTI will be needed to evaluate these and other non-motor manifestations of Vo et al. 14 the disorder.

Pathological and adaptive network configurations in dystonia
The changes in functional connectivity seen in dystonia were not limited to preferential gain and loss of links in specific H-DytRP subspaces. However, by quantifying a small number of prespecified graph metrics, we obtained valuable information on the connectivity patterns that characterized the various groups. As with gain and loss, the changes in network architecture denoted by the metrics were most pronounced in the core zone. In MAN, degree centrality and characteristic path length were significantly increased in this module, which was consistent with the substantial gain in connections that was observed. That said, the clustering coefficient, an index of parallel processing, was reduced in the core zone of MAN patients. By the same token, core small-worldness is also reduced in this group, denoting an imbalance between segregation and integration of neural signal and inefficient information transfer through this module. In addition, the core zone in MAN exhibited high assortativity, i.e., the tendency for connections to form between nodes of similar degree centrality (34)(35)(36). This gives rise to relatively homogeneous nodal interactions, a feature of unstable, pathological connectivity responses in disease networks (20). The H-DytRP core also exhibited a rich club, i.e., a set of mutually interconnected high degree nodes, in MAN but not in the other groups. Rich clubs are an important means of facilitating functional integration in normal brain networks (37)(38)(39)(40). They may, however, also signify subspaces with enhanced vulnerability to pathology in neurodegenerative disorders (41,42). The current data suggest that rich club organization may also play a pathological role in neurodevelopmental disorders such as hereditary dystonia.
The pattern of functional connectivity recorded in the H-DytRP core was substantially Vo et al. 15 different in NM carriers. In contrast to MAN, core degree centrality and characteristic path length were normal in this group, while the clustering coefficient and small-worldness were both increased in this subgraph. These findings, along with reduced assortativity, i.e., heterogeneity of nodal interactions, are consistent with enhanced adaptive capacity of a network in response to perturbation (43,44). In the context of disease networks, this connectivity pattern suggests a beneficial adaptation, analogous to that observed in PD patients with the slowly progressive LRRK2-G2019S genotype or sporadic PD patients following network rewiring as a consequence of subthalamic gene therapy (19)(20)(21).
The distinctive network configurations observed in the H-DytRP core in MAN and NM can be further considered in the context of self-organization criticality, the process in which network performance is optimized to a critical point or a range of configurations (45,46). Selforganized criticality is inherently non-linear (47): network behavior is chaotic and unstable for configurations sampled above the critical point (supercritical range). Below this point (subcritical range), however, the network behavior is quiescent, but performance is suboptimal. Between these extremes lies a narrow range of configurations (critical region), in which network performance is optimized for stable and efficient information processing (48).
Recent studies of experimental systems suggest that networks can assume different configurations as criticality develops. Indeed, in dissociated cortical cultures from the early developmental period, self-organized criticality was reached after a stereotyped sequence of maturational steps (46,49,50). According to this model, network development begins in a quiescent state in the subcritical regime. Subsequently, dysregulated supercritical behavior emerges as increasing numbers of connections are formed that give rise to homogeneous nodal interactions. A critical point is reached later, however, after a period of synaptic pruning in Vo et al. 16 which optimal network efficiency is achieved by removing excess connections. One might speculate that this process is incomplete in DYT1/DYT6 carriers: network maturation is arrested in the initial quiescent phase in NM, and in the later uncontrolled phase in MAN. Given that the process of self-organization is itself non-linear (47), transitions from one network configuration to another are unlikely to occur spontaneously, particularly after the system has matured to a certain level. That said, arrest of network development in the supercritical phase, as proposed for MAN, may involve additional features associated with the disease process. For example, the appearance of a set of tightly interconnected high degree nodes in the core, referred to as a "richcore" (39), may enhance vulnerability to extrinsic insults occurring later in brain development (51). Self-organized criticality is also influenced by factors distinct from network mesostructure, such as GABAergic inhibition and plasticity of connections in key modules (46,52,53). The role of each of these variables in an individual patient may depend on genotype as well as environmental factors.

Hereditary and sporadic dystonia networks as disease markers
The connectivity changes in the H-DytRP space seen in MAN and NM were not necessarily generalizable to SPOR. Firstly, while H-DytRP and S-DytRP topographies are similar, they are by no means identical. The same ICs were selected by the algorithm when it was applied to the MAN patients or to combined MAN and SPOR sample, suggesting that the two networks share a common topographic signature. Even so, contributions from the cerebellum (IC18) were accentuated in the former network and relatively reduced in the latter, while those from sensorimotor and occipital association cortex (IC10) were in the opposite direction. It is likely that the SPOR group, which included patients with focal, segmental, and generalized dystonia, Vo et al. 17 were too heterogeneous for reliable IC selection by our algorithm. It is possible that an S-DytRP identified in a larger, more clinically homogeneous SPOR sample would exhibit a somewhat different topography than currently observed. By the same token, the mesostructure of the H-DytRP and S-DytRP also cannot be viewed as equivalent, in that the nodes that comprise the core and periphery may differ for the two networks. Indeed, in SPOR, the absence of abnormalities in some graph metrics may be explained by a shift in the boundaries of the core zone for the two dystonia networks. Despite these limitations, generating S-DytRP from the ICs identified in the MAN data proved to be a reasonable strategy, given the close relationship of H-DytRP and S-DytRP expression values measured in the two patient groups. Importantly, the presence or absence of motor manifestations in DYT1/DYT6 carriers was not explained by differences in H-DytRP or S-DytRP expression. That said, functional connections linking the thalamus and Rolandic gyrus were seen in MAN, but not in NM. This observation mirrors the earlier DTI studies, which points to penetrance-related differences in the microstructure of cerebellothalamocortical (CbTC) motor pathway in dystonia gene carriers (22,30). Whereas cerebellothalamic fiber tracks, which constitute the proximal segment of the pathway, were reduced in both MAN and NM, thalamocortical projections along the distal segment were reduced only in the latter group. These data accord with the current functional connectivity studies, suggesting that transmission of aberrant cerebellar output to the motor cortex -by way of relatively intact thalamocortical pathways -is essential for clinical penetrance in this population. Despite the microstructural changes in both CbTC segments observed in NM, the cerebellar vermis and the precentral and paracentral gyri can still interact through a set of functional connections that are present in these individuals but not in HC subjects (Fig. 4B, middle; Table S3A). Given that these abnormal connections were identified in Vo et al. 18 NM but not in MAN, they may play an adaptive role in the disorder.
The network data may also help explain differences in the severity of motor symptoms among dystonia patients. In an earlier DTI study, we found that clinical dystonia ratings for affected upper limbs were worse in those with more intact structural connectivity (higher FA) in the distal (thalamocortical) segment, in the CbTC pathway measured in a fixed subrolandic white matter (SWM) volume (30) (Fig. S3C). In aggregate, the findings suggest that individual patient differences in dystonia severity can be explained by a combination of dystonia network activity measured over the whole brain and local microstructural integrity measured in the distal CbTC segment. As a functional marker, however, dystonia network expression is likely to be more sensitive to the effects of treatment than DTI measurements. Further studies are needed, however, to validate these measures, singly or in combination, as potential biomarkers of primary dystonia.

Limitations
Several limitations should be considered in the interpretation of these data. Despite the power of the dystonia networks to detect significant group differences in the testing data, we acknowledge the need for further validation of the respective topographies in larger, multisite training sets.

Vo et al. 19
While H-and S-DytRP are topographically similar, we cannot exclude specific subnetworks that may be unique to SPOR patients. In this regard, network mesostructures may also vary across groups, with different core-periphery arrangements and connectivity patterns for MAN-and SPOR-related topographies. As stated above, limiting comparisons to similar clinical phenotypes may also be helpful in this regard. Another limitation was that gain and loss of functional connectivity in the H-DytRP space was not systematically compared with DTI measures of structural integrity along the same pathways. The focus here was on the CbTC motor system, which had previously been associated with penetrance and phenotypic variation, as well as clinical disability. A deeper multimodal analysis is needed to map connectivity changes in nonmotor pathways (31), and to compare network metrics for graphs defined by structural versus functional nodal interactions. Finally, the rs-fMRI approach to network identification can potentially be applied to individual scan data to map the network on a single subject basis. While potentially useful for customized interventions, such individualized reconstructions can be limited by noise. Methods are currently being developed to improve the accuracy of this approach for use in individual patients. Even so, the utility of the general approach in two separate brain disease, Parkinson's disease, which is neurodegenerative, and primary dystonia, which is likely neurodevelopmental, supports its potential value in the study of other CNS disorders.
Ten of the patients were manifesting carriers (MAN) of the DYT1 (n=6) or DYT6 (n=4), and 10 Vo et al. 20 had sporadic dystonia (SPOR). We additionally studied 10 age-matched non-manifesting (NM) mutation carriers (4 M/6 F, age 46.4±16.9 years). Twenty healthy volunteers (9 M/11 F, age 47.6 ± 7.8 years) served as controls. The healthy control (HC) subjects were divided into two groups of 10 (HC1 and HC2) based on age and gender. The HC1 subjects (4M/6F, age 46.7 ± 9.0 years) had similar demographics to the dystonia patients; their scans were used in the network identification procedure (see below). Scans from the HC2 subjects (5M/5F, age 48.5 ± 6.6 years) were used for testing. The clinical and demographic features of these patient and control groups are presented in Table S1.
All subjects underwent T1-weighted imaging, diffusion tensor imaging (DTI) and rs- Ethical permission for these studies was obtained from the Institutional Review Board of Northwell Health. Written consent was obtained from each subject following detailed explanation of the procedures.

Resting-state functional magnetic resonance imaging
Subjects were scanned in an awake resting state; no specific task was performed in this condition. Subjects receiving botulinum toxin treatment were scanned 3-4 months after the last  (55). The preprocessing included motion correction, brain extraction, spatial smoothing (kernel=8 mm; full width at half maximum) and temporal high-pass filtering (cutoff frequency=1/150 Hz).
To evaluate the effects of motion during rs-fMRI, we examined absolute and relative displacements during the scan. There was less than 1.4 mm absolute displacement and less than 0.35 mm relative displacement, in all scans without significant differences between groups (absolute: p>0.50; relative: p>0.43; Student's t-tests) ( Table S1). The resulting fMRI volumes were registered to the individual subject's structural T1 image and then to the standard Montreal Neurological Institute (MNI) 152 template. The rs-fMRI data were intensity normalized to reduce variability and improve the reliability of ICA maps.

Identification and validation of dystonia-related patterns
Two approaches were employed to identify and validate dystonia-related networks in the rs-fMRI data. In the first, we identified a hereditary dystonia-related pattern (H-DytRP) in scans from manifesting (MAN) dystonia gene carriers and age/gender-matched healthy control subjects (HC1). We then prospectively measured individual H-DytRP expression levels (subject scores) in the sporadic dystonia (SPOR) patients, the non-manifesting (NM) mutation carriers, and in the remaining healthy subjects (HC2). In the second, we identified a sporadic dystonia-related pattern (S-DytRP) in scans from SPOR and HC1 subjects, and prospectively measured expression levels for this network in the MAN, NM, and HC2 groups. The specific computational steps employed for network identification and validation are outlined in Fig. S1 and presented in detail elsewhere (15,17).

Vo et al. 22
Step 1: The rs-fMRI data from all 50 participants were analyzed using group-wise spatial ICA (56) with GIFT software (http://mialab.mrn.org/); 50 group independent components (ICs), and associated time courses and power spectra, were obtained. Of these, 10 ICs were discarded because of motion artifacts and/or scanner-related or physiological noise on visual inspection.
The procedures for network identification and testing that follow were applied to the remaining 40 ICs.
Step 2: Subject spatial maps and temporal dynamics were estimated using dual regression (57).
Step 3: Subject scores, representing pattern expressions in individual rs-fMRI scans, were computed for the 40 ICs in each of the patients and control subjects. The subject score of component in subject (denoted ! (#) ) is the scalar projection of the group spatial map for component onto the subject spatial map for component in subject .
Step 4: The resulting measures for the training dataset were entered into a logistic regression model with bootstrap resampling (1000 iterations) to identify a subset of ICs whose subject scores best discriminated between the two groups. The weight (coefficient) for each of the selected ICs was estimated through an additional bootstrap resampling procedure (1000 iterations) (15).
Step 5: Using the coefficients estimated in Step 4, we linearly combined the selected components to form a composite disease network associated with hereditary or sporadic dystonia, i.e., H-DytRP or S-DytRP. The same coefficients are applied to the corresponding subject scores to quantify expression levels for the respective networks.
Steps 1-5 were performed for network identification. The derived H-DytRP and S-DytRP topographies were compared by computing voxel-wise correlation coefficients, which were

Clinical correlations of network expression
To assess the relationship of dystonia network expression and the severity of clinical manifestations of the disorder, we separately correlated H-DytRP and S-DytRP subject scores with BFMDRS motor ratings in the combined MAN and SPOR patient sample. This was done by computing Pearson product-moment correlation coefficients, which were considered significant for p<0.05.

Vo et al. 24
In an earlier study using magnetic resonance DTI in primary dystonia, we found that BFMDRS motor ratings correlated with the microstructural integrity of the distal, i.e., thalamocortical, segment of CbTC motor pathway (30). carriers and HC subjects (Fig. S2) (22,30). This volume was defined by significant bilateral FA reductions in NM compared to either MAN or HC (22). Given that SWM FA was also found to correlate with BFMDRS motor ratings (30), we entered these values, averaged across hemispheres, into linear regression models as an additional predictor of clinical severity. The best quality model, i.e., that with the greatest generalizability as defined by the lowest out-ofsample prediction error, was selected based upon the Akaike Information Criterion (AIC) (61).
The relationship between each predictor and motor ratings, independent of variation in the other predictor, was displayed using partial correlation leverage plots (Statistics and Machine Learning Vo et al. 25 Toolbox, MATLAB R2020a), and were considered significant for p<0.05.

Connectivity maps in the network space
To document the abnormal node-to-node interactions that underlie the dystonia-related functional topographies, we normalized the 3D voxel-wise network displays to MNI space and thresholded the images at T=4.8 (p<0.001). Using the AAL atlas (62), we parcellated the brain into 95 regions-of-interest (ROIs) and identified those corresponding spatially to significant network clusters (18)(19)(20)(21). For each node, we computed mean time series for rs-fMRI scans in the By this scheme, the magnitude of the correlation (|r|) provided a measure of connectivity between network nodes for each group. For a given pair of nodes, group differences in connectivity were described by the absolute difference (|dr|) in the two correlation coefficients (19,21). For the gain of a connection in a dystonia group (i.e., MAN or SPOR) or mutation carriers (i.e., NM) to be significant, we required that the magnitude of the correlation coefficient (|r|) that defined the graphical edge be greater than or equal to 0.6 (p<0.05; Pearson correlation) in that group but not in HC, and that the corresponding absolute difference (|dr|) from HC be greater than 0.4 (p<0.05; permutation test, 1000 iterations). The latter threshold was determined Vo et al. 26 using the HC graph and permuting the regional labels 1000 times to create a set of pseudorandom correlations for each iteration as described previously (19,21). Likewise, the loss of a normal connection in MAN, NM, and/or SPOR was significant if |r| for the graphical edge was greater than or equal to 0.6 in HC but not in the particular disease group, and that the corresponding absolute difference from HC was greater than 0.4. For each edge, the connections that satisfied these criteria were confirmed by bootstrap resampling (100 iterations) using the Statistics and Machine Learning Toolbox in MATLAB R2020a. The resulting connections were considered for further analysis only if: (1) they also conformed to known anatomical pathways determined by diffusion tractography and/or postmortem analysis, and (2) the two nodes were separated by no more than two sequential hops along the graph.
Additionally, connections were mapped within and between nodal communities defined according to a graph partitioning algorithm (19). This was done based on a modularity maximization algorithm as described elsewhere (63)(64)(65), utilizing the "core_periphery_dir" function in the Brain Connectivity Toolbox (66). According to this method, the optimal core/periphery subdivision is a partition of the network into two non-overlapping zones, such that the number/weight of within-core edges is maximized while the number/weight of withinperiphery edges is minimized (63,67,68). In this study, we used this approach to partition the dystonia network into core and periphery subgraphs for further analysis. Specifically, we counted the total number of connections gained relative to HC in the network space for the MAN, NM, and SPOR groups. Network connections were further classified according to edge type: corecore, periphery-periphery, or core-periphery. Bootstrap resampling (100 iterations) was used to estimate the mean and standard error of this measure for each group and edge category.

Network metrics
To evaluate connectivity patterns in network subspaces, we computed the following graph metrics for each group on weighted undirected graphical links as described elsewhere (19)(20)(21)69): 1. Degree centrality: the number of connections (edges) within the network divided by the total number of nodes in the same space. This is a measure of overall connectivity of nodes within a given network or subgraph.
2. Normalized clustering coefficient: the likelihood that the nearest neighbors of a node will also be connected. This measure provides an index of parallel processing in locally interconnected (closed triples) in a network or subgraph.
3. Normalized characteristic path length: the shortest path length between two nodes averaged over all pairs of nodes in a given network (34,66). This measure was also normalized to the corresponding value from an equivalent random graph.
4. Small-worldness: the ratio of clustering coefficient to characteristic path length, normalized to the corresponding value from an equivalent random graph (70). This measure quantifies the ratio of segregation to integration of information sources in the network space.

5.
Assortativity coefficient: the correlation coefficient between the degrees of all nodes on two opposite ends of a link (35,36,71). For a given network, the coefficient is described as assortative for positive values, neutral for ≈ 0, and disassortative for negative values.
6. Rich-club coefficient: the fraction of edges that connect a set of nodes of degree k or greater relative to the maximum number of edges possible in the subnetwork (72). The measure assays the ability of the "rich club" nodes to influence global network function by virtue of their dense interconnections (38,39,73).

Vo et al. 28
These parameters were computed using the Brain Connectivity Toolbox (66) and an inhouse Matlab script (MATLAB R2020a). We note that group differences may be difficult to discern because of inclusion of random, non-specific links at low thresholds (r<0.3; graph density >60%), and graph disconnection at high thresholds (r>0.6; graph density <25%) (19)(20)(21).
Network metrics are therefore presented for a range of connectivity thresholds (r=0.3 to 0.6, at 0.05 increments) corresponding to graph densities between 25% and 60%. By plotting the results over the range, we show that group differences for a given metric are valid over multiple adjacent levels.
To visualize group differences in connectivity patterns, relevant network subspaces were displayed at minimum threshold (Level 1, r=0.30). For specific metrics such as the assortativity coefficient, exemplars in the range of the mean value ±0.5 SD were selected from the bootstrap samples described above. Individual graph configurations were represented as force-directed displays (74) using an in-house Matlab script (mathematics/graph and network algorithms toolbox, MATLAB R2020a).

Statistical analysis
Expression values for H-DytRP and S-DytRP were computed in the MAN, NM, SPOR, and HC groups and compared using ANOVA with the Tukey-Kramer HSD adjustment for multiple comparisons. The gain and loss in H-DytRP connections in MAN, NM, and SPOR (i.e., those present in these groups but not in HC (gain) or those present in HC group but not in these groups (loss)) was computed as a percentage of the total within and between core and periphery subgraphs, as well as across the whole network. The difference between the connections gained and lost inside (core-core) versus outside (core-periphery and periphery-periphery) the H- Vo et al. 29 DytRP core for each dystonia group (100 bootstrap iterations) was also computed, and group differences were evaluated as described above. These analyses were considered significant for p<0.05, corrected for multiple comparisons.
For graph analysis, the bootstrapped data were used to assess group differences in the network metrics for the relevant graph or subgraph. Group differences in each of the metrics were evaluated using the general linear model (GLM) with post-hoc Tukey-Kramer HSD tests across graph thresholds. These analyses were performed using IBM SPSS Statistics for Windows, version 21 (IBM Corp., Armonk, N.Y., USA). Results were considered significant for p<0.05.

Data availability
Deidentified data will be made available on reasonable request from interested investigators for the purpose of replicating results.  (75). c AAL = Anatomical Automatic Labeling atlas (62). The number given for each significant region denotes the standardized region-of-interest (ROI) from the atlas that was used in the graph theory analysis (see text). BA = Brodmann area. Bold font indicates core nodes and regular font indicates periphery nodes.