Age Trajectories of the Structural Connectome in Child and Adolescent Offspring of Individuals With Bipolar Disorder or Schizophrenia

Background Offspring of parents with severe mental illness (e.g., bipolar disorder or schizophrenia) are at elevated risk of developing psychiatric illness owing to both genetic predisposition and increased burden of environmental stress. Emerging evidence indicates a disruption of brain network connectivity in young offspring of patients with bipolar disorder and schizophrenia, but the age trajectories of these brain networks in this high-familial-risk population remain to be elucidated. Methods A total of 271 T1-weighted and diffusion-weighted scans were obtained from 174 offspring of at least 1 parent diagnosed with bipolar disorder (n = 74) or schizophrenia (n = 51) and offspring of parents without severe mental illness (n = 49). The age range was 8 to 23 years; 97 offspring underwent 2 scans. Anatomical brain networks were reconstructed into structural connectivity matrices. Network analysis was performed to investigate anatomical brain connectivity. Results Offspring of parents with schizophrenia had differential trajectories of connectivity strength and clustering compared with offspring of parents with bipolar disorder and parents without severe mental illness, of global efficiency compared with offspring of parents without severe mental illness, and of local connectivity compared with offspring of parents with bipolar disorder. Conclusions The findings of this study suggest that familial high risk of schizophrenia is related to deviations in age trajectories of global structural connectome properties and local connectivity strength.

entire sample (N = 286), to determine whether there are group differences in how many times each of the two scanners was used, and whether the effect of age on which scanner was used differs between groups, we ran a generalized linear model as follows: 'Scanner ~ Group + Age + Group x Age', which showed that there are no significant (age-related) group differences in which of the two scanners was used (all p's > 0.117) (Supplemental Table S5).

Recruitment process and exclusion criteria
Offspring at high familial risk were recruited via adult psychiatrists and child-and-adolescent psychiatrists in the Netherlands, advertisements or hear-say.Control families were recruited via advertisements on schools and leisure clubs or hear-say.Exclusion criteria for all participants were severe physical illness, neurological problems (e.g., open or closed head injury, epilepsy) or an IQ below 70 and, for controls only, a first-degree relative with a severe mood or psychotic disorder.

Instruments
At both time points, participants were psychiatrically evaluated using the Schedule for Affective Disorders and Schizophrenia for School-Age Children-Present and Lifetime Version (1) by interviewing participants and their parents separately.The summary information of this face-to-face semi-structured interview was used to establish the presence and severity of symptoms and current and lifetime DSM-IV axis I diagnoses.IQ was estimated using four subtests (block design, picture completion, information and vocabulary) of the Wechsler Intelligence Scale for Children-III (≤16 years) (2) or the Wechsler Adult Intelligence Scale-III (>16 years) (3).

MRI preprocessing
All imaging data were coded to ensure blinding for participant identification and diagnoses during image processing.Imaging data were (pre)processed and analysed on SURFsara's Snellius compute cluster.Before migration to the compute cluster, T1-weighted images were defaced using mri_deface version 1.22 (4).Visual quality control was performed by S.P. to determine whether defacing was performed successfully (i.e., the participant's face was made unrecognisable with the brain left intact).
FreeSurfer (v7.1.1)(5) and the FMRIB Software Library (FSL, v6.0.6) (6) were used to process and segment the T1-weighted image and to perform preprocessing on the DWI images.T1-weighted image processing included automated Talairach transformation, intensity normalisation, removal of non-brain tissue, segmentation of subcortical white matter and grey matter, tessellation of the grey/white matter boundary, and automated topology correction (7)(8)(9)(10).Visual quality control was done according to the ENIGMA procedures (http://enigma.ini.usc.edu/).DWI data were corrected for susceptibility-induced distortions (11), eddy current distortions and motion artefacts (12), outliers (13), intra-volume movement (14) and susceptibility-by-movement (15).The b-vectors were updated to adjust for DWI corrections (16).Visual quality control of the processed DWI data was performed by S.P. and M.B., where each researcher thoroughly inspected each volume of half of the scans, checking residual motion, scanner and metal-induced artefacts, and brain coverage.In case of uncertainty about whether a scan should be excluded or not, the other researcher was included in the process and a decision was reached through consensus.Visual quality control was always done blind to the group the participant was in.After our visual data quality control, we re-evaluated the scan quality of the scans where participants showed a steep change in network metrics.These scans were all of acceptable quality and the steep changes could not be explained by changes in scan quality (e.g., having one borderline scan and one of really good quality).We also did not see group differences in how often participants show a steep change (defined as two standard deviations away from the mean) in connectivity strength (Co: n = 1; BDo: n = 1; SZo: n =1).Therefore, we did not see a reason to exclude these data.
The Connectivity Analysis TOolbox (CATO) (v3.2.1) (17) was used to reconstruct structural connectivity from the processed DWI data.A reference image was computed based on the corrected b=0 volumes.Next, the registration matrix between the DWI reference image and the T1 image was created using bbregister (18).The FreeSurfer segmentation files were registered to the DWI reference image for anatomical alignment.Diffusion directions were determined in each voxel using the informed RESTORE algorithm to reduce the impact of physiological noise artefacts on the DTI modelling (19,20).Deterministic tractography was performed to reconstruct the white matter pathways, using an extended version of the Fiber Assignment by Continuous Tracking algorithm (21).
While advanced reconstruction methods such as those based on crossing-fiber models and probabilistic tractography allow for more complete reconstructions of complex fibre pathways and high sensitivity (22), deterministic tractography was chosen to obtain an adequate false negative to false positive fibre reconstruction ratio, which has been shown to have a substantial impact on connectome analyses (23)(24)(25).Fibre reconstruction was initiated from eight seeds in each cerebral white matter voxel.For each seed, a tractography streamline was propagated in the mean diffusion direction from voxel to voxel.
Fibre tracking was terminated when (1) FA<0.1, (2) a>45˚ turn was made, (3) a stop region was entered, (4) a forbidden region was about to be entered, (5) the current or previously visited voxel was about to be revisited or (6) the brain mask was exited.

Structural connectome topology
A selection of global characteristic graph metrics was computed to investigate possible differences in the development of overall connectome topology between the three groups (26) (Figure 1E).The connectivity strength of a brain network was quantified as the sum of the weights of all connections in that network.Global efficiency, representing the network's overall communication capacity, was computed as the harmonic mean of the shortest path length (i.e., fewest number of edges that must be traversed for one node to connect to another) between all pairs of nodes in the network.As a measure of segregation, the clustering coefficient was calculated as the average likelihood that a node's neighbours are also connected among themselves.The modularity represents the degree to which a network can be broken down into community structures (i.e., cliques of nodes with the highest possible number of within-group connections and the lowest possible number of between-groups connections) (27).

Demographic and clinical characteristics
Group differences in demographic and clinical characteristics were evaluated at both time points using Fisher's exact test for categorical variables (i.e., clinical diagnosis, psychotropic medication use and sex) and analyses of variance for continuous variables (i.e., age, IQ and scan interval) followed by Tukey's test in case of significant group effects.

Confounding variables
Several potentially confounding variables are added to our (sensitivity) analyses.Sex was included because there are known differences between boys and girls in age-trajectories of structural brain development (28)(29)(30)(31).IQ, medication use and psychopathology are confounders of the familial risk of mental illness: IQ is known to be lower, and medication use and psychopathology are known to be higher in the high-familial-risk groups (32)(33)(34).Therefore, we applied sensitivity analyses with these variables, as it allows us to examine whether they explain the group differences, i.e., it enables us to compare group differences before and after adding the covariate.

Statistical assumptions
Normality of the model residuals was visually inspected using histograms and Q-Q plots.Additionally, model residuals were tested for skewness, kurtosis and normality using the D'Agostino, Anscombe-Glynn test and Shapiro-Wilk tests, respectively.Based on these plots and statistics, we concluded that the model residuals are sufficiently normally distributed.

Time point 1 connectome analyses
Analyses to compare the graph theoretical metrics between groups at time point 1, as was done previously in a subgroup of the current sample (35), were repeated on the metrics obtained from the NOS-weighted networks of the time point 1 sample.Specifically, a linear mixed model analysis was performed with group, age and sex as fixed effects and family as a random effect.In contrast to our previous analyses, in which participants were only scanned in one scanner, scanner was now added as a fixed effect.Also, we now selected rich club nodes a priori instead of defining them for each subject individually.The resulting formula was as follows: 'lmer(NetworkMetric ~ Group + Age + Sex + Scanner + (1|FamilyID))'.

Effects of head motion
In addition to comprehensive preprocessing and rigorous visual quality control of the DWI and T1weighted scans, we investigated potential effects of head motion statistically.For each DWI dataset, we used the average of the root mean square (RMS) volume-to-volume movement values, accounting for the eddy current component, taken from the my_eddy_output.eddy_restricted_movements_rmsoutput file (see: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide).To see whether there was an effect of age and whether this differed per group, we ran an analysis of variance (ANOVA) as follows: 'AverageRMSmovement ~ Group + Age + Group x Age'.There was no group difference (Group: p = 0.238) and there was a strong effect of age (Age: p < .001),but this did not differ significantly between groups (Group x Age: p = 0.457).
For each T1-weighted scan, we calculated the Euler number for each hemisphere using the surface holes given by FreeSurfer's recon-all, and then took the sum of the Euler numbers of both hemispheres for a total Euler number per scan.To see whether there was an effect of age and whether this differed per group, we ran an ANOVA as follows: 'Euler ~ Group + Age + Group x Age'.There was no significant group difference (Group: p = 0.691) and there was a strong effect of age (Age: p < .001),but this did not differ significantly between groups (Group x Age: p = 0.857).
Next, we added both the AverageRMSmovement variable and the Euler number to the main analyses' lmer() models as fixed effects to evaluate their impact on each of the seven outcome measures, resulting in the following formula: 'lmer(NetworkMetric ~ Group + Age + Group x Age + Sex + Scanner + AverageRMSmovement + Euler + (1|FamilyID/SubjectID))'.Adding these two parameters did not change the nature of our findings (Supplemental Table S7).

Rich club organization
Rich club organization was investigated by normalising the weighted rich club coefficient (36) relative to a set of 10.000 comparable random networks (26,37,38).Specifically, a weighted group-averaged network was reconstructed by constructing a binary network of all connections identified in at least 50% of the participants (participants with scans at both time points and therefore two network matrices had their network matrices averaged first).The weight of each connection in the group-averaged network was averaged across subjects.Next, a weighted rich club curve was computed for the empirical group-network and compared to the average rich club curve of a set of 10.000 comparable randomised networks comprising the same connection weights (26,37,38).Then, a normalised rich club curve was calculated as the ratio between the empirical and random rich club curve (35).The presence of rich club organization is indicated where the normalised rich club curve significantly exceeded 1 over a range of degree k.
To confirm consistent rich club membership across groups, a separate group-averaged network was computed for each group as well.The top 20 high-degree nodes of each group's group-averaged network were compared.Additionally, hubs with a degree of 19 or higher of each group's groupaveraged brain network were compared.

Time point 1 connectome analyses
Repeating analyses on only those scanned at time point 1 largely aligned with our previous findings (36,37), despite now including individuals scanned at two different scanners and selecting rich club nodes a priori instead of defining them for each subject individually.We find no significant effects between groups on global network metrics (i.e., connectivity strength, clustering coefficient, global efficiency or modularity), local connectivity or feeder connectivity.Now, SZo showed reduced rich club connectivity compared to Co and BDo, but these differences did not reach statistical significance (p = 0.290 and p = 0.052, respectively) (Supplemental Table S6), while Collin et al. (33) reported significantly reduced rich club connectivity in SZo as compared with Co (p = 0.018) and BDo (p = 0.010).

Rich Club Organization
Rich club organization was verified in the current cohort at k=15, k=16 and over range k=19 to k=23 (Supplemental Figure S3).The top 20 high-degree nodes and the nodes with a degree of 19 or higher of each group's group-averaged network were examined to compare rich club membership between groups, revealing it to be highly consistent between groups (Supplemental Table S11 and S12).Across groups (with edges present in at least 50% of all subjects in the respective group, after averaging the network matrices of those with scans at both time points), the top 20 high-degree nodes included parts of the bilateral inferior parietal gyrus, posterior cingulate gyrus, precuneus, rostral middle frontal gyrus and superior parietal gyrus, the left superior frontal gyrus and superior temporal gyrus and the right insula, isthmus of the cingulate gyrus and postcentral gyrus.In the group connectome based on all participants (with edges present in at least 50% of all subjects, after averaging the network matrices of those with scans at both time points), the top 20 high-degree nodes comprised all of these regions and, additionally, subregions of the left lateral orbitofrontal gyrus and the right lateral occipital gyrus and superior temporal gyrus (Supplemental Table S11).Rich club, feeder and local connectivity analyses with hub selection based on each group's averaged brain network are shown in Supplemental Table S13.Running these analyses with hub selection based on the group connectome of all participants yielded highly similar results (Supplemental Table S13).

FA-weighted networks
Running the analyses on network measures computed on FA-weighted networks yielded a significant age effect in Co in connectivity strength (p<0.001) and global efficiency (p=0.001),but not in clustering (p=0.066) and modularity (p=0.302).There were no significant group differences in age effects in the global metrics computed on the FA-weighted networks (p≥0.072,Supplemental Table S8).
Running the analyses on network measures computed on FA-weighted networks showed a significant age effect in Co in rich club, feeder and local connectivity (all p's<0.001).There were no significant group differences in age effects between the three measures (all p's<0.138).S1 S11.Top 20 high-degree nodes per group based on each group's average network.Top 20 high-degree nodes in the group-averaged network of all participants are in italic.Supplemental Table S13 (34,35).

Supplemental Figure S1 .
Correlation matrices of the seven investigated NOS-weighted network metrics.(top) Entire dataset (values of participants who were scanned at both time points were averaged first); (bottom) Datasets at time point 1 (left) and time point 2 (right).Supplemental Figure S3.Normalised rich club curve (relative to random rich club curve based on average 10000 random networks) for the weighted group-averaged network of all participants.A normalised coefficient Φ greater than 1 over a range of k suggests the existence of rich club organization in a network

.
Exclusion reason per group per time point.
a Reasons include braces, piercings or metal in the body.

Table S4 .
Number of offspring per group scanned with each scanner at each time point.Generalized linear model were run with scanner as outcome variable and group, age and group*age as predictor variables.Analyses run with bipolar disorder offspring as reference group.
BDo = bipolar disorder offspring, Co = control offspring, SZo = schizophrenia offspring.a These numbers concern the 41 BDo, 23 SZo and 33 Co that were scanned at both time points.a b Analyses run with control offspring as reference group.c

Table S6 .
Time point 1 connectome analyses for NOS-weighted global network metrics.Co = controls, NOS = number of streamlines, SZo = schizophrenia offspring.a Linear mixed model analyses were run for each network metric with group, age, sex and scanner as fixed effects, and family as random effect.Analyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group. b

Table S7 .
Main effect of Age (representing the effect of age on network metric in Co) and Group*Age interaction effects for NOS-weighted global network metrics, corrected for head motion during the DWI and T1-weighted scans.Linear mixed model analyses were run for each network metric with group, age, group*age, sex, scanner, average RMS movement and Euler number as fixed effects, and subject and family as random effects.
a b Analyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group.

Table S8 .
Main effect of Age (representing the effect of age on network metric in Co) and Group*Age interaction effects for FA-weighted global network metrics.Linear mixed model analyses were run for each network metric with group, age, group*age, sex and scanner as fixed effects, and subject and family as random effects.
a b Analyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group.

Table S9 .
Principal component analysis loadings (correlations between the z-scores of the NOSweighted variables and the principal components) and results.

with those with two scans averaged first, n = 174)
Supplemental

Table S10 .
Effects of group, sex and scanner in the main linear mixed-effects model analyses for NOSweighted global network metrics.Co = controls, NOS = number of streamlines, SZo = schizophrenia offspring.a Linear mixed model analyses were run for each network metric with group, age, group*age, sex and scanner as fixed effects, and subject and family as random effects.
bAnalyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group.

.
Main effect of Age (representing the effect of age on network metric in Co) and Group*Age interaction effects for NOS-weighted rich club, feeder and local connectivity with top 20 highdegree nodes selected separately for each group based on their group-averaged network, or based on the total sample's group-averaged network (same selection for each group).

20 nodes based on the group-averaged network of the total sample
Linear mixed model analyses were run for each network metric with group, age, group*age, sex and scanner as fixed effects, and subject and family as random effects.b Analyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group. a

Table S15 .
Main effect of Age (representing the effect of age on network metric in Co) and Group*Age interaction effects for NOSweighted global network metrics with IQ added as covariate BDo = bipolar disorder offspring, Co = controls, NOS = number of streamlines, SZo = schizophrenia offspring.a Linear mixed model analyses were run for each network metric with group, age, group*age, sex, scanner and IQ as fixed effects, and subject and family as random effects.
bAnalyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group.

Table S16 .
Main effect of Age (representing the effect of age on network metric in Co) and Group*Age interaction effects for NOSweighted global network metrics with psychotropic medication use (yes/no) added as covariate disorder offspring, Co = controls, NOS = number of streamlines, SZo = schizophrenia offspring.a Linear mixed model analyses were run for each network metric with group, age, group*age, sex, scanner and psychotropic medication use (yes/no) as fixed effects, and subject and family as random effects.
bAnalyses run with control offspring as reference group.c Analyses run with bipolar disorder offspring as reference group.