We construct a Brain Network Model (BNM) to study the role of symmetry-breaking due to structural connectivity. The local dynamics in the network model are simulated by a 2D Neural Mass Model (NMM), derived from an ensemble of all-to-all coupled QIF neurons, and the brain states correspond to the stable fixed points of the system, encompassing its dynamical repertoire, as also confirmed by our extensive numerical simulations. We build a low dimensional representation of these brain states from the Laplacian of the connectivity matrix and their distribution is mapped in the framework of parameter-degeneracy in the system. A noise-driven exploration of the brain states is employed to validate the contours of degeneracy and the consequential flow on the low dimensional manifold. We also provide a geometric interpretation of the degeneracy and a possible expansion to include symmetry-breaking by regional heterogeneity.

Fig. 3

Eigenvectors of the Laplacian matrix. The four leading eigenvectors of the Laplacian matrix (L) are shown as surface plots. (A) The leading eigenvector \((e_0)\) corresponding to the zero eigenvalue. The coefficients of \((e_0)\) are proportional to the in-strengths of the regions. Since the network is connected, \(\lambda _0 = 0\) is unique and the coefficients of \(e_0\) are all of the same sign, indicating that the largest structural organization in the network is the network itself. (B) The second eigenvector \((e_1)\), corresponding to the first non-zero eigenvalue \(\lambda _1 = 0.1016\). \(e_1\) extracts the largest partition in the connectome, that is the separation of the hemispheres. The left and right hemispheres are separated by positive and negative coefficients of \(e_1\) respectively. We ignore \(e_2\) as it only highlights the cerebellum and gives no further information. (C) The next eigenvector \(e_3\), corresponding to \(\lambda _3 = 0.1974\), brings out the anterior-posterior separation. The anterior and posterior regions are partitioned, respectively, by negative and positive coefficients of \(e_3\). (D) The next eigenvector \(e_4\) corresponds to \(\lambda _4 = 0.2836\), and it extracts the hemispheric and anterior-posterior partitions simultaneously. Further eigenvectors bring out increasingly smaller structural partitions

3.1 Brain network model

The BNM is constructed as described in Methods-1. The local dynamics are modeled by the NMM (Montbrió et al., 2015), which describes the evolutions of the firing rate (r) and membrane potential (v). This model has already been used in virtual brains to explain the emergence of various dynamical features relevant for neuroimaging at different time-scales (Rabuffo et al., 2021; Fousek et al., 2024), as well as for explaining the compensation effects during aging (Lavanga et al., 2023) and the change in brain’s dynamical working point during different stages of consciousness (Breyton et al., 2025). Moreover, the model was used to explore different strategies for tackling the non-identifiability in the pathology-linked degradation of brain’s connections (Hashemi et al., 2024). The parameters of the single network node are chosen such that the system exhibits bistability, with a stable focus for a high firing rate and a stable node for a low firing rate. These are used to define our up and down states respectively. The network dynamics arises from interacting bistable nodes, where the oscillatory dynamics of the up-state provide the fast oscillations and set the characteristic frequency of the system. The slow oscillations in the network arise from noise-induced transitions between the up and down states (via the unstable fixed point) of the network nodes. A representative quantifier of the system behavior thus becomes the occupancy of the up-state of the network. In Fig. 1A, we show the nullclines (in red), the basin of attraction (in black), and the phase-flows (as arrows). Sample trajectories of the system, in the vicinity of both the stable focus and node, are also shown in blue.

3.2 Laplacian eigenbasis

The Laplacian form of the structural connectivity matrix plays a key role in the definition of a low-dimensional representation of brain states. Decomposing the symmetric normalized form of the Laplacian matrix into its eigen-basis extracts the hierarchy of spatial organizations, that are manifested as structural partitions. Surface plots of the leading vectors are plotted and discussed in Fig. 3. We exploit this hierarchy to define a low-dimensional embedding and extract the resting state manifold. The symmetric normalized form also ensures that the eigenvalues are real, non-negative and bound in the interval [0, 2], expanding the comparability of results across various connectomes.

Fig. 4

Defining the RS manifold. (A) Comparison of the distributions of compositions of an ensemble of 50000 initial conditions and stable network fixed points for \(g = 0.55\). (B) Plot of composition versus stability of the stable fixed points, showing the decrease in stability at higher compositions. (C) Projection of the numerically calculated stable states, relative to the initial conditions, onto the two leading eigenvectors of L. The initial conditions and stable fixed point states are colored black and red respectively. (D) Organization of the network states, by composition, on the resting state manifold

3.3 Fixed points of the connectome

The stable fixed points of the network \((\ = \)\), that represent the brain states, are numerically calculated. The details of this calculation are given in the methods section. While each network state is a 2N-dimensional vector, there is redundancy in the sense that, for each region, both variables exhibit qualitatively similar behavior. Hence we limit our quantification of system dynamics to the dynamics of the membrane potential without any loss of generality.

Based on this, the network states are quantified using two scalar metrics, namely composition and (in)stability. The composition of a network state is defined as the number of regions in the up-state and takes integer values between 0 and N, Fig. 4(A). The instability \((\lambda )\) of a state is defined as the magnitude of the real-part of the largest eigenvalue of the Jacobian of the system, evaluated at the fixed point, Fig. 4(B). A negative value implies that the state is fully stable while a positive value indicates that the state is partially unstable in at least one component. The sampling of the possible solutions is implemented using Newton-Raphson algorithm. In short, from a small region around the symmetric (uncoupled) case fixed points, we uniformly sample random initial conditions (50 000 in total), black dots in Fig. 4(C), from which we calculate the convergence to fixed points (red dots). The results show that for significant values of coupling strength, states with higher compositions easily lose stability, Fig. 4(A,B). The global coupling strength (g), which regulates the cumulative effect of the connectome, quantifies the possible set of dynamical states, which in the symmetric case (\(g=0\)) spans \(2^N\) states. An increased role of the structure constrains the dynamical repertoire of the BNM. The contraction of the sample-space can be explained by understanding the association between composition and stability of brain states.

3.4 Structural asymmetry constrains the dynamical repertoire of the BNM

Consider a network state with composition C, that has a corresponding instability \(\lambda _C\). When the network is at a fixed point, each region (component) can be treated independently from all other regions, which implies that one unstable region translates to an unstable brain state. The dynamics the fixed point of the \(i^\) region in a brain state are given by Eq. (3). Assuming that the network has an average in-strength k, the interaction term in Eq. (3) can be written as \(g k \sum _j r^*_j\). The assumption of an average in-strength is valid because the in-strengths are quite narrowly distributed around an average value. The fixed-point states of individual regions \((r^*_j)\) can be coarse-grained into up (\(r^*_\)) and down (\(r^*_\)) states. Then, for a state with composition C, the interaction term becomes \(g k C r^*_ (N-C)r^*_\). This can be further simplified to gkC because the values of \(r^*_ \approx 0\) and \(r^*_ \approx 1\), as seen in Fig. 5 (A). The interaction term, for each region, can now be treated as an input current term that depends directly on the composition C of the global state. Going back to Fig. 5 (B), we now have the relationship between the composition of a brain state and the dispersion of individual regions from their respective symmetric-case fixed-point values. More specifically, the deviations grow proportionally with the composition, and this is driven by the global coupling strength, Fig. 5 (C).

Next, if we consider the stability of the \(j^\) region in a brain state, the real-part of the eigenvalue of the Jacobian is equal to \(2v_j^*\). From Fig. 5(A), we observe that in the up-state, \(v^* \approx 0\) but negative; while in the down-state, \(v^* \approx -2\). Therefore, between the two, the up-state is significantly less stable than the down-state. And Fig. 5(B) indicates that an increasing input pushes \(v_j^*\) closer towards 0, especially in the up-state. Hence, for any given region, a higher composition of the network state leads to higher regional inputs, thus resulting in lower stability. Between the up and down states, the up-state is more susceptible to this effect.

Therefore, if a global state has a higher composition, then there are more regions in the up-state that receive a higher input, thus lowering the regional stability more quickly compared to a state with a lower composition. This explains the observations for \(g>0\), in Fig. 5(B), where the states that quickly lose stability are those with higher compositions. For small values of g, the above effect is limited to states with very high compositions. But as g is increased, regions in the down-state get pushed to the up-state and the same process continues, and this explains the observations in Fig. 5(C).

Fig. 5

Calculation of fixed-point states. (A) Components (pairs of mean membrane voltage and firing rate) of an initial condition with randomly chosen composition and the corresponding fixed point. (B) Change in the distribution of compositions of stable fixed points of the system for increasing values of global coupling strength. (C) Decrease in the fraction of stable fixed points as a function of global coupling strength

3.5 Resting state manifold (RSM)

A low dimensional representation of brain states, that we refer to as the RSM, is constructed by projecting the mean membrane voltage component of the fixed points of the system onto the leading eigenvectors of the Laplacian matrix. The details of this process are given Methods-3. This projection nicely separates the network states and spans the manifold. The projection is not unique and others exist Runfola et al. (2025b), but it is convenient. Based on our knowledge of the eigenvectors of L (Fig. 3), we observe along the horizontal axis, a smooth left-right gradient in the composition of states (Fig. 4D) and on the vertical axis, is the separation of hemispheres.

In this new representation, it is easy to visualize how the stable fixed points of the system define a sample space or set of possible configurations for the network. In response to any perturbations (noise or stimulating currents), the system traces a trajectory in this space, as the perturbation induces a switching behavior which in turn is the consequence of symmetry-breaking by the network structure.

3.6 Degeneracy

As mentioned in the previous section, given a system of N identical nodes connected by a network W, the complex topology introduces a symmetry-breaking in the system. In the symmetric case, where the regions are either isolated or identically connected, all regions in the up-state reach the same fixed point (\(r^*_\), \(v^*_\)) and in the down-state, they all go to \((r^*_, v^*_)\). For \(g>0\), this symmetry is broken as each region receives a unique input, that is a function of the topology (in-strength) and the existing configuration of the network. Regions that receive different inputs go to slightly different fixed points and this excursion from the symmetric fixed point is directly proportional to the input. The manifestation of symmetry-breaking and the connection to the framework of degeneracy are depicted in Fig. 6, with the plot (A) showing the values of the mean membrane potential for each of the network states, and the plot (B) showing the overall distribution of these network states.

As seen from Fig. 6C, if a region in the down-state receives high-input, the new location of the fixed point moves closer to the up-state and vice-versa. This implies that each region has a different likelihood of transitioning to the other state, depending on the input at that instant. This effect of symmetry-breaking, that results in a differential response to the same perturbation, is transformed and interpreted in the framework of the degeneracy of components of the network state. In the remainder of this section, we will discuss the definition, calculation, and results pertaining to the degeneracy at the regional and global level.

Fig. 6

Leading up to Degeneracy. (A) Scatter plot of all numerically calculated network fixed points (only v variable). (B) Distribution of membrane potentials of fixed-point components of all network states. (C) Deviation of the membrane potentials of fixed-point components from the symmetric-case fixed point as a function of input. The deviations in the up/down state correspond to curves with lower/higher slopes

Degeneracy is defined at the regional level and for a specified parameter of the system and is subsequently extended to network states. Consider parameter \(\eta\), for region i. The degeneracy of i w.r.t \(\eta\) is defined as the absolute maximum excursion, relative to its initial value, that is required to switch the region from its local state. If \(\eta _\) and \(\eta _\) are the initial and final values of the parameter respectively, then the degeneracy is given by

$$\begin Degeneracy (\eta _i) = \delta \eta _i = \frac - \eta _ |}} \end$$

(5)

Although the degeneracy is defined as the absolute difference, in certain sections of this work, the sign is taken into account. This is primarily to differentiate the roles of the local states. Further, if a region is in the up-state, the value of the parameter is decreased in order to induce a transition to the down-state and vice versa. This results in upper/lower-bounds of the parameter for regions in the down/up-states. Henceforth, upper/lower-bounds and down/up-states will be used synonymously. The procedure for the calculation of degeneracy is discussed in the methods section.

Fig. 7

Degeneracy at the regional level for the chosen working point of the brain. After calculating the degeneracy w.r.t all three parameters, we study the role of local-states and the effect of in-strength. (A - C) Distributions of \(\delta \eta _i\), \(\delta J_i\) and \(\delta \Delta _i\), respectively, and the distinction based on local states of the regions. (D - F) The dependence of average degeneracy on in-strength, together with the separation of local states

3.6.1 Role of connectivity

As a first step in understanding the distributions of degeneracies, only the role of the local state is considered. In Fig. 7 (A - C), the distributions of \(\delta \eta _i\), \(\delta J_i\), and \(\delta \Delta _i\) are shown, distinguished only by the local states and independent of regional properties and global states. From Fig. 7(A), it is immediately evident that the distributions for \(\delta \eta _i\) are distinctly different for up and down states. In the down-state, regions are more likely to have higher values of \(\delta \eta _i\), whereas they are more likely to have lower degeneracy w.r.t \(\eta\) in the up-state. However, in both cases, the range of \(\delta \eta _i\) is approximately the same. The distributions of \(\delta J_i\) are qualitatively similar to those of \(\delta \eta _i\). The actual range of values are, however, well-separated for up and down states. For regions in the down-state, the distribution peaks at much higher values of \(\delta J_i\). This can be attributed to the scaling by the local state (\(Jr_i\) terms) in the network and is discussed in greater detail in the section on geometric interpretation. For degeneracy w.r.t \(\Delta\), regions in the down-states are again more likely to have higher degeneracy, while regions in the up-state have no visible distribution and are concentrated in a narrow range on the lower end.

Since the degeneracy is a manifestation of the node-to-node variability introduced by network structure, or the connectome induced asymmetry as we have referred to in the earlier works (Fousek et al., 2024; Runfola et al., 2025b; Jirsa, 2020; Woodman & Jirsa, 2013), it is important to understand exactly how it is influenced by the structural properties. The local states, together with the network structure, define the perturbations to the locations of fixed points of individual regions, relative to the symmetric (\(g = 0\)) case. These deviations significantly affect the response to perturbations, a.k.a degeneracy, and consequently the switching behavior. We focus on the relationship between average degeneracy and the in-strength of the regions. In order to do this, the values of \(\delta \eta _i\), \(\delta J_i\) and \(\delta \Delta _i\) are each averaged over all global states for all i and plotted as functions of in-strengths of the respective regions and shown in Fig. 7(D - F).

In the down-state, the mean degeneracy decreases linearly with in-strength for all 3 parameters. And for regions in the up-state, it increases linearly with increasing in-strength, except in the case of \(< \delta \Delta>\), where it takes a low value and remains independent of in-strength. Therefore, for regions with low connectivity, the local state plays an important role in defining their degeneracy and hence also their response to perturbation. Regions in the up-state are less degenerate and therefore respond more easily to perturbation while in the down-state, they have higher resilience. For regions with intermediate connectivity, the local state plays an important role in the case of \(< \delta J>\) and \(< \delta \Delta>\). For \(< \delta \eta>\), the response to perturbation cannot be predicted based on the local state as the degeneracy takes the same range of values. For regions with high connectivity, there is a greater spread in degeneracy as is observed from the increasing size of variance. Due to this overlapping of degeneracies, it is hard to make the case for a distinct role of local states, for such brain regions. In the next section, we will see how this relationship with in-strength affects the switching behavior as a response to white noise.

Fig. 8

Degeneracy at the global level. Numerically calculated stable fixed-point states are plotted on the resting state manifold and represented by each of the two definitions of global degeneracy. In A - C, each global state is represented by \(min(|\delta \eta _i|), min(|\delta J_i|)\) and \(min(|\delta \Delta _i|)\) respectively. And in D - F, the states are represented by the mean values of degeneracy w.r.t each parameter respectively. In all cases, the values are placed on the blue to red color scheme. Negative values of minimum/mean degeneracy indicate that the main contribution to the global value comes from a region(s) in the up-state and positive values correspond to contributions of regions in the down-state

3.6.2 Degeneracy on the manifold

So far, we have limited our discussion of degeneracy to the regional level. Each region has been considered in isolation, without concern for the global configuration of the network. Now, we develop two simple approaches to extend the measure of degeneracy to the global state. In both cases, we look for a scalar representation of the network. In this manner, we can study the distributions of degeneracy across the manifold. This allows us to identify subsets of stable and not-so-stable states and their locations, thus laying the foundation for predicting regions with high occupancy and possible transitions between different parts of the manifold in response to perturbations.

In the first approach, the degeneracy (w.r.t each parameter) of the global state is defined as the minimum of absolute values of the degeneracy of all regions (\(min(|\delta \eta _i|, min(|\delta J_i|,\)\(min(|\delta \Delta _i|), \forall i\)). Once the minimum value has been identified, the local state of the selected region is accounted for by assigning a \(+/-\) for down/up states (Fig. 8(A-C)). By assigning the minimum value, we identify, for each global state, the most vulnerable component that is most likely to flip its state, hence changing the global state itself. From an alternate perspective, it provides a lower limit for the amplitude of perturbation, that would be required to induce switching.

Starting with the distribution of \(min(|\delta \eta _i|)\), on either ends of the manifold, the least degenerate regions are in the up-state while in the center of the manifold, it is regions in the down-state that have the smallest degeneracy. As we start from the left-end of the manifold and move towards the center, we find that the absolute value of degeneracy increases very slowly until regions in the down-state start to appear. This is marked by the appearance of the blue-points. Over here, the regions already have higher degeneracy and this begins to decrease as we move further to the right, at which point the minimum degeneracy reaches zero. And further to the other end of the manifold, the values increase again. In case of \(min(|\delta J_i|)\) and \(min(|\delta \Delta _i|)\), the least degenerate states are all largely associated with regions in the up-state, with very little contribution from the down-states.

A second approach is to assign to the global-state, the mean value calculated over all regional degeneracies. This provides an average measure of the expected response to perturbation and is particularly relevant when there is a stochastic element to the perturbation. As we observe from Fig. 8(D-F), w.r.t all three parameters, there is a very similar pattern of degeneracy. Along the horizontal axis, the mean degeneracy is high on one side and towards the center of the manifold, it goes to zero and then starts to increase again towards the opposite end. However, on the left-side of the manifold, it is positive due to the dominance of low composition states. Similarly, the opposite side of the manifold consists of high composition states and therefore, we have high degeneracy but with negative values. For each parameter, the balance between degeneracies of up/down states occurs on different sections of the manifold.

In both approaches, it is immediately evident that the gradient of degeneracy is primarily along the axis of the leading eigenvector. Looking back on Fig. 4(D), we see that composition also follows the same pattern, indicating a possible relationship between the two. Therefore, in Fig. 9, we associate the minimum/mean degeneracy to the composition and study the connection between the two global properties. The benefit of this approach is that all discussions concerning the response to a perturbation can now be framed in terms of transitions between states of different compositions.

Fig. 9

Degeneracy versus Composition at the global level. For each global state, the minimum/mean degeneracy is plotted against the composition of the state. This is shown for \(min(|\delta \eta _i|), min(|\delta J_i|)\) and \(min(|\delta \Delta _i|)\), respectively from A - C. The unique values of degeneracy for each state are shown in blue while the average values of degeneracy of all states with a given composition are shown in black. From, D - F, the mean degeneracy w.r.t each parameter is plotted against the composition of the state. Here, we study three cases, where we consider only the regions in the up-state (red), only regions in the down-state (blue), and the entire global state (black)

From Fig. 9(A - C), we see how the minimum value of degeneracy changes with the composition of the state. On average, for low and high compositions, the degeneracy is relatively high, with the mean degeneracy in the down-states decreasing with the composition for every parameter, while showing generally opposite trends for the up-states, Fig. 9(D - F). For the latter, however, the degeneracy of the synaptic strength J and of the heterogeneity of the external currents \(\Delta\) is much less dependent on the network’s composition, than it is the case for the mean external currents, \(\eta\). For states with intermediate compositions, there is a fairly sharp dip (although it is not very pronounced in case of J). On the lower end, the average minimum degeneracy is almost constant, pointing to a limited role of perturbations in this range. If regions in a low composition (typically \(<25\)) state are perturbed, they switch to the respective opposite states, thus increasing/decreasing the composition on either side. However, the new state also has similar degeneracy, and therefore there is no net preferred direction for the change of composition. At a slightly higher range, the average minimum degeneracy decreases with increasing composition. In this case, if the composition increases, the new state has lower degeneracy and is more susceptible to perturbation. However, if the new state has a lower composition, it is more resilient and therefore, there is a net preference for decreasing the composition. Towards higher compositions, the degeneracy increases rapidly with composition. Following the earlier argument, we should expect a preferential movement towards higher compositions. Note that the mean degeneracy dependence on the composition, Fig. 9(D - F), follows similar trends as the local degeneracy dependence on the in-strength, Fig. 7(D - F), in agreement with the strong dependence of the composition on the in-strength as seen in Fig. 4(D), where the leading Laplacian eigenvector captures the in-strength.

So far, we have studied the patterns of degeneracy w.r.t all three parameters of the model. We have looked at the role of local states, the relationship to structural properties and the gradients on the manifold. This distribution of degeneracy defines the contours of stability in the sample-space of fixed-point states of the network. In the following sections, we will narrow down the choice of parameters using a geometric interpretation of degeneracy, with the ultimate goal of understanding and predicting the occupancy and likely transitions on the manifold in a stochastic environment.

3.6.3 Geometric interpretation of degeneracy

In this section, we briefly steer away from the calculation and consequences of degeneracy and focus on the underlying changes in the phase-space of individual regions. This geometric interpretation allows us to take an analytic approach to degeneracy and opens it up for more general interpretation. Starting from Eq. (1), the fixed points of the network \(\\) are calculated for a given set of parameters \((\eta _0, J_0, \Delta _0)\), by solving the following the system Eq. (3). This system has homogeneous distributions for all 3 parameters and each region goes to a unique fixed point value depending on the respective initial condition and network input. However, the local fixed-point states can be coarse-grained into \(r^*_\) and \(v^*_\). This allows us to categorize each region into up and down states, which then makes it possible to bring degeneracy into the discussion. If we consider parameter \(\eta _0\) and let \(\delta \eta _i\) be the degeneracy in \(\eta\) for region i, then Eq. (3) can be rewritten to include the degeneracy.

$$\begin \begin 0&= \frac + 2 r^*_i v^*_i \\ 0&= ^2 - \pi ^2 ^2 + J_0 r^*_i + \eta _0 + \delta \eta _i + g \sum _^N W_r^*_j \end \end$$

(6)

When the network is in a fixed-point state, each region can be treated as an isolated node with a network input term \(I^_i = g \sum _ W_r^*_j\). And the entire set of regions can be treated as a set of independent nodes receiving a heterogeneous distribution of input currents \(I_i = \delta \eta _i + I_i^\). Therefore, for the \(i^\) region, the fixed-point states can be calculated from

$$\begin \begin 0&= \frac + 2 r^*_i v^*_i \\ 0&= ^2 - \pi ^2 ^2 + J_0 r^*_i + \eta _0 + I_i \end \end$$

(7)

Similarly, if we are interested in \(J_o\), and \(\delta J_i\) is the corresponding degeneracy of the \(i^\) node, then the system of equations can be recast as

$$\begin \begin 0&= \frac + 2 r^*_i v^*_i \\ 0&= ^2 - \pi ^2 ^2 + J_0 r^*_i + \eta _0 + \delta J_i r_i^* + g \sum _^N W_r^*_j \end \end$$

(8)

where we can now analogously define an input term as \(I_i = \delta J_i r_i^* + g \sum _^N W_r^*_j\). Recasting the equations for degeneracy in \(\Delta\) results in a more complicated form, but the principle remains just as valid.

Fig. 10

Geometric interpretation of degeneracy. The dynamics of the network in the bistable regime, are plotted in the state space, for different values of input current \(I_i\). For each value of \(I_i\), we plot the nullclines (red), the basin of attraction (black) and 2 sample trajectories, starting from different initial conditions, to indicate the choice of stable fixed point

Based on Eq. (7), we can study the effect of the \(I_i\) term, as shown in Fig. 10. When \(I_i = 0\), any degeneracy w.r.t a parameter is offset by the input received from the neighborhood. If this is taken as a reference, then the node has a basin of attraction of a particular size, and the up/down states are well-separated. When \(I_i\) becomes negative, the size of the basin of attraction of the up-state shrinks rapidly, bringing the stable focus closer to the unstable fixed point. As a result, such regions are constrained to the down-state. On the other hand, when \(I_i> 0\), the size of the basin increases, and brings the stable node closer to the unstable fixed point. In this case, the regions are more likely to be locked in the up-state.

If the symmetry-breaking due to the connectome is to be recast as the degeneracy of the parameters, then \(I_i = 0\) can be treated as the case where the role of the connectome is compensated by the degeneracy of the parameter. Therefore, if a region is in the down-state and has high connectivity, then it is more likely to transition to the up-state than a region with low connectivity. And a region in the up-state, with high in-strength is less likely to transition to the down-state than a region with low in-strength. Alternatively, regions in the down-state with high in-strength have low degeneracy and those with low in-strength have higher degeneracy. And well-connected regions in the up-state have higher degeneracy than those with low connectivity.

Another point worth noting is the relevance of parameter \(\eta\). Comparing Eqs. (6) and (8), we see that the degeneracy w.r.t \(\eta\) is not dependent on the local-state, whereas, to relate \(\delta J_i\) to the in-strength, the local-state has to be distinctly accounted for. For this reason, in the following results, we will depend on \(\delta \eta\) to connect the role of the connectome with the role of degeneracy and study the stochastic exploration of the manifold.

3.7 Noise-driven BNM explores the contours of degeneracy on the RSM

In this section, we study the noise-driven dynamics of the BNM, as reflected in the simulated BOLD fMRI activity. The description of the simulation of the BNM is detailed in the methods section, while the BOLD signals were obtained by converting the simulated neural activity with the Balloon-Windkessel model (Friston et al., 2003) of The Virtual Brain (Sanz-Leon et al., 2015) and the resulting time-series were then downsampled at 0.7s, to correspond to the empirical BOLD signals from HCP (Van Essen et al., 2013).

First, we analyze the dynamics at the regional level. Then, we combine this information with our knowledge of the distribution of degeneracy on the manifold, to understand and explain the broader dynamic patterns of the manifold, specifically its occupancy. The low-frequency time-averaged signal is used to study the switching properties at the regional level and the high-resolution source signals are used to study the occupancy and exploration of the manifold. In Fig. 11, for each region, we plot the fraction of time spent in the up-state and the total number of switches versus the in-strength of the nodes. This is done for three cases: low noise, medium noise, and high noise, where the amplitude of noise is controlled by \(\sigma _\).

As seen from Fig. 11, in the low noise scenario, there is only a quasi-static effect. A majority of the regions remain their respective initial states. Noise-induced switches are few and far in-between and there is no significant time spent in the up-state. This is because the average dispersion due to noise is much lower than the average separation between the up and down states. Hence, low noise has little to no effect on the switching dynamics. Only regions with very low in-strengths show some likelihood for switching. This is further limited to those regions that are initialized in the up-state, which then switch down due to the relatively low stability/degeneracy in the up-state.

In case of high noise, the fraction of time spent in the up-state grows proportionally with in-strength. As for the number of switches, we observe an overall increase. Now, the separation between states is comparable to the amplitude of noise. Hence, switching can occur with relative ease and increases linearly for regions with low in-strength but on the higher end, it remains quite constant. Overall, low noise does not generate much activity expressed by the switching of the states, while high noise leads to activity that is much less constrained and predictable by the network asymmetry, making the activity of the single nodes more random.

Fig. 11

Switching Dynamics. The role of local connectivity on the excitability of regions is studied as a function of in-strength and calculated for one fixed point using the simulated BOLD signal. (A) The fraction of time spent by a region in the up-state is plotted as function of the in-strength of the regions. (B) Plot of number of switches as a function of the in-strength of regions. Both these metrics are plotted for three levels of noise. These two metrics, together, provide the information required to understand the likelihood of individual regions to undergo switching, which in turn lays the groundwork to understand the transition between global states

When the amplitude of the noise is at an intermediate/optimal level, we observe the most interesting behavior. The fraction of time spent in the up-state shows a sigmoidal behavior with regions with low in-strength spending very little time in the up-state while well-connected regions spent the majority of their time in the up-state. This implies that regions with low connectivity are mostly locked in the down-state while those with high connectivity are locked in the up-state. This is consistent with the results from Fig.