Large-scale exploration of whole-brain structural connectivity in anorexia nervosa: alterations in the connectivity of frontal and subcortical networks

Background: Anorexia nervosa (AN) is characterized by disturbances in cognition and behavior surrounding eating and weight. The severity of AN combined with the absence of localized brain abnormalities suggests distributed, systemic underpinnings that may be identified using diffusion-weighted MRI (dMRI) and tractography to reconstruct white matter pathways. Methods: dMRI data acquired from female patients with AN (n = 147) and female healthy controls (HC; n = 119), aged 12–40 years, were combined across five studies. Probabilistic tractography was completed, and full cortex connectomes describing streamline counts between 84 brain regions generated and harmonized. Graph theory methods were used to describe alterations in network organization in AN. The network-based statistic tested between-group differences in brain subnetwork connectivity. The metrics strength and efficiency indexed the connectivity of brain regions (network nodes), and were compared between groups using multiple linear regression. Results: Individuals with AN, relative to HC, had reduced connectivity in a network comprising subcortical regions and greater connectivity between frontal cortical regions (p < 0.05, FWE corrected). Node-based analyses indicated reduced connectivity of the left hippocampus in patients relative to HC (p < 0.05, permutation corrected). Severity of illness, assessed by BMI, was associated with subcortical connectivity (p < 0.05, uncorrected). Conclusions: Analyses identified reduced structural connectivity of subcortical networks and regions, and stronger cortical network connectivity, amongst individuals with AN relative to HC. These findings are consistent with alterations in feeding, emotion and executive control circuits in AN, and may direct hypothesis-driven research into mechanisms of persistent restrictive eating behavior.


Abstract:
Anorexia nervosa (AN) is characterised by disturbances in cognition and behaviour surrounding eating and weight, which may relate to the structural connectivity of the brain that supports effective information processing and transfer. In the current investigation, diffusion-weighted MRI data acquired from female patients with AN (n = 148) and female healthy controls (HC; n = 119), aged 12-40 years, were combined across five cross-sectional studies. Probabilistic tractography was completed, and full cortex connectomes describing streamline counts between 84 brain regions generated and harmonised. The network-based statistic tested between-group differences in connectivity strength of brain subnetworks. Whole-brain connectivity of brain regions was indexed using graph theory tools, and compared between groups using multiple linear regression. Associations between structural connectivity variables that differed between groups, and illness severity markers, were explored amongst AN patients using multiple linear regression. Statistical models included age, motion, and study as covariates. The network-based statistic indicated AN patients, relative to HC, had reduced connectivity in a network comprising subcortical regions and greater connectivity between frontal cortical regions (p < 0.05, FWE corrected). Graph theory analyses supported reduced connectivity of subcortical regions, and greater connectivity of left occipital cortex, in patients relative to HC (p < 0.05, permutation corrected). Reduced subcortical network connectivity was associated with lower BMI among the AN group. Collectively, the findings suggest that structural differences in subcortical and cortical networks are present in AN, and may reflect illness mechanisms.

Background
Anorexia nervosa (AN) is characterised by abnormal behaviour and cognition surrounding eating and weight, including persistent restriction of food intake, fear of gaining weight or of becoming fat, and distorted self-image (1). The illness has high rates of mortality (2), physical and psychiatric morbidity (3,4), and severe effects on quality of life (5). The efficacy of current treatments for AN is limited, and relapse rates are high (6). Improvement in outcome is likely to depend on an improved understanding of the neural underpinnings of illness, including understanding the structural architecture of the brain in AN. To address existing gaps in knowledge, this study examined the structural (white matter) connectivity across the whole brain in AN. Diffusion-weighted MRI (dMRI) data from five studies were pooled to compare a large sample of individuals with AN and healthy control participants (HC). The brain operates as a large-scale system comprising densely connected and interacting brain regions. Co-ordinated patterns of activation, or functional connectivity, across spatially diffuse brain regions facilitates the information processing and integration that is necessary for behaviour and cognition across a range of domains, including visual memory and attention (7), language (8), and cognitive/inhibitory control (9). Structural, or white matter, connections support functional connectivity across the brain (10,11). Analysis of structural connectivity can be a fruitful approach to identifying functionally significant biological correlates of illness in AN.
White matter connectivity may be studied in vivo using dMRI (12). To date, most dMRI studies of AN have focused on the microstructural properties of white matter (13), which describe the quality of nerve cell axons but do not provide information about connectivity between regions of the brain. Studies that have compared connections between brain regions in individuals with AN relative to HC have typically focused on a small number of white matter tracts or connections (13,14). This approach could mean that important between-group differences are missed because they were not included in a priori hypotheses. As the neuropathological basis of AN is relatively unclear, additional data-driven approaches are warranted.
Whole-brain approaches explore structural connectivity across the entire brain. Pairwise structural connections between all regions of the brain can be represented in the structural connectome, a construct derived via mathematical graph theory. The structural connectome (15) contains "nodes" (brain regions) and "edges" (the white matter connections between brain regions). Using mathematical graph theory tools, regional, subnetwork, and whole-brain, connectivity properties can be quantified, to describe characteristics of the brain system that affect information transfer. There are few structural connectome studies of AN. Those that do exist are generally small (16), though the single large study (including 96 AN and 96 HC; (17)) reported increased white matter integrity within parietal and occipital subnetworks. Existing structural connectome investigations in AN populations have typically used a tract reconstruction method that estimates a single fibre orientation within each voxel. This approach can lead to misrepresentation of the white matter tracts, since white matter often contains complex fibre architectures such as bending or crossing fibres (18,19). More accurate reconstructions are obtained by modelling multiple fibres per voxel, as probabilistic tractography algorithms do (20).
Probabilistic methods have the additional advantage of accounting for the noise inherent in dMRI data, and consequent uncertainty in the tract reconstruction process.
The aim of the current study was to leverage advances in dMRI methods, including the application of a probabilistic algorithm for tract reconstruction, to explore brain-wide structural connectivity in a large sample of patients with AN and HC. In addition to considering the strength of white matter connections between brain regions, we probed differences in how connected brain regions are with the rest of the brain (i.e., node connectivity). This approach offers a window into brain organisation differences that may contribute to psychopathology in AN.

Participants and Methods
Analysis of dMRI data collected from 267 participants in five different cross-sectional studies (14,(21)(22)(23)(24) was undertaken. Data from a subsample of participants from one of these studies (22 AN, 18 HC) were included in a prior publication comparing specific frontostriatal tracts between AN and HC (14); other dMRI data (for 122 AN, 101 HC) are unpublished.

Participants
Participants were 148 females who met DSM-5 criteria for AN, and 119 female HC. All participants were between the ages of 12 and 40 years. Exclusion criteria across all studies were: presence of psychotic disorder, current substance use disorder, history of a neurological disorder/injury, estimated IQ less than 80, and MRI contraindications. Additional exclusion criteria for HC included current/past psychiatric illness, significant medical conditions, and body mass index (BMI) outside of the 18.5-25 kg/m 2 range. Aside from 11 participants, all were free from psychotropic medication use at the time of MRI scanning. Medication included selective serotonin reuptake inhibitors (n = 1), serotonin and noradrenaline reuptake inhibitors (n = 1), and tricyclic antidepressants (n = 1). Eating disorder diagnoses were made by the Eating Disorder Examination (25) or Eating Disorder Diagnostic Assessment for DSM-5 (26). Co-occurring diagnoses were assessed (for adults) by Structured Clinical Interview for DSM-IV or 5 (27,28), and (for adolescents) by the Schedule for Affective Disorders and Schizophrenia for School-Aged Children (29). Height and weight were measured on the day of the MRI scan for calculation of BMI. Each study was approved by the New York State Psychiatric Institute (NYSPI) Institutional Review Board. Researchers provided a complete description of the study to participants, following which adult participants provided written informed consent, and adolescents provided written assent with parental consent.
Amongst the AN patients, 85 (57.4%) were diagnosed with the restricting subtype. All adult AN participants, and a subset of adolescent AN participants, received inpatient or outpatient treatment at NYSPI. HC were recruited from the community.

Scanning parameters
The MRI scanner and DWI scanning parameters were study-specific and are described in the Supplement (Table S1).

MRI Image processing
Anatomical images Structural images were processed using the automated Freesurfer pipeline (30,31), which includes motion correction, brain extraction, Talairach transformation and intensity normalisation. Grey matter regions were parsed into 84 distinct regions based on the Desikan-Killany atlas (32).
The parcellated image was used to generate a five-tissue type image, in which white matter, cortical and subcortical grey matter, cerebral spinal fluid (CSF), and pathological tissue are distinguished. Anatomical images were warped into diffusion space using a rigid body transformation performed with FSL's flirt method (50, 51), followed by a non-linear normalization using ANTS (52,53), and resampled such that the voxel size was 1.25 mm 3 .

Diffusion images
Each scan was visually checked for artifacts prior to image pre-processing, and volumes containing artifacts were removed. When more than 25% of the total volumes were removed, the participant was excluded (n = 3 AN, 0 HC). A fully automated image pre-processing pipeline implemented in MRtrix3 (33) was applied to the data, including: denoising, gibbs-artifact removal, B0 susceptibility-induced distortion correction, eddy current and motion correction (using the FSL eddy tool; (34)), and B1 field inhomogeneity correction (using the ANTS N4 algorithm; (35)). When a pair of B0 images with opposite phase encoding was available, the inhomogeneity field was estimated and applied to the DWI data using the eddy topup and apply topup tools (36) during the eddy current/motion correction phase. Where fieldmap images were available but B0 image pairs were not, distortion correction was applied prior to the processing pipeline by applying a representative affine transformation to all DWI volumes using FSL's epi_reg command (37,38). The preprocessed dMRI data were up-sampled to an isotropic voxel size of 1.25 mm 3 . Motion was quantified with average framewise displacement across the dMRI scan and did not differ between AN and HC participants: HC=0. 16 (Table S2).

Probabilistic Tractography
Anatomically-constrained probabilistic tractography was completed through application of a standard MRtrix3 pipeline (33). First, response functions of white matter, grey matter and CSF were derived with the dhollander algorithm (39), using the dMRI signal and the five-tissue type image generated from the structural MRI scan. Next, the response functions of white matter fibres were estimated in the voxels with the 300 greatest fractional anisotropy values (indicative of more directed diffusion). Using the dMRI signal and the calculated response functions, the underlying distribution of fibres in each voxel was recovered using multi-shell multi-tissue constrained spherical deconvolution (40), with default maximum spherical harmonic degree parameters. This approach allows dMRI signal resulting from white matter and CSF components to be distinguished, preventing CSF contaminating estimates of white matter fibre orientation distributions. Consequently, the influence of partial volume effects (e.g., due to brain atrophy in AN) on tract reconstruction was minimised. White matter fibre orientation distributions were normalised, to remove residual inhomogeneity in intensities. Seed points for streamlines were determined dynamically. Seeding continued until 20,000,000 streamlines between 10 and 250 mm were generated. Streamline step size was 0.625 mm; the maximum angle between successive steps was 22.5 degrees.
The iFOD2 tracking algorithm (41) determines the probability of candidate streamline paths by sampling fibre orientation distributions along the entire path. More probable, and thus reconstructed, streamlines are generally those where fibre orientation distribution amplitudes along the path are large. The tractography was anatomically constrained, with the five-tissue type image used as a biological prior to guide streamline termination and acceptance (42). Streamlines were re-tracked if a poor termination was encountered, and cropped at the grey matter white matter interface.

Graph construction
Graphs can provide a mathematical representation of the brain network, or structural connectome. To create brain graphs, a network node was assigned to each of the 84 regions from the anatomical parcellation. For each participant, the number of streamlines between all possible pairs of regions (Ni, Nj) was counted to provide edge values, generating an undirected, weighted, graph that quantified connectivity strength within each tract of the brain. Streamline endpoints were mapped to nodes within a 4mm radius of the endpoint. When no streamlines connected a pair or regions, the edge value was set to 0.
Individual graphs (or connectomes) were harmonised using the R implementation of neuroCombat (43,44), to remove effects of differences in scanner/scan parameters across and within studies. To retain effects of diagnosis, age and BMI, these variables were specified as covariates in the harmonisation model. In exploratory analyses probing subtype differences in structural connectivity, the harmonisation model was applied to data from AN patients only, and included subtype as an additional covariate. Following recent recommendations graphs were not thresholded prior to completing network and graph theory analyses (45).

Data analysis
Streamline counts between brain regions comprised the measure of connectivity in the current study.
Network-based statistic (46) Edge-wise connectivity was compared between AN and HC groups using the network-based statistic method implemented in Matlab (https://sites.google.com/site/bctnet/comparison/nbs; (47)). To complete the network-based statistic analysis, a t-test was performed to test group differences at each connection, or edge, to create an 84 x 84 matrix of t-scores. The t-score matrix was thresholded to identify connections with a between-group difference t-score greater than 3.1, which corresponds to a p value of 0.001. This relatively stringent threshold reduces the risk of detecting false positive effects. Sets of supra-threshold edges that are connected (components) are identified. Permutation testing derives p-values for a test of the hypothesis that the size of components differing in connectivity strength between groups is larger than what would be expected if there was no between group difference: for each permutation, the group to which each participant belongs is randomly exchanged, the test statistic computed for group differences at each edge, and the size of the maximal supra-threshold component stored. The null distribution of the maximum component size can be then empirically estimated and used to compute a p-value for the observed components (i.e., those identified in the unpermuted data).
10,000 permutations were specified, and identified components were significant at a family-wise error corrected rate of p < 0.05. The network-based statistical analysis was completed for the contrasts HC > AN, and AN > HC; p-values were Bonferroni corrected for two tests. Models were adjusted for age, motion and study. Study was included as a covariate given the potential for residual batch effects to remain following harmonisation.

Node-based analyses
Two measures of regional network topology, nodal strength and efficiency, were computed for each network node (or brain region), using the BCT toolbox (https://sites.google.com/site/bctnet/) in Matlab (47,48). Strength is the total number of connections (streamlines) of each node, and was calculated as the sum of the node's streamline connections. Efficiency describes the ease at which information can move from one region to another within a network of the node's neighbours ("ease" here refers to fewer intermediary nodes). To calculate efficiency, the average inverse of the shortest path length between all nodes directly connected to the node in question was determined. Path length was defined as the inverse of streamline count. The contribution of a given path to the efficiency calculations was proportionate to the streamline count between the path nodes and the node for which efficiency was calculated (i.e. the average was weighted). Global efficiency, the ease at which information can move from one region to any other region on average across the whole brain, was also calculated.
Node characteristics were compared between AN and HC in linear regression models adjusted for age, motion and study. Permutation tests established statistical significance and corrected for multiple comparisons across the 84 nodes (with the minP method), using the R package Permuco (49). 10,000 permutations for each test were specified; the Freedman-Lane procedure managed effects of nuisance variables (age, motion and study). A linear regression model adjusted for age, motion and study compared AN and HC groups on the global efficiency metric.
The tractography, connectome generation, connectome metric calculation and analysis procedures are detailed in Figure 1.

Figure 1: Neuroimaging and analysis methods of the current study
For each participant, a diffusion MRI and structural MRI scan were acquired (panel A). White matter tracts in the brain were reconstructed using a probabilistic tracking algorithm applied to diffusion MRI data, and the structural brain image was parcellated into 84 regions (panel B).

Combining the tractography output with the brain parcellation allowed the creation of a structural connectome for each participant, which described the number of streamlines between each pair of regions (panel C). This information was represented in a weighted adjacency matrix, which was used to calculate node metrics (panel D). Anorexia nervosa and healthy control groups were compared on the strength of structural connections and node metrics, across the whole brain (panel E).
The main analyses probing group differences in graph theory metrics and subnetwork connectivity were repeated using the unharmonised data. This was to verify whether findings were consistent when different methods were used to combine data across the five studies, which would support the robustness of study conclusions.

Association of structural connectome metrics and AN severity
Network components and node metrics that significantly differed between groups (p < 0.05, corrected) were tested for associations with BMI and illness duration within the AN group only.
For each component reaching statistical significance, the mean number of streamlines across edges of these components was calculated. To combine BMI for adults and BMI percentile for adolescents, scaled scores (mean centred) were created in each group, using BMI for individuals over 19, and BMI percentiles for individuals 19 and under. Component connectivity metrics were regressed onto BMI and illness duration, in separate models, in the patient sample. Using the same method, node metrics that differed between groups were tested for their association with clinical variables amongst individuals with AN. All statistical models were adjusted for motion, age and study. Models involving the regression of illness duration were additionally adjusted for the scaled BMI variable. Model predictors and outcomes (aside from the already scaled BMI variable) were mean centred and scaled, for comparison of effect estimates across exploratory analyses.

Subtype differences in structural connectome metrics
Subtype differences in component connectivity and node metrics were explored by regressing outcomes that differed between AN and HC groups onto subtype, in linear regression models adjusted for motion, age, study, and the scaled BMI variable.

Code and data availability
The code for dMRI data preprocessing, tractography, and analysis, steps is available on github: https://github.com/CaitlinLloyd/Anorexia.Nervosa_Tractography. Data for between group differences in streamline counts between all pairs of regions are available as Supplementary material (see Supplementary Tables file), for harmonised and unharmonised data.

Participant characteristics
Participant characteristics are displayed in Table 1. For a summary of results by study and diagnosis, see the Supplement (Table S2).      (Tables S4 and S5).

Unharmonised data analyses
Analyses with unharmonised data yielded outcomes that were consistent (in direction, and generally size and significance) with outcomes of harmonised data analyses. Full results are reported in the Supplement (Tables S7-S9).

Association of structural connectome metrics and AN severity
Results of the exploratory analyses are presented in Table 3. Amongst AN, the scaled BMI  Table 2: Outcomes of linear regression analyses exploring associations of illness severity markers with identified structural connectome differences amongst individuals with AN

Subtype differences
Individuals with AN binge-purge subtype had greater strength of the right cerebellum relative to individuals with AN restricting subtype. There were no other significant differences between AN subtypes in terms of structural connectivity metrics (Supplementary material, Table S10).

Discussion
The current study harmonised data across multiple studies to reconstruct anatomical connections of the brain using probabilistic tractography, in a large sample of individuals with AN and HC.
Analyses yielded a comprehensive description of the structural connectome in AN, including subnetwork and regional connectivity. The network-based statistic indicated that, relative to HC, individuals with AN had reduced connectivity within a component comprising connections between structures located subcortically in the brain. In contrast, patients displayed enhanced connectivity in a cluster of tracts connecting regions of the frontal cortex, as compared to HC.
The node-based analysis findings were consistent with network-based statistic outcomes: individuals with AN had reduced whole-brain connectivity of regions within the subcortical component. Amongst individuals with AN, lower BMI was associated with reduced connectivity within the subcortical component, highlighting an association between structural brain connectivity and illness severity in AN.
While this study did not measure behaviour directly, the finding of heightened connectivity in the frontal lobe along with weaker connectivity between subcortical regions is consistent with several proposed neurobiological models of AN. Heightened cortical connectivity, particularly between inferior frontal, lateral orbitofrontal, and anterior cingulate regions implicated in cognitive control, may relate to the highly inflexible behaviour observed in patients (50).
Weakening within subcortical circuits may relate to decreased appetitive drive (51). These alterations may together support the persistent restriction of food intake that is the central behavioural characteristic of AN (52). The structural differences present amongst individuals with AN may underpin the alterations reported in functional MRI studies. In particular, reduced connectivity between subcortical regions (e.g., amygdala, thalamus), and altered connectivity within networks comprising frontal and cingulate regions, has been reported amongst individuals with AN, relative to HC, at rest (53). Altered activation within frontal brain regions during cognitive tasks is also reported amongst individuals with AN, as compared to HC (54).
Our findings are consistent with some earlier dMRI studies probing white matter properties that have reported individuals with AN, as compared with HC, to have greater integrity within frontal portions of the inferior fronto-occipital fasciculus (55,56), and reduced integrity of white matter in the fornix (13). The inferior fronto-occipital fasciculus traverses the inferior frontal and ventromedial regions in which elevated connectivity was seen amongst individuals with AN in the current study. The fornix is the major output tract of the hippocampus and encompasses white matter within the subcortical component in which we identified reduced connectivity amongst individuals with AN (e.g., connections between hippocampus, thalamus, parahippocampal cortex and entorhinal cortex).
The current findings do not precisely replicate prior investigations of region-to-region white matter connectivity in AN. Several investigations identified stronger connectivity in intrahemispheric frontostriatal tracts amongst individuals with AN, both between the nucleus accumbens and OFC (14,57), and between the putamen and motor areas (58). Studies observing elevated frontostriatal connectivity in AN patients adopted a region of interest approach, whereby the unit of analysis was single connections between brain regions. The current study tested differences at the level of subnetworks (or constellations of connected white matter connections), meaning inferences about differences within certain tracts cannot be made.
Different analytical approaches could explain inconsistencies in findings relating to connectivity alterations within frontostriatal circuits amongst individuals with AN.
One large-scale tractography study of AN applied a brain-wide approach and probed networklevel characteristics. This study reported greater integrity of white matter in occipital-parietal regions amongst individuals with AN, though no between-group differences in streamline counts (17). The earlier study used a tract reconstruction method based on the presence of a single fibre direction per voxel, which may have contributed to the discrepancy between outcomes of the earlier study and outcomes of the current study. Single-fibre methods produce different (and less biologically accurate) reconstructions of white matter pathways relative to the multi-fibre probabilistic approach implemented herein (18)(19)(20). Replication attempts that implement similar tract reconstruction and analytical methods to those of the current study are necessary to evaluate the reliability and validity of inferences concerning brain organisation in AN that arise from this investigation.
Several limitations of this study should be considered. First, the data are cross-sectional, meaning causes and consequences of starvation cannot be parsed. Longitudinal studies may provide insight into whether the structural connectivity associated with AN in this study predicts illness onset/maintenance, or normalises with weight gain. Second, dMRI data were not optimal for constrained spherical deconvolution-based tractography. Ideally, dMRI scans would have been acquired with higher b-values (e.g., 3000), for more precise estimation of fibre orientation distributions within a given voxel. However even at low b-values, constrained spherical deconvolution-based tractography outperforms methods assuming a single-fibre orientation within each voxel, in terms of producing an anatomically accurate reconstruction of white matter pathways (59). Third, data harmonisation does not completely remove batch effects, and collecting all dMRI data in one study would have been ideal.
In summary, this large-scale investigation of whole-brain structural connectivity in patients with AN implemented advances in diffusion tractography and graph theory to yield information about the neural architecture associated with AN. Individuals with AN had reduced connectivity of and between subcortical regions, and elevated connectivity within tracts connecting frontal regions.
The identified connectivity profile hints at weaker processing within motivational networks, and greater processing within networks supporting cognitive control, the combination of which may drive restrictive eating behaviour. Future investigations using improved dMRI data acquisition parameters should seek to confirm our findings, and, using longitudinal designs, elucidate the nature of association between structural connectivity features and AN pathology.

Author Statements
Funding and disclosures Funding The study was supported by funding from the following sources: Global Foundation for Eating Disorders; Klarman Family Foundation; Translating Duke Health Initiative; NIMH (MH099388, MH076195, MH110445, MH105452, MH079397, MH113737). Disclosures Dr. Steinglass reports receiving royalties from UpTo Date. Dr. Posner has received support from Takeda (formerly Shire), Aevi Genomics, and Innovation Sciences. Dr Lloyd, Dr Foerde, Ms Aw, Mr Semanek, and Dr Muratore have no conflicts of interest to report.

Author contributions
ECL developed the analysis plan, and completed MRI preprocessing, diffusion tractography, data aggregation/harmonisation and data analysis, and drafted the manuscript. KEF refined the analysis plan and manuscript drafts. NA developed the diffusion imaging preprocessing and tractography pipelines for the current study and assisted with diffusion MRI preprocessing and diffusion tractography. DS assisted with preprocessing of the structural MRI data. AFM contributed to data aggregation and harmonisation. JES and JP designed the study and refined the analysis plan and manuscript drafts. All authors approved the final submission.  White matter tracts in the brain were reconstructed using a probabilistic tracking algorithm applied to diffusion MRI data, and the structural brain image was parcellated into 84 regions (panel B). Combining the tractography output with the brain parcellation allowed the creation of a structural connectome for each participant, which described the number of streamlines between each pair of regions (panel C). This information was represented in a weighted adjacency matrix, which was used to calculate node metrics (panel D). Anorexia nervosa and healthy control groups were compared on the strength of structural connections and node metrics, across the whole brain (panel E).