The unbalanced reorganization of weaker functional connections induces the altered brain network topology in schizophrenia

Network neuroscience shed some light on the functional and structural modifications occurring to the brain associated with the phenomenology of schizophrenia. In particular, resting-state functional networks have helped our understanding of the illness by highlighting the global and local alterations within the cerebral organization. We investigated the robustness of the brain functional architecture in 44 medicated schizophrenic patients and 40 healthy comparators through an advanced network analysis of resting-state functional magnetic resonance imaging data. The networks in patients showed more resistance to disconnection than in healthy controls, with an evident discrepancy between the two groups in the node degree distribution computed along a percolation process. Despite a substantial similarity of the basal functional organization between the two groups, the expected hierarchy of healthy brains' modular organization is crumbled in schizophrenia, showing a peculiar arrangement of the functional connections, characterized by several topologically equivalent backbones. Thus, the manifold nature of the functional organization’s basal scheme, together with its altered hierarchical modularity, may be crucial in the pathogenesis of schizophrenia. This result fits the disconnection hypothesis that describes schizophrenia as a brain disorder characterized by an abnormal functional integration among brain regions.

specific macroscopic networks (e.g. default mode, salience, somatosensory, etc.) or the role of single brain regions as possible hubs in the whole brain web and their coordination with the first neighbours.
In order to overwhelm such limitations, in this paper, we investigated the functional connectivity of brain regions by addressing their hierarchical integration at multiple scales (from the single region level to clusters and entire subnetworks level) within the whole functional network 39 in a cohort of patients diagnosed with SCZ compared to healthy controls (HC). Based on a varied multilevel pipeline, the approach we propose has significant advantages over previously used methods of data-driven brain functional connectivity investigation, in terms of sensitivity of small changes due to small synchronization shifts. We employed recently developed 39 advanced methods of network neuroscience 40,41 to investigate the cerebral region-to-region interaction both in SCZ and HC and the possible dysfunctional reshuffling of connections across the whole network. Briefly, to test the network resilience to disconnection, we ran a percolation analysis by following the disconnection process from the progressive removal of links from the network and the maximum spanning tree (MST) filtration by keeping the most robust links discarding all the others. Finally, we applied an allometric scaling to the MST to quantify its divergence from a linear tree-like organization.
We hypothesize that the altered topology of functional brain networks in SCZ is a consequence of an unbalanced distribution of large and small connectivity strength values. This abnormal connectivity organization includes areas belonging to systems highly engaged in the same processing roles (i.e., areas functionally segregated) and others with different processing assignments (i.e., areas globally integrated) 28 . A possible mechanism underpinning SCZ may arise from a not optimized topological arrangement of the functional network due to the altered connectivity distribution.
We further hypothesize that the cerebral region-to-region interaction is more resistant to disconnection in patients than in healthy subjects due to a dysfunctional reshuffling of significant connections across the whole network. Consequently, the loss of higher cortical hierarchies, responsible for predicting representations, should be a prodromal condition for the reduced specialization of the patients' brain functional connectivity patterns.

Methods
Participants. Forty-four patients diagnosed with SCZ according to the DSM-5 42 were recruited. The clinicians who had been treating the patients and knew their clinical history made the preliminary diagnosis. A senior research psychiatrist (GS) confirmed all prior diagnoses using the Structured Clinical Interview for DSM-5-research version (SCID-5 for DSM-5, Research Version; SCID-5-RV) 43 . Other inclusion criteria were: (1) age between 18 and 65 years; (2) at least eight years of education; (3) no dementia or cognitive deterioration according to the DSM-5 criteria, and a Mini-Mental State Examination score 44 higher than 24 and (4) suitability for a Magnetic Resonance Imaging (MRI) scan. Exclusion criteria were: (1) a history of alcohol or drug dependence or abuse in the last 2 years; (2) traumatic head injury; (3) any past or present major medical or neurological illness; (4) any other psychiatric disorder or mental retardation diagnosis and (5) MRI evidence of focal parenchymal abnormalities or cerebrovascular diseases. All patients were in a stable clinical compensation phase and were receiving stable oral doses of one or more atypical antipsychotic drugs. Forty HC were also recruited. They did not differ with patients in terms of age, education, and gender. HC were screened for a current or lifetime history of DSM-5 psychiatric and personality disorders using the SCID-5-RV 43 and SCID-5-PD 45 , were right-handed, and with normal or corrected-to-normal vision. Table 1 shows HC and SCZ patients' principal demographic characteristics and the associated statistics (two-sample t test).
Participants gave written informed consent to partake in the study after the procedures had been fully explained. The Ethics Committee of the IRCCS Santa Lucia Foundation approved the study protocol which was conducted in accordance with the Helsinki Declaration 46 .
We also acquired high-resolution T1-weighted whole-brain structural scans (1 × 1 × 1 mm voxels). Subjects were asked to keep their eyes open, and their cardiac and respiratory cycles were recorded using the scanner's built-in photoplethysmograph and a pneumatic chest belt, respectively. For each subject and each time-series, we removed possible sources of physiological variance: time-lock cardiac and respiratory artifacts through linear regression 47 (i.e., two cardiac and respiratory harmonics, respectively, together with four interaction terms). We also looked for the effect of low-frequency respiratory and heart rates [48][49][50] . The pre-processing of fMRI data consisted of: head-motion and slice timing corrections plus the discard of voxels not belonging to the brain (with FSL: FMRIB's Software Library, http:// www. fmrib. ox. ac. uk/ fsl). We used head motion parameters estimation to obtain the Framewise Displacement (FD). The subjects were included in this study only if they presented a www.nature.com/scientificreports/ Framewise Displacement smaller than 0.5 mm; however, time points with high FD (FD > 0.2 mm) were replaced using a least-squares spectral decomposition following 51 . Then, we detrended, demeaned, and band-pass filtered (frequency range 0.01-0.1 Hz) data using custom Matlab algorithms. We performed a 2-step registration in line with group-analysis.
We first transformed fMRI data from the functional space to the structural space with FLIRT (FMRIB's Linear Registration Tool) for each subject. Using Advanced Normalization Tools (ANTs; Penn Image Computing and Science Lab, http:// www. picsl. upenn. edu/ ANTS/), the data were subsequently non-linearly sent to the standard space (Montreal Neurological Institute MNI152 standard map).
A spatial smooth was applied to the data (5 × 5 × 5 mm full-width half-maximum Gaussian kernel).
Network analysis. The whole network analysis was characterized by several steps hereby summarized (all the steps are graphically reported as a flowchart in Fig. 1): first, we built the brain functional network by computing the correlation matrices among brain areas; then, we calculated the variability of the correlation matrices across subjects in order to evaluate the inter-subject variability; subsequently, we ran a percolation analysis and applied a MST filtration in order to test network resilience. Briefly, percolation analyses can elucidate network redundancy or robustness where, for example, apparently redundant nodes provide alternative dispersal pathways should any nodes be lost or impacted. MST filtration was computed by keeping the most robust links discarding all the others. Finally, we applied an allometric scaling to the MST in order to study its topological properties.
Network construction. We used an AAL mask 52 to parcellate the human brain in 58 anatomical regions for each hemisphere. The fMRI signals were extracted from each voxel and averaged in each region of interest (ROI), ending with 116 BOLD time-series. Then, Pearson's correlation coefficient was used to estimate timeseries pairwise similarity, which led to a full and symmetric correlation matrix. The debate on the meaning of negative correlation values in fMRI studies is still open 53 without a consensus supporting their inclusion or exclusion in the analysis. Since they do not bring any relevant information to the actual study's scope, each correlation coefficient was transformed in Z-values 54 and squared, returning a weighted adjacency matrix associated with a weighted fully connected graph.
Inter-subject variability. To evaluate the variability of the correlation matrices across subjects, for each diagnostic group (i.e., HC and SCZ), the variation coefficient was computed as the standard deviation ratio to the mean across subjects of the correlation values between each pair of brain areas. The coefficient of variation distribution across ROIs was fitted to a log-normal function. The Mann-Whitney-Wilcoxon test 55 was used to address their statistical difference, with the number of degrees of freedom equal to the number of links in the network (N(N − 1)/2, N = 116).

Percolation analysis.
A subject-wise percolation analysis 39 for the two groups was performed. Given the individual weighted adjacency matrix, its entries were ranked in ascending order. One at a time, the links corresponding to the ordered list values (i.e., the corresponding weights) were removed by deleting them from the www.nature.com/scientificreports/ rank. The number of connected components and the size of the largest one (Giant Component) were evaluated progressively after each link was removed (i.e., at each percolation threshold). The percolation curves, intended as the number of connected components as a function of the weight removed, were averaged across subjects.
Results are reported together with the 95% confidence interval. The size of the giant component across subjects and the node degree variation across ROIs were investigated to quantify the difference between the two percolation processes. We estimated the lines of reduction as the boundaries over which the giant component's size resulted reduced by 5-50% of its original size (namely 116). For each line of reduction, the area under the curve (AUC) was calculated. The difference between the areas corresponding to two consecutive curves, ∆AUC, normalized to the AUC associated with the maximum giant component size (AUC MaxGCSize ), was computed in the two groups. The rate of change of the two trends was established by fitting them to the exponential decay.
As far as the comparison between HC and SCZ ROIs degree at different percolation thresholds, we calculated the square root of the sum of the squared differences between consecutive degree values across the thresholds (Degree Variation Coefficient): where t is the threshold index, and R is the ROI index. The distribution density of these values has been fitted to a log-normal curve. The statistical difference between the two distributions (HC and SCZ) was established through the Mann-Whitney-Wilcoxon test 55 , with the number of degrees of freedom equal to the number of links in each network (N(N − 1)/2, N = 116).

Maximum Spanning Tree.
We computed the Maximum Spanning Tree (MST) for each brain network, keeping each node's maximum weighted link and discarding all the others. The resulting network components were connected through only one connection: the strongest one not forming cycles 39 .
Allometric Scale. Once the MST was computed for each node, two quantities were estimated: (i) A i , the number of nodes forming the subtree having node i as root (including i); and (ii) C i = k A k , where k runs over all nodes in the subtree having root i (including i). It has been observed that in many cases, there exists a power-law relation between these two quantities, C ∝ A η , with the exponent universally identified for food-webs, river, and vascular networks 56,57 . In general, the exponent η, which characterizes the allometric relation of a planar tree, ranges between 1 and 2, where η → 1 for a star-like topology and η = 2 for a one-dimensional chainlike structure. Thus, the eventual observation of η close to 1 (2) is the signature of a global star-like (chain-like) structure of the whole tree. The MSTs associated with the human (correlation) functional brain networks are undirected by construction 39,58 . Thus, to compute the quantities mentioned above (A i , C i ) for each node i, we have artificially introduced a directionality by choosing one ROI of the undirected tree as the root, determining the different paths towards the remaining ROIs. In Fig. 2b-d, we show three examples of such directionality on a toy model (Fig. 2a) composed of ten nodes related to three possible choices of the root: nodes number 5, 7, and 1. Furthermore, in Fig. 2e, we give a pictorial representation of the directed version of the MST of the human functional brain network once the choice of a particular node (ROI 51) as root was made. We considered all the possible directed versions, considering each node as root and associating to the MST the induced direction. In other words, we built 116 replicas of the MSTs related to each subject. Such MSTs have the same list of edges and weights, with different links direction.
Since this outcome makes meaningless the analysis of an "average brain" for HC and SCZ, we performed a subject-wise study, showing the average and the 95% confidence interval of the analysed quantities.
Percolation analysis. For each subject separately, we evaluated the percolation process in the functional network 39,58 . Figure 4 shows the average percolation curve computed for the two groups (HC and SCZ) together with the 95% confidence interval. A net separation emerges between the average percolation curves of the two groups. Specifically, the functional brain network in SCZ appears more resistant to a disgregation process induced by the removal of weak links. Accordingly, the number of connected components remains smaller than the number observed in HC in the range [0.2, 0.8] of thresholds. It suggests a stronger resistance to the decomposition in separate components of the global network architecture for the SCZ group.
The size of the giant component within each network (HC and SCZ) was estimated across thresholds (Fig. 5a,b) to estimate the source of such difference in the percolation processes quantitatively. Firstly, the area under the curve associated with the maximum giant component size (AUC MaxGCSize ) normalized to the total area shows that the SCZ group is more resistant to the initial fragmentation than the HC one (AUC SCZ /AUC We also checked if the differences between the two percolation processes can be ascribed to a different distribution of the weights within each network. Figure 6 shows the correlation values of the networks of all subjects pooled together (squared values in the inset). Even if the distribution of the weights collected from SCZ is wider than the one coming from HC, not significant differences were found according to a Mann-Whitney-Wilcoxon test. On the contrary, by comparing how each ROI degree (averaged across subjects) changes at each threshold value (Fig. 7a,b), we discovered relevant differences between the two groups. Specifically, the distributions of the Degree Variation Coefficient for HC and SCZ (Fig. 7c) have been found significantly different (Z = 51, p < 0.001) according to a Mann-Whitney-Wilcoxon test.
Accordingly, in the case of SCZ, the node degree decrease appears almost homogeneous for all ROIs. In other terms, our results reveal that during the percolation, weak links are removed almost 'randomly' from the functional brain network of schizophrenic patients, with the consequence of affecting almost in the same way the variation of the degree of each ROI. Simultaneously, for the HC group, nodes disconnect rapidly from the  www.nature.com/scientificreports/ rest of the network and become soon isolated. It suggests that weaker links are uniformly distributed in patients' functional brain networks instead of being concentrated around the same node or groups of nodes as for the healthy subjects.

Maximum Spanning Tree.
In order to investigate more on the differential functional organization in HC and SCZ, we filtered the correlation network of each individual by extracting its Maximum Spanning Tree (MST) 39 . Both HC and SCZ showed a chain-like arrangement of their MSTs. The divergence of the MSTs structure from a linear organization was quantified according to the analysis of the allometric exponent. Once the MST from each individual was computed, and the root induced the directions along the tree, the exponent η of the allometric scaling law was computed for each individual. This procedure was separately repeated for 116 trees (one for each ROI chosen as root) for both HC and SCZ subjects. Figure 8 shows the η values for each ROI, averaged across subjects, reported with the 95% confidence interval. The mean allometric exponents for the two groups were the same across roots (Fig. 8), revealing on average a similar organization of the functional brain network backbone for HC and SCZ. www.nature.com/scientificreports/ To investigate further the possible origin of the differences showed by the percolation analysis, we extended the allometric analysis beyond the correlation network backbone. Specifically, we defined higher rank MSTs as follows: first, we rename the correlation network's MST as the first rank MST. Then, we identify a second rank MST as the one obtained from the original correlation network by eliminating the links belonging to the first rank MST. Accordingly, we can define the nth rank MST as the MST of the correlation network from which we have already removed all the previous n − 1 MSTs. We computed the first four MSTs for each ROI used as root and for each of the subjects. Figure 9 shows that, starting from a similar condition, allometric exponents tend to reduce as the MST rank increases. However, the different speed of reduction (SR) determines an evident and growing difference between the two groups of subjects, with the allometric exponents decreasing more slowly for SCZ than HC. In fact, an exponential fit (η(Rank) = (η 0 − η ∞ ) * exp(− Rank/SR) + η ∞ ) of the average η across the ROIs returned SR HC = (0.610 ± 0.001) and SR SCZ = (0.802 ± 0.001), for HC and SCZ respectively.

Discussion
Here we report an advanced network-based analysis of fMRI data recorded at rest in SCZ and HC. We aimed to assess possible alterations in the brain's basal functional network organization of SCZ through percolation and MST analyses 39,58 . Our findings demonstrate a global change of the connectivity strength distribution in the functional networks of SCZ as compared to HC, consisting of an increased homogeneity of the distribution of the weighted links across the whole network.
Specifically, from the inspection and analysis of the percolation curves, we found that the functional networks in SCZ are characterized by a region-to-region interaction which is more resistant to disconnection than in HC. Indeed, by scrolling the threshold (below which links are erased) from low to high connectivity strength values, the number of separate clusters is systematically smaller in SCZ. Similarly, we found that, when comparing the size of the giant component 59 as a function of the threshold, the size reduction in HC is always faster than that observed in SCZ. A wider distribution of the connectivity strengths in SCZ cannot explain such differences in the percolated networks, mainly because the two distributions become very similar when the squared correlation values are considered.
On the other hand, analysing the distribution of the node degree computed for each percolated network, a clear discrepancy was found between the two groups with consequent delayed emergence of a second connected component in the percolating network. This outcome sheds light on the existence of a more homogeneous distribution of connections in the functional architecture of brains of patients diagnosed with SCZ compared to HC. Since the modular structure of a network is determined by the unbalance or predominance of inward and outward links of individual communities, rather than by the average total distribution of edge weights, an effect of the re-arrangement of the weights may result in a progressive loss of hierarchy of functional connections, which, in turns, leads to an alteration or disruption of modularity.
Modularity structure modifications in SCZ were reported as due to changes in community participation of nodes in the somatosensory, subcortical, auditory, default mode, and salience networks 35 and as a specific fragmentation of somatosensory cortices 36 . Nonetheless, the local specificity of the aberrant cerebral functioning in SCZ is still under debate. Hypo-connectivity was reported in the olfactory cortex, temporal pole, angular gyrus, parahippocampus, amygdala, caudate, and pallidum 28,32 . Conversely, hyper-connectivity was found within the default mode network 60 , the bilateral striatum 61 , and between the default mode and visual and motor regions 62 . These findings suggest a widespread modulation of connectivity in SCZ that involves both directions (i.e., increases and decreases) and are in line with our results: functional brain networks in SCZ are characterized by a more homogeneous distribution of weights, where strong correlation patterns and weak ones share the same topology. Analogous conclusions were reached by authors 28 , who pointed to the altered graph topology in SCZ as responsible for changes in the brain's activity and connectivity complexity.
Conversely, the MST analysis of the functional network does not differ between the two groups, thus implying that the backbone of the functional organization seems to be very similar between groups 39 . A quantitative www.nature.com/scientificreports/ confirmation of this result is given by the allometric exponent calculated for each subject. MSTs derived from functional brain networks can, in principle, range between two different topologies: the linear configuration (allometric exponent equal to 2), the least efficient with the lowest cost of functioning, and the star-like organization (allometric exponent that tends to 1), the most efficient but with the highest cost of functioning 63 . The configuration of the ones we derived from functional brain networks is associated with an allometric exponent approximately equal to 1.5, a configuration that balances efficiency and cost. In the context of functional brain connectivity, which involves interaction between regions not necessarily close one to the other, and where the proximity is just a consequence of their synchronization, efficiency and cost are intended as the ability to organize information transfer via a synchronous coordination 64 and the provision of energetic resources that sustain such organization 65 respectively. Since we observed a widespread change in the distribution of connectivity strength in brain networks of SCZ along with an almost unaltered MST topology, we explored the allometric properties  www.nature.com/scientificreports/ of the MSTs at increasing rank. Our findings show that by progressively removing MSTs from a functional network, the value of the allometric exponent decreases both in SCZ and in HC; however, the rate of reduction of the exponent is faster in HC, leading to a net separation of the two groups at the third rank MST. It suggests that not only the hierarchy and modular structure of the functional network are reduced by a homogeneous distribution of the connectivity strengths but also that in terms of the global integration of brain regions in the network, weaker links in SCZ guarantee the same topology shown by the stronger ones 28 . The higher topological similarity of the MSTs of different ranks, observed in SCZ as compared to HC, suggests that while a given stimulus can engage a single functional connectivity path 39 in patients, it determines the simultaneous involvement of different ones. A possible explanation of this conclusion can be found in the disconnection hypothesis 8,66 , which relates the molecular and the neuronal pathophysiology to SCZ symptoms. It states that this psychiatric disorder may stem from an abnormal response of the NMDA receptor to specific neuromodulatory receptor activation. A failure of such mechanism may lead to an inability to modulate the congruity of sensory evidence, corresponding to the congruity of beliefs about the causes of sensory cues, and consequently to false inference (e.g., hallucinations and delusions) 66 . Accordingly, our findings show that a reshuffling of connectivity strengths with the realization of alternative connectivity patterns occurs in the illness that disrupts the local modular organization (i.e., the somatosensory community). A consequence of this multi-choice configuration is the loss of higher levels of cortical hierarchies that generate predictions of representations in lower levels, jeopardizing in SCZ the brain's ability to process sensory information by optimizing explanations for its sensations [67][68][69] . It is worth further commenting on the possible physiological causes of the heterogeneity reduction, and the loss of hierarchy showed by the correlation strength observed in SCZ. In a computational model 70 , it was proposed that the best way to obtain functional networks with topological properties matching those reported experimentally in SCZ would be the decrease of the strength of excitatory synaptic inputs as a consequence of the disruption of synaptic mechanisms. At the neurophysiological level, the disconnection hypothesis accounts for such altered mechanisms 9 . Nonetheless, an excess of synaptic refinement (enhanced pruning) has also been hypothesized to underlie the neuropathology of SCZ 71,72 . Our results suggest that the brain activity in SCZ is characterized by a subtle change of the global functional architecture, which is not random but involves both an increase and a decrease of the local connectivity strengths. It probably occurs due to the brain's attempt to compensate for an imbalance of the local homeostatic signaling 73 . Such untargeted compensation may homogenize the patterns of information spread in the brain, which is expected to reach all areas without the whole system's involvement in the activity 74 .
In summary, this study indicates how previously reported fragmentation of the modular structure of functional connectivity in medicated patients diagnosed with SCZ might be due to a redistribution, and consequently a homogenization, of the connectivity strengths among all the regions of the brain. Our findings support the theory that aberrant connectivity may induce deficits that propagate to higher functions through a bottom-up process 66 . Moreover, we report several equivalent basal functional skeletons in SCZ (MSTs), which implies the lack of specificity in the brain's functional connectivity organization that is actually expected to be hierarchically set up to account for the correct integration for the proper functioning.  Allometric exponent for higher-order MSTs. Mean allometric exponent η, separately averaged over the two groups HC and SCZ, with 95% confidence interval for the second, third and fourth rank MSTs. The speed of reduction (SR) of the allometric exponent is larger for SCZ than for HC. An exponential fit (η(rank) = (η 0 − η ∞ ) * exp(− rank/SR) + η ∞ ) of the average η across the ROIs returned SR HC = (0.61 ± 0.08) SR SCZ = (0.80 ± 0.10). Each ROI was chosen as a directionality root.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.