Abstract
Isolated focal dystonia is a debilitating movement disorder of unknown pathophysiology. Early studies in focal dystonias have pointed to segregated changes in brain activity and connectivity. Only recently has the notion that dystonia pathophysiology may lie in abnormalities of large-scale brain networks appeared in the literature. Here, we outline a novel concept of functional connectome-wide alterations that are linked to dystonia phenotype and genotype. Using a neural community detection strategy and graph theoretical analysis of functional MRI data in human patients with the laryngeal form of dystonia (LD) and healthy controls (both males and females), we identified an abnormally widespread hub formation in LD, which particularly affected the primary sensorimotor and parietal cortices and thalamus. Left thalamic regions formed a delineated functional community that highlighted differences in network topology between LD patients with and without family history of dystonia. Conversely, marked differences in the topological organization of parietal regions were found between phenotypically different forms of LD. The interface between sporadic genotype and adductor phenotype of LD yielded four functional communities that were primarily governed by intramodular hub regions. Conversely, the interface between familial genotype and abductor phenotype was associated with numerous long-range hub nodes and an abnormal integration of left thalamus and basal ganglia. Our findings provide the first comprehensive atlas of functional topology across different phenotypes and genotypes of focal dystonia. As such, this study constitutes an important step toward defining dystonia as a large-scale network disorder, understanding its causative pathophysiology, and identifying disorder-specific markers.
SIGNIFICANCE STATEMENT The architecture of the functional connectome in focal dystonia was analyzed in a large population of patients with laryngeal dystonia. Breaking with the empirical concept of dystonia as a basal ganglia disorder, we discovered large-scale alterations of neural communities that are significantly influenced by the disorder's clinical phenotype and genotype.
Introduction
Laryngeal dystonia (LD), or spasmodic dysphonia, is an isolated task-specific focal dystonia that selectively affects the production of speech due to impaired voluntary control of vocal fold movements caused by involuntary spasms in the laryngeal muscles. As with other focal dystonias, the exact causative pathophysiology of this debilitating movement disorder remains unclear. However, several neuroimaging studies have demonstrated disorder-characteristic functional and structural abnormalities in the laryngeal/orofacial sensorimotor and parietal cortices, basal ganglia, thalamus, and cerebellum, as well as altered striatal dopaminergic and GABAergic function in LD patients compared with healthy individuals (Haslinger et al., 2005; Ali et al., 2006; Simonyan et al., 2008; Simonyan and Ludlow, 2010, 2012; Simonyan et al., 2013; Samargia et al., 2016). More recent reports of the presence of dystonia in up to 16% of patient families (Kirke et al., 2015), threefold penetrance of dystonia mediational endophenotype in familial compared with sporadic LD cases (Termsarasab et al., 2016), and the association of brain abnormalities with gene mutations (i.e., DYT 25) and the polygenic risk (Putzel et al., 2016a, 2016b) have laid the foundation for a consideration of the interplay between potential etiological genetic factors and pathophysiological mechanisms that may lead to distinct clinical phenotypes and putative genotypes of LD.
Previously, we described LD phenotype- and genotype-specific alterations of functional activity within the sensorimotor and frontoparietal networks, specifically involving primary sensorimotor, premotor, and inferior parietal cortices (Battistella et al., 2016). LD phenotype-specific structural abnormalities have been reported in the primary and associative motor regions, whereas genotype-related structural changes have been found in the brain regions regulating phonological and sensory processing (Bianchi et al., 2017). Recent experimental evidence has also suggested that focal dystonias in general, including LD, may represent a disorder of large-scale functional networks as opposed to a basal ganglia pathology alone (Battistella et al., 2017).
In support of this emerging multifaceted and novel view on the causative pathophysiology of LD, we conducted detailed qualitative and quantitative investigations of the architecture of large-scale functional brain networks in a uniquely large population of LD patients compared with healthy subjects. Taking into account the potential for distinct pathophysiological mechanisms underlying different phenotypes and genotypes of this disorder, we examined the large-scale network topology across LD genotypes (sporadic and familial forms) and its clinical phenotypes (adductor and abductor forms). Based on the clinical history of the disorder, a patient without any familial history of LD and/or other isolated dystonias was considered to have the sporadic form, whereas a patient with at least one other blood relative affected with LD and/or other isolated dystonias was considered to have familial LD. Phenotype classifications were based on a patient's clinical diagnosis: adductor LD (ADLD) was characterized by involuntary vocal fold closure causing voice breaks and strained voice quality on vowels, whereas abductor LD (ABLD) showed intermittent breathy voice breaks extending voiceless consonants.
Using a graph theoretical approach (Bullmore and Sporns, 2009), we quantified large-scale network architecture by investigating the influence of a single region on network communication pathways (i.e., local network features) and the formation of neural communities of dense regional connectivity profiles (i.e., global network features). This approach allowed us to assess regional and network-wide aberrations related to LD, permitting a detailed investigation of whole-brain dystonic activity not limited to regions traditionally focused on in studies of LD pathophysiology such as basal ganglia and cerebellum. Based on previous knowledge, we hypothesized that the functional embedding of central motor control areas would be abnormally weak in LD accompanied by an overrepresentation of basal ganglia and thalamus within the network. We further hypothesized that our approach to a tiered, bottom-up analysis of network differences in the functional footprint of LD phenotypes and genotypes would allow us to identify genotype- and phenotype-specific differences in the dystonic cascade, which may ultimately underlie LD pathophysiology.
Materials and Methods
Subjects
A total of 90 LD patients (54.5 ± 12.9 years; 73 females/17 males) and 32 healthy individuals (50.5 ± 10.2 years; 21 females/12 males) participated in this study (see detailed demographics in Table 1). All subjects were right-handed, native English speakers and none had any history of neurological (other than LD in the patient group), psychiatric, or laryngeal problems. All participants underwent neuroradiological evaluation, which found normal brain anatomy. Those LD patients who received botulinum toxin injections as part of their clinical symptom management participated in the study at least 3 months after the last treatment to ensure that they were fully symptomatic at the time of study participation. Genetic screening confirmed that none of the participants were carriers of TOR1A (DYT1), THAP1 (DYT6), TUBB4A (DYT4), or GNAL (DYT25) mutations.
Participant demographics
Patients and healthy individuals were assigned to experimental groups as follows. Group 1 was a comparison between the LD and normal states and included 32 LD patients (49.8 ± 11.1 years; 22 females/10 males) and 33 healthy individuals (50.5 ± 10.2 years; 21 females/12 males). To allow for a balanced representation of different forms of LD, this patient group comprised 16 sporadic (8 ADLD/8 ABLD) and 16 familial LD cases (9 ADLD/7 ABLD). Group 2 was a comparison between putative LD genotypes and included 27 sporadic LD patients without any past or present history of dystonia in family members (57.6 ± 13.0 years; 22 females/5 males; 20 ADLD/7 ABLD) and 27 familial LD patients with a history of LD and/or other isolated dystonia in one or more family members (56.6 ± 14.9 years; 22 females/5 males; 20 ADLD/7 ABLD). As a next step, to assess potential phenotype-related influences on LD genotype, we assembled more homogeneous patient groups with only the ADLD phenotype (20 sporadic ADLD: 58.25 ± 11.7, 17 females/3 males vs 20 familial ADLD: 57.25 ± 15.0 years, 17 females/3 males). Because our study included only seven familial ABLD patients due to the extreme rarity of this particular form of the disorder, this cohort was not powered enough for a comparison with sporadic ABLD patients. This constituted a limitation of the current study, which should be addressed in a future study with a larger sample of familial ABLD patients. Group 3 was a comparison between different LD phenotypes and included 32 sporadic ADLD patients (55.1 ± 11.2 years; 24 females/8 males) and 31 sporadic ABLD patients (52.2 ± 12.5 years; 27 females/4 males). Due to the low number of familial ABLD patients, we did not compare this group with familial ADLD patients. However, we assessed the influence of putative genotype on ADLD and ABLD phenotypes by assembling more heterogeneous groups of mixed-genotype ADLD patients (55.1 ± 12.7 years; 29 females/10 males; 32 sporadic/7 familial ADLD) and mixed-genotype ABLD patients (52.6 ± 13.0 years; 32 females/6 males; 31 sporadic/7 familial ABLD).
Written informed consent was obtained from all subjects before study participation and the study was approved by the Institutional Review Board of the Icahn School of Medicine at Mount Sinai. Some subjects participated in other studies conducted around the same time that assessed the resting-state and structural brain organization and neural correlates of abnormal sensory discrimination in LD (Battistella et al., 2016; Termsarasab et al., 2016; Bianchi et al., 2017).
Experimental design
Whole-brain functional images were acquired with a gradient-weighted echo-planar imaging (EPI) pulse sequence using an event-related sparse-sampling design (TR = 2 s per volume and 8.6 s between volumes with total TR = 10.6 s, TE = 30 ms, FA = 90, FOV = 240 mm, voxel size = 3.75 × 3.75 mm, 36 slices with 4 mm slice thickness). Four functional runs, including 24 sentence and 16 resting conditions, were acquired per subject. Subjects were instructed to listen to a meaningful grammatically correct English sentence presented via MR-compatible headphones within a 3.6-s interval, then repeat the same sentence within a 5-s period, after which a functional image was acquired within a 2-s interval while subjects silently fixated on a cross presented on a screen. Subjects listened to 10 sample sentences containing a high load of vowels (e.g., “Are the olives large?”; “Jack sat on a tack”) to predominantly elicit ADLD symptoms and a high load of voiceless consonants (e.g., “He is hiding behind the house”; “Who says a mahogany highboy isn't heavy?”) to predominantly elicit ABLD symptoms. All sentences were presented in a random order one at a time. Resting trials, during which the subjects were asked to silently fixate on a cross, were acquired as a baseline and pseudorandomized with the sentence production stimuli. To minimize motion artifacts, each subject's head was tightly cushioned within the coil, and all participants were monitored constantly for possible movements during the scanning session.
A high-resolution sagittal T1-weighted image was acquired using a magnetization prepared rapid acquisition gradient echo (MPRAGE) sequence (172 contiguous slices, 1 mm3 voxel-size, TR = 7.5 ms, TE = 3.5 ms, FOV = 210 mm) for registration of the functional images and to rule out any structural lesions.
Functional data analysis was performed using AFNI software (Cox, 1996) following a standard image preprocessing pipeline as described previously (Simonyan and Ludlow, 2012; Fuertinger et al., 2015; Termsarasab et al., 2016). Briefly, the first four volumes of all EPI datasets were removed to account for potential T1 stabilization effects and subsequently registered to the volume temporally closest to the high-resolution anatomical image. After spatial smoothing with a 4 mm Gaussian kernel, we performed multiple linear regression with the task regressor convolved with a canonical hemodynamic response function, six motion parameters (3D translation, pitch, roll, yaw rotation) to account for residual motion and three polynomials to control for low-frequency components such as scanner drift (Perrachione and Ghosh, 2013; Battistella et al., 2016; Termsarasab et al., 2016).
Functional connectivity analysis: construction of functional networks
To construct whole-brain functional networks of regional interactions, we used 212 regions of interest (ROIs) based on the cytoarchitectonic maximum probability maps and macrolabel atlas (Eickhoff et al., 2005; Fuertinger et al., 2015), including 142 cortical, 36 subcortical, and 34 cerebellar areas. In all subjects, pairwise interactions between regionally averaged time series during speech production were estimated by computing normalized mutual information (NMI) coefficients (Strehl and Ghosh, 2002; Fuertinger and Simonyan, 2016). The NMI was obtained by dividing the classic mutual information (Cover and Thomas, 1991) by the geometric mean of the associated Shannon entropies (Shannon, 1948), yielding a statistical dependence measure with values between zero (statistical independence) and one (mutual dependence). Therefore, in contrast to classic Pearson correlations, NMI edges are always non-negative, circumventing the need to account for negative weights in the computation of nodal influence metrics or to control for the emergence of negative cycles particularly affecting path-based measures (Fakcharoenphol and Rao, 2006). Because NMI coefficients reflect linear correlations and dependencies in higher moments (Brillinger, 2002), a nonzero Pearson correlation necessarily implies a positive NMI coefficient. Therefore, by definition, the nonzero structure of a pairwise Pearson correlation matrix is preserved in an NMI matrix, and NMI is not antagonistic but rather complementary to linear correlation (Kinney and Atwal, 2014).
In each subject, the constructed 212 × 212 NMI matrix was interpreted as the connectivity matrix of a weighted undirected graph, with ROIs representing the graph's nodes and NMI values serving as edge weights. The first estimate of each graph's connectivity profile was obtained by computing the network density, which was represented by the number of present to maximally possible edges in a graph. All subject networks were fully connected with a density of 100% and were thus subjected to thresholding (Fuertinger et al., 2015) following the rationale that very dense networks (graphs with a density >50%) tend to exhibit random graph characteristics (Humphries et al., 2006; Lynall et al., 2010). The density of the networks was reduced to 50% by removing edges with a weight less than an iteratively incremented percentage of the maximum weight in the graph (starting with 1%) while simultaneously ensuring that the network did not fragment by fully disconnecting nodes. After the thresholding of each individual graph, group-averaged networks were computed within each experimental group.
Graph theoretical analysis and statistical analysis
Network processing and visualization codes were written in Python using the open-source packages NumPy, SciPy (van der Walt et al., 2011), and Matplotlib (Hunter, 2007). Optimal modular decompositions were calculated in MATLAB (The MathWorks) using the Brian Connectivity Toolbox (Rubinov and Sporns, 2010). 3D renderings of the networks embedded in reference brain models were generated using the BrainNet Viewer (Xia et al., 2013).
Network metrics.
Functional networks were analyzed to assess network segregation, integration, and nodal influence. Network segregation was estimated by computing the nodal clustering coefficient as the geometric mean of weights in triangles around each node (Onnela et al., 2005). In this manner, the presence of functional cliques in a node's local neighborhood was quantified. Network integration was approximated using the measure of global efficiency estimating a network's predilection for long-range communication (Sporns, 2011). Global efficiency was calculated as the average inverse shortest path length in the network (Latora and Marchiori, 2001) with connection lengths based on inverse edge weights. The influence of a single node within the graph was quantified by computing the nodal degree (i.e., the number of connected edges) and the nodal strength (i.e., the sum of connected edge weights; Rubinov and Sporns, 2010). Although the information provided by clustering coefficient, local efficiency, and modular organization of a network partially overlaps, the combination of these metrics provides a complete picture of local and global network properties (Sporns, 2011). The statistical significance of differences in network metrics across groups was assessed by using a two-tailed permutation t test using 10,000 Monte Carlo randomizations at p ≤ 0.05 adjusted for multiple comparisons by a Tmax approach for the maximal test statistic (Pesarin, 2001).
Community-based network organization.
To examine the global reconfiguration of functional networks in normal and diseased states, we performed a graph community analysis by estimating the optimal modular decomposition of a network, which divides a graph into nonoverlapping, connected groups of nodes (i.e., modules) with a maximal number of within-module edges and minimal number of between-module links (Fortunato, 2010; Fuertinger et al., 2015, 2016). A heuristic optimization strategy based on the Kernighan–Lin algorithm (Sun et al., 2009) was used to maximize the Newman modularity (Newman, 2006) to approximate each graph's optimal modular decomposition. An initial artificial community assignment, in which every node formed a module, was iteratively refined by solving the modularity maximization problem 100 times per network to account for the random nature of the used optimization algorithm, as described previously (Fuertinger et al., 2015). The final community affiliation was computed based on the average nodal module assignment, with node #1 serving as a reference to account for random module numbering. As a first step to estimate the similarity of community structures between groups, we calculated their respective partition distances pd, which represented the normalized variation of information between the respective community affiliation vectors (Meila, 2007). In a follow-up analysis, we examined differences in spatial community structure and hub formation between the groups.
Hub formation.
Nodes exhibiting both degree and strength values at least 1 SD larger than the respective average degree and strength in the network were defined as bivariate hubs. Nodes with both degree and strength values in the top 25% of the network but not satisfying the hub criterion were classified as high-influence nodes. To examine hub formation based on the number of connected intermodule vs intramodule edges, the computed optimal modular decomposition of each network was used to calculate the nodal participation coefficient pci (i = 1, …, 212; van den Heuvel and Sporns, 2011; Fuertinger et al., 2015). In a network with m ≥ 2 modules, the maximum value of pci is 1 − 1/m. If a hub exhibited a pci value within 10% of this maximum value, then it was classified as a connector hub linking two or more communities. If a hub showed a lower pci value, it was defined as a provincial hub connecting nodes within the same community.
Results
Group 1: overall characterization of the LD functional connectome
Comparisons of functional networks between healthy individuals and LD patients (including ADLD, ABLD, sporadic, and familial forms) showed significantly increased network segregation (estimated by the clustering coefficient) in the patient group (controls: 0.14 ± 0.02; LD: 0.15 ± 0.02, p < 1e-3) but similar network integration (estimated by the global efficiency) between the two groups (controls: 0.27 ± 0.04; LD: 0.27 ± 0.03, p = 0.9). Nodal influence quantified by degree and strength failed to exhibit pronounced differences between the healthy and patient groups (both p ≥ 0.9; Table 2).
Summary statistics of graph metrics in the group-averaged networks
As hypothesized, the topological architecture of the LD functional connectome demonstrated widespread alterations of large-scale networks in this disorder. The optimal modular decomposition of the group-averaged network in healthy individuals revealed the following community formation: two frontotemporoparietal communities (modules I and III) each with clear hemispheric prevalence (Fig. 1A, red and purple), a bilateral basal ganglia/cerebellar community (module IV, orange), and a bilateral prefrontal community (module II, green). Compared with healthy controls, LD patients showed pronounced changes in network community structure (Fig. 1A) and a pd of 0.24. Specifically, LD formed an additional fifth functional community (module V, blue) by fragmenting the solitary basal ganglia/cerebellar community (module IV) seen in healthy subjects into two communities in patients. The additional module V in LD was left lateralized to the basal ganglia and thalamus (Fig. 1A), pointing to abnormal functional integration of these areas in the large-scale network.
Functional architecture of the group-averaged LD network compared with healthy individuals. A, Spatial topology of detected functional communities in the group-averaged healthy control (left) and LD patient network (right) is shown on 3D brain renderings in the axial and sagittal views, where regional module affiliation is illustrated by nodal color. B, Bivariate (with respect to strength and degree) provincial (yellow) and connector hubs (red) and their connectivity profile with high-influence nodes (purple; nodes with degree and strength in the top 25% of the network) in the group-averaged healthy volunteer network (left) and in the mean LD graph (right). C, Visualization of variations in degree (light gray) and strength (dark gray) across provincial (yellow) and connector hubs (red) in the group-averaged healthy control and LD patient networks (left and right, respectively). White circles indicate nodes that did not qualify as hubs in the respective graph. 3b/6/17, Areas 3b/6/17; 4a, anterior part of area 4; 5L/5M, subdivisions of area 5; 7A/7M/7P/7PC, subdivisions of area 7; CbI-I/IV/Cbl-V/Cbl-VI/Cbl-VIIIav/Cbl-VIIbv, cerebellar lobules I/IV/V/VI/VIIIa vermis; Cu/PCu, cuneus/precuneus; FG, fusiform gyrus; IOG, inferior occipital gyrus; LG, lingual gyrus; MCC/PCC, middle and posterior cingulate cortex; PF, area PF in the inferior parietal cortex; SOG, superior occipital gyrus; Tp/Tpf/Tt, parietal/prefronta/temporal subdivisions of the thalamus; hOC5v, ventral part of area hOC5; L, left; R, right.
Within the detected modules, alterations in hub formation between the mean network of healthy controls and LD patients were observed on several levels, including: (1) degradation of connector to provincial hubs, (2) dropout of regions from hub level to high-influence nodes or below, and (3) formation of new hubs in the LD network that were not present in the healthy network (Fig. 1B). Specifically, the functional connectome of healthy individuals was characterized by an approximately balanced representation of provincial and connector hubs (12/17, respectively), whereas hub formation in the LD network was shaped by a dominant presence of provincial hubs with only few connectors (26/6, respectively). These hubs established strong connections with high-influence nodes, which were underrepresented in LD patients compared with healthy individuals. The brain regions that were most prone to these changes included the primary sensorimotor and parietal cortices, cerebellum, and thalamus (Fig. 1C).
Most notably, left primary somatosensory cortex (area 3b, connector hub), primary motor cortex (area 4a, provincial hub), and premotor cortex (area 6, connector hub) were hubs in healthy subjects, but were not even among high-influence nodes in LD patients, suggesting a pathophysiological functional degradation of these regions in LD. Particularly, area 3b showed the most pronounced drop in nodal influence (degree: −31%, strength: −29%) that was indicative of marked impairment of sensorimotor processing and integration in LD patients. Similarly, the functional impact of bilateral superior parietal cortex, which was found to host hub regions in normal controls (bilateral area 7A, left area 5L: connectors; right area 5L: provincial), was drastically lower in LD patients illustrated by a decrease in degree and strength of up to 16% and 14%, respectively. Therefore, bilateral parietal area 5L failed to show even elevated nodal influence values in the LD network and left and right parietal areas 7A were downgraded to a high-influence node and a provincial hub, respectively.
Conversely, left posterior parietal cortex exhibited abnormally increased functional influence (up to 22% and 30% in degree and strength, respectively) in the LD connectome, upgrading these regions to provincial hubs, albeit with decreased connectivity in the right hemisphere compared with healthy individuals (areas 7P/7PC; degree: −22%, strength: −24%). A similar pattern was found for the cerebellum (right lobules VIIIa/b, vermis) and thalamus, the left parietal and bilateral temporal subdivisions of which played the role of connector and provincial hubs, respectively, in the LD connectome, but were not classified as hubs in healthy volunteers. Finally, nodal influence metrics of the left extrastriate cortex (ventral part of area hOC5), which has been implicated in lateralization of language perception (Selpien et al., 2015) and was a network-wide relay in the healthy speech connectome, dropped by 20% (degree: −20%, strength: −21%) in LD patients, rendering its importance below high-influence nodes.
In summary, the loss of network-wide hub relays in left premotor, primary motor, somatosensory, and parietal regions, together with the emergence of connector hubs in the thalamus as part of the formation of a segregated nodal community in LD points to deficient global functional coupling of principle motor control areas accompanied by abnormal overrepresentation of thalamic activity in the LD functional network.
Group 2: influence of genotype on the LD functional connectome
A comparative analysis of functional networks constructed based on familial versus sporadic LD patients yielded no statistically significant differences in either nodal influence (all p > 0.7) or global efficiency (all p > 0.2) and no pronounced deviations in clustering coefficients (all p > 0.8; Table 2). However, a more detailed analysis of graph segregation and integration based on the community structure of the respective group-averaged networks revealed visible alterations in network topologies across genotypes (Figs. 2, 3; sporadic vs familial cases with mixed heterogeneous ADLD + ABLD phenotype: pd = 0.22 or ADLD-only homogeneous phenotype: pd = 0.25).
Functional architecture of the group-averaged familial LD and sporadic LD networks with mixed phenotype. A, Spatial topology of detected functional communities in the group-averaged mixed-phenotype familial LD (left) and sporadic LD networks (right) is shown on 3D brain renderings in the axial and sagittal views, where regional module affiliation is illustrated by nodal color. B, Bivariate (with respect to strength and degree) provincial (yellow) and connector hubs (red) and their connectivity profile with high-influence nodes (purple; nodes with degree and strength in the top 25% of the network) in the group-averaged mixed-phenotype familial LD (left) and in the mean mixed-phenotype sporadic LD graph (right). C, Visualization of variations in degree (light gray) and strength (dark gray) across provincial (yellow) and connector hubs (red) in the group-averaged mixed phenotype familial and sporadic LD networks (left and right, respectively). White circles indicate nodes that did not qualify as hubs in the respective graph. 6/17/18, areas 6/17/18; 4p, posterior part of area 4; 5L/5M/7P, subdivisions of parietal areas 5 and 7; CbI-I/IV/Cbl-V/Cbl-VI/Cbl-VIIIav/Cbl-VIIa/Cr1, cerebellar lobules I/IV/V/VI/VIIIa vermis/VIIa crus I; Cu/PCu, cuneus/precuneus; FG, fusiform gyrus; IOG/MOG/SOG, inferior/middle/superior orbital gyrus; LG, lingual gyrus; MCC/PCC, middle/posterior cingulate cortex; OP1, subdivision 1 of the operculum; PF/PFcm, areas PF/PFcm in the inferior parietal cortex; Tm/Tp/Tpf/Ts/Tt, motor/parietal/prefrontal/somatosensory/temporal subdivisions of the thalamus; hOC3v/hOC4v, ventral parts of areas hOC3/hOC4; L, left; R, right.
Functional architecture of the group-averaged familial LD and sporadic LD networks with adductor phenotype. A, Spatial topology of detected functional communities in the group-averaged adductor familial LD (left) and adductor sporadic LD network (right) is shown on 3D brain renderings in the axial and sagittal views, where regional module affiliation is illustrated by nodal color. B, Bivariate (with respect to strength and degree) provincial (yellow) and connector hubs (red) and their connectivity profile with high-influence nodes (purple; nodes with degree and strength in the top 25% of the network) in the group-averaged adductor familial LD (left) and in the mean adductor sporadic LD graph (right). C, Visualization of variations in degree (light gray) and strength (dark gray) across provincial (yellow) and connector hubs (red) in the group-averaged adductor familial and adductor sporadic LD networks (left and right, respectively). White circles indicate nodes that did not qualify as hubs in the respective graph. 6/17/18, areas 6/17/18; 4a, anterior part of area 4; 5L/5M/7P, subdivisions of parietal areas 5/7; CbI-I/IV/Cbl-IXv/Cbl-V/Cbl-VI/Cbl-VIIIav/Cbl-VIIIb, cerebellar lobules I/IV/IX vermis/V/VI/VIIIa vermis/VIIIb; Cu/PCu, cuneus/precuneus; FG, fusiform gyrus; IOG/MOG/SOG, inferior/middle/superior occipital gyrus; Ig1, insula subdivision Ig1; LG, lingual gyrus; MCC/PCC, middle/posterior cingulate cortex; OP1, part 1 of the operculum; PF, parietal area PF; TE1.1, area TE1.1; Tm/Tp/Tpf/Ts/Tt, motor/parietal/prefrontal/somatosensory/temporal subdivisions of the thalamus; hOC3v/hOC4v, ventral parts of areas hOC3/hOC4; L, left; R, right.
The partitions of both sporadic networks showed a similar spatial distribution across phenotypes and a partition distance of pd = 0.17 (heterogeneous vs homogeneous). Both heterogeneous and homogeneous sporadic LD networks comprised four communities (Figs. 2A, 3A) with a spatial topology comparable to the modules found in healthy controls (Fig. 1A). Sporadic LD networks were characterized by a large central module that contained 38% and 40% of all nodes in the graph across phenotypes and comprised all cerebellar and most basal ganglia and thalamic regions (Figs. 2A, 3A, orange).
The spatial architecture of communities in familial networks was also consistent across phenotypes and exhibited a partition distance of pd = 0.14 (heterogeneous vs homogeneous). Both heterogeneous and homogeneous familial LD networks were decomposed into five functional communities (Figs. 2A, 3A), which resembled the partitioning of the composite LD patient network (Fig. 1A). The left basal ganglia thalamic and cerebellar regions formed a separate community (module V, blue) in familial LD networks independent of phenotype. These findings suggest different mechanisms by which subcortical areas are functionally integrated in familial versus sporadic LD neural networks.
Influence of phenotype associations on the functional connectome of LD genotypes
Hub connectivity in the familial network with heterogeneous phenotype was characterized by a high number of connectors vs provincial relays (22/7, respectively), which was contrasted by only three connectors out of a total of 25 hubs in the respective sporadic network (Fig. 2B). Similarly, the mixed-phenotype familial network recruited 19 high-influence nodes, whereas the mixed-phenotype sporadic graph only contained eight high-influence regions. Despite the fact that neither strength nor degree values showed pronounced differences between the two putative LD genotypes (with heterogeneous phenotype), the formation of hubs exhibited distinct deviations between the networks. Specifically, module V (Fig. 2B, blue) was one of the smallest nodal communities detected within the familial graph (comprising 9% of all nodes), but it contained five connector hubs. In contrast, module II (green) represented a community of the same size but did not contain any hubs. Notably, four of five connectors found in module V were left thalamic regions (motor, parietal, somatosensory, and temporal subdivisions), three of which (motor, parietal, and somatosensory thalamic subdivisions) failed to show elevated nodal influence levels in the sporadic network (Fig. 2C).
In contrast to the mixed-phenotype familial network, the familial ADLD-only group showed a prevalence of provincial hubs with only a few connectors (connector/provincial: 4/25). This effect was also seen in the sporadic ADLD-only group, which exhibited only two connectors but 24 provincial hubs. Remarkably, module V (Fig. 3B, blue) played again a key role in the functional architecture of familial ADLD, containing two of four connectors (motor and parietal subdivisions of the left thalamus) and one provincial hub (somatosensory part of the left thalamus), while being the smallest nodal community in the network. Analogously to the mixed-phenotype networks, these three thalamic areas were not even considered high-influence nodes in the homogeneous sporadic group. Similarly, the familial ADLD network showed more high-influence nodes than the corresponding sporadic graph. These findings suggest that abnormal left thalamic function crucially shaped the familial LD connectome. Independent of phenotype, both familial networks showed an increase in high-influence nodes of >36% compared with the corresponding sporadic graphs (Figs. 2B, 3B).
Group 3: influence of phenotype on the LD functional connectome
Network segregation estimated by the clustering coefficient showed pronounced differences between the group-averaged ABLD and ADLD networks independent of LD genotype (mixed heterogeneous sporadic/familial genotypes vs homogeneous sporadic-only genotype: all p < 0.002). Deviations in network integration between ADLD and ABLD graphs, which was approximated by the global efficiency, were not statistically significant for the networks with mixed-genotype (p = 0.16), but exhibited a trend for sporadic-only graphs (p = 0.07). Both used nodal influence metrics were similar in the LD phenotype networks across examined genotypes (all p > 0.3; Table 2).
These findings suggest that the two groups were characterized by distinctly different functional network architectures, which was confirmed by the computed optimal modular decompositions. Although sporadic ABLD and ADLD networks gave rise to four functional communities, the spatial distribution of these modules showed pronounced differences and a partition distance of pd = 0.26. Notably, the modular architecture of the mixed-genotype ABLD network was characterized by the emergence of an additional fifth module and showed a higher degree of similarity to the mixed-genotype ADLD graph with a partition distance of pd = 0.2. The difference in ABLD network architecture across genotypes was initially quantified by calculating the partition distance between mixed-genotype and sporadic-only ABLD networks (pd = 0.2). In contrast, the modular makeup of the mixed vs homogeneous ADLD networks showed a high degree of spatial consistency and a partition distance of pd = 0.11.
The modular structure of both networks exhibited characteristics similar to the partitioning of healthy individuals with two lateralized nodal communities (modules I, red, and III, purple), a large central group (module IV, orange), and a smaller frontal component (module II, green) (Figs. 4A, 5A). In contrast, the partitioning of the sporadic ABLD network showed little similarity to any other group. Most notably, module I (red) in the sporadic ABLD network was not restricted to the left hemisphere, but also included bilateral parietal regions, and module II (green) did not only comprise bilateral frontal, but also included left basal ganglia and temporal areas (Fig. 4A). The remaining two communities of the sporadic ABLD network (modules III, purple, and IV, orange) were similar to the corresponding nodal groups in sporadic ADLD patients. However, many of the nodes that were part of module III (purple) in the sporadic ADLD network were members of module I (red) in the sporadic ABLD graph, thus shrinking this right-lateralized community in sporadic ABLD compared with sporadic ADLD.
Functional architecture of the group-averaged abductor LD and adductor LD networks with sporadic genotype. A, Spatial topology of detected functional communities in the group-averaged abductor LD (left) and adductor LD networks with sporadic genotype (right) is shown on 3D brain renderings in the axial and sagittal views, where regional module affiliation is illustrated by nodal color. B, Bivariate (with respect to strength and degree) provincial (yellow) and connector hubs (red) and their connectivity profile with high-influence nodes (purple; nodes with degree and strength in the top 25% of the network) in the group-averaged abductor LD (left) and in the mean adductor LD graph with sporadic genotype (right). C, Visualization of variations in degree (light gray) and strength (dark gray) across provincial (yellow) and connector hubs (red) in the group-averaged sporadic abductor and adductor LD networks (left and right, respectively). White circles indicate nodes that did not qualify as hubs in the respective graph. 6/17, area 6/17; 7A/7M/7P, subdivisions of area 7; CbI-I/IV/Cbl-V/Cbl-VI/Cbl-VIIIav/Cbl-VIv, cerebellar lobules I/IV/V/VI/VIIIa vermis/VI vermis; Cu/PCu, cuneus/precuneus; FG, fusiform gyrus; Ig1, part Ig1 of the insula; LG, lingual gyrus; MCC/PCC, middle/posterior cingulate cortex; PF, parietal area PF; SOG, superior orbital gyrus; Tp/Tpf/Ts/Tt, parietal/prefrontal/somatosensory/temporal subdivisions of the thalamus; hOC4v, ventral part of area hOC4; L, left; R, right.
Functional architecture of the group-averaged abductor LD and adductor LD networks with mixed genotype. A, Spatial topology of detected functional communities in the group-averaged abductor LD (left) and adductor LD networks with mixed genotype (right) is shown on 3D brain renderings in the axial and sagittal views, where regional module affiliation is illustrated by nodal color. B, Bivariate (with respect to strength and degree) provincial (yellow) and connector hubs (red) and their connectivity profile with high-influence nodes (purple; nodes with degree and strength in the top 25% of the network) in the group-averaged abductor LD (left) and in the mean adductor LD graph with mixed genotype (right). C, Visualization of variations in degree (light gray) and strength (dark gray) across provincial (yellow) and connector hubs (red) in the group-averaged mixed-genotype abductor and adductor LD networks (left and right, respectively). White circles indicate nodes that did not qualify as hubs in the respective graph. 6/17, area 6/17; 7A/7P, subdivsions of area 7; CbI-I/IV/Cbl-V/Cbl-VI, cerebellar lobules I/IV/V/VI; Cu/PCu, cuneus/precuneus; FG, fusiform gyrus; Ig1, part Ig1 of the insula; LG, lingual gyrus; MCC/PCC, middle/posterior cingulate cortex; SOG, superior orbital gyrus; Tp/Tpf/Ts/Tt, parietal/prefrontal/somatosensory/temporal subdivisions of the thalamus; hOC4v, ventral part of area hOC4; L, left; R, right.
Conversely, the most notable characteristic of ABLD with mixed genotype was the emergence of a fifth functional module (blue) that was centered on left thalamic and basal ganglia regions (Fig. 5A). Although this module was comparable to module V in the composite LD and mixed-phenotype familial LD networks (Figs. 1A, 2A), it was smaller in size and restricted to the left hemisphere. Interestingly, this left-lateralized community almost exclusively comprised areas that were part of module II in the sporadic ABLD graph (Fig. 4A), with the exception of a single node (fascia dentata subdivision of the left hippocampus). Aside from this thalamic/basal ganglia community, the mixed-genotype ABLD graph showed a modular architecture similar to mixed-genotype ADLD that was shaped by two large clearly lateralized frontotemporoparietal communities (modules I, red, and III, purple; Fig. 5A), a bilateral prefrontal module (module II, green), and a bilateral basal ganglia/cerebellar community (module IV, orange). However, most basal ganglia regions that were part of module IV in the mixed-genotype ADLD network migrated to modules V and III in the left and right hemisphere, respectively, in mixed-genotype ABLD. Therefore, module IV was much smaller in the ABLD network with mixed genotype, comprising only 29% of nodes in contrast to 37% in mixed-genotype ADLD.
In summary, the observed variations in the modular architecture of ABLD networks across genotypes suggest an interplay of ABLD phenotype with familial genotype, which appeared to be a crucial driver for community formation in ABLD but not ADLD patients. Conversely, the latter showed a remarkably stable community architecture across genotypes.
Influence of genotype associations on the functional connectome of LD phenotypes
Hub formation varied markedly between sporadic ADLD and sporadic ABLD networks. Sporadic ADLD patients exhibited predominantly provincial hubs (provincial/connector: 24/3), whereas ABLD patients were characterized by a prevalence of connectors (provincial/connector: 4/27) (Fig. 4B). Similarly, the number of high-influence nodes was very low in the sporadic ABLD compared with the corresponding ADLD network (ABLD/ADLD: 4/12). Most notably, the left insula (subdivision Ig1) and the right premotor cortex (area 6) represented two out of a total of three connector hubs in the sporadic ADLD network, but failed to show any pronounced increase in nodal influence in the corresponding ABLD network (Fig. 4B).
Comparing LD phenotype networks with mixed genotypes revealed a similar pattern. Mixed-genotype ABLD showed twice as many connector hubs as the corresponding ADLD network (8 vs 4), whereas the number of high-influence nodes was much larger in mixed-genotype ADLD than in the ABLD graph (15 vs 2; Fig. 5B). Similarly, the same areas (left insular subdivision Ig1 and right area 6) represented provincial and connector hubs, respectively, in the mixed-genotype ADLD network but were not even classified as a high-influence node (insula) or downgraded to a provincial hub (area 6) for ABLD with mixed genotype (Fig. 5C).
Together, these findings point to a markedly different functional connectome of ABLD compared with ADLD with the ABLD network showing a distinct functional architecture dominated by a strongly coupled connector hub network spanning several nodal communities.
Discussion
In this study, we defined the functional connectome of the laryngeal form of isolated focal dystonia by means of an in-depth graph theoretical analysis of functional networks during symptomatic speech production. In Group 1, we identified that, in general, the LD functional connectome is shaped by abnormal functional integration of basal ganglia, cerebellum, and thalamus, as well as by pronounced functional deficiency of parietal and primary sensorimotor cortices. In Group 2, we demonstrated a driving role of left thalamus in forging the functional architecture of familial LD. In Group 3, we showed that the hub connectome of ABLD was dominated by connectors, whereas a prevalence of provincial hubs was a characteristic feature of ADLD independent of LD genotype. The interplay of familial LD genotype with ADLD phenotype consistently resulted in five nodal communities by segregating left thalamic/basal ganglia regions, whereas the combination of sporadic genotype and ADLD phenotype always yielded four communities. Therefore, at the level of functional networks, the association of LD genotypes and phenotypes was reflected by the interplay of community formation and intramodule vs intermodule hub configuration.
Group 1: Overall characterization of the LD functional connectome
The connectivity profile in LD patients differed from that in healthy individuals on two levels. First, a clearly delineated left thalamic/basal ganglia community played a crucial role in the functional architecture of LD network; second, the majority of LD hubs were local relays. The additional fifth functional community in LD was centered on left basal ganglia and thalamus and only comprised 12% of network nodes, although it contained two of the total six connector hubs, which were both left thalamic regions. This indicates that, in a dystonic network, the left thalamus acted as a transmitter that relayed network-wide information to left basal ganglia forming a closely interacting community. This finding is consistent with the results of earlier study that illustrated the segregation of basal ganglia, thalamic and cerebellar resting-state communities across different focal dystonia (Battistella et al., 2017). Considering that both basal ganglia and cerebellar networks converge in the motor cortex (Bostan and Strick, 2010), our results suggest that pathophysiological integration of basal ganglia structures in the motor control network may propagate alterations via left thalamic structures throughout the entire network and likely contribute to the dystonic cascade. Together, these findings not only support previous reports of functional alterations in these regions in dystonia (Galardi et al., 1996; Neychev et al., 2008; Simonyan and Ludlow, 2012), but also reveal important details pertaining to the highly abnormal functional embedding of these key areas in the global large-scale brain network in LD.
A detailed assessment of hub structure showed a marked prevalence of local relays in the LD network, which was in stark contrast to a balanced distribution of hubs participating in long- and short-range edges in the healthy network. LD connectivity patterns exhibited a shift toward nodal cliques with strong intragroup connections but minimal network-wide relays, pointing to an abnormally segregated network architecture compared with healthy individuals. The latter was also reflected by significantly increased global clustering coefficient values in LD patients. Specifically, left primary motor and somatosensory cortices (areas 4a and 3b), which are known to drive the functional connectome of speech in healthy individuals (Fuertinger et al., 2015), failed to show significant nodal influence in the LD network. The loss of these areas as provincial and connector hubs points to severely impaired cortical sensorimotor processing and execution in LD. This also provides an explanation for abnormal sensorimotor alterations, particularly in task-specific focal dystonias (Bara-Jimenez et al., 2000; Fiorio et al., 2003; Tinazzi et al., 2006; Termsarasab et al., 2016). Similarly, left premotor cortex (area 6) was downgraded from a connector hub in the healthy network to a provincial relay in LD. Given its importance for sensorimotor integration, motor learning, and planning during speech production (Rauschecker, 2011; Jenson et al., 2014; Hardwick et al., 2015), altered function of this hub may additionally contribute to abnormal execution of motor commands in LD. Profound alterations of sensorimotor network components observed during selective impairment of speech production as a complex voluntary behavior point to a hierarchical mechanism underlying dystonic symptoms.
This hierarchy also pertains to the loss of pivotal regions of higher-order sensorimotor processing in the bilateral superior parietal cortex (areas 5L and 7A; Culham and Valyear, 2006; Brownsett and Wise, 2010; Shum et al., 2011; Sereno and Huang, 2014). Particularly, functional degradation of these areas may have further deteriorating effects, not only on general aspects of motor planning and execution, but also on coding of task-specific aspects of complex voluntary movement patterns (Rathelot et al., 2017) such as speech production, presumably via a direct access to the laryngeal motor cortex (Simonyan and Jürgens, 2002; Kumar et al., 2016). Conversely, increased functional influence of other parietal regions (left areas 7M and 7P, right area 5M), which were local functional relays in the dystonic network but were not even among the high-influence nodes in the healthy network, may point to a compensatory and abnormal upregulation of sensorimotor processing. Because recent studies suggested a potentially important role for parietal cortex in the pathophysiology of dystonia (Battistella et al., 2016; Gallea et al., 2016; Putzel et al., 2016a; Bianchi et al., 2017), our current findings provide more precise insights into abnormal functional embedding and specialized influences of the parietal cortex within the large-scale LD network.
Group 2: influence of genotype on the LD functional connectome
The functional architecture of familial LD networks was characterized by the formation of five nodal communities. Among these, a segregated community was centered on the basal ganglia and thalamus, where left thalamus acted as a global network relay station, functionally embedding the basal ganglia into the large-scale brain network. Notably, the parietal thalamic subdivision was found to be a shared connector hub across all familial networks independent of LD phenotype. This points to further abnormal integration of the parietal cortex via its subcortical connectivity.
In contrast to familial LD, basal ganglia and thalamus did not form a separate community in sporadic LD, but were rather integrated in a larger subcortical–cerebellar community. This network topology resembled that of the healthy network, possibly suggesting that the cerebellothalamocortical circuit was less impaired in the sporadic than familial LD. As another characteristic feature of the sporadic LD connectome, the hub distribution in both ADLD-only and mixed-phenotype sporadic LD networks was shaped by numerous provincial hubs and only a few connectors. This suggests that the pathophysiology underlying the combination of ADLD phenotype with sporadic genotype was characterized by a loss of network-wide hub relays in key areas of central motor control and sensorimotor processing.
Group 3: influence of phenotype on the LD functional connectome
Independent of genotype, the ADLD network exhibited a community structure similar to sporadic networks, whereas the ABLD network was shaped by highly distinct modular characteristics. Specifically, sporadic ABLD network was characterized by a large functional community that comprised 30% of all network nodes, including 10 of 27 connector hubs, and spanned bilateral superior parietal cortex. The latter was abnormally coupled with premotor, primary sensorimotor, and middle and posterior cingulate cortices. Altered functional connectivity profile of this nodal community suggests multilevel aberrations of sensory and motor transformations in the pathophysiology of dystonia due in part to substantial involvement of the parietal cortex as an important region of higher-level sensorimotor processing and integration.
Four connector hubs of the sporadic ABLD-only network, all of which were left thalamic areas (parietal, prefrontal, somatosensory, and temporal subdivisions), were embedded with basal ganglia and bilateral frontal regions in the same nodal community. This may suggest that the left thalamus plays a key role in sporadic ABLD hub connectome by establishing network-wide connections as part of frontothalamic nodal community. This network topology may underlie alterations of executive processing in LD such as abnormal sensory temporal discrimination (Termsarasab et al., 2016). It is also notable that the same thalamic regions (except for the prefrontal subdivision) were found to be connector hubs located in a segregated subcortical community of the mixed-genotype ABLD network, once again pointing to the function of this structure as a global relay station.
While the hub formation of sporadic ABLD network was dominated by connectors, sporadic ADLD was characterized by numerous provincial hub relays. The cerebellar frontoparietal community was the largest nodal group in the sporadic ADLD network that comprised 40% of all network nodes. It encompassed thalamus, cerebellum, cingulate, and parietal and frontal cortices, containing 16 of 24 provincial hubs. The formation of this large community, together with an overrepresentation of local intramodular hubs, might be associated with dystonia-related cortical hyperexcitability because the cerebellum has been shown to establish projections with association cortex in prefrontal, parietal, and superior temporal areas (Kelly and Strick, 2003). Although the mixed-genotype ABLD network exhibited more connectors than corresponding ADLD network, its hub connectome was primarily shaped by provincial relays. Therefore, abnormally high connector hub formation was not just driven by the ABLD phenotype; the interplay of sporadic genotype with ABLD phenotype yielded excessive connector hub recruitment.
Summary
On a network level, LD genotype and phenotype associations were reflected by the interplay of community organization and the formation of intramodule vs intermodule relay nodes. LD was characterized by abnormal global functional coupling of sensorimotor cortical areas, particularly primary sensorimotor and parietal cortices, as well as cerebellum, basal ganglia, and thalamus. Our results point to left thalamus as a key driver of functional network aberrations and abnormal connector hub recruitment, as well as to altered connectivity of parietal cortex as an important hallmark of LD connectome. The combination of sporadic genotype with ADLD phenotype was associated with a connectome that was mainly populated by provincial hubs, whereas familial genotype with ABLD phenotype yielded numerous long-range connector hubs spread out over functional communities. The connectome-wide approach taken in this study may serve as a blueprint for similar studies in other forms of focal dystonia for the discovery of pathophysiological mechanisms underlying this disorder that are both form-specific and shared across dystonias.
Footnotes
This work was supported by the National Institute on Deafness and Other Communication Disorders–National Institutes of Health (Grants R01DC01805 and R01NS088160 to K.S.). We thank Amanda Pechman, Melissa Choy, Ian Farwell, Estee Rubien-Thomas, and Heather Alexander for help with subject recruitment and imaging data acquisition; Andrew Blitzer, Miodrag Velickovic, and Steven Fruch for patient referrals, Laurie Ozelius for initial genetic testing, and Giovanni Battistella for help with image processing.
The authors declare no competing financial interests.
- Correspondence should be addressed to Kristina Simonyan, M.D., Ph.D., Department of Otolaryngology, Massachusetts Eye and Ear Infirmary, Harvard Medical School, 243 Charles St., Boston, MA 02114. dystonia.research{at}gmail.com
