Introduction

Human brain is an inherently complex network even at rest, comprising nearly 10 billion neurons connected by about 100 trillion synapses1. Yet, the brain owes its beautiful complexity and the consequent cognitive capacities not solely to the size but also the interconnected interactions between its constituent elements. In the last two decades, graph theory has provided a valuable framework to study the structure of the complex brain networks2. The C. elegance connectome was modeled as a binary graph, as one of the first contact points between neuroscience and modern network science3. The tract-tracing data of cat and macaque monkey was examined in early 20004 and later graph theory started to be well applied to the human neuroimaging data5. Along with the progress in neuroimaging techniques, there has also been growing interest in investigating patterns of associations between time series extracted from the brain regions, i.e., functional connectivity. Altogether, graph-theoretic studies are an important foundation for computational modeling of the complex brain networks, and have revealed fundamental properties of their organization and function, in addition to their alterations in brain disorders6,7. While all these advances and the potential of graph theory’s perspective is undisputed, a vital question have yet to be asked: Suppose in a signed brain network x, y and z regions are connected, what is the inevitable impact of xz and yz interactions on the sign and weight of the interaction between x and y? Is it realistic to study the xy interaction independent from the triadic interconnection (xyz) it lives in? Specially in case of heterogeneous brain disorders such as ASD, a complex neurodevelopmental disorder with substantial heterogeneity not only in multiple causes and courses but also in the onset8,9,10,11. Recently, this concern has been well addressed by the long-standing structural balance theory (SBT) in variety of scenarios, and triadic interactions are accepted to play key role in organization of real-world signed networks12,13,14,15,16. Here, we apply concepts from SBT to the weighted signed rs-fMRI networks of individuals with ASD compared to healthy controls in order to answer the following questions: Considering recent verification of SBT by many real-world signed networks, do the brain networks display structural balance as well? Regarding controversies on strong and weak notions of SBT, which will be confirmed on the brain networks? Knowing the previously proposed hypo (hyper) functional connectivity in ASD, if we take into account the interplay between positive and negative links, which types of triadic interactions (in terms of both the frequency and energy distribution) are atypical in ASD? Is the consequent structural balance altered for the whole brain ASD network in course of development? What about the functional sub-networks? These questions in essence, are the cornerstone of this study.

Initially proposed by Heider17 and mathematically formulated by Cartwright and Harary18, SBT has been a promising approach in arguing why the structures of many real signed social12,13,14, political15, ecological16 and biological19 networks are the way they are. The distinguishing feature of SBT is that it highlights the role that balanced (unbalanced) triadic interactions play in forming the global structure of a network. This theory has led to better understanding of how a tendency to reduce overall stress directs and arranges the organization of signed networks. Specifically, SBT argues that a link between two nodes being positive (negative) is hugely affected by a third node’s presence within a triadic structure. Accordingly, four types of triads are defined in this context (Eq. 1, Fig. 6), namely, strongly balanced \(T_3: (+ + +)\), weakly balanced \(T_1: (+{-}{-})\), strongly unbalanced \(T_2: (+ + -)\), and weakly unbalanced \(T_0: (- - -)\), details of which is provided in methods. Frequent building blocks of two to five vertices, known as the structural and functional motifs, have been identified and investigated in the brain networks through applying framework’s other than SBT as well20,21. However, there are two main differences between triads as motifs from one hand and triads as has been defined in SBT from the other hand: (1) Sign of links are not of concern in motifs, yet they are at the heart of triads’ definition in the context of SBT, (2) Motifs are mostly directed structures while triads in SBT are not originally. Moreover, measures of segregation such as the clustering coefficient or transitivity are also based on the frequency of triads, yet again they measure properties of a network regardless of signs22. This is while, it has been recently confirmed on the real-world networks that studying the interplay between positive and negative links within triadic interactions, instead of focusing exclusively on positive links, opens new doors to existing challenges12,13.

Signed networks not only represent the structure but also embody additional information on the state of relationships between nodes, and has been long of interest to network scientists in different domains from social science to biology. In the last decade, many social scientists have reported that in analyzing large-scale social networks considering the content of interactions, positive (negative) links representing friendship (enmity), contains much promise. In this regard, the key role that triadic interactions play in forming the global structure of signed social networks have been strongly confirmed23,24,25. Similarly in biological and biochemical signed networks the positive (negative) ties correspond to activating, correlating (inhibiting, anti-correlating) interactions26. There has been various evidences of anti-correlated functional sub-networks in the intrinsic brain networks of both human27 and animal28 as well. These studies imply that negative associations are meaningful in the resting-state brain networks and not only the unwelcome result of preprocessing steps, such as the global signal regression29. Although it has been shown that regressing the global signal introduces anti-correlations to the network30, thus here we did not regressed the global signal to make sure existing negative links are in a way reflections of neural activities. To explore real-world signed networks, SBT as both a general theory as well as a practical framework for conducting experiments has been very reassuring. Since its proposal, there has been two main directions in the literature of SBT: (1) Theoretically extending its analytical aspects31,32,33,34,35, (2) Empirically verifying it on real signed networks13,16, or both12,14. Moreover, a recent study has investigated the brain evidence for SBT from psychological viewpoint36, this is while we investigate the network evidence for the theory through examining Blood Oxygen Level-Dependent (BOLD) signals. Furthermore, besides SBT which investigates the structural balance of complex networks through studying triadic interactions, there has been valuable research to explore higher-order interactions in the brain networks based on topological properties as well37,38,39.

In this study, we explore the weighted signed rs-fMRI networks of individuals with ASD compared to healthy participants in the context of SBT. We conduct our analysis in three age ranges, namely, \(1{\mathrm{st}}\) childhood (6–9 years old), \(2{\mathrm{nd}}\) childhood (9–13 years old) and adolescence (13–18 years old). First, according to previous empirical studies, we expect the brain networks to display structural balance as well—hypothesis #1, meaning that balanced triads are overrepresented while unbalanced triads are underrepresented compare to chance. Although, due to controversies regarding the over (under) representation of \(T_0\), it is nontrivial whether the strong or weak notion of balance would be confirmed. Next, we study the frequency of pair and triadic interactions in the brain networks. Considering previous results on hypo (hyper) functional connectivity in the ASD brain network, it seems that taking into account sign of links and exploring triadic interactions provides us with a more contextual insight into atypical functional connectivity in the ASD brain network—hypothesis #2. Moreover, we illustrate the energy distributions of triads and suppose that there are probably of note differences regarding energy distributions of triads between ASD and CON networks—hypothesis #3. However, it is not clear which types of triads are different in terms of energy and during which age ranges. Last but not least, we examine the overall energy of the whole-brain network and 17 Yeo sub-networks and expect to observe alterations in the energy of the ASD brain (sub) networks, possibly suffering from lower or more limited energy scales in some of the key (sub) networks—hypothesis #4. Our study, while appreciating standard approaches in analyzing functional connectivity, proposes a new perspective in exploring rs-fMRI brain networks, which can be specially of interest in investigating complex brain disorders such as ASD.

Results

Empirical evidence of SBT on the brain networks

First, we have investigated the notion of structural balance in the brain networks. According to the strong version of this notion, real-world networks evolve in a way that eventually unbalanced triads are underrepresented while balanced triads are overrepresented. The weak notion however, is less strict and argues that \(T_0\) triads, although unbalanced yet may be overrepresented as well. To this aim, we have applied a method proposed by Leskovec et al.12 which states that for a triad \(T_i\) if \(p(T_i) > p_0(T_i)\) then \(T_i\) is overrepresented, and if \(p(T_i) < p_0(T_i)\) then \(T_i\) is underrepresented. Here \(p_0(T_i)\) is the fraction of triads of type \(T_i\) after shuffling, details of which are provided in methods. Furthermore, to investigate how big these over (under) representations are, we have applied the concept of surprise, i.e, \(s(T_i)\), which on the order of tens would be significant, due to the distribution of \(s(T_i)\) being standard normal.

As our results show, on average in the brain networks of both CON (Table 1A) and ASD (Table 1B) groups and in all three age ranges both balanced triads are overrepresented relative to chance, that is for \(T_3\) and \(T_1\) we have \(p(T_i) > p_0(T_i)\). On the contrary, on average both unbalanced triads are underrepresented compared to the null model, that is, we have observed \(p(T_i) < p_0(T_i)\) for both \(T_2\) and \(T_0\) triads. Additionally, all the surprises have been significant. It should be mentioned that these results are group-level, that is, averaged over all participants in each group. Results on the level of each participant is the same for \(T_3\), \(T_2\) and \(T_1\) triads, meaning that for each and every participant in all groups we have overrepresentation regarding \(T_3\) and \(T_1\) while underrepresentation for \(T_2\) triads. However, for \(T_0\) triads a small percentage of participant’s networks have more \(T_0\) triads relative to chance as follows: for ASD and CON groups during \(1{\mathrm{st}}\) childhood the percentage of networks with overrepresentation of \(T_0\) triads were \(22.2 \%\) and \(10 \%\), respectively. During \(2{\mathrm{nd}}\) childhood this percentages were \(3.8 \%\) in ASD group and \(5.7 \%\) in CON group. Finally, during adolescence in ASD group \(3.4\%\) and in CON group only \(1.7\%\) of networks had more \(T_0\) triads relative to chance.

Table 1 Number and probability of triads in the brain networks compared to the null model. On average, in both CON and ASD groups, both types of balanced triads, \(T_3\) and \(T_1\), are overrepresented (positive \(s (T_i)\)), while both types of unbalanced triads, \(T_2\) and \(T_0\), are underrepresented (negative \(s (T_i)\)).
Full size table

Why considering triadic interactions in studying the brain networks?

After verifying SBT on the brain networks, we have explored the value of examining ternary interactions (triads) in the brain network along with the pair interactions (links). In other words, why consider triadic interactions while analyzing brain networks? To study this question, we have first compared mean differences in the frequency of positive and negative links between ASD and CON groups in three age ranges. Then, we have conducted the same analysis for the mean frequency of different types of triads. As results show, although exploring links and triads are both valuable in revealing statistically significant differences between groups, acknowledging triads provides results with bigger effect sizes, that is, practical significance, details of which are as follows.

First, to compare the mean differences in the frequency of links between ASD and CON groups during development, we have conducted two two-way ANCOVAs, one with the frequency of positive and the other with negative links as a dependent variable (Table 2A). We have considered group and age as independent variables while controlling for FIQ, medication, mean frame-wise displacement as head motion parameter and site information. We are interested in the main effect of group and the interaction between group and age. While the former provides an overall difference between ASD and CON groups, the latter allows us to explore the differences during development. Then, we have conducted the same analysis for triads, that is, we have defined four two-way ANCOVAs one for each type of triads as has been defined in Eq. 1 (Table 2B). Another option was to conduct two MANCOVAs, one for links and the other for triads, however due to multicollinearity between the two types of links and the four types of triads we chose ANCOVA over MANCOVA. As results show (Table 2), the main effect of group is neither significant on the mean frequency of links nor triads. However, while the effect of interaction between group and age is only statistically significant for positive links (a medium effect), it is significant for all types of triads. More interestingly, it is only practically significant on the mean frequency of \(T_2\) triads \((F(1, 251)=95.17 \ p<0.001,\ partial \eta ^2=0.46\)).

Table 2 Analysis of covariance for the frequency of links versus triads. (A) The effect of group is not significant on the mean frequency of links. The effect of interaction between group and age is only significant for the positive links, yet not practically significant. (B) For triads regarding the effect of interaction between group and age, not only all the differences are statistically significant but the effect size (partial \(\eta ^2\)) is interestingly considerable for \(T_2\).
Full size table

To further explore in which age levels the differences in mean frequency of links (Fig. 1A) and triads (Fig. 1B) have occurred between ASD and CON groups, we have conducted post-hoc tests using Mann-Whitney U test. As results depict, during \(1{\mathrm{st}}\) childhood the frequency of positive links from one hand \((U(20,18)=102.00,\ z=-2.28, p=0.02)\), and \(T_1\) triads from the other hand \((U(20,18)=101.00,\ z=-2.31, p=0.02)\) are both practically significant \((\eta ^2=0.14)\). During \(2{\mathrm{nd}}\) childhood, while there is a medium difference \((\eta ^2=0.10)\) in the mean frequency of positive links \((U(52,52)=1853.00,\ z=3.25, p=0.001)\), there is a large difference \((\eta ^2=0.39)\) in the mean frequency of \(T_2\) triads \((U(52,52)=2338.00,\ z=6.41, p<0.001)\). Additionally, there is also a medium difference \((\eta ^2=0.10)\) in the mean frequency of \(T_3\) triads between ASD and CON groups during \(2{\mathrm{nd}}\) childhood \((U(52,52)=1868.00,\ z=3.35, p=0.001)\). Similarly, for adolescents although there is a significant difference for both types of links (for positive links: \(U(58,58)=1181.50,\ z=-2.76, p=0.006\) and for negative links: \(U(58,58)=1010.00,\ z=-3.71, p<0.001)\), yet both of these differences are medium, \(\eta ^2=0.06\) and \(\eta ^2=0.12\) for positive and negative links, respectively. However, surprisingly for \(T_2\) triads the difference is far more bigger than both types of links \((U(58,58)=162.00,\ z=-8.39, p<0.001, \eta ^2=0.60)\). There is also a big difference \((\eta ^2=0.23)\) in the mean frequency of \(T_1\) triads \((U(58,58)=745.00,\ z=-5.17, p<0.001)\). Additionally, there is a medium difference \((\eta ^2=0.09)\) in the mean frequency of \(T_0\) triads \((U(58,58)=1102.00,\ z=-3.20, p=0.001)\), and a small difference \((\eta ^2=0.04)\) in the mean frequency of \(T_3\) triads \((U(58,58)=1270.00,\ z=-2.27, p=0.02)\).

Figure 1
figure 1

Mann-Whitney U test on the frequency of links and triads. (A) Difference between ASD and CON groups is only practically significant \((\eta ^2 \ge 0.14)\) for the mean frequency of positive links during \(1{\mathrm{st}}\) childhood. (B) Yet, considering triads provides us with practical differences not only during \(1{\mathrm{st}}\) childhood but also during \(2{\mathrm{nd}}\) childhood and adolescence. \(\eta ^2\), the effect size; CON, control; ASD, autism spectrum disorder (Color Online).

Full size image

Energy distribution of triads

After studying the frequency of triads in the brain networks of ASD and CON groups, we have explored the energy distributions of different types of triads. To this aim, we have calculated the energy of triads, \(U(T_i)\), as shown in Fig. 6, and then for each group in three age ranges outlined the corresponding distributions. As results indicate (Supplementary Fig. S1), for all types of triads the brain networks of ASD and CON groups have many triads with small energies and a few with considerable energies. Furthermore, results of comparisons between energy distributions of triads in ASD and CON groups show that during \(1{\mathrm{st}}\) childhood the pattern of energy distributions differs between the two groups for \(T_1\) and \(T_0\) triads (Fig. 2A). That is, looking at the energy distributions of these two triads in ASD group (purple/dark triangles) compare to CON group (green/ light squares) it is clear that the energy distributions in ASD group lag behind distributions as they are in CON group (\(KL= 0.03\), for both types). For other types of triads and during \(2{\mathrm{nd}}\) childhood and adolescence we have not observed such a lag (Supplementary Fig. S1).

To better clarify this lag, for each \(T_i\) we have investigated node participation which for any given node specifies: (1) In how many \(T_i\) this node has participated? and (2) How big are the energies of those \(T_i\)s? In Fig. 2B and Supplementary Fig. S2, we have addressed the former question with size of a given node and the latter with its color, which is a color-map from blue to red that respectively corresponds to the minimum and maximum levels of triads’ energy. For example, a big blue node is present in many \(T_is\) which are small in terms of energy. On the other hand, a small red node although lives in just a few \(T_i\)s but in those that are important in terms of energy levels. As can be seen, moving from the minimum to maximum level of energy (Fig. 2B.1,B.2; left to right), while node participation decreases sharply to zero in ASD group, the speed of this decrease is slower in CON group. In other words, we have observed a threshold in the energy distributions of both \(T_1\) and \(T_0\) triads above which node participation is zero for ASD group yet nonzero for CON group. To be specific, in ASD group no regions of interest participate in creating \(T_0\) triads that have energies higher than 0.04, while in CON group many nodes are participating so that the resulting network has \(T_0\) triads with energies higher than 0.04. This pattern holds for the case of \(T_1\) triads. That is, while in ASD group node participation is zero for \(T_1\) triads with \(\mid Energy \mid \ \ge 0.32\), the network of CON group during \(1{\mathrm{st}}\) childhood enjoys the presence of \(T_1\) triads with \(\mid Energy \mid \ \ge 0.32\) (Supplementary Fig. S2).

Figure 2
figure 2

Energy distributions and node participation in \(T_1\) and \(T_0\) triads during \(1{\mathrm{st}}\) childhood. (A) The energy distributions of ASD group (purple/dark triangles) lag behind CON group (green/light squares) for \(T_1\) (A.1) and \(T_0\) (A.2) triads during \(1{\mathrm{st}}\) childhood. (B) For both triads we have observed a threshold above which node participation is zero for the brain network of ASD group, yet nonzero for CON group (B.1.3 for \(T_1\) and B.2.3 for \(T_0\)). A detailed version of B.1.3 and B.2.3 with regions of interest’s labels is provided in the Supplementary Fig. S2. BrainNet Viewer 1.6340 (https://www.nitrc.org/projects/bnv) has been used for visualization of the brain networks. \(|T_i|\), total number of \(T_i\); \(U(T_i)\), the energy of \(T_i\); KL, the Kullback–Leiber divergence; CON, control; ASD, autism spectrum disorder (Color Online).

Full size image

Structural balance of the whole brain network and 17 Yeo sub-networks

Next, we have analyzed energy levels of the brain networks in the context of SBT. In this regard, first we have calculated the energy levels of each participant’s network using Eq. 2. For each individual, we have computed 18 different energies, one for the whole brain network and the rest 17 each corresponding to 17 Yeo sub-networks. Then, we have explored the effect of group and its interaction with age on these energy levels controlling for FIQ, medication, mean frame-wise displacement as head motion parameter and site information. We have performed this by conducting 18 ANCOVAs each regarding one of the energy levels as mentioned above, results of which are as follows:

Figure 3
figure 3

Energy of the brain (sub) networks during development. (A) For both groups, energy of the whole brain network increases monotonically with age, yet ANCOVA results show that the pattern of this increase is different between the two groups. (B,C) The effect of age on energy of the SN (A) and DMN (B) is significant, with ASD having less energies in both sub-networks. Results of Mann-Whitney U tests with corresponding effect sizes (\(\eta ^2\)) are shown above each of the box plots. (D,E) Regions of interest related to the SN (A) and DMN (B), as has been defined in the Schaefer atlas. BrainNet Viewer 1.6340 (https://www.nitrc.org/projects/bnv) has been used for visualization of the brain networks. SN, the salience/ventral attention network; DMN, the default mode network; Ins, the insula; PFC, the prefrontal cortex; Temp, the temporal cortex; FrOper, the frontal opercular; ParOper, the parietal opercular; FrMed, the frontal medial; IPL, the inferior parietal lobule; CON, control; ASD, autism spectrum disorder (Color Online).

Full size image

  • The effect of group on the whole brain energy is statistically significant \((F(1, 251)=12.08,\ p<0.001,\ partial \eta ^2=0.05)\). Similarly, the effect of interaction between group and age on the whole brain energy is significant as well \((F(2, 251)=4.57,\ p=0.01, \eta ^2=0.04)\).

  • The effect of group is significant on the energy levels of the SN (A) \((F(1, 251)=12.52,\ p<0.001,\ partial \eta ^2=0.05)\) and the DMN (B) \((F(1, 251)=6.73,\ p=0.01,\ partial \eta ^2=0.03)\).

Moreover, to identify the groups to which these differences are related to we have conducted post-hoc tests using Mann-Whitney U test. Results for the energy levels of the whole brain networks are as follows (Fig. 3A):

  • During \(1{\mathrm{st}}\) childhood, energy of the whole brain network is significantly greater for ASD than CON group \((U(20,18)=261.00,\ z=2.36, p=0.01)\) and it is not only statistically but also practically significant \((\eta ^2=0.15)\).

  • For both the ASD and CON brain networks, there is a practically significant difference between energy levels during \(1{\mathrm{st}}\) childhood from one hand, with \(2{\mathrm{nd}}\) childhood and adolescence from the other hand. Results of comparison between \(1{\mathrm{st}}\) childhood and \(2{\mathrm{nd}}\) childhood for ASD group is \(U(18,52)=854.00,\ z=5.18, p<0.001, \eta ^2=0.38\), and for CON group it is, \(U(20,52)=1018.00,\ z=6.26, p<0.001, \eta ^2=0.55\). Additionally, between \(1{\mathrm{st}}\) childhood and adolescence results are as follows for ASD and CON groups, respectively: \(U(18,58)=938.00,\ z=5.08, p<0.001, \eta ^2=0.34\) and \(U(20,58)=1152.00,\ z=6.54, p<0.001, \eta ^2=0.55.\)

  • More interestingly, there is a significant difference \((\eta ^2=0.04)\) between \(2{\mathrm{nd}}\) childhood and adolescence in CON group \((U(52,58)=1876.00,\ z=2.20, p=0.02)\), yet for ASD group there is no such a change in energy level of the whole brain network from \(2{\mathrm{nd}}\) childhood to adolescence, that is, energy of the whole brain network remains statistically unchanged after \(1{\mathrm{st}}\) childhood.

For the Yeo sub-networks, during \(1{\mathrm{st}}\) childhood a Mann-Whitney U test have indicated that energy of the SN (A) is greater for CON group than ASD group \((U(20,18)=104.00,\ z=-2.22, p=0.02, \eta ^2=0.13)\). Later during adolescence, there is also a significant difference on the energy level of the SN (A) between ASD and CON groups \((U(58,58)=1280.00,\ z= -2.22, p=0.02, \eta ^2=0.04\); Fig. 3B). Additionally, there is a practically significant difference \((\eta ^2=0.14)\) on the energy level of the DMN (B) between ASD and CON groups during \(2{\mathrm{nd}}\) childhood \((U(52,52)=758.00,\ z= -3.86, p<0.001\); Fig. 3C).

Figure 4
figure 4

Pearson’s correlation coefficients between network measures and behavioral scores. (A–C) During \(1{\mathrm{st}}\) childhood, ADOS-G severity and total scores show negative associations with the whole brain energy. As like, energy of the SN (A) is negatively associated with ADOS stereotype behaviors. (D,E). During \(2{\mathrm{nd}}\) childhood, energy of the whole brain network positively correlates with ADOS-G social scores and negatively with ADI restricted behaviors. (F) For adolescents energy of the whole brain network is positively correlated with ADI restricted behaviors. SN, the salience/ventral attention network; ADI, the autism diagnostic interview; ADOS, the autism diagnostic observation schedule; ADOS-G, the autism diagnostic observation schedule-generic.

Full size image

The brain-behavior relationship: correlation of energy levels with behavioral scores

Ultimately, to investigate whether energy of the whole brain network and Yeo sub-networks are associated with the clinical symptoms of ASD we have measured the correlation between the two (Fig. 4). On the x axis, we have scores from modules defined in three well-known instruments for diagnosing and assessing ASD, that is, the autism diagnostic interview (ADI), the autism diagnostic observation schedule (ADOS) and the autism diagnostic observation schedule-generic (ADOS-G). From these tests, score from some of the modules were available and among them those modules have been chosen that had less missing values. On the y axis, we have energy levels of the brain (sub) networks.

As results of two sided Pearson’s correlation between behavioral scores and energy levels indicate, during \(1{\mathrm{st}}\) childhood there was a significant negative association between energy of the whole brain network and ADOS-G severity \((r(18)=-0.49, p=0.03)\) (Fig. 4A), this holds for the ADOS-G total scores as well \((r(18)=-0.59, p=0.009)\) (Fig. 4B). Moreover, during \(1{\mathrm{st}}\) childhood energy of the SN (A) was negatively correlated with ADOS stereotype behaviors \((r(13)=-0.66, p=0.01)\) (Fig. 4C). Furthermore, during \(2{\mathrm{nd}}\) childhood and adolescence there were significant correlations between energy of the whole brain network and the following behavioral modules:

  • There was a positive correlation with ADOS-G social scores \((r(51)=0.30, p=0.03)\) during \(2{\mathrm{nd}}\) childhood (Fig. 4D), yet for adolescence the correlation with ADOS communication was significantly negative \((r(32)=-0.32, p=0.04)\). It is worth mentioning that the ADOS communication scores for adolescents were from both module-3 and -4. Moreover, the correlations of each module separately with the whole brain energy did not reach significance, possibly due to missing data. Thus, results regarding this correlation are reported in Supplementary Fig. S3.

  • There was a negative correlation during \(2{\mathrm{nd}}\) childhood with ADI restricted and repetitive behaviors \((r(49)=-0.30, p=0.03)\) (Fig. 4E), which was a significantly positive association during adolescence \((r(34)=0.30, p=0.04)\) (Fig. 4F).

Discussion

The current study has analyzed the weighted signed brain networks of ASD and CON groups in the context of structural balance theory (SBT) during development. Although analyzing pair interactions between brain regions has revealed fundamental network properties in the last decades, yet in exploring a complex system such as the brain questioning the unavoidable impact of the interactions each of the two regions has with a third region within triadic interrelationships seems plausible. In other words, the crucial role that weighted signed triadic interactions play in the organization of real-world complex networks which have been widely accepted, is worth considering. Accordingly, the current study has demonstrated the following results:

Evidence for the strong formulation of Heider’s balance theory in the brain networks (Table 1). Our findings have provided first-time empirical evidence for Heider’s balance theory on the ASD and healthy brain networks during development. To be specific, we have found that on average for both ASD and CON groups, from \(1{\mathrm{st}}\) childhood all the way to adolescence, balanced triads are heavily overrepresented in the brain networks while unbalanced triads are underrepresented—hypothesis #1. Our result on the group-level interestingly supports the strong notion of structural balance theory, which states that both unbalanced triads, \(T_2\) and \(T_0\), are underrepresented in real networks. This is while in social networks, the weak formulation of structural balance has been reported, which states that only unbalanced triads of type \(T_2\) are underrepresented in real signed networks12,13. While unbalanced triads of type \(T_0\) are in some cases (such as the Epinions and Wikipedia) overrepresented and in some although underrepresented yet to a lesser degree than the other three types (such as the Slashdot)12. These results are consistent with the fact that on the level of each participant, we have found that there is a small percentage of individuals within each group in the networks of which \(T_0\) triads are overrepresented. However, these percentages are small. Thus, the overall underrepresentation of \(T_0\) triads seems to be the case in the brain networks.

Lower frequency of \({{T}}_{{1}}\) during \({{1}}{{{\mathrm{st}}}}\) childhood, higher \({{{T}}_{{2}}}\) during \(2{\mathrm{nd}}\) childhood and lower \({{T}}_{{{1}}}\) and \({T}_{{{2}}}\) during adolescence in ASD (Fig. 1). Regarding practical differences (\(\eta ^2 \ge 0.14\)) in the mean frequency of different types of triads during development, we have observed the following results: (1) During \(1{\mathrm{st}}\) childhood, the mean frequency of \(T_1\) triads is significantly greater in CON group. (2) During \(2{\mathrm{nd}}\) childhood, the mean frequency of \(T_2\) triads is significantly greater in ASD group. (3) During adolescence, the mean frequency of both \(T_1\) and \(T_2\) triads are greater in CON group. As results depict, the frequency of \(T_3\) triads in which all three links are positive, and \(T_0\) triads that are consisted of three negative links, are not of much practical concern when studying differences between the brain network of ASD compared to healthy controls. This can pave the way for a candidate explanation that for the ASD network compared to the CON network regarding links/triads frequency, it is the interplay of negative and positive links that gives rise to the practical differences between ASD and healthy individuals—hypothesis #2, and not solely the hypo (hyper) connectivity of one type of links independent from the other type and detached from the setting they belong to. Furthermore, the frequency of triads mentioned above in ASD group being lower during \(1{\mathrm{st}}\) childhood (mean age 8) and adolescence (mean age 16), yet higher during \(2{\mathrm{nd}}\) childhood (mean age 11) reflects the hypo (hyper) functional connectivity in the backbone of networks which is pretty much in line with previous studies besides many different controversies in this regard41,42,43. These changes in the pairwise functional connectivity are as well consistent with studies that have suggested critical developmental shift during the time of puberty, which is typically between 9 and 13 years42,44.

Lower participation of brain regions in \({{T}}_{{{1}}}\) and \({{{T}}_{{0}}}\) triads with high energies in ASD during \({1}{{{\mathrm{st}}}}\) childhood (Fig. 2). As results indicate, the energy distributions of \(T_1\) and \(T_0\) triads in ASD group lag behind CON group during \(1{\mathrm{st}}\) childhood. That is, we have observed a threshold in the energy distributions of \(T_1\) and \(T_0\) triads above which node participation is zero for ASD group, yet nonzero for CON group. In other words, the ASD network lacks high energy \(T_1\) and \(T_0\) triads compared to healthy individuals during \(1{\mathrm{st}}\) childhood. This can be interpreted according to the role of \(T_1\) and \(T_0\) triads in the organization of networks from the perspective of SBT, which is to provide connected modularity. Specifically, when a network has \(T_1\) or \(T_0\) triads on its stable state, we expect it to have different groups of nodes (the so-called modules) with negative links connecting them to each other. This may go quite consistently with the theory of functional segregation and integration in the brain networks2. According to this theory, functional segregation corresponds to the presence of specialized modules or clusters within the brain network. In comparison, functional integration in the brain is the potential to combine specific information from local distributed brain regions. The Absence of high energy \(T_1\) and \(T_0\) triads in the ASD network (compared to healthy individuals) may provide us with an alternative explanation for the reduced functional integration and segregation in ASD—hypothesis #3. This finding is well consistent with previous studies that have proposed the reduced integration and segregation of information within the large-scale brain networks in ASD43.

Furthermore, as illustrated in Supplementary Fig. S2, different regions from various Yeo sub-networks are present in high energy \(T_1\) and \(T_0\) triads in the CON network yet absent in the ASD network, among which regions from three sub-networks are dominant: the DMN, the SN, and the central executive network (CEN). Most involved regions in high energy \(T_1\) and \(T_0\) triads from the SN are the insula, the frontal opercular (FrOper) and the parietal opercular (ParOper), from the DMN the prefrontal cortex (PFC), the precuneus and posterior cingulate cortex (PCUN/PCC) as well as the inferior parietal lobule (IPL), and from the CEN mostly regions from the intraparietal sulcus (IPS). Interaction between these three sub-networks, known as a triple network model of the brain, has been recently proposed to underlie a wide range of disorders, including ASD45. Zero participation of important regions of the DMN, the SN and the CEN in high energy \(T_1\) and \(T_0\) triads in the ASD network compered to the CON network, provides evidence for the triple network model from a perspective of SBT. That is to say, high energy \(T_1\) and \(T_0\) triads provide a needed structure for a consistent and reliable three-way interaction between these three sub-networks which seems to be missing in ASD group.

Narrow scale of change in energy of the whole brain network in ASD during development (Fig. 3A). Our results have revealed that from \(1{\mathrm{st}}\) childhood to adolescence, there is an overall increase in energy of the whole-brain networks in both ASD and CON groups. We can thus hypothesize that increase of the whole-brain energy with age provides networks with the necessary structure to accommodate more dynamism needed for higher cognitive abilities, which are also increasing with age. Yet as results have revealed, the pattern of this increase from \(1{\mathrm{st}}\) childhood to \(2{\mathrm{nd}}\) childhood and from there to adolescence is significantly different between ASD and CON groups. That is, while networks in CON group start from being more balance during \(1{\mathrm{st}}\) childhood and gradually gain more energy during development, networks in ASD group start with higher mean energy during \(1{\mathrm{st}}\) childhood and then seem to be frozen after \(2{\mathrm{nd}}\) childhood. In other words, in ASD group the increasing change in energy of the whole-brain network is missed from \(2{\mathrm{nd}}\) childhood to adolescence, while it is significant in CON group.

It is worth mentioning that the whole-brain ASD network during \(1{\mathrm{st}}\) childhood having significantly higher energy than the CON network is due to the ASD network having practically less \(T_1\) triads while the frequency of unbalanced triads is nearly the same in both networks (Fig. 1B). As both groups have nearly equal number of unbalanced triads, the higher energy of the ASD network during \(1{\mathrm{st}}\) childhood cannot be assigned to more dynamism, but to the less functional integration and segregation. This is because, as discussed above, \(T_1\) triads provide needed structure for the integration between local segregated modules in networks. The severity and total ASD symptoms, as measured with ADOS-G, seem to confirm these results as well. That is, less energy during \(1{\mathrm{st}}\) childhood is associated with the severity of ASD. Altogether, our results have shown that changes in the energy of the whole-brain ASD networks during development are limited to a narrower band compared to CON group—hypothesis #4.1. Moreover, during \(2{\mathrm{nd}}\) childhood and adolescence, although the total energy levels are not different between the ASD and CON networks, the underlying triadic interactions that give birth to these final energies are different (as has been discussed previously). Finally, during \(2{\mathrm{nd}}\) childhood the whole brain energy, as measured in the context of SBT, is positively associated with ADOS-G social scores and negatively with restricted and repetitive behaviors. While during adolescence, the whole-brain energy is positively correlated with restricted and repetitive behaviors as measured using ADI.

Less energy of the SN and the DMN in ASD and its association with stereotype or repetitive behaviors (Fig. 3B,C). We have observed that generally during development, the SN (A) and the DMN (B) are more balanced, have less energy, in the ASD than CON networks. In other words, in these two sub-networks in ASD group unbalanced triads are in minority compared to balanced triads. This is while unbalanced triads are known to be sources of dynamism through injecting energy into networks. That is to say, unbalanced triads are structures that excite a network towards change, unlike balanced triads that drive a network back to its stable states, that is, the minimum level of energy. Thus, unbalanced triads although underrepresented compared to chance, are playing a crucial role in healthy networks. As Heider himself stated, in a healthy community there may be a tendency to leave the balanced comfortable equilibrium to seek new horizons17. In other words, unbalanced triads are necessary if a community is to leave its comfort zone towards new experiments, that is, growth. As our results indicate, the mean energy of the SN (A) is nearly practically greater in the CON network compared to the ASD network during \(1{\mathrm{st}}\) childhood. The SN (A) in the Schaefer atlas46 includes eight regions from the insula (Fig. 3D), namely, insula left 1 \((-38, 2, -4)\), insula left 2 \((-40, -14, -2)\), insula left 3 \((-32, 18, 8)\), insula left 4 \((-36, 4, 10)\), insula right 1 \((40, 6, -16)\), insula right 2 \((40, 8, -2)\), insula right 3 \((40, -10, -4)\) and insula right 4 \((40, -2, 6)\). Classically, the insula has been considered a limbic region, yet recent network neuroscience studies have suggested its vital role in detecting novel salient stimuli across multiple modalities. The SN itself is known to be involved in attentional processes and dynamic switching between task-positive (CEN) and task-negative (DMN) processes47. Moreover, the ability to detect and attend from one event to the other flexibly has been associated with the SN’s well functioning. It is noteworthy that our finding shows a negative association between energy of the SN (A) and stereotype behaviors during \(1{\mathrm{st}}\) childhood. To be specific, greater energy of the SN (A) is associated with reduced stereotype behaviors during \(1{\mathrm{st}}\) childhood. Altogether, it seems that the SN (A) having less energy in the ASD network is reflecting the difficulty in dynamic switching, which is manifested in form of increased repetitive and restricted behaviors—hypothesis #4.2.

Furthermore, our results have unveiled that energy of the DMN (B) is practically greater in CON group than in ASD group during \(2{\mathrm{nd}}\) childhood. The DMN (B) in the Schaefer atlas46 heavily relies on the regions from the prefrontal cortex (PFC) on both hemispheres, such as the ventral PFC, the dorsal PFC as well as the lateral PFC (Fig. 3E). The role of the PFC on both social impairment and restricted behaviors in ASD has been suggested, specifically the proper level of dopamine in the PFC seems to be necessary for jumping out of repetitive behaviors48. For example, when an antagonist of dopamine was injected into the PFC of rats, it induced ipsiversive circling49. From a network point of view, it is known that systems with more energy are more probable to explore different possible states, while for networks with less energy the chance to be trapped in one of the local minima is higher. Less energy of the DMN (B), which includes several regions from the PFC, seems to be a candidate explanation for a deficiency to move from one local minimum to another in the ASD network, which may be expressing itself as increased repetitive behaviors—hypothesis #4.3. The general pattern of association between energy and restricted behaviors, the more the energy of the DMN (B) the less repetitive behaviors, nearly holds for the DMN (B) during \(2{\mathrm{nd}}\) childhood as well. Although it did not reach significance \((p=0.10, r = -0.23)\) due to missing values, thus further investigations can be helpful in this regard.

All things considered, the current study proposes SBT as a promising context to understand alterations found in the brain networks of individuals with ASD compared to healthy controls. A limitation of this study is that the data analyzed here are cross-sectional due to the small sample size of open longitudinal fMRI dataset. Thus, further investigations based on longitudinal data can sure be insightful. According to the results obtained here, studying triadic interactions and the consequent structural balance of the weighted signed brain networks, and the role that their disruption may play on the complex neurodevelopmental disorders seems to be of value. Furthermore, in addition to SBT for studying higher-order interactions in the brain networks methods such as the topological data analysis can be insightful.

Methods

Participants

T1-weighted MRI and resting-state fMRI along with the corresponding demographic data of 311 individuals, 152 with ASD, and 159 healthy volunteers (CON), aggregated from multiple sites in the publicly available Autism Brain Imaging Data Exchange (ABIDE I Preprocessed repository)50, were processed and analyzed in this study (Fig. 5A). Inclusion criteria for sites were to have (1) Similar repetition time (\(TR=2\ ms\)), in order to limit the multi-site variability and (2) Instructed participants to relax with their eyes open while a cross was projected on a screen. The reason for this choice was that, as has been reported previously51, the reliability of resting-state analysis is higher when subjects are instructed to rest with their eyes open compared to the eyes closed, due to the possibility of falling sleep when eyes closed. These two criteria resulted in five sites, namely, New York University Langone Medical Center (NYU), San Diego State University (SDSU), University of Michigan (both samples: UM-1, UM-2), University of Utah School of Medicine (USM) and Yale Child Study Center (Yale). Furthermore, it is worth mentioning that in the NYU dataset 14% of subjects were eyes closed, which have been excluded prior to preprocessing.

Table 3 Demographic data. p values are according to paired t-test for full/verbal/perform IQ, Mann-Whitney U test for age, mean FD and Chi-square test for gender and handedness.
Full size table

Inclusion as a participant in ASD group required receiving ASD diagnosis based on the autism diagnostic observation schedule-generic (ADOS-G) and an expert clinical opinion upon DSM-IV criteria. Individuals in CON group must have no history of psychiatric disorders in themselves or their first-degree relatives. Participants in both groups should have no prior or concurrent diagnosis of neurological disorders (e.g., epilepsy, meningitis, encephalitis, head trauma, or seizures). Being fully verbal and IQ > 70, as has been measured via WASI and/or WASI-IV, were also required for all. Moreover, to met the assumptions of ANCOVA we have discarded 53 participants due to being outliers (Fig. 5D), details of which are provided in statistical analysis. Altogether, the final number of participants has been 258 individuals, 128 with ASD and 130 healthy controls. We have further divided these individuals based on age, that is, 6–13 years old as children and 13–18 years old as adolescents. Additionally, as middle childhood is a crucial period during development, we have considered the opening years of middle childhood and the closing years separately. We have referred to the former as \(1{\mathrm{st}}\) childhood (6–9 years old) and the latter as \(2{\mathrm{nd}}\) childhood (9–13 years old), that are age, gender, handedness, IQ and head motion matched (Table 3).

Approval for human experiments

For each of the included sites following Institutional Review Boards (IRB) and ethics committees have approved the experiments, and data acquisition procedures were in accordance to their guidelines and regulations: (1) IRB Operations at NYU, (2) San Diego State University’s Human Research Protection Program (3) IRB of the University of Michigan Medical School (4) The University of Utah’s IRB (5) Yale University’s IRB. Furthermore, all experiments were in accordance with HIPAA guidelines and 1000 Functional Connectomics Project / INDI protocols, that is, all data were fully anonymized with no protected health information and legal guardians of all participants have signed informed consent.

Figure 5
figure 5

Graphical abstract (A) After preprocessing, 311 individuals (159 and 152 in CON and ASD groups, respectively) have been included. (B) Brain construction was based on Schaefer-400 atlas as nodes and significant Pearson’s correlations as edges. (C,E) Procedures as has been defined in the context of SBT. D Statistical analysis for comparing group means. For head motion and site information covariates, frame-wise displacement and site codes have been used, respectively. As a results of applying ANCOVA’s assumption 53 individuals were discarded. After excluding these individuals, network’s edges have been redefined based on the remaining 258 participants. Step C has been reapplied on the final networks. BrainNet Viewer 1.6340 (nitrc.org/projects/bnv) has been used for visualization of the brain networks. n, number of individuals; \(|T_i|\), total number of \(T_i\); \(p(T_i)\), fraction of \(T_i\); \(p_0(T_i)\), fraction of \(T_i\) in the null model; \(s(T_i)\), the amount of surprise; SBT, structural balance theory; NYU, New York University; SDSU, San Diego State University; UM, University of Michigan; USM, University of Utah School of Medicine; YALE, Yale Child Study Center; CON, Control; ASD, autism spectrum disorder.

Full size image

Data preprocessing

We have used the preprocessed version of ABIDE I neuro-imaging data that is provided by the Preprocessed Connectomes Project (PCP) and is public at https://preprocessed-connectomes-project.org/abide/. Among the four available preprocessing pipelines in ABIDE I preprocessed, we have used data from Configurable Pipeline for the Analysis of Connectomes (CPAC). The following preprocessing steps were performed: (1) Basic processing that includes: Realignment, slice timing correction, registration (Boundary-based rigid body (BBR) for functional to anatomical and ANTs for anatomical to standard template), intensity normalization (4D Global mean = 1000). (2) Nuisance signal removal to clean confounding variations resulting from head motion, physiological processes (such as heart beat and respiration) and scanner drifts that have included regressing: 24 head motion parameters (six standard parameters R= [ X Y Z pitch yaw roll], along with their squares \(R^2\) and temporal derivatives \(R'\)), tissue signals using component-based noise correction method (CompCor52), linear and quadratic scanner’s low frequency drifts. (3) As it has been reported53, even after regressing head motion parameters resting-state functional connectivity may still be affected by motion thus scans with frame-wise displacement \(> 0.5\) mm and global BOLD signal changes \(> 3\) SD (the conservative option in the default preprocessing steps of Conn toolbox52) were flagged and scrubbed. Participants with at least 150 valid scans (\(\sim 5\) min or more) have been included. Moreover, results of the Mann-Whitney U test have confirmed that the difference between ASD and CON groups was not statistically significant regarding head motion as has been quantified by mean frame-wise displacement (Table 3). (4) We have not regressed the global signal out, that is, the global mean signal was not included with nuisance variable regression. Moreover, spatial smoothing using a Gaussian kernel of 6 mm FWHM has been conducted. We have visually inspected the de-noised data using Conn QA plots52. Band-pass filtering (Slow-4: 0.027–0.073 HZ) was applied after preprocessing and during the de-noising step because it has been suggested that resting-state brain networks derived in this frequency exhibit grater reliability54. (5) Finally, to harmonized data across the five sites MatLab implementation of the ComBat method55 have been applied through following steps: (1) For each \(i = 1\) to 258 subjects we have created the \(tril_i\) \((79{,}800 \times 1)\) matrix with the lower triangle elements of its \((400 \times 400)\) functional connectivity matrix. Then, we have designed the dat input matrix by concatenating tril matrices column-wise, resulting in a \((79{,}800 \times 258)\) matrix with columns regarding subjects and rows corresponding to their functional connectivities. (2) We have coded different sites from 1 to 5 corresponding to NYU, SDSU, UM, USM and YALE, respectively. Then, we have set the batch input to be a \((258 \times 1)\) matrix, whose element (i, 1) corresponds to subject’s i site code. (3) For the mod input, which aims to preserve the original biological variations within data while harmonizing the effect of site, we have designed a \((258 \times 2)\) matrix with first and second columns being age and group (0 for ASD, 1 for CON), respectively. (4) We have conducted a series of Kruskal–Wallis tests (FDR corrected p values at \(5\%\) significance) to investigate site effect on functional connectivities between all pair regions of interest before and after applying ComBat harmonization. As results have shown before the harmonization \(1.37 \%\) of all functional connectivities were significantly different across the five sites, which have decreased to zero after harmonization.

Construction of the brain networks from resting-state fMRI data

In this study, we have defined the whole brain networks’ nodes based on 400 regions of interest as introduced by the 2-mm Schaefer-400 atlas46. Moreover, for the analysis of brain sub-networks, 17 Yeo parcellation was applied as defined in the Schaefer-400 atlas as well. Specifically, on this atlas the DMN (B) and the SN (A) have 32 and 34 regions of interest, respectively. For links of the whole brain network, we have assessed functional connectivity using Pearson’s correlation between all pairs of 400 regions of interest for each participant across the full length of the resting-state Blood Oxygen Level-Dependent (BOLD) time series, creating 311 weighted correlation matrices \(r(400 \times 400)\). The same calculations have been carried for the 17 Yeo functional sub-networks. Next, as Fisher’s z-transformation stabilizes the variance for better use in statistical testings, we have converted Pearson’s r to the normally distributed z variables. Then, we have performed two-tailed one sample t-tests on Fisher’s Z-transformed correlation coefficients (z), to check whether correlation coefficients are spurious or significantly different from zero1. Specifically, for each edge (i, j) we have performed a one-sample t-test based on the distribution of z values between i and j across the group. Then, for every (i, j) we have set it to zero if it did not pass the \(5\%\) significance. Of course, a multiple comparison correction was necessary to account for multiple comparisons, thus Bonferroni corrected p values were used. We applied no further manual threshold on networks to avoid its inevitable effects on network parameters. This procedure has been conducted on MatLab and Statistics Toolbox Release 2017b and the brain networks have been visualized using BrainNet Viewer version 1.6340 (Fig. 5B).

Structural balance theory

Our approach to studying the brain networks is structural balance theory (SBT) that provides a framework to go beyond the assumption that pair interactions are independent from each other in signed networks, through analyzing triadic interactions18,32. Similar to the graph-theoretical framework, SBT has also been developed to investigate the organizational properties of complex networks31. However, unlike graph theory, which has been formulated initially to model dyadic relations between information units, SBT goes beyond pairwise interactions and study triadic interactions. Specifically, SBT proposes that a network evolves in a direction that leads to the minimum level of tension and frustration between triadic interrelationships56,57. When applying graph theoretical analysis on the brain networks, one defines a graph G(V, E) with a set of nodes \(V=\{v_n; \ n=1,2,\ldots ,N\}\) and estimates a measure of association between each pair regions of interest as a set of weighted links \(E=\{w_{xy} \mid w_{xy}>0\) if \(v_x\) is correlated with \(v_y\), and \({w_{xy}}<0\) if \(v_x\) is anti-correlated with \(v_y \}\). Yet, in the context of SBT as network’s primary building blocks one defines sets of triads as follows14(Fig. 6):

$$\begin{aligned} T_i=\{\Delta _{xyz} \ | \ v_x, \ v_y\text {\ and } \ v_z \ \text {are (anti) associated} \}, \text {where} \ i= {\left\{ \begin{array}{ll} \text {3: strongly balanced,} &{} \quad \text {if } w_{xy},w_{yz}, w_{zx}>0 \\ \text {2: strongly unbalanced,} &{} \quad \text {if } w_{xy},w_{yz}>0 \ \text {and} \ w_{zx}<0 \\ \text {1: weakly balanced,} &{} \quad \text {if } w_{xy},w_{yz}<0 \ \text {and} \ w_{zx}>0 \\ \text {0: weakly unbalanced,} &{} \quad \text {if } w_{xy},w_{yz}, w_{zx}<0 \\ \end{array}\right. } \end{aligned}$$
(1)
Figure 6
figure 6

Four types of triads as has been defined in SBT. The number of negative links is even for balanced triads, and odd for unbalanced triads. Furthermore, being strong or weak refers to how much frustration a presence of a triad imposes on a network, with \(T_2\) triads injecting more frustration than \(T_0\) triads, and \(T_3\) being more stable than \(T_1\) triads14. w, weights of links; \(U(T_i)\), energy of \(T_i\).

Full size image

This way, information about the organization of network that cannot be detected on the level of pair connections would get a chance to be revealed. A well-known analogy for this definition is that positive (negative) links are considered friendship (enmity) relations, respectively. Then, balanced triads are those that fulfill following axioms: (1) A friend of my friend is my friend, that is, strongly balanced triad, \(T_3\) (2) An enemy of my enemy is my friend, a friend of my enemy is my enemy and an enemy of my friend is my enemy, that is, weakly balanced triads, \(T_1\), otherwise we have unbalanced triads, \(T_2\) and \(T_0\). A network satisfies the structural balance property in two ways: (1) Either all its triads are balanced (known as heaven), or (2) It is divided into sub-networks such that within each sub-network positive links are present, while there are negative links between sub-networks (known as bipolar). In this study, we have taken the correlations (anti-correlations) between each two regions of interest in the brain network as positive (negative) relations, respectively.

It can be seen from Eq. 1 that, for a balanced interaction the product of its edges is a positive number, whereas for an unbalanced interaction this product is negative. As has been previously proposed58, if one sums the negative of these products and divides it by the total number of ternary interactions, a quantity U(N) or energy of a network would be obtained. Energy represents the extent to which a network is structurally balanced. Explicitly,

$$\begin{aligned} U(N)= - \frac{1}{\Delta } \sum _{x<y<z} w_{xy} w_{yz} w_{zx} \end{aligned}$$
(2)

in which N is the number of nodes and \(\Delta\) is the total number of triads in the network. Applying this measure on a fully balanced network results in a lowest U(N), that is, \(-1\), while a fully unbalanced network obtains the highest possible U(N), that is, \(+1\). In this study, for any given network we have:

  • Counted the number (frequency) of \(T_i\) (\(|T_i|\)), then calculated the energy of \(T_i\), i.e., \(U(T_i) = w_{xy} w_{yz} w_{zx}\), where x, y and z are \(T_i\)’s nodes and \(w_{xy}, w_{yz}\) and \(w_{zx}\) are link’s weights. Afterwards, we have estimated the distribution of \(U(T_i)\) during \(1{\mathrm{nd}}\) childhood, \(2{\mathrm{nd}}\) childhood and adolescence for both ASD and CON groups. Finally, we have computed total energy of the network, U(N) (Fig. 5C).

  • As has been previously proposed12, we have created a null model, which is a network with the exact fraction of positive (negative) signs which are randomly assigned to existing links. The purpose of this null model is to compare the empirical frequencies of \(T_i\), as in the real brain networks, with the corresponding frequencies if signs of links were generated randomly from the same distribution of positive (negative) signs. In other words, the null model represents a condition where no underlying structure (organization) directs the placement of signs, rather it is random. Thus, after generating this shuffled version of a given network we have computed \(p_0 (T_i)\), fraction of triad \(T_i\) in the null model, and compared it to \(p (T_i)\), fraction of triad \(T_i\) in the real brain network. If \(p(T_i) > p_0(T_i)\) then \(T_i\) is overrepresented, and if \(p(T_i) < p_0(T_i)\) then \(T_i\) is underrepresented. Furthermore, we have calculated the value of surprise, \(s (T_i) = {(T_i -E[T_i]) }/{\sqrt{\Delta p_0 (T_i) (1-p_0 (T_i)) }}\), in which \(E[T_i] = p_0 (T_i)\Delta\) is the expected number of triads \(T_i\) and \(\Delta\) is the total number of triads in the network (Fig. 5E).

Statistical analysis

Throughout this study, to determine the group mean differences we have conducted two-way Analysis of Covariance (ANCOVA) with group and age as independent variables while controlling for FIQ, medication, mean frame-wise displacement as head motion parameter and site information. In total, we had 24 dependent variables, namely, the frequency of positive and negative links (Table 2A), the frequency of each types of triads (Table 2B), energy of the whole-brain network (Fig. 3A), energy of each of the 17 Yeo sub-networks among which results of ANCOVA were significant for SN (A) and DMN (B) (Fig. 3B,C). For each of these 24 models, first we have checked the assumptions of ANCOVA regarding the dependent variable of that model as follows: (1) We have explored if the dependent variable have outliers across different groups of independent variables using box plots, which in total have resulted in discarding 53 out of 311 participants. (2) Results of the Kolmogorov–Smirnov tests have indicated that except for the frequency of \(T_2\) triads in ASD group during \(2{\mathrm{nd}}\) childhood and adolescence, as well as the frequency of \(T_0\) triads in ASD group during adolescence, which are only moderately skewed (\(-1 \le skewness \le +1)\), all other dependent variables in this study are normally distributed. (3) Homogeneity of variances was investigated through Levene’s Test of Equality of Error Variances. (4) The linear relation between the dependent variable and covariates was studied through scatter matrices. Additionally, when results of ANCOVA were significant we have applied the Mann-Whitney U test as the post-hoc test to determine the specific groups that are different (as shown using the asterisks above the box plots in Fig. 1 for the frequency of links and triads as dependent variable, and Fig. 3A–C, for energy of the whole brain and sub-networks as dependent variables). In case of ANCOVA, effect sizes are reported as has been readily estimated by IBM SPSS Statistics version 26, that is, \(partial \eta ^2 = {SS_{effect}}/{SS_{effect}+SS_{error}}\). For the interpretation, Cohen’s guideline was applied, i.e., \(\eta ^2\) at least 0.01:  small, 0.06:  medium and 0.14:  large effects59). Likewise, for the post-hoc tests we have calculated \(\eta ^2= {Z^2}/{n-1}\), using the Z statistics from the Mann-Whitney U test. Finally, the Kullback–Leiber Divergence between distributions P and Q in Fig. 2 has been calculated through \(\mathbb K(P||Q) + \mathbb K(Q||P) = \sum _{i}\log _2 (p_i/q_i)p_i + \sum _{i}\log _2 (q_i/p_i)q_i\).