Determinants of collective failure in excitable networks

Sustained oscillatory rhythms are a hallmark of biological systems and manifest themselves at different timescales and levels of biological organization. Intrinsic rhythmic activity typically arises from many nonlinear regulatory mechanisms and feedback processes and is crucial for ensuring physiological functions.¹ On the tissue level, the oscillatory behavior is often regulated by large networks of interacting cells. Intercellular communication is the basis for the occurrence of synchronization or phase-locked dynamics, which, together with the responses to environmental perturbations, facilitates the collective rhythm that is essential for normal function^2–4 However, under certain conditions, such as during pathogenesis, individual cells are subjected to malfunction and become inactive. This leads to an interesting dynamical phenomenon known as aging transition.⁵ Investigating how a networked system operates dynamically, when segments of the system malfunction, is not only theoretically appealing but also of high relevance for the evaluation of the robustness of the biological network to perturbations and malfunctions.

The phenomenon of the aging transition was originally reported by Daido and Nakanishi for a system of globally coupled limit-cycle oscillators in which the fraction of inactive oscillators was progressively increased.^6,7 The notion of a non-trivial transition toward a macroscopically inactive state due to the deactivation of individual microscopic constituents has stimulated further research and studies on the aging transition have been further elaborated to locally coupled^8,9 and heterogeneous ensembles of oscillators,^10–13 as well as to specific coupling schemes^14,15 that can include time-delays^16,17 or repulsive interactions.¹⁸ Furthermore, along with the advances in network science, the scope of studying aging transitions has been extended to various network structures. In complex networks, the degrees of oscillator nodes, i.e., the number of connections to other oscillators, are heterogeneous, which not only highly affects the transition dynamics of the aging network but also puts forward different strategies with which the oscillators get inactivated.¹⁹ While, in general, the collective activity of heterogeneous networks is highly vulnerable to hub removal,²⁰ under specific circumstances such as weakly weighted networks with weak interactions, the macroscopic activity is highly fragile to targeted inactivation of low-degree nodes.^21,22 This result does not only highlight a perhaps counterintuitive discrepancy between structural and dynamical robustness of a network but also indicate that aging transition on complex networks is a multifactorial phenomenon. Specifically, further research has shown that edge weights,^20,23 the topology of the network,^24,25 heterogeneity,¹² as well as a multilayered structure of interactions^26,27 profoundly affect the deterioration of macroscopic oscillatory behavior.

While the majority of earlier works was focused solely on limit-cycle oscillators, the scope has later shifted also to populations composed of oscillatory and excitable elements. In this scenario, the silenced self-oscillatory units are poised in an excitable state so that, for example, the oscillatory behavior is lost in a saddle-node bifurcation^28–30 Such a composition of oscillatory and excitable units has many possible applications, particularly, in biology and neuroscience. In this vein, Ratas and Pyragas studied a globally coupled network of integrate-and-fire neurons, which are canonical examples of class I neurons. They reported that this setup offers rich dynamical behavior,³¹ including non-intuitive phenomena at certain parameter regions, where a decrease of active neurons can enhance macroscopic oscillations.³² Moreover, Tanaka et al.¹⁰ have shown that in networks of coupled oscillators, population heterogeneity enhances dynamical robustness of the neuronal dynamics in a network of Morris–Lecar models coupled through electrical synapses, irrespective of the type of bifurcation. Very recently, Biswas and Gupta studied aging transitions in a random network of Rulkov neurons. They have shown that the interplay between the coupling strength and the network density determines how the neuronal dynamics deteriorates as the fraction of inactive neurons increases and that under certain conditions, the transition can be abrupt.³³ Further research has revealed that a heterogeneous distribution of connections strengths between neurons,²³ the presence of mixed attractive and repulsive coupling,¹⁵ or multilayered interactions governed by ephaptic coupling³⁴ are also important determinants in shaping the dynamical robustness of the macroscopic neuronal activity.

Nonetheless, there are still many open questions and unexplored issues regarding the dynamical robustness of coupled excitable oscillators. In the present study, we, therefore, focus on studying the degradation of macroscopic neuronal activity in biologically realistic networks of heterogeneous excitable elements. To simulate interactions between oscillators, we utilize a spatially embedded network model that yields a heterogeneous architecture with small-world characteristics, as such topological structures are omnipresent in various biological systems including ensembles of excitable cells.^35,36 The dynamics of individual nodes are driven by the FitzHugh–Nagumo model as it represents a paradigmatic example of excitable dynamics. In our simulations, the oscillators are progressively placed from the oscillatory to the excitable state. We differentiate between three scenarios of aging, random and by degree, where either the high-degree nodes or the low-degree nodes are first to be made inactive. We explore how the interplay between the ratio of inactive oscillators, the coupling strength, and the distance from the bifurcation point affect the aging transition and demonstrate that the route to diminished macroscopic oscillations is characterized by non-trivial and interesting behavior. Our findings, therefore, provide new insights into how different physiological determinants can affect the robustness and fragility of realistic networks of excitable elements.

In Eq. (1), $x i$ mimics the membrane potential and $y i$ is a recovery variable, $J i$ plays the role of external stimulus, the variables $a i$ and $b i$ are proportional, respectively, to the ratio between inductance and capacitance and to the cell membrane, K is the coupling strength, and $S i$ defines the subset containing neighbors of the ith oscillator. A single uncoupled oscillator can be characterized with a threshold value $ε i= 3 a i 2 − 2 a i 2 b i − b i 2 3 a i 3 a i 2 − b i$ for parameter $J i$⁠.³⁹ For values of $J i∈[ − ε i , ε i]$⁠, the equilibrium point is unstable, and the system is in an oscillatory regime. For outside values of $J i$⁠, the equilibrium point is stable and the system falls into a fixed point. $J i<− ε i$ corresponds to a constant negative value of variable $x i$ and the system is said to be in an excitable state, while for $J i> ε i$⁠, the system is in excitation block with a constant positive value of $x i$⁠.⁴⁰ In our simulations, $J i$ is initially fixed for all oscillators and set to 0, whereas during aging, the parameter is set to the excitable regime, i.e., $J i<− ε i$⁠, as explained further in continuation. To mimic heterogeneity of oscillators, cell-to-cell variability is introduced by randomly distributing parameter $a i$ for each oscillator from values 3 to 5, whereby higher values of $a i$ correspond to longer oscillation periods and a slight increase in the peak amplitude. The range of $a i$ was selected so that the oscillation frequencies range in approximately ±15% from the mean. Parameter $b i$ is constant and set to 1. The system of equations was numerically solved with the Runge–Kutta method of second order and a time step of $dt=0.01$⁠. The initial conditions were chosen randomly from the interval [1.25, 1.75].

The oscillators are arranged in two subsets, inactive and active, with their sizes being $p∗N$ and $( 1 − p)∗N$⁠, respectively, where p is the ratio of inactive oscillators. By increasing the ratio of inactive oscillators by setting the external stimulus J outside the oscillatory regime, the system then undergoes an aging transition. An individual oscillator is considered inactive if its parameters are changed so that its solution, if the oscillator was uncoupled, becomes a fixed point. We define $J i=−σ ϵ i$ if the $i$th oscillator is in the inactive set and $J i=0$ if it is in the active set. A parameter $σ$ is introduced to specify how far from the bifurcation point (oscillatory regime) the oscillator is placed. Higher values of σ indicate a larger distance from the bifurcation and thereby a stronger degree of inactivation. To quantify the global oscillatory activity of the network M, we compute the root mean square amplitude on the variable x. M is then defined as $M= ⟨ ( X − ⟨ X ⟩ ) 2 ⟩,$ where $X= ∑ i = 1 N x i / N$ is the sum over all oscillators at time t and the angle brackets signify an average over time T. To minimize the effect of any transient phenomena on amplitude M, the traces for first t = 200 were discarded. The results presented in the paper were obtained by averaging M over up to 80 different realizations of network configurations and initial conditions for each set of parameters.

In our study, we systematically investigated how the global amplitude M changes as a function of the ratio of inactive oscillators p, coupling strength K, and aging strength $σ$⁠. Additionally, in our simulations, we considered and compared three aging scenarios. Random inactivation of nodes and two scenarios where the oscillators were made inactive progressively in either decreasing (hub nodes first) or increasing (peripheral nodes first) order of their degree ranks. We start by showing in Fig. 2 how the global oscillatory activity deteriorates as a function of inactivated oscillators for different sets of parameters and strategies. For weak coupling (K = 0.02), M decreases rather continuously, whereby the fastest deactivation is attained when high-degree nodes are targeted first and the slowest when the low-degree nodes are preferentially placed before the bifurcation. These differences become more pronounced if the distance to the oscillatory regime is larger (σ = 8). Interestingly, as the coupling increases (K = 0.056), the network remains intact for a higher fraction of inactivated oscillators if hubs are targeted first, but only if the oscillators are not placed far from the bifurcation. Remarkably, when the coupling is further increased (K = 0.19), the sequence changes completely so that the excitable network loses its dynamic integrity the fastest when the peripheral nodes are targeted first, but again only for smaller displacements from the oscillatory regime (σ = 2.5). Finally, if the coupling is very strong (K = 0.64), targeting the high-degree nodes is the most effective way to abolish global network activity, irrespective of the distance from the bifurcation point. Apparently, with increased coupling, the transitions do not only become steeper but can also display non-trivial behavior, particularly for intermediate values of coupling, which will be further and more systematically investigated in continuation.

In Fig. 3, we present the color-coded values of fractions p at which the global network amplitude M falls below 20% (upper row) and 80% (lower row) of its initial value (p = 0) as a function of the coupling strength (K) and the distance from the bifurcation point (σ). It can be noticed that the network's dynamical activity strongly depends on the strategy with which individual components are deteriorated and that particularly the coupling strength affects the aging transition in an interesting manner. Specifically, for intermediate coupling strengths, macroscopic oscillations in the excitable network vanish at higher fractions of inactivated elements than in the cases of a weak or strong coupling regime, which is most pronounced when the hub nodes are inactivated first. In these circumstances, the network seems to be dynamically more robust to targeted inactivation also when compared to random inactivation, or when peripheral nodes are put below the bifurcation point first. This has been noticed already in Fig. 2 and will be further explored in continuation. Moreover, the parameter σ reflecting the distance from bifurcation has in principle a more monotonous role, but the analysis reveals that the most interesting behavior can be noticed only if the oscillators are not put too deep into the excitable regime.

To assess how the strategy of inactivations with regard to the coupling strength and the distance from the bifurcation affects the excitable network dynamics in further detail, we show in Fig. 4 the color-coded differences between the critical p values presented in Fig. 3. Specifically, we subtracted the values of fractions of inactivated oscillators at which the amplitude of macroscopic oscillations falls below 80% or 20% for the three different combinations: (i) random inactivation minus targeted inactivation of high-degree nodes (left column), (ii) random inactivation minus targeted inactivation of low-degree nodes (middle column), and (iii) the difference between targeted inactivation of low-degree and high-degree nodes (right column). Negative values in these plots indicate that the subtracted strategy is more favorable in terms of the dynamical robustness of the excitable network. Please note that the color-scales differ in each panel. The interesting patterns in these panels indicate that the examined parameters indeed affect the aging transition in a non-trivial way, especially when the oscillators are not kicked too far from the oscillatory regime. This in combination with intermediate values of coupling strengths results in counterintuitive behavior and specifies thereby the parameter regions where the differences between the structural and dynamical robustness are the most pronounced.

Finally, to gain more mechanistic insight, we assess how for specified sets of parameters, the amplitudes of individual oscillators $M i$ in the network change when they are targeted at random or when the high-degree or low-degree nodes are inactivated first. Results for two intermediate values of coupling strength (K = 0.056 and K = 0.19) and for selected sets of parameters that correspond to the second and third panels in Fig. 2 are shown in Fig. 5. When peripheral nodes are put in the excitable regime first (right column, upper row), their amplitude immediately drops close to zero. In contrast, when hub nodes are targeted first, their amplitude remains nearly normal, as they are coupled to many other oscillators that are normally active. Moreover, when oscillators are disabled by random, most often peripheral nodes are affected, which, in turn, leads to the deterioration of their activity. For that reason, this strategy is less advantageous than targeted inactivation of high-degree oscillators and only slightly more favorable than targeted inactivation of low-degree oscillators (see Fig. 2, the second panel in the upper row). Conceptually, very similar behavior is noticed at the higher coupling value (lower row, K = 0.19), except that due to stronger interaction the amplitude of individual oscillators drops close to zero at higher fractions of inactivated units. The results presented in Fig. 5 explain the counterintuitive notion that the network's dynamical robustness can be the greatest when preferentially hub nodes are inactivated.

The loss of the macroscopic oscillatory activity in different classes of networks due to the failure of its microscopic constituents is a vibrant topic in the scientific community and has many practical applications in physical, biological, and engineering systems. It is, therefore, of paramount importance to investigate how different circumstances affect the macroscopic activity when some of the microscopic dynamical units lose their self-oscillatory behavior due to breakdowns. This issue is particularly relevant in biological systems research, where the collective function often arises from the interactions of individual oscillatory units that can progressively malfunction under pathological conditions. In this paper, we aimed to further investigate aging transitions in a realistic biological network constituted by heterogeneous excitable elements. This setting can be used to describe a plethora of biological systems, such as neuronal ensembles, populations of cardiac and other smooth muscle cells, multicellular syncytium formed by pancreatic endocrine cells, social systems, etc.^45–50 In our study, we numerically investigated how different physiological determinants affect the transition to a macroscopically inactive state. Specifically, the interplay between the distance from bifurcation, coupling strength, and the strategy of targeting the oscillators in the network was found to profoundly affect the aging transition, whereby under specific circumstances even a counterintuitive and non-trivial behavior could be noted. Perhaps, the most striking observation is that for some intermediate coupling strengths, the network remains its global activity for the highest fraction of inactivated units when the high-degree nodes are targeted first. This goes in hand with previous observations that oscillatory networks can be highly fragile to targeted inactivation of low-degree nodes, but only under specific circumstances, such as weak coupling.^21,22 We studied this phenomenon in further detail and revealed that the most efficient strategy to inactivate the network does not only non-monotonically depend on the coupling strength but is also crucially affected by how far from the bifurcation point individual oscillators are poised. Notably, we also found that within reasonable limits, the observed behavior does not depend on the network size.

The distance from a bifurcation point and the coupling strength have both important implications for real excitable systems. The distance from the oscillatory regime affects the switching behavior of the system and the sensitivity to perturbations and can, therefore, play a critical role in determining the stability of neuronal, cardiac, or hormonal rhythms.⁵¹ Experimentally, the degradation of the collective activity in ensembles of excitable cells is often studied by means of optogenetic silencing.^52,53 By these means, the degree of hyperpolarization, which resembles how far from the active state a cell is poised, can be tuned, representing thereby a relevant factor to consider when studying such perturbations of the global activity in tissues. Moreover, the coupling strength between excitable cells can change under various physiological and pathophysiological conditions. Changes in the extracellular conditions can, for example, regulate the expression or activity of gap junctions, which can, in turn, modulate intercellular interactions on different timescales.^54,55 Changes in coupling between excitable cells also occur during pathogenesis, such as in neurodegenerative diseases, epileptogenic activity, and diabetes^56–58 As a result, the nature of the aging transition in excitable tissues during aging, or disease progression, when individual cells lose their self-oscillatory behavior due to malfunction can also be influenced by the associated changes in cell-to-cell interactions. The strategy how individual cells are targeted also has a strong biological meaning. Random targeting refers to malfunctions without any specific criteria, such as random mutations or damage of cells. On the other hand, preferential targeting involves selecting a subset of cells with specific characteristics. Hub cells often have higher metabolic demands and can be more susceptible to the accumulation of oxidative damage. In addition to their high degree of connectivity and their involvement in multiple cellular pathways, this can make them more vulnerable to stressors, insults, or apoptosis. In this regard, theoretical as well as experimental studies on networks of excitable beta cells have shown that targeting certain subpopulations of cells can much more affect the function of the multicellular syncytium than others.^52,59,60 Furthermore, in certain physiological processes, such as cell competition or neuronal development, the elimination of peripheral elements is favorized before the central ones. Apparently, targeting of different types of cells in the multicellular networks may have distinct consequences for tissue function during the onset of age-related or other diseases.^52,61,62 Finally, heterogeneity is another genuine characteristic of biological networks. In real excitable networks, individual components are not identical but heterogeneous in their activities and responses.^63,64 This aspect has also been included in our computational model by introducing randomness to one of the parameters, which has led to variability in intrinsic oscillation frequencies.

In sum, in our study, we have shown that aging transition in a network of coupled excitable oscillators is a dynamical phenomenon that is influenced by many factors. Our systematic analyses indicate that the degradation of macroscopic activity is governed in a non-trivial manner by the interplay between the coupling strength, the distance from the bifurcation point, and the type of strategy used to place the oscillators from the oscillatory to the excitable regime. Our findings thus provide novel insights into the understanding of the dynamical robustness in specific biological networks and, in addition, point out possible future experimental studies in excitable tissues in which these concepts could be further investigated by optogenetic silencing or photoablation.

The authors acknowledge the support from the Slovenian Research Agency (Research Core Funding Nos. P3-0396, P1-0403, and I0-0029, as well as Research Projects Nos. J3-3077 and J1-2457).

The authors have no conflicts to disclose.

Uroš Barać: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Software (lead); Visualization (equal); Writing – review & editing (equal). Matjaž Perc: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Writing – review & editing (equal). Marko Gosak: Conceptualization (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal).

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