Introduction

Data-driven connectivity analysis of resting-state functional magnetic resonance imaging (rs-fMRI)1,2 has revealed abnormalities in brain network topology in several mental disorders, particularly schizophrenia (SCZ)3,4,5,6. Often described in terms of circumscribed alterations in neural circuits7, SCZ has been explained in terms of aberrant interactions among brain regions leading to an integration of distributed brain networks, referred to as “misconnection” or “dysconnection” syndrome8,9. According to this hypothesis, physiological interactions between distributed neuronal ensembles are critical for the effective coordination of information processing and, consequently, the production of coherent action and cognition10,11,12,13,14,15.

Synaptic pruning has been hypothesized to underlie the neuropathology of SCZ16. In particular, neural networks in the development result from pruning processes, including expansive growth of synapses followed by activity-dependent elimination17,18. A dysfunctional synaptic pruning generally implies normal synapses formation, followed by an unbalanced process of elimination19.

Brain functional connectivity has been proven abnormal in specific brain circuits in SCZ20,21,22 and associated with alterations in perception, thoughts, mood, and behavior, that are typical of this neuropsychiatric illness23,24,25,26,27.

Most of the whole-brain approaches to the study of brain functional and structural connectedness showed widespread disturbances in the network dynamics28,29,30,31,32,33,34 and alterations of the whole cerebral functional architecture's modular organization in SCZ32,35,36,37,38. Nonetheless, the alterations in the global functional integration and the brain's local functional connectedness reported in literature need to be clarified. In fact, previous research addressed the local and global brain connectedness disjointly, characterizing either the collective behavior of 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 network39 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 developed39 advanced methods of network neuroscience40,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-542 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 score44 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-RV43 and SCID-5-PD45, 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).

Table 1 Principal demographic characteristics of the studied subjects.

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 Declaration46.

Data acquisition and pre-processing

FMRI data were collected at 3 T (Philips Achieva equipped with a thirty-two channel receive-only head coil) using gradient-echo echo-planar imaging for the registration of the blood oxygen level-dependent (BOLD) signal (TR = 3 s, TE = 30 ms, matrix = 80 × 80, FOV = 224 × 224, slice thickness = 3 mm, flip angle = 90 Å, 50 slices, 240 vol).

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 regression47 (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 rates48,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 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 following51.

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.

Figure 1
figure 1

Flowchart of the network analysis. The temporal progression of the network analysis can be divided into three main parts: the network construction from BOLD time-series averaged in the regions of interest of the AAL atlas, the percolation analysis that highlights the integration process of the network's nodes into modules according to the strength of their connectivity and finally the unpacking of the whole brain graph in terms of Maximum Spanning Trees (MSTs) that allows one to identify deviations from the hierarchical organization of the brain.

Network construction

We used an AAL mask52 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 time-series 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 open53 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-values54 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 test55 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 analysis39 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 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 (AUCMaxGCSize), 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):

$${\Delta d}^{{\text{R}}} = \left[ {\mathop \sum \limits_{{{\text{t}} = 1}}^{{\text{N}}} \left( {d_{t + 1}^{R} - d_{t}^{R} } \right)^{2} } \right]^{\frac{1}{2}} ,$$

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 test55, 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 cycles39.

Allometric Scale

Once the MST was computed for each node, two quantities were estimated: (i) Ai, the number of nodes forming the subtree having node i as root (including i); and (ii) \(C_{i} = \sum\nolimits_{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 networks56,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 chain-like 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 construction39,58.

Thus, to compute the quantities mentioned above (Ai, Ci) 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.

Figure 2
figure 2

Allometric scale. (a) A toy model of 10 nodes. According to the different nodes chosen as roots, it is possible to obtain ten directed versions of the undirected tree (a). (b)–(d) Three directed versions of the tree in (a) with roots equal to, respectively, nodes 5, 7, and 1. Numbers in brackets represent the quantities (Ai, Ci) associated with each node. (e) Maximum spanning tree of the human functional brain network of an individual chosen at random in healthy subjects with directionality induced by choosing node fifty-one as root. Colors in (e) represent anatomical regions according to the grouping of AAL parcellation shown on the right bottom: Dark blue circles: Frontal Lobe; Orange circles: Insula; Light blue circles: Cingulate; Green circles: Temporal Lobe; Red circles: Occipital Lobe; Yellow circles: Parietal Lobe; Grey circles: Deep Grey Matter; Magenta circles: Cerebellum.

Results

Inter-subject variability

Figure 3 shows the coefficient of variation for HC and SCZ. According to the analysis of the two cohorts of subjects, HC (Fig. 3a) returned a level of homogeneity much larger than the one observed for SCZ (Fig. 3b). Specifically, the distribution of the coefficients of variation across ROIs was fitted to a log–normal function both for HC and SCZ, with the following mean values and standard deviations with the relative errors of curve fitting: [(μHTH ± ∆μHTH) = (0.74 ± 0.01), (σHTH ± ∆σHTH) = (0.21 ± 0.01)] and [(μSCZ ± ∆μSCZ) = (1.41 ± 0.01), (σSCZ ± ∆σSCZ) = (0.85 ± 0.01)] (Fig. 3c). According to a Mann–Whitney–Wilcoxon test, the two distributions were found significantly different (Z = 57, p < 0.001).

Figure 3
figure 3

Inter-subjects variability. (a) Coefficients of variation matrix for HC, (b) coefficients of variation matrix for SCZ, and (c) density of the distribution of the pairwise coefficients of variation for HC (blue) and SCZ (red). The coefficient of variation was calculated as the standard deviation ratio to the mean across subjects of the correlation values between each pair of brain areas. Data were fitted to a log-normal distribution with mean μ and variance σ2. The curve fitting results are reported for both HC and SCZ with the errors of the fitting procedure. According to a Mann–Whitney–Wilcoxon test, the two distributions resulted significantly different (Z = 57, p < 0.001).

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 network39,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.

Figure 4
figure 4

Percolation analysis. The number of connected components of the percolated network versus the related correlation threshold for HC (blue) and SCZ (red). The two curves represent the average of the individual percolation curves and are reported with a 95% confidence interval.

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 (AUCMaxGCSize) normalized to the total area shows that the SCZ group is more resistant to the initial fragmentation than the HC one (AUCSCZ/AUCtotal = 0.30 and AUCHC/AUCtotal = 0.16). It implies a delay at the beginning of the disgregation process between SCZ and HC. Moreover, the estimation of the disconnection rate (Fig. 5c) demonstrates that, even if the process starts later in SCZ, its progression is faster, pairing the two groups when the giant component size is halved. The difference between the areas under two consecutive curves, ∆AUC, normalized to the area under the curve associated with the maximum giant component size, AUCMaxGCSize, at different values of the giant component size reduction was fitted to the exponential decay. The speed of decay (SoD) was found to be SoDHC = (5.03 ± 0.02) and SoDSCZ = (3.11 ± 0.01).

Figure 5
figure 5

Percolation and giant component analysis. (a) Giant component size variation for HC, (b) giant component size variation for SCZ, and (c) the difference between the areas under two consecutive curves, ∆AUC, normalized to the area under the curve associated with the maximum giant component size, AUCMaxGCSize (yellow line), for HC (blue) and SCZ (red). ∆AUC/AUCMaxGCSize, for different values of the giant component size reduction, was fitted to the exponential decay. The speed of decay (SoD) was found to be SoDHC = (5.03 ± 0.02) and SoDSCZ = (3.11 ± 0.01) for HC and SCZ, respectively. Dotted lines represent 95% of prediction bounds.

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.

Figure 6
figure 6

Weights distribution comparison. The density of the distributions of the human functional brain networks' correlation values for healthy subjects (blue) and schizophrenic patients (red). Inset: distribution of the squared correlation values.

Figure 7
figure 7

Percolation and node degree analysis. (a) Variation of the node degree averaged across subjects for HC. (b) Variation of the node degree averaged across subjects for SCZ, (c) density of the distribution of the Degree Variation Coefficients for HC (blue) and SCZ (red). Data were fitted to a log-normal distribution with mean μ and variance σ2. The curve fitting results are reported for both HC and SCZ with the errors of the fitting procedure. According to a Mann–Whitney–Wilcoxon test, the two distributions resulted significantly different (Z = 51, p < 0.001).

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 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.

Figure 8
figure 8

Allometric exponent for the first rank MST. Mean allometric exponent η, separately averaged over the two groups HC and SCZ, with 95% confidence interval for the first rank MST. Each ROI was chosen as a directionality root.

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 SRHC = (0.610 ± 0.001) and SRSCZ = (0.802 ± 0.001), for HC and SCZ respectively.

Figure 9
figure 9

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 SRHC = (0.61 ± 0.08) SRSCZ = (0.80 ± 0.10). Each ROI was chosen as a directionality root.

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 analyses39,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 component59 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 networks35 and as a specific fragmentation of somatosensory cortices36. 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 pallidum28,32. Conversely, hyper-connectivity was found within the default mode network60, the bilateral striatum61, and between the default mode and visual and motor regions62. 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 authors28, 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 groups39. A quantitative 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 functioning63. 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 coordination64 and the provision of energetic resources that sustain such organization65 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 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 ones28. 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 path39 in patients, it determines the simultaneous involvement of different ones. A possible explanation of this conclusion can be found in the disconnection hypothesis8,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 sensations67,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 model70, 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 mechanisms9. Nonetheless, an excess of synaptic refinement (enhanced pruning) has also been hypothesized to underlie the neuropathology of SCZ71,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 signaling73. 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 activity74.

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 process66. 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.