Disrupted structural network of inferomedial temporal regions in relapsing–remitting multiple sclerosis compared with neuromyelitis optica spectrum disorder

Multiple sclerosis (MS) and neuromyelitis optica spectrum disorder (NMOSD) are two representative chronic inflammatory demyelinating disorders of the central nervous system. We aimed to determine and compare the alterations of white matter (WM) connectivity between MS, NMOSD, and healthy controls (HC). This study included 68 patients with relapsing–remitting MS, 50 with NMOSD, and 26 HC. A network-based statistics method was used to assess disrupted patterns in WM networks. Topological characteristics of the three groups were compared and their associations with clinical parameters were examined. WM network analysis indicated that the MS and NMOSD groups had lower total strength, clustering coefficient, global efficiency, and local efficiency and had longer characteristic path length than HC, but there were no differences between the MS and NMOSD groups. At the nodal level, the MS group had more brain regions with altered network topologies than did the NMOSD group when compared with the HC group. Network alterations were correlated with Expanded Disability Status Scale score and disease duration in both MS and NMOSD groups. Two distinct subnetworks that characterized the disease groups were also identified. When compared with NMOSD, the most discriminative connectivity changes in MS were located between the thalamus, hippocampus, parahippocampal gyrus, amygdala, fusiform gyrus, and inferior and superior temporal gyri. In conclusion, MS patients had greater network dysfunction compared to NMOSD and altered short connections within the thalamus and inferomedial temporal regions were relatively spared in NMOSD compared with MS.


Results
summarizes the demographics and clinical features of RRMS (N = 68), NMOSD (N = 50), and 26 HC. The NMOSD group was older than the MS and HC groups (p < 0.001 and p = 0.004, respectively), but there were no differences in sex. The disabilities of patients were also more severe in the NMOSD group than in the MS group (p < 0.001). Brain attacks were more prevalent in MS than in NMOSD (cerebral hemisphere: 100% vs 44%, p < 0.001; brainstem/cerebellum: 65% vs. 46%, p = 0.043). Neither patient group had gadolinium enhancement in their cerebral lesions. Each MS patient had brain lesions in more than one MS-typical locations, such as adjacent to lateral ventricles body, in the inferior temporal lobes or subcortical U-fibers. Among patients with NMOSD, 56% (N = 28) of patients had more than one NMOSD-typical brain lesions. The imaging findings of NMOSD patients in our study were summarized in Supplementary Table S1. Supplementary Figure S1 shows the brain lesion probability maps for the MS and NMOSD groups. MS lesions were more often immediately adjacent to Table 1. Demographics and clinical characteristics of study subjects. MS multiple sclerosis, NMOSD neuromyelitis optica spectrum disorder, HC healthy control, EDSS expanded disability status scale, IQR interquartile range, NA not applicable. a p value from chi-squared test. b p value from ANOVA, posthoc tests results: MS vs. NMOSD (p < 0.001), HC vs. NMOSD (p = 0.004), and MS vs. HC (p > 0.05). c Symptomatic involvement. d Taken at the time of brain MRI. e Interferon β-1b (n = 24), interferon β-1a (n = 11), teriflunomide (n = 8), azathioprine (n = 5), fingolimod (n = 3), glatiramer acetate (n = 3), dimethyl fumarate (n = 2), mitoxantrone (n = 1). f Azathioprine (n = 29), mycophenolate mofetil (n = 12), rituximab (n = 2), cyclosporine (n = 1), methotrexate (n = 1). Disrupted brain network topology in MS and NMOSD. Overall, the MS and NMOSD groups had lower total strength, clustering coefficient, global efficiency, and local efficiency and had longer characteristic path length (CPL) than HC, but there were no differences between the MS and NMOSD groups (Table 2). Disease duration or Expanded Disability Status Scale (EDSS) score in the MS and NMOSD groups had positive correlation with CPL and negative correlation with the other global properties of network (Table 3). Since the longer CPL is coincident with the lower global efficiency, all showed the same trend. When compared with the HC group, brain regions with altered nodal characteristics were more widespread in the MS group than in the NMOSD group. The left medial orbital region of superior frontal gyrus had lower nodal degree and nodal strength and the left superior parietal gyrus, and precuneus had lower nodal strength in the MS and NMOSD groups compared with the HC group. The left fusiform gyrus and the right superior occipital gyrus in MS also had lower nodal strength compared with HC. Furthermore, significant differences were indicated among the nodal measures of several brain regions between the MS and NMOSD groups. The local efficiencies of the left hippocampus, left parahippocampal gyrus, left superior temporal gyrus, and right Heschl's gyrus, and the regional efficiency of the left middle frontal gyrus were significantly lower in MS than in NMOSD. The brain regions with significant between-group differences within nodal measures are listed in Supplementary Table S2. The correlation results between nodal measures of the brain regions and clinical parameters were listed in Supplementary Table S3.
Disrupted subnetworks in MS and NMOSD. NBS identified two disrupted subnetworks (threshold = 3.0) among the MS, NMOSD, and HC groups (Table 4, Fig. 1A, Supplementary Fig. S2). The first subnetwork (subnetwork 1, p < 0.001) consisted of regions mostly located in the left temporo-parieto-occipital lobes, although the precuneus, cuneus, superior occipital gyrus, and calcarine cortex had bilateral involvement. Posthoc tests indicated that compared with HC, the MS group had significant decreases in most edges of subnetwork 1 (28 out of 29 edges); however, in NMOSD only about half of the edges (13 out of 29 edges) were significantly decreased (Fig. 1B,C). Some edges were more disrupted in the MS group than in the NMOSD group, which were the connections between the left hippocampus and fusiform gyrus, between the fusiform gyrus and amygdala, between the superior temporal and inferior temporal gyri, and between the right precuneus and paracentral lobule (Fig. 1D). The second subnetwork (subnetwork 2, p = 0.028) was only disrupted in the MS group and was generally located in the inferomedial temporal region and consisted of five connections between the hip-  www.nature.com/scientificreports/ pocampus, parahippocampal gyrus, fusiform gyrus, and thalamus in the right hemisphere. Compared with the NMOSD group, most edges (four out of five edges) were significantly disrupted in the MS group. We also identified that several edges of the disrupted subnetworks were associated with disease duration or EDSS score in the MS and NMOSD groups (Supplementary Table S4).

Discussion
Disrupted global topological organization of WM networks was demonstrated in both MS and NMOSD, with no differences between these two conditions. However, at the nodal level, the MS group had more brain regions with altered metrics than did the NMOSD group when compared with the HC group. Furthermore, the left hippocampus, left parahippocampal gyrus, left superior temporal gyrus, and right Heschl's gyrus of the MS group had greater decreases in local efficiencies compared with the NMOSD group. Two distinct subnetworks were identified to characterize the disease groups. Disconnected edges were more widespread in the MS group compared with the NMOSD group, and discriminative connectivity changes were mostly found in the thalamus and inferomedial temporal regions. Network alterations were associated with EDSS score and disease duration in both groups. www.nature.com/scientificreports/  www.nature.com/scientificreports/ Disorganized WM network in MS and NMOSD. The whole-brain WM networks of MS and NMOSD were characterized by reduced total strength, global efficiency, and local efficiency compared with HC. These changes were also correlated with higher EDSS scores and longer disease durations in both disease groups, which is consistent with previous studies 11,15 . Similar patterns of change in structural connectivity have been reported in MS 11,16,17 and NMOSD 12,18 . Our previous study reported that patients with NMOSD only had reduced total strength compared with HC, which was likely due to including fewer patients with cerebral lesions (14% vs 44% in this study) 10 .
More disrupted regions with altered network topology in MS than NMOSD. Brain regions with altered nodal topologies were more widely distributed in MS than in NMOSD. Decreased nodal degree or nodal strength with connection losses was indicated for MS in the left fusiform gyrus and several other brain hub regions, such as the left medial orbital part of the superior frontal gyrus, precuneus, and superior parietal gyrus, and the right superior occipital gyrus 19 . Decreased nodal degree or strength was also indicated for NMOSD in the left medial orbital part of the superior frontal gyrus, precuneus, and the right superior occipital gyrus. Many of the connections to and from hub regions were long-range association fibers, which are more susceptible to WM damage and disconnection from brain disorders 14 . Medial orbital part of the superior frontal gyrus receives various sensory input and is associated with reward, mood and emotion; the lower connectivity with parahippocampal gyrus/medial temporal lobe might contribute to reduced happy memories and anhedonia 20 . Precuneus, superior parietal and superior occipital gyri are functional and structural hubs important in global communication and associated with memory, visuospatial function, or object recognition 21 . The fusiform gyrus is the largest component of the ventral temporal cortex and is also a conduit for long association fibers such as the inferior longitudinal fasciculus and inferior frontal occipital fasciculus 22 . Dysfunction and atrophy of the fusiform gyrus occurred in MS 23,24 . Another study using graph theoretical analysis and diffusion MRI found that MS was associated with a larger number of nodes with reduced nodal strength (26 out of 84 nodes) compared with HC, more than what the present study indicated (5 out of 90 nodes) 17 . This may be due to our patients having shorter disease durations with less brain damage, or differences in analysis methods (e.g., construction of WM tractography and the number of comparison groups). Similar to the present study, previous studies have also found that MS was associated with a larger number of disrupted regions than NMOSD 11,12 and similar brain regions with decreased efficiency were found in MS 11 and NMOSD 12 compared with HC. However, these studies did not directly compare structural networks between MS and NMOSD. Therefore, it was particularly interesting that the differences between the MS and NMOSD patients of the present study were primarily in the local efficiencies of several temporal regions; the left hippocampus, left parahippocampal gyrus, left superior temporal gyrus, and right Heschl's gyrus exhibited significantly smaller local efficiencies in MS than in NMOSD. These regions had a role in memory, visuospatial/auditory processing, language and social cognition 25,26 and were susceptible to microstructural damage or atrophy in MS even compared with NMOSD 27,28 Previously, MS patients exhibited severe cognitive impairment, especially verbal and visual memory, compared with NMOSD patients 29 , which might be associated with left superior temporal gyrus volume loss 30 .
Disrupted subnetworks within the thalamus and inferomedial temporal regions differentiating MS from NMOSD. Two subnetworks that characterized the structural connectivity abnormalities in the MS and NMOSD groups were also identified. Subnetwork 1 consisted of broad regions in the temporoparieto-occipital lobes (mostly on the left side) and encompassed a part of default-mode, visual/visuospatial, and memory systems. Disrupted edges included decreased strength nodes, which were the left fusiform gyrus, left precuneus, left superior parietal gyrus, and the right superior occipital gyrus. Several disrupted edges, including the connection between the right and left precuneus, were significantly associated with EDSS score or disease duration in both disease groups. Subnetwork 2 consisted of five edges between the hippocampus, parahippocampal gyrus, fusiform gyrus, or thalamus in the right hemisphere, which were the right-side counterparts of those in subnetwork 1. The weights of most edges (four out of five) were negatively associated with the duration of MS. A previous study of WM structural networks in MS found decreased efficiency in brain regions related with default-mode, visual, memory, and language function 11 . Within these subnetworks, all connections except for one were disrupted in the MS group. However, in the NMOSD group there was relatively less disruption of the short connection between the thalamus and inferomedial temporal regions, including edges between the thalamus, hippocampus, parahippocampal gyrus, amygdala, fusiform gyrus, and inferior and superior temporal gyri. These regions overlapped with brain regions that had decreased local efficiencies in MS than in NMOSD. These findings suggest that the pathological alterations in the corresponding GM regions and their surroundings were greater in MS than NMOSD, which was supported by a previous report that the GM volume of the left parahippocampal, superior temporal gyri, and the right hippocampus was smaller in MS than in NMOSD, while GM atrophy was diffuse in MS compared to HC 28 . In RRMS patients, both cortical and deep GM atrophies were reportedly associated with the disrupted integrity of the connected WM tracts, with the associations being strongest in the temporal lobes such as the hippocampus, parahippocampal gyrus, fusiform gyrus, and lateral temporal regions 31 . The present study supports previous reports that demyelination and atrophy of the thalamus and hippocampus was evident in MS patients 1 , and lesions in the inferior temporal lobe that were common in MS were helpful to be distinguished from NMOSD 3 .
There were several limitations in this study. First, brain network alterations in NMOSD might vary between patients since brain involvement does not occur in all patients. Second, the duration of MS could affect the network comparison results since edge weights between the hippocampus, parahippocampal gyrus, fusiform gyrus, and thalamus-which differentiated MS from NMOSD-were negatively correlated with MS duration. Third, inclusion of patients within 3 months after myelitis relapse (11% and 28% of MS and NMOSD patients, www.nature.com/scientificreports/ respectively) may affect the correlation between EDSS score and network measures. Fourth, the HC group was not exactly age-matched with the patient groups, especially the NMOSD group, which may confound the comparison results. Fifth, the tractography resulted in a longer z-direction, which may have affected its overall quality since the voxel dimensions of our DTI protocol were anisotropic (1.72 mm × 1.72 mm × 3 mm). Finally, deterministic tractography is subject to a major inherent limitation of the crossing-branching problem. The present study compared connectome-level differences in structural networks between MS, NMOSD, and HC. WM network disruption was more widespread and severe in MS than in NMOSD. Significant disconnections within inferomedial temporal regions might be a differentiated feature in MS from NMOSD.

Methods
From the prospective cohort of CNS demyelinating disease at Samsung Medical Center in Seoul, South Korea, patients with relapsing-remitting MS (RRMS, N = 68) and NMOSD (N = 50), who visited the clinic between January 2014 and December 2018, were enrolled by retrospective chart review. Inclusion criteria were patients who (1) over the age of 18 years, (2) met the revised 2017 McDonald criteria for RRMS 32 or revised 2015 NMOSD diagnostic criteria 33 , and (3) had analyzable brain MRI performed at least 3 months after relapse in the brain. The NMOSD patients included 45 (90%) with AQP4-IgG, which was measured with a cell-based indirect immunofluorescence assay as previously described 34 . AQP-IgG-negative NMOSD patients were also seronegative for myelin oligodendrocyte glycoprotein-IgG, which was determined by in-house live cell-based immunofluorescence assay using an anti-human IgG1-Fc secondary antibody, as previously described 35 . At the time of selection, 8 (11%) MS patients and 14 (28%) NMOSD patients had experienced a spinal cord attack within the previous 3 months, while the rest were all in remission. We also included 26 HC who had no history of medical, neurological, or psychiatric disorders. Of those, 21 people were recruited from a previous study 10 and the rest were prospectively enrolled. This study was approved by the Institutional Review Board of Samsung Medical Center (IRB No. 2020-04-119), and written informed consent was obtained from all participants. All procedures were performed in accordance with relevant guidelines and regulations.
Image acquisition. All participants underwent a three-dimensional volumetric brain MRI scan. An Achieva 3.0-Tesla MRI scanner (Philips, Best, the Netherlands) was used to acquire 3D T1 Turbo Field Echo MRI data using a sagittal slice thickness of 1.0 mm, over contiguous slices with 50% overlap and no gap, a repetition time (TR) of 9.9 ms, an echo time (TE) of 4.6 ms, a flip angle of 8° and matrix size of 240 × 240 pixels reconstructed to 480 × 480 over a field of view of 240 mm. 3D fluid attenuated inversion recovery (FLAIR) MRI data were acquired in the axial plane with the following parameters: axial slice thickness of 1 mm, no gap; TR 11,000 ms; TE 125 ms; flip angle 90°; and matrix size of 512 × 512 pixels. In the whole-brain DTI, sets of axial diffusion-weighted single-shot echo-planar images were collected with the following parameters: 128 × 128 acquisition matrix; 1.72 × 1.72 × mm 3 voxels; 70 axial slices; 220 × 220 mm 2 field of view; TR 7696 ms, TE 60 ms; flip angle 90°; slice gap 0 mm; b-factor of 600 s/mm 2 . With the baseline image without diffusion weighting (the reference volume), DTI were acquired from 45 different directions. All axial sections were acquired parallel to the anterior commissure-posterior commissure line and perpendicular to the mid-sagittal plane.
Image preprocessing and network construction. We defined brain regions as network nodes using the automated anatomical labeling atlas 36 , which consists of 78 cortical and 12 subcortical regions, excluding cerebellum regions. We registered them onto the DTI space of each subject using Statistical Parametric Mapping software (version 12, SPM12) 37 . The overall procedure is briefly introduced as follows. We first obtained deformation fields between the DTI space and the tissue probability map (TPM) space through segmentation procedure. In the segmentation procedure, we used average sized TPM template and mutual information-based affine registration. Using the obtained deformation fields, we normalized the DTI of each subject onto the TPM space. After we co-registered the AAL atlas to the TPM space, the AAL atlas, transformed in the TPM space, was projected onto the standard space by applying the inverse of the deformation fields obtained in the segmentation procedure; we matched its voxel size with the original DTI by reslicing it. In the normalization and coregistration procedure of the AAL atlas, we used the nearest neighborhood interpolation method that led clearer boundaries of regions. As a result, we obtained the registered AAL atlas in each subject's DTI.
We defined network edges by obtaining the averages of fractional anisotropy (FA) values of all voxels on the streamlines connecting any two brain regions. We first performed eddy-current correction of FSL's Diffusion Toolkit (version 3.0) to adjust unwanted movements by registering all volumes with gradient directions to the reference volume. The gradient directions were appropriately rotated during this alignment procedure 38 . We performed tractography using the Fiber Assignment by Continuous Tracking algorithm with 45 degrees as the angular threshold determined using Diffusion Toolkit with TrackVis 39,40 on the movement adjusted DTI 38 . We then collected all the voxels on the streamline connection paths and obtained the averages of their FA values. WM voxels where FA values of the seed voxels exceeded 0.2 were included and streamlines shorter than 20 mm were excluded. As a result, we developed a 90-by-90 connectivity matrix consisting of 90 brain regions whose edge weights were the mean FA values between all pairs of brain regions.
Network topology. We quantified global and nodal properties of network using the Brain Connectivity Toolbox (http:// sites. google. com/ site/ bctnet) 41 . We measured degree, strength, clustering coefficient, characteristic path length (CPL), local efficiency, global efficiency, and regional efficiency. www.nature.com/scientificreports/ Degree. Degree is the number of neighboring nodes linked to a node. When an edge between the ith and jth node exists, a ij = 1 ; otherwise, a ij = 0. The degree of the ith node, K i , is computed by summing all values of a ij between the ith node and others, where N is the number of the nodes in network.
Strength. Strength at the nodal level, nodal strength, is computed by summing all edge weights linked to a node. The nodal strength of the ith node, S i , is defined by the summation of w ij , where w ij is the edge weight between the ith and jth node. Strength at the global level, total strength, is computed by summing all edge weights in whole network.
Clustering coefficient. Clustering coefficient measures how well nodes are clustered in a network using the number of triangles between those nodes 42 . In a weighted network, the total intensity of triangles is used instead of the number of triangles. The total intensity of triangles of the ith node, t i , is computed by the summation of the cubic root of products of all edge weights for all connected triangles, where j and k are the indices of its neighboring nodes 43 .
where w ij , w jk , and w ki indicated edge weights between two neighboring nodes linked to the ith node. The clustering coefficient at the nodal level, nodal clustering coefficient is defined by dividing the total intensity of the ith node by the number of possible triangles linked to ith node, that is, The clustering coefficient at the global level is defined by averaging of the nodal clustering coefficients over all nodes.
Characteristic path length (CPL). CPL is defined as the average of shortest path lengths for all pair of nodes in a network 42 . We transformed path length between two nodes as reciprocal of the edge weight connecting them, because the edge weight in our network represents the strength of connection between two nodes. Then, we used the Dijkstra algorithm to find shortest path length between any pair of nodes. If the edge between two nodes is disconnected, because the reciprocal of edge weight is infinite, we ignored them in averaging.
where d ij is the shortest path length between the ith and jth nodes, and D is the number of finite d ij .
Local and global efficiency. Local efficiency is defined as the average efficiency of local subgraphs where the local subgraph is the set of neighboring nodes linked to a certain node excluding the centered node 44 . It captures the efficient communication between neighbors of a certain node. We defined local efficiency of nodal level by averaging them 44 . The local efficiency at the global level is defined as the average of the nodal local efficiency of all the nodes 44 .
Global efficiency is defined by the average of the reciprocal of shortest path lengths between all pairs of nodes 44 .
Regional efficiency. Regional efficiency defined as the average of reciprocal of the shortest path length between a certain node and all the other nodes in a network. It shows how well a node communicates with all the others.
In weighted network, path length between two nodes represents a reciprocal of edge weight between them. Edge screening. Since the connectivity matrices obtained through tractography may include false-positive and false-negative connections, their effects were controlled using a method proposed by de Reus and van den Heuvel 45 . We first computed prevalence rate of each connection. The prevalence rate is a ratio of the number of subjects for which a certain connection exists within the HC group. As an example, if all the subjects in the HC group have a certain connection, the prevalence rate of the connection is one. We identified whether a certain connection actually exists using a group threshold which is set with a value between 0 and 1. As an example, if the prevalence rate of a certain connection is bigger than a chosen group threshold, the connection is considered as an existing connection. Otherwise, it is considered as a non-existing connection. After deciding which connection is an existing connection, we defined the false positive and false negative as follows: if a non-existing edge exists in a subject, it is a false positive; if an existing edge does not exist in a subject, it is a false negative. As a chosen group threshold increases, since the number of existing connections increases, the number of false positives decreases, and that of the false negatives increases. We counted the number of false positives and false negatives for various group thresholds and chose a balancing point which minimizes both the number of the false positives and false negatives as the optimal group threshold value, which was 0.66 in our data set. Connections that were weaker than the optimal group threshold were not considered connections.

Network based statistics.
To identify a subnetwork that consisted of significant differences in WM connections between the MS, NMOSD, and HC groups, NBS 46 was performed with analysis of covariance (ANCOVA) while controlling for the effects of age and sex. Within NBS, cluster-based thresholding of statistical maps was performed to control the familywise error rate during mass univariate testing of every network connection. We thresholded the test statistics from ANCOVA, and the connected edges formed subnetworks (clusters). NBS also used permutation testing to estimate the significance levels of group differences according to the number of larger subnetworks when compared with randomly formed subnetworks. NBS therefore extracted subnetworks with significantly different connections between the three groups. Thresholds were selected when stable clusters form over multiple runs. The threshold was 3.0 and 10,000 permutations were performed.
Lesion probability map. Lesions from the FLAIR images were automatically segmented using the lesion prediction algorithm 47 implemented using the Lesion Segmentation Toolbox (LST) (version 3.0.0, www. stati stical-model ing. de/ lst. html) for SPM12 (https:// www. fil. ion. ucl. ac. uk/ spm/). This estimated the lesion probability of each voxel using a logistic regression model that was trained using the data of 53 MS patients with severe lesion patterns. We applied this model to each voxel within the FLAIR images to segment the lesions and estimate the lesion probability. Voxels with lesion probabilities below 50% were discarded and binarized to develop a lesion map image for every subject. We then registered the lesion map images to the Montreal Neurological Institute standard-space template. The averages of each patient group were calculated to help develop the group lesion probability map.
Statistical analyses. Network measures were compared using permutation-based ANCOVA 48,49 which controlled for the effects of age and sex. Because 90 nodes received group comparisons for nodal measures, multiple comparison corrections were performed using the false discovery rate (FDR) 50,51 .
Since ANCOVA does not determine differences between specific groups, post-hoc tests were performed using three pairwise permutation-based ANCOVA, and then FDR was used to perform multiple comparison corrections on these three pairwise comparisons. Identical post-hoc tests were also performed on the subnetwork edges obtained from NBS. N was set as 10,000 for permutation-based ANCOVA and number of permutations.
We also investigated the association between topological network characteristics and clinical variables using the generalized linear model. We controlled for the effects of age and sex: where each beta value is the regression coefficient of each term. Significance levels of the regression coefficient of network measures ( β 1 ) were determined through ANCOVA.