# Fluctuation of Dynamical Robustness in a Networked Oscillators System

## Abstract

In this work, we study the dynamical robustness in a system consisting of both active and inactive oscillators. We analytically show that the dynamical robustness of such system is determined by the cross link density between active and inactive subpopulations, which depends on the specific process of inactivation. It is the multi-valued dependence of the cross link density on the control parameter, i.e., the ratio of inactive oscillators in the system, that leads to the fluctuation of the critical points. We further investigate how different network topologies and inactivation strategies affect the fluctuation. Our results explain why the fluctuation is more obvious in heterogeneous networks than in homogeneous ones, and why the low-degree nodes are crucial in terms of dynamical robustness. The analytical results are supported by numerical verifications.

###### pacs:

89.75.Hc,05.45.Xt,89.20.-aRobustness is an important aspect of the nature of a system. In the past years, topological robustness has been deep studied, but the dynamical robustness did not catch much attention. In this paper, we study the dynamical robustness of a system comprised of Stuart-Landau oscillators. The oscillators are divided into two kinds: active and inactive. When the ratio of inactive oscillators increase, the status of the system will from oscillation change to amplitude death and the critical point can be used to measure the Robustness of the sysrem. Using numerical simulation and theoretical analysis, we find that the cross link density between active and inactive subpopulations determins the dynamical robustness of the system.

## I Introduction

Dynamical systems in nature and human labs usually comprise a large number of interacting individual elements, such as synchronizing fireflies (1), neurons in human brain (2) , cardiac pacemaker cells (3), power grids (4); (5), and Josephson junction arrays (6), just to name a few, which can be modeled by networked oscillators (7). One important issue in the study of such dynamical systems is the robustness, i.e., the ability to maintain basic structure and function under attacks or disfunctions (8); (9); (10); (11); (12). It has been shown that the huge blackouts, which inevitably cause tremendous economic loss, are related to the cascading failure of power-grids (8). Such structure robustness involving network connectedness has been intensively investigated in previous works (8); (9); (10). On the other hand, networked systems typically carry dynamics, for example, the circadian rhythms of mammals, the synchronization of cardiac cells, and spatio-temporal patterns in brain (2); (3); (2). Of equal importance is the dynamical robustness, i.e., the ability of a networked system to maintain its normal dynamical activity when the topology or the local dynamics of network components are subjected to changes.

In Ref. (13), Daido et al. investigated the dynamical robustness of a globally coupled system which simultaneously consists of active and inactive oscillators. With the increase of the ratio of inactive oscillators, the dynamical activity in the system decreases until finally it totally vanishes at certain critical point. Such phenomenon is termed as aging transition which characterizes the dynamical robustness of such networked systems (13); (11). Recently, Tanaka et al. further studied this phenomenon in various complex network topologies (12). Surprisingly, it is found that the low-degree nodes, rather than the hubs, play a dominant role in terms of dynamical robustness. This finding uncovered the great differences between dynamical robustness and structure robustness.

For such dynamical systems, we always expect the dynamical robustness of them could be controllable. While sometimes we wish a strong robust system to maintain its functionality when suffering attacks, sometimes we prefer a fragile system because it is harmful and we need to destroy it. Tanaka’s work reported that coupled oscillators network with smaller mean degree appears more robust, which give us a hint that we could increase the dynamical robustness of the system by reducing connections between elements. However, in this paper, we report the fluctuation of dynamical robustness in the above mentioned networked system. We find that even the network topology and the coupling strength in the dynamical system are fixed, there is still a fluctuation for the critical point of aging transition measured by the ratio of inactive oscillators in the system, which characterizes the dynamical robustness of the system. Moreover, the fluctuation could be large enough that totally changes the robustness property of the system. So it is important to understand the fluctuation for the potential applications of dynamical robustness. In addition, this fluctuation occurs more obviously in heterogeneous networks than in homogeneous ones. By applying mean-field approximation and linear stability analysis, we show that the cross link density between subpopulations of active and inactive oscillators plays a dominant role. In principle, different inactivation realizations lead to different link densities, which finally result in fluctuations of the critical points. Based on this analysis, we can explain why the fluctuations in heterogeneous networks are severer than that in homogeneous networks. Our results are also helpful to understand the crucial role of the low-degree nodes in determining the dynamical robustness of networked oscillators reported in Ref. (12). The rest of this paper is organized as follows: section II introduces the dynamical model; theoretical analyses and numerical verifications are presented in section III; finally, conclusions are presented in section IV.

## Ii Model

In this work, we investigate a dynamical model of networked oscillators. It has following features: (1) the dynamical system is described by coupled oscillators; (2) there are two types of oscillators in the system, i.e., active and inactive, representing oscillating state and non-oscillating state, respectively; (3) the coupling among oscillators forms a complex network. The general form of this model can be written as:

(1) |

Here, is the index of oscillator or node. is the state vector describing the dynamics of oscillators. The first term at the RHS of Eq. (1) describes the local dynamics of an oscillator, and the second term is the diffusive coupling that denotes interactions among different oscillators. is the coupling strength. is the element in the adjacent matrix of coupling network, which equals to 1 if nodes and are connected, and 0 otherwise. It should be pointed out that similar models have been studied previously. For example, aging transition in such models was investigated in fully coupled network in Ref. (13), and later in regular ring topology in Refs (14); (17). Recently, it has been extended to different complex topologies (11); (12).

In the present work, we mainly choose the Stuart-Landau (SL) oscillators as the local dynamics. Specifically, the networked SL oscillators system can be described by the following coupled ODEs:

(2) |

where and are the complex amplitude and the inherent frequency of the th SL oscillator, respectively. is the control parameter denoting the distance from a Hopf bifurcation point. When , the oscillator is a limit cycle with an amplitude . However, it settles down to a fixed point when . Thus generally the oscillator is active when and inactive when . The oscillator will lose its activity as its value converts from positive to negative. This can be used to model the two distinct dynamical states of oscillators. We define parameter as the ratio of inactive oscillators in the network. Reasonably, the global activity of the networked system can be characterized by the normalized order parameter , defined as with . As increases to a critical value , the networked system will totally loses its global activity until , i.e., an aging transition occurs (13). Because the ratio is the largest ratio with which the dynamical system can maintain certain activity, it measures the dynamical robustness of this networked system. The larger the , the better the dynamical robustness.

For simplicity, we will not consider the case in which the frequencies of oscillators have a distribution. Throughout this paper, we set . Since the model simultaneously contains both active and inactive local dynamics characterized by positive and negative , respectively, we first examine how the configuration of values affect the global dynamics. To this end, we consider three typical distributions of . The first and the simplest case is for all active oscillators, while for all inactive ones, i.e., is binary. In the second and the third cases, satisfies uniform distribution and Gaussian distribution, respectively, i.e., varies in a range. The ratio of the inactive oscillators can be tuned by gradually changing the mean of the distribution from positive to negative. Fig. 1 illustrates the aging transition for the three types of settings. Interestingly, it is found that although the values of satisfy different distributions, what really determines the transition point is the the ratio of the inactive oscillators in the system. Therefore, we can always equivalently simplify the distribution of into binary case by coarse-graining process. Without loss of generality, in the following we only consider the situation where takes binary values.

Besides SL oscillators, we also consider the following networked Rössler oscillators in this paper:

(3) |

where are the state variables of Rössler oscillator, and are parameters. In this system, there are also two types of oscillators by choosing different parameters, i.e, for active oscillators and for inactive oscillators. To quantify the global activity, the order parameter can be defined as , where is the centroid and the bracket means the long time average. Numerical results thoughout this work were obtained for random initial conditions by means of the fourth order Runge-Kutta method with time step 0.01.

## Iii Analysis and results

Recently, in our study we find an interesting phenomenon, i.e., even for fixed network topology and coupling strength, the critical point has fluctuations depending on specific network realizations. Typical examples are shown in Fig. 2. In our study, in order to conveniently compare different network topologies, we propose a simple but effective method to continuously change a heterogeneous network into a homogeneous one. The detail of the algorithm is explained in Appendix. As shown in Fig. 2, fluctuations of critical point are observed in both heterogeneous and homogeneous networks. To quantify this feature, in Fig. 3, we plot the variance of when the network topology continuously changes from heterogeneous to homogeneous with the increase of rewinding probability . It is found that this phenomenon is more obvious in heterogeneous networks. In the following, we present both theoretical analyses and numerical verifications to understand this phenomenon. We notice that some non-mean-field methods, e.g., in Ref. (18), have been proposed to deal with synchronization of Kuramoto-type phase oscillators on complex networks. However, a general theoretical framework to treat networked periodical and chaotic oscillators involving amplitudes is still not available. Basically, the analyses in the present study are based on mean-field approximation.

### iii.1 Linear stability analysis

We start from Eq. (2). When , i.e., the system contains only identical active oscillators. In this case, it will easily achieve global synchronization as the coupling strength increases. When we randomly choose some nodes and change the oscillators into inactive states, i.e., the proportion of inactive nodes , the system simultaneously contains both active and inactive oscillators. Such a process is called random inactivation. Throughout this paper, it is adopted by default unless otherwise stated. Usually, the non-trivial aging transition would only occur when the coupling strength is large enough (13), so before the transition all the active oscillators already well synchronize into a cluster, and so does the inactive oscillators. This has been numerically verified and one such example is illustrated in Fig. 4. Although the synchronization is not rigorous, it is good enough for us to approximately reduce the high-dimensional system into two identical clusters of oscillators, denoted as and for active and inactive oscillators, respectively. We use and to represent the dynamics of active and inactive subpopulations, respectively, as treated in Ref. (13). Then the dynamics of the networked system, i.e., Eq. (2), can be written as:

(4) | |||

(5) |

Here, means the degree of oscillator , and is the proportion of inactive oscillator among the neighbors of oscillator . Summing the equations of all active and inactive oscillators, we obtain

(6) | |||

(7) |

In the above equations, the summing terms at the RHS, i.e., and , both represent the total number of cross links, denoted by , between the subpopulations of active oscillators and inactive ones. Actually they are the same thing in different ways. To avoid the influence of network size, we define the density of cross links as . Then Eqs. (6) and (7) become

(8) | |||

(9) |

From these equations, we immediately find that the network topology actually affects the global dynamics through parameter . It is the cross link density between active and inactive subpopulations that plays a dominant role in this networked system. Now we analytically study how the dynamical robustness of such system can be affected by parameter . With the increase of control parameter , the dynamics of the networked system will gradually lose global activity. This process can be characterized by when . At the transition point, the networked system loses its stability and in the mean time the trivial fixed point becomes stable. By a linear stability analysis, we obtain the critical point as:

(10) |

with . For typical dense complex network, . Therefore, is a small quantity. We apply Taylor expansion to in Eq. (10) and keep the linear term. It finally becomes

(11) |

### iii.2 Analysis to fluctuation of critical point

Eq. (11) shows that only has functional relation with given that the parameters of the local dynamics and the coupling strength are fixed. Apparently, there is a maximal value when . With the increase of , will monotonically decrease. When , approaches the minima , as shown in Fig. 5. This result provides us insight that the dynamical robustness of such system is determined by the density of cross links . Physically, this can be reasonably understood. Since there are both active and inactive local states in the network, the interaction or influence between the two subpopulations are crucial for the global activity. The density of cross links just characterizes this interaction. For example, when it is large, there exists strong interaction between the two subpopulations of active and inactive oscillators. As a result, a small critical value can be expected.

Now we explain why there exists fluctuation of critical point even when the network topology is fixed. Let us analyze the phase diagram on the parameter panel of -. As shown in Fig. 5 (a), Eq. (10) actually defines the curve of bifurcation, i.e., the boundary of two areas corresponding to distinct dynamical phases of the system. In the upper right area, the networked system is in the quenching state losing global activity; while in the bottom left area, it is in active state oscillating to some extent. The active state loses its stability when the system passes through the curve Eq. (10). For an inactivation process, i.e,. in a specific realization, is a function of . Actually, is a unimodal function satisfying , as shown in Fig. 5 (a). The curve intersects the bifurcation curve Eq. (10), and the crosspoint just determines the critical point . The key point here is that is a multi-valued function of even when the network topology is fixed. For different numerical realizations, i.e., in different specific inactivation processes, we will get different curves , which lead to different crosspoints in the bifurcation curve Eq. (10). This is the origin of fluctuation of critical point observed in our study.

Naturally, one may want to know which inactivation process will give the maximal or the minimal , which directly limits the fluctuation scope of given the network topology. Interestingly, this problem can be mapped into the ground state problem of anti-ferromagnetic Ising model or the MAX-CUT problem in combinatorial optimization (19); (20); (21). It has been proved that the solution of these problems in general networks is a NP-Complete problem, i.e., we cannot find an effective algorithm to determinate it in polynomial time (19). Since analytical method is not available to get the fluctuation scope of , we turn to numerical study. To this end, we numerically plot the possible curves on the parameter panel of -. The two crosspoints that define the largest range on the axis of roughly gives the fluctuation range of . As shown in Fig. 5(b), it is found that for a heterogeneous network has a much larger fluctuation area compared to a homogeneous network. Therefore, more obvious fluctuation of is observed in the former case. Physically, this result can be heuristically understood. We know that the fluctuation of is caused by the multi-valued , which depends on specific inactivation processes. Obviously, an inactivation process is more sensitive to the topology in heterogeneous networks than in homogeneous ones because there are huge differences among the node degrees in the former case. This leads to the fluctuation of to a greater degree, which finally contributes to the fluctuation of the critical point .

As one application of the above theory, we now consider the case of targeted inactivation rather than random activation in network. In Ref. (12), It has been reported that under targeted inactivation the dynamical robustness of the system depends more on the low-degree nodes rather than the hubs. This important finding reveals the crucial role of the low-degree nodes in the context of dynamical robustness. Based on our analytical treatment, we can provide an explaination here. In our study, we apply three typical strategies of inactivation: (1) Inactivation goes from the node with the maximal degree to the one with the minimal degree; (2) Inactivation takes the inverse order of (1); (3) Random inactivation. Numerical results for both networked SL oscillators and networked Rössler oscillators are shown in Fig. 6. In both cases, under strategy 1 is always greater than that under strategy 2, regardless of the heterogeneity/homogeneity of the network topology. To understand the above results, we plot the curves on the - parameter plane. As shown in Fig. 7, all curves are unimodal with ==0. Considering the specific processes of the three strategies, it is not difficult to figure out that with the increase of , curve 1 increases faster than curve 2 to the maximum, while the decrease from maximum to zero is just the opposite. Strategies 1 and 2 are two extreme cases, and strategy 3 is between them. As shown in Fig. 7, the crosspoints of curves with the bifurcation curve defined by Eq. (10) determine that . Therefore, our theory successfully explain why the low-degree nodes plays more important role than the hubs in terms of dynamical robustness.

### iii.3 Random inactivation

In the above analysis, the only limitation for the network topology is that it should be dense enough so that the cross link density , and there is no requirement for the strategy of inactivation. Thus the above result holds for general inactivation processes. In fact, in our study, we find for typical random inactivation, Eqs. (8) and (9) can take another forms, and the theoretical treatment can be simplified.

We notice that in principle for different oscillator , and are not necessarily the same. However, they must satisfy the constraint and over all the nodes when the mean degree of the network is fixed. Let , where is the deviation of from its mean value. Normally, a node with larger degree would have smaller deviation. In a dense network, it is reasonable to expect that and distributes symmetrically around its mean value. For random inactivation, the active and inactive oscillators are mixed evenly. Thus the distributions of in active and inactive subpopulations will approximately the same as in the whole system, i.e., we can approximately assume for all . Then we have , and . Substitute these relations into Eqs. (8) and (9), we get

(12) | |||

(13) |

where and are the mean degrees of active and inactive subpopulations, respectively. Similarly, by applying a linear stability analysis, we can analytically obtain the critical point as:

(14) |

From this equation, we can immediately find that the dynamical robustness of the system are determined by the mean degrees of active and inactive subpopulations in the case of random inactivation. The point is, even though the mean degree of the whole network is fixed, there is still some degree of freedom for and as long as they satisfy the following constraint:

(15) |

On the parameter panel of -, only certain area can satisfy this condition, as shown in Fig. 8. In particular, parameters and not only are related to the network topology, but also to the specific strategy of inactivation. Usually, each different realization results in different and , causing the small fluctuations of observed. We can see in Fig. 8 that our theoretical result is qualitatively consistent with the numerical verifications in both networked SL system and Rössler system. Extensive numerical results have shown that Eq. (14) is valid as long as the mean degree is large enough, e.g., .

Furthermore, a trivial solution always exists for Eq. (15), i.e., . In this circumstance, Eq. (15) degenerates as:

(16) |

In strictly sense, this result only holds for absolutely homogeneous topologies, i.e., the regular networks. However, for random inactivation in very homogeneous network, the above conclusion can be expected to hold approximately. In this case the critical point only involves the mean degree and is almost independent with the strategy of inactivation. Interestingly, this coincides with one situation investigated in Ref. (12), where Eq. (16) is derived by a different approach. Moreover, Eqs. (14) and (16) can recover Eq. (4) in Ref. (13) when the topology is globally connected; and under the strong coupling limit, Eq. (16) can reproduce the results in Refs. (16); (15). Therefore, all these studies provides insights from different perspectives into the dynamical robustness of such networked system.

## Iv Conclusion

In this work, we have studied the dynamical robustness in a networked system, in which both active and inactive oscillators coexist. With the increase of inactive oscillators in such system, it will gradually lose its global activity. The critical point of the transition characterizes the dynamical robustness of such networked system. Interestingly, we found that the critical point fluctuates even thought the coupling strength and the network topology are fixed. By analytical treatment and numerical experiments, we successfully explain the origin of this phenomenon. Mainly, our study provides the following results: (1) The fluctuation of dynamical robustness in such networked system is caused by the variation of the cross link density. It is the multi-valued dependence of the cross link density on the control parameter, i.e., the ratio of inactive oscillators in the system, that leads to the fluctuation of critical points; (2) Since the inactivation process affects the cross link density more in heterogeneous networks than in homogeneous ones, the fluctuation turns out to be more obvious in the former case; (3) Fluctuations under three typical strategies of inactivation can be predicted based on our theory and verified by numerical simulations. This helps understand the importance of the low-degree nodes in terms of dynamical robustness recently reported in (12); (4) For random inactivation, the critical point is determined by the mean degrees of subpopulations of both active and inactive oscillators; and in very homogeneous networks, it is only related to the mean degree of the network. This paper is helpful to understand the collective behaviors of networked dynamical systems, and shed light on certain practical applications as well.

SGG is sponsored by: The Science and Technology Commission of Shanghai Municipality grant No. 10PJ1403300; The Innovation Program of Shanghai Municipal Education Commission grant No. 12ZZ043; and The NSFC grant Nos. 11075056 and 11135001. ZHL is sponsored by The NSFC grant Nos. 10975053 and 11135001.

## Appendix A Method to tune the network topology

Previously, there are several methods to change network topology from homogeneous to heterogeneous (22); (23). For example, the degree distribution of network can be tuned between exponential and power-law. However, in the present work, we expect degree distribution to be Poisson form or even Delta function. So we propose a simple but effective method for this purpose. The main idea is to gradually rewind edges from an initially heterogeneous network. Here are the main steps: (1) Generate a heterogeneous network, e.g., the BA network characterized by a power law degree distribution; (2) Choose an arbitrary edge and compare the degrees of both its ends; (3) Cut the edge from the node with higher degree, and randomly rewind it to a node in the network. We define the rewinding probability as the number of edges rewound normalized by the total number of edges in the network. When varies from 0 to 1, the network topology will change from heterogeneous to homogeneous as verified by numerical experiments. Therefore, is a parameter that can characterize the heterogeneous extent of a network. In our simulations, we have started from scale-free networks with power law exponents between 2 and 3, the results are qualitatively the same.

### References

- Buck J 1988 Q. Rev. Bio. 63 265
- Izhikevich E 2007 Dynamical systems in neuroscience: the geometry of excitability and bursting (Massachusetts: MIT Press)
- Glass L 2001 Nature 410 277
- Filatrella G, Nielsen A and Pedersen N 2008 Eur. Phys. J. B 61 485
- Rohden M, Sorge A, Timme M and Witthaut D 2012 Phys. Rev. Lett. 109 064101
- Wiesenfeld K, Colet P and Strogatz S H 1996 Phys. Rev. Lett. 76 404
- Pikvosky A, Rosenblum M and Kurths J 2003 Synchronization: A universal concept in nonlinear sciences (Cambridge: Cambridge University Press)
- Buldyrev S, Parshani R, Paul G, Stanley H and Havlin S 2010 Nature 464 1025
- Callaway D S, Newman M E J, Strogatz S H and Watts D J 2000 Phys. Rev. Lett. 85 5468
- Vespignani A 2010 Nature 464 984
- Morino K, Tanaka G and Aihara K 2011 Phys. Rev. E 83 056208
- Tanaka G, Morino K and Aihara K 2012 Scientific Reports 2 232
- Daido H and Nakanishi K 2004 Phys. Rev. Lett. 93 104101
- Daido H 2008 Europhys. Lett. 84 10002
- Daido H 2011 Phys. Rev. E 83 026209
- Pazó D and Montbrió E 2006 Phys. Rev. E 73 055202
- Daido H 2011 Phys. Rev. E 84 016215
- Restrepo J G, Ott E and Hunt B R 2005 Phys. Rev. E 71 036151
- Garey R and Johnson D 1979 Computers and Intractability: A Guide to the Theory of NP-Completeness (New York: WH Freeman and Company)
- Hyafil L and Rivest R 1976 Info. Proc. Lett. 5 15
- Zhou H 2005 Phys. Rev. Lett. 94 217203
- Liu Z, Lai Y C, Ye N and Dasgupta P 2002 Phys. Lett. A 303 337
- Gómez-Gardeñes J and Moreno Y 2006 Phys. Rev. E 73 056124