Stability threshold approach for complex dynamical systems

A new measure to characterize stability of complex dynamical systems against large perturbation is suggested, the stability threshold (ST). It quantifies the magnitude of the weakest perturbation capable to disrupt the system and switch it to an undesired dynamical regime. In the phase space, the stability threshold corresponds to the"thinnest site"of the attraction basin and therefore indicates the most"dangerous"direction of perturbations. We introduce a computational algorithm for quantification of the stability threshold and demonstrate that the suggested approach is effective and provides important insights. The generality of the obtained results defines their vast potential for application in such fields as engineering, neuroscience, power grids, Earth science and many others where robustness of complex systems is studied.


Introduction
Complex systems science is strongly based on linear stability analysis considering small perturbations of dynamical systems. In a seminal paper [1] this concept was extended even to the stability of synchronization in complex networks leading to the efficient master stability formalism. However, for various applications often the influence of large perturbations is also of crucial importance. Typical examples are climatological systems, in particular ocean circulations. It is well accepted that the Atlantic Meridional Overturning Circulation may be sensitive to changes in the freshwater balance of the northern North Atlantic. When an anomalous freshwater flux is applied in the subpolar North Atlantic, this circulation collapses in many ocean-climate models [2]. Another example is power grids which are networks of connected generators and consumers of electrical power. For proper function of such networks synchronization between all the nodes is essential. Local failures, overloads or lines breaks may cause desynchronization of nodes and lead to large-scale blackouts [3,4].
The study of a system's stability against large perturbations implies treating the following challenging problem: definition of the class of 'safe', or admissible perturbations after which the system returns back to the initial regime. In contrast, 'unsafe' perturbations switch the system to another, often unwanted, dynamical regime. The definition of the class of safe perturbations of a nonlinear system is very complicated and basically different from the linear stability analysis. The reason is that for large perturbations linearization is inadequate and the perturbed dynamics is governed by nonlinear equations whose analytical study is impossible in general. Some analytical methods do exist, for example the method of Lyapunov functions [5]. However, this method has serious limitations since a Lyapunov function for a particular dynamical system is often not constructive. The 'safe' perturbation class can also be analytically estimated for some specific systems, e.g. networks of spiking neurons [6]. Nevertheless, an important task is to develop numerical methods of defining and describing the class of safe perturbations [7].
From the viewpoint of nonlinear dynamics, established dynamical regimes of the system correspond to attractors in the phase space. The class of safe perturbations is equal to the attraction basin, i.e. the set of the points which converge to the attractor. A perturbation is safe if it brings the system to a point within the attraction basin. The first attempt to characterize attraction basins in complex networks was undertaken in [8] where the concept of basin stability was introduced. The basin stability equals is the density of the perturbed states, and x c ( ) equals one if the point x converges to the attractor and zero otherwise. The value S 0; 1 B Î ( ]expresses the likelihood that the perturbed system returns to the attractor. An important advantage of this measure is that it can be easily calculated by the Monte-Carlo method. Namely, one should just take a large number of random points in the phase space and check how many of them converge to the attractor.
. Basin stability is an important characteristic extending the concept of linear stability for the case of large perturbations. However, many real dynamical systems, especially complex networks, possess highly dimensional phase space with complicated structure. This makes it problematic to characterize an attraction basin by just a single scalar value. Moreover, basin stability depends on the perturbation class Q which should be chosen a priori.
In this paper we suggest a new measure to characterize stability against large perturbation, the stability threshold (ST). We were inspired by the observation that for real systems it is often important to know the maximal magnitude of perturbation which the system is guaranteed to withstand, like the maximal voltage jump for a stabilizer or the maximal bullet energy for a bulletproof vest. In the following, in section 1 we introduce the ST in detail and explain how to calculate it. In sections 2 and 3 its potential is demonstrated for two paradigmatic model systems. In section 4 we relate the two measures, the ST and the basin stability. In section 5 we briefly summarize and discuss our results.

Definition and quantification of stability threshold
We define ST as the minimal magnitude of a perturbation capable of disrupting the established dynamical regime, i.e. to push the system out of the attraction basin. In the phase space, ST is the minimal distance between the attractor  and the border  d of its attraction basin, i.e.
where dist , (· ·) is the distance between the points in the phase space. Further, we use the Euclidean metric to calculate the distance. However, any other metrics can be used as well.
To better understand the physical meaning of ST consider the system settled to the attractor  as depicted in . Consider now a perturbation of the system which results in a quick change of its state. Namely, let the system state change from the unperturbed one x u to the perturbed one x p so that the perturbation x x x p u D = has the magnitude q x = D | |. If q s < , the perturbation can never kick the system out of the attraction basin ( x 1 D in figure 1). But if q s > and the system is near the point a just before the perturbation, it may be kicked out of the basin if the direction of the vector x D is close to the vector D b a = -( x 2 D in figure 1). The above reasoning shows that besides the value of σ, the direction of the corresponding vector D is critical. This vector corresponds to the most 'dangerous' direction of perturbations in which the distance to the basin border is the shortest. Let us now come to the question of quantification of the ST. For this purpose we suggest a two-stage algorithm whose basic principles are described below and also illustrated in figure 1.
(i) At the first stage we identify some point K 1 on the border of the attraction basin. For this purpose we choose an arbitrary point K 0 in the vicinity of the attractor and start to move from the attractor until leaving the basin. The point K 1 is found then by the bisection method.
(ii) At the second stage we move along the basin border. On each step we draw a tangential hyperplane to the border at the current point K n . In the hyperplane we find the point closest to the attractor  and make a step towards this point and so obtain the new point K n 1 + . Such steps bring us closer and closer to the attractor and finally converge to the point M on the border with the minimal distance to the attractor.
The second stage relies on smoothness of the border for only in this case it can be approximated by a tangential hyperplane. Further we assume the border is smooth (defined by a function of class C 1 ) which is the case for many realistic systems.
The suggested algorithm allows us to determine the local minima of the distance between the attractor and the basin border, which we call further 'local threshold' (LOCT) points. Starting from different initial points we get different LOCT points M M M , , , m ) .
This brute-force search to obtain all the local minima does not seem to be a very effective strategy. However, the effectiveness of the method is essentially improved in a parametric study, i.e. when the system properties are studied versus its parameters. Note that such tasks are typical since all realistic systems depend on parameters and usually one wants to know what happens if they are varied. Suppose that for a certain parameter value . In a robust system, the phase space structure changes continuously when p is changed. Thus, the coordinates of LOCT points depend continuously on p. So, when p is changed by a small value p p p 0 D =one should start the algorithm from the points M p j 0 ( ). Since the actual positions of M j (p) are close, the algorithm converges to them quickly. In this manner one can effectively trace the positions of LOCT points over the parameter value.
One can still argue that to find all the LOCT points even for one parameter value p 0 may be a practically impossible task for complex high-dimension systems. However, it is important to emphasize that this task is significantly simplified if the system is a network, i.e. it is composed of many low-dimensional subsystems. In this case the coupling strength is a natural choice for the parameter p one can trace the system over. For the case of no coupling (p=0) the high-dimension phase space of the whole system is just a direct product of the lowdimension subspaces. In each of these subspaces the LOCT points can be found relatively easily which allows to spot them in the full space as well. Further one can gradually increase the coupling and trace the points.
Another issue is the possible emergence of new LOCT points. We have an effective method to trace once found points over the parameter, but how can one assure that no other points have emerged closer to the attractor? This problem is typical for the global extremum seeking in the case when full sampling of the space is impossible [9]. Although there is no way to guarantee that the found ST point is indeed the global minimum, there is a way to make the probability of that arbitrarily close to one. For this sake one has to check sufficiently many random initial points and make sure that none of them gives a better result. We revisit this issue at the end of section 4 where we provide an estimate of the number of trials one should run to decrease the mistake probability below a given level.

Stability threshold of pendulum
Now we show how the ST approach can be applied to study some paradigmatic dynamical systems. First we consider a classic pendulum under an external force P: Here, θ and ω are the deviation angle and the angular velocity, and α describes friction. Noteworthy models similar to (3) are often used to describe the dynamics of nodes of power grids, i.e. generators or consumers [4,10]. The phase space of the model is a cylinder S R 1 1 and includes two attractors: a stable steady state O P arcsin , 0 ( )and a stable limit cycle L ( figure 2(a)). In the context of power grids, the steady state corresponds to the state when the generator operates in synchrony with the grid, and the limit cycle corresponds to an undesired asynchronous regime.
Next we use the concept of ST to study the attraction basin of the steady state O. The identified LOCT points are depicted by red dots in figure 2(a). The most important ones are M 1 corresponding to positive perturbations and M 2 corresponding to negative ones. Because of the complex shape of the attraction basin other LOCT points exist further from the attractor, e.g. M .

È =
). When P increases, the basin stability for all classes of perturbations decreases as well as the ST. Thus, both measures indicate that the system becomes less robust. However, basin stability fails to detect which perturbations are more dangerous: S B2 is sufficiently larger than S B1 for all values of P.
We also checked that the efficiency of our algorithm is essentially improved by tracing LOCT points over the parameter. For higher-dimensional systems the improvement is even much higher. Note also that T c in the third setup increases by less than 30% with respect to the second setup, although the number of datapoints is twice as large. The reason is that with a smaller parameter step the positions of the LOCT points change less and they are found faster.

Stability threshold of complex networks
The second example is a network of coupled one-dimensional maps. We chose maps for two reasons: first, because of simpler implementation, and second, to demonstrate the generality of our approach. The network on N nodes is governed as follows: Here, a 0 1 < < is the system parameter, κ stands for the global coupling coefficient and c ij are the elements of the coupling matrix. Coupling between two nodes i and j equals c ij k . The network has the only attractor, the stable fixed point O 0, 0, , 0 ¼ ( ) . However, after a large perturbation the system trajectories may go to infinity. For network (4), a natural way to find LOCT points is to trace them over the coupling coefficient κ.
to negative ones. We start from these points for 0 k = , then gradually increase κ and trace their positions. We also periodically check for emergence of new LOCT points, but failed to detect any.
We study various networks with N 2 100   and different types of topology: all-to-all, random [11], small-world [12], scale-free [13], and cluster networks [14]. In all the cases, the behavior of LOCT points is quite similar. When κ increases, the positions of the points change, so that the coordinates  perturbations have a lower threshold than negative ones, and this threshold increases with .
i k This finding leads to an easy and intuitively clear rule: the stronger the node is connected to the network the harder it is to tear it off. However, too strong coupling ( i *  k k ) is undesirable, since it increases susceptibility to negative perturbations. The global ST of the network is defined by the lowest LOCT. Figure 3(b) illustrates a typical dependence of the ST on κ.

Stability threshold and basin stability
It is interesting to compare the two measures of stability against large perturbations, the ST and the basin stability for the same network ( figure 3(d)). As the perturbation class Q we use a hypersphere of radius q=0.8 with constant density ρ which means that we consider perturbations of amplitude q and random direction. [15] One may see that S 1 B = when the ST exceeds q for 0.26 q k k > » . This confirms that the ST indeed characterizes the weakest perturbation that can disrupt the network.
From figure 3(d) one may acquire the wrong impression that the basin stability reaches unity much earlier than κ reaches q k . The reason is that S B approaches unity very quickly when σ approaches q. This can be seen in the inset of figure 3(d) which has a logarithmic scale. This feature seems to be typical for high-dimensional dynamical systems. Indeed, consider an arbitrary dynamical system in the N-dimensional phase space settled into the attractor  with the ST .
s Consider the perturbation class Q consisting of perturbations with the amplitude q. For q s < , the set Q resides inside the attraction basin , therefore S 1 B = . For q s = , the set Q contacts the border of the basin  d . For q s > some part of the set Q gets out of the basin  and S B becomes smaller than one ( figure 4(a)). The probability of the perturbed state to be out of the basin is proportional to the surface area s of the protrusive part (gray in the figure), so The corresponding slope is given by the blue line in the inset of figure 4(d) and agrees with the numerical results. The scaling law (5) suggests that for high-dimensional systems it is very unlikely that the system will be disrupted by a random perturbation of magnitude just above the ST. From the other side, a wisely designed perturbation can disrupt the system even being slightly above the ST.
From computational prospective, the estimate (5) shows that attempts to estimate the ST from basin stability is inefficient for high-dimensional systems. Indeed, on the example from figure 4(d) one can see that it is very complicated to detect the exact point where the basin stability reaches exactly unity. The probability to randomly hit a point which is ε above the ST is of the order of N 1 2 e -( ) , so the time to determine the threshold with the accuracy ε by the brute-force strategy grows inversely. In contrast, the method suggested herein allows us to reach the LOCT point and determine the ST in quite a few steps.
The estimate ( In figure 4(b) the posterior probability is plotted versus the prior probability for several values of α. Note that if α is small enough and P(B) is not close to one (7) can be approximated as P B M P B 1 a » ( | ) ( ), which means that each trial decreases the probability by the factor α. Thus, after M trials the posterior probability of B can be estimated as The estimate (8) gives the probability that after M trials we have still missed a LOCT point whose distance to the attractor is by ds lower than 1 s . By sufficient increasing of M this mistake probability can be made arbitrarily small.

Conclusions and discussion
To conclude, we have introduced a novel measure to describe stability of dynamical systems against external perturbations. This is the ST which equals the magnitude of the weakest perturbation capable of disrupting the established dynamical regime. The ST provides important information, since it guarantees the system to withstand any perturbation of smaller magnitude. In the phase space, the ST is the minimal distance between the system's attractor and the border of its attraction basin. From this prospective, the ST defines the 'thinnest site' of the basin. And as the saying goes, where something is thin, that is where it tears: the direction corresponding to ST is the most dangerous for the system.
For dynamical networks, different directions in the multidimensional phase space are associated with different nodes. To this end, the ST approach allows us to determine the nodes which are mostly susceptible to perturbations. Applying external perturbations to these nodes, one may disrupt the network comparatively easily. However, sometimes the ST is associated with perturbations involving several nodes. An example of such a situation is depicted in figure 3(a). Under such circumstances, it is easier to disrupt the network by simultaneous perturbation of several nodes rather than by perturbing just one of them.
We have also suggested an algorithm to calculate the ST for arbitrary dynamical systems and demonstrated its effectiveness. The generality of the ST-based approach defines its vast potential for applications. Possible fields include engineering, neuroscience, power grids, Earth science and many others where robustness of complex systems against large perturbations is important.
Besides application to the study of particular systems, further development and extension of the approach per se are of great interest. To this end, several directions are seen. The most obvious modification of the suggested approach is to apply it to problems opposite to the one considered herein. Particularly, in many applications the task is to change the behavior of the system by an external action [16]. From the nonlinear dynamics prospective, this means pushing the system from one attractor to the basin of another. In this situation, quantification of the ST immediately provides the information about the optimal perturbation to induce the switching.
Another possibility is to extend the concept of the ST to multistable systems and reversible transitions between different attractors. In many applications, not a single, but several different stable regimes play essential roles, and the system may switch these regimes. The examples are switching between various patterns of activity in neural circuits [14,17], and episodic emergence of El Niño events [18], transitions between free flow and congestion in traffic [19] or between failure states and active states in self-recovering networks [20]. In all these cases, transitions in both directions are of interest. To study such transitions from one attractor to another and backwards, two STs can be introduced associated with each of the transitions. The values of these thresholds may provide important information about the transition rates when the system is externally perturbed.
Finally, the numerical algorithm can be further developed and extended for a broader class of systems. In particular, the current version of the algorithm relies on smoothness of attraction basins borders, which is a limitation of the method. It is known that basin boundaries for some dynamical systems can be fractal [21]. In this case the border can not be locally approximated by a tangential hyperplane, which is crucial to determine the direction in which to move along it. However, it is still possible to move along the border if the steps are taken in random direction. Each step is accepted if it brings the point closer to the attractor and declined otherwise. This modification of the method requires additional investigation to the define the conditions of its applicability and estimate its convergence.