Front. Netw. Physiol. Frontiers in Network Physiology Front. Netw. Physiol. 2674-0109 Frontiers Media S.A. 838142 10.3389/fnetp.2022.838142 Network Physiology Original Research Towards a Data-Driven Estimation of Resilience in Networked Dynamical Systems: Designing a Versatile Testbed Fischer et al. Data-Driven Estimation of Resilience Fischer Tobias 1 2 * Rings Thorsten 1 2 Rahimi Tabar M. Reza 3 4 Lehnertz Klaus 1 2 5 * 1 Department of Epileptology, University of Bonn Medical Centre, Bonn, Germany 2 Helmholtz Institute for Radiation and Nuclear Physics, University of Bonn, Bonn, Germany 3 Department of Physics, Sharif University of Technology, Tehran, Iran 4 Institute of Physics, Carl von Ossietzky University of Oldenburg, Oldenburg, Germany 5 Interdisciplinary Center for Complex Systems, University of Bonn, Bonn, Germany

Edited by: Eckehard Schöll, Technical University of Berlin, Germany

Reviewed by: Alexander N. Pisarchik, Polytechnic University of Madrid, Spain

Ralph G. Andrzejak, Pompeu Fabra University, Spain

 *Correspondence: Tobias Fischer, tobias.fischer@uni-bonn.de; Klaus Lehnertz, Klaus.Lehnertz@ukbonn.de

This article was submitted to Networks of Dynamical Systems, a section of the journal Frontiers in Network Physiology

 17 03 2022 2022 2 838142 17 12 2021 03 02 2022 Copyright © 2022 Fischer, Rings, Rahimi Tabar and Lehnertz. 2022 Fischer, Rings, Rahimi Tabar and Lehnertz

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Estimating resilience of adaptive, networked dynamical systems remains a challenge. Resilience refers to a system’s capacity “to absorb exogenous and/or endogenous disturbances and to reorganize while undergoing change so as to still retain essentially the same functioning, structure, and feedbacks.” The majority of approaches to estimate resilience requires exact knowledge of the underlying equations of motion; the few data-driven approaches so far either lack appropriate strategies to verify their suitability or remain subject of considerable debate. We develop a testbed that allows one to modify resilience of a multistable networked dynamical system in a controlled manner. The testbed also enables generation of multivariate time series of system observables to evaluate the suitability of data-driven estimators of resilience. We report first findings for such an estimator.

 resilience coupled oscillators complex networks time series analysis testbed FitzHugh-Nagumo 1 Introduction

The term resilience is commonly used to describe the ability of a system to return to its normal condition after its state has been perturbed. It is closely related to the notion of local stability (Lyapunov, 1892; Holling and Goldberg, 1971). When dealing with adaptive dynamical systems, the nonlocal stability concept of ecological resilience is increasingly employed throughout a number of scientific disciplines. With this concept, resilience refers to the system’s capacity “to absorb exogenous and/or endogenous disturbances and to reorganize while undergoing change so as to still retain essentially the same functioning, structure, and feedbacks” (Walker et al., 2004; Barabási and Posfai, 2016; Schoenmakers and Feudel, 2021). Despite the wide use of this concept, there is by now no commonly accepted method to assess resilience. Rather, a plethora of different indicators of resilience were proposed which are often not directly comparable [for an overview, see Meyer (2016) and Quinlan et al. (2016)]. Moreover, the vast majority of indicators of resilience require precise knowledge of the governing equations of motion and are thus of only limited value for the analysis of empirical data, such as for example time series of physiological observables of the human organism (Lehnertz et al., 2020; Romero-Ortuño et al., 2021).

Among the few data-driven indicators of resilience (or loss thereof) are the ones related to the concept of critical slowing down [see for example Scheffer et al. (2009); Dai et al. (2012); Lenton et al. (2012); Dakos et al. (2015); Scheffer et al. (2018)], namely lag-1 autocorrelation and variance (or other higher-order statistical indicators) estimated from time series of appropriate system observables. The suitability of these indicators has been extensively investigated with a variety of domain-specific models (see for example Weinans et al. (2021) and references therein). Nevertheless, a number of field applications indicate that these indicators are not always reliable (Ditlevsen and Johnsen, 2010; Boettiger and Hastings, 2013; Clements et al., 2019; Diks et al., 2019; Wilkat et al., 2019; Marconi et al., 2020; Hagemann et al., 2021; van der Bolt et al., 2021). In part, this can be attributed to the fact that the assumed mechanism behind critical slowing down (bifurcation-induced tipping) may be too simplistic for many natural systems which calls for data-driven indicators of resilience related to other transition scenarios (Kuehn, 2011; Ashwin et al., 2012; Kuehn et al., 2015; Ritchie and Sieber, 2017; Vanselow et al., 2019; O’Keeffe and Wieczorek, 2020; Hernández-Navarro et al., 2021; Schoenmakers and Feudel, 2021).

Another, more recently proposed, fully data-driven indicator of resilience – dynamical resistance (Rings et al., 2019) – takes into account that the dynamics of some elementary unit of a N-dimensional networked dynamical system can be described by its self-dynamics as well as by interactions with other units: x ̇ i = f x i + σ ∑ j = 1 N A i j h x j ; x i , where f ( x i ) determines the self-dynamics of unit i. The coupling between units i and j is defined by a coupling strength σ, a coupling matrix A (binary adjacency matrix), and a coupling function h , each of which can be time-dependent. With the ansatz of Rings et al. (2019), it is assumed that perturbations predominantly affect the dynamical coupling structure (second term on r.h.s of Eq. 1), so that a possible influence of a unit’s self-dynamics can be neglected. The authors proposed to approximate this structure with bivariate time-series-analysis techniques (Pikovsky et al., 2001; Kantz and Schreiber, 2003; Hlaváčková-Schindler et al., 2007; Marwan et al., 2007; Stankovski et al., 2017) in a time-resolved manner, and at the example of epileptic seizures in human brains, they demonstrated their approach to efficiently monitor dynamical resistance of a complex networked system prone to extreme events.

It is, however, not clear whether the assumptions underlying data-driven indicators of resilience are fully justified and whether indicators are generally applicable to estimate resilience of any real world system. An important step towards answering these questions is the development of a versatile testbed that would allow one to verify the reliability of data-driven estimators of resilience of networked dynamical systems. Here, we develop such a testbed that allows one to modify resilience of a multistable networked dynamical system in a controlled manner and to generate multivariate time series of system observables. We report first findings on the suitability of dynamical resistance as a data-driven indicator of resilience.

 2 Setting up the Testbed

Before going into details, let us first list some basic requirements for a testbed to be versatile:

1 in order to simulate a multistable system, a testbed should allow for an adjustable number of system states but more than two;

2 in order to simulate normal (desired) and aberrant (undesired) states, a testbed should generate distinguishable dynamics for each state;

3 in order to allow data-driven indicators of resilience to reliably characterize different system states (including those with fragmented boundaries), waiting times within each state should be sufficiently long;

4 in order to allow for a number of modifications of the system’s resilience, a testbed should have sufficiently many control parameters.

2.1 Dynamics: From Oscillators to Networks of Oscillator Networks

For our testbed, we consider one of the most simple and widely studied excitable systems, namely the FitzHugh–Nagumo (FHN) oscillator (also known as Bonhoeffer–van der Pol oscillator) (van der Pol and van der Mark, 1928; Bonhoeffer, 1948; FitzHugh, 1961; Nagumo et al., 1962; Rocşoreanu et al., 2000), which is a prototypical model for excitable behavior in neural and cardiac nonlinear activities (Glass et al., 1991; Koch, 1999). The equations of motion for the FHN oscillator read: x ̇ = x a − x x − 1 − y + I y ̇ = b x − c y , with x and y denoting the (fast) excitatory and (slow) inhibitory dynamical state variables. Here a, b, and c are dimensionless control parameters: a and c are scaling parameters, and b controls the emergence of various dynamical regimes (such as tonic and phasic spiking or sub-threshold oscillations). I is the magnitude of stimulus current.

Networks of coupled, non-identical FHN oscillators can exhibit much richer dynamics depending on the coupling and the coupling topology. Apart from various synchronization phenomena (Chernihovskyi and Lehnertz, 2007; Omelchenko et al., 2015; Plotnikov et al., 2016; Masoliver et al., 2017; Ramlow et al., 2019; Gerster et al., 2020), such networks are capable of generating extreme-event-like phenomena (Ansmann et al., 2013; Karnatak et al., 2014; Rings et al., 2017; Saha and Feudel, 2017) and self-induced switching between multiple space-time pattern (Ansmann et al., 2016). The complexity of dynamics can further be enhanced, if one considers networks of networks of (diffusively coupled) FHN oscillators (Rydin Gorjão et al., 2018). Here the equations of motion for the i-th oscillator ( i ∈ 1 , … , N o ; N o denotes the number of oscillators) in the k-th sub-network ( ( k , l ) ∈ 1 , … , N n ; N n denotes the number of sub-networks) read: x ̇ i k = x i k a i − x i k x i k − 1 − y i k + C w k N o − 1 ∑ j = 1 N o S i j x j k − x i k + 1 N ∑ l = 1 N n C b k , l B k l ∑ j = 1 N o x j l − x i k y ̇ i k = b i x i k − c i y i k , where C w ( k ) and C b ( k , l ) denote the global coupling strengths within and between sub-networks k and l. For a given sub-network, S ∈ 0,1 N o × N o denotes the symmetric adjacency matrix (S ij = S ji = 1, if oscillators i and j are coupled, else S ij = S ji = 0). The symmetric adjacency matrix B ∈ 0,1 N n × N n characterizes the coupling structure between sub-networks. For the numerical integration of such large systems of differential equations, we use the python module jitcode (Ansmann, 2018).

For our study, we consider N n = 2 fully connected sub-networks, each consisting of N o = 25 diffusively coupled non-identical FHN oscillators (FHN-NoN). The control parameters a and c are identical for all oscillators: a i = a = − 0.027 6 ∀i and c i = c = − 0.02 ∀i; parameter b i is linearly distributed on 0.006 , 0.014 to prevent immediate synchronization among oscillators for non-zero coupling strengths (cf. Ansmann et al. (2013) and Rydin Gorjão et al. (2018)). The set of coupling strengths β c = ( C w ( 1 ) , C w ( 2 ) , C b ) is adjustable (note that we dropped the superscripts k and k| for the sake of readability). As an example, we show in Figure 1 time series of the averaged dynamical variables x ̄ ( k ) ( t ) = ∑ i = 1 N o x i ( k ) ( t ) and y ̄ ( k ) ( t ) = ∑ i = 1 N o y i ( k ) ( t ) of sub-network 1 for various β c , following Ansmann et al. (2013); Karnatak et al. (2014); Ansmann et al. (2016), and Rydin Gorjão et al. (2018). For the single FHN oscillator, we obtain the two Lyapunov exponents as λ 1 = − 3.45 ± 6.14 ⋅ 10−6 and λ 2 = − 2.89 ± 0.01 ⋅ 10−1. For two coupled FHN, the three largest Lyapunov exponents amount to λ 1 = 4.07 ± 0.14 ⋅ 10−3, λ 2 = 5.04 ± 12.06 ⋅ 10−7, and λ 3 = − 7.00 ± 0.09 ⋅ 10−2. For the fully connected network with 25 FHN oscillators, we obtain for the three largest Lyapunov exponents λ 1 = 5.91 ± 0.23 ⋅ 10−3, λ 2 = 4.81 ± 13.01 ⋅ 10−6, λ 3 = − 2.17 ± 0.05 ⋅ 10−2 and for the FHN-NoN λ 1 = 6.42 ± 0.16 ⋅ 10−3, λ 2 = 5.40 ± 8.85 ⋅ 10−6, λ 3 = − 5.94 ± 0.78 ⋅ 10−4. These Lyapunov exponents (Benettin et al., 1980a,b) were derived from 20 realizations of the systems with different initial conditions.

Exemplary time series and corresponding phase space representations of the dynamics of a single FHN oscillator (β c = (0, 0, 0)), two coupled FHN (β c = (0.130 8, 0, 0)), a fully connected network with 25 FHN oscillators (β c = (0.115, 0, 0)), and a network of two networks with 50 FHN oscillators (β c = (0.115, 0.215, 9.5 ⋅ 10−5)).

 2.2 Modeling a Multistable System

The next building block of our testbed is the modeling of a multistable system, in which transitions between more than two states (dynamical regimes) are not induced by a change of control parameters but occur in a self-induced manner. To this end, we consider a potential landscape with transitions between states that are mediated by the rare recurring, high-amplitude phenomena generated by the network of networks of FHN oscillators. We approximate the potential landscape L(z) by a succession of N s potential wells (modeled as inverted Gaussian functions; z is the dynamical variable that describes the motion within the potential landscape), that mimic the basins of attraction of the N s states of our multistable system (see upper part of Figure 2): L z = ∑ n = 1 N s a n G 2 π σ n G 2 e x p − z − μ n G 2 2 σ n G 2 .

(A): exemplary potential landscape with four system states. The potential landscape is driven by the mean field of FHN-NoN and, together, the system exhibits self-induced transitions. (B): Exemplary phase space representations of the averaged dynamical variables for each state. (C,D): Exemplary time series of the averaged dynamical variable x ̄ ( 1 ) ( t ) of the FHN-NoN and of the dynamical variable z(t) of the potential landscape. In (C), Γ(t) is a random sequence of + 1s and − 1s and in (D), Γ(t) is chosen state-dependent such that the system evolves alternating from the first to the last state (S1 → S4: Γ(t) = + 1) and back (S4 → S1: Γ(t) = − 1). Settings of other parameters and conditions: random initial conditions for x i ( k ) × y i ( k ) ∈ [ 0,1 ] × [ 0,1 ] and for z ∈ [ − 5, 5]; a G = [5.5, 5.5, 5.5, 5.5]; μ G = [ − 6, − 2, 2, 6]; ( σ G ) 2 = [ 1.8 , 1.8 , 1.8 , 1.8 ] ; for all states: C w ( 1 ) = 0.115 , for state S1: C w ( 2 ) = 0.116 and C b = 1.045 ⋅ 10−5; for state S2: C w ( 2 ) = 0.116 and C b = 1.425 ⋅ 10−5; for state S3: C w ( 2 ) = 0.125 and C b = 1.045 ⋅ 10−5; for state S4: C w ( 2 ) = 0.125 and C b = 1.5 ⋅ 10−5. The initial transient (105 data points) was dismissed.

Here a n G , μ n G , and σ n G denote amplitude, mean and standard deviation of the n-th Gaussian function. Together with the distance Δz between potential wells, these control parameters allow one to modify resilience of a system (see Mitra et al. (2015) and references therein). We couple the FHN-NoN’s mean field M(t) to the potential landscape and derive the equation which governs the motion of z: z ̇ t = L ̇ z z t + ζ t M t , where M t = Γ t N o ∑ i = 1 N o x i 1 C w 1 , C b , t + x i 2 C w 2 , C b , t and ζ(t) = (0.6 − 0.01z(t)) is a scaling factor. The second term on the r.h.s. of Eq. 5 enables the switching of the dynamics into different potential wells mediated by the high-amplitude phenomena generated by the FHN-NoN. This self-induced switching thus solely depends on the amplitude of the FHN-NoN’s mean field M(t) and is not mediated by a change of control parameters of the potential landscape. In order to guarantee distinguishable dynamics for the different potential wells, we change the coupling strengths C w ( k ) and C b once the mean field M(t) drives z(t) to overcome the local barrier between the respective wells. The asymmetric amplitude distribution (cf. Figure 1) of the mean field dominates the sequence in which states are switched. As an example, setting Γ(t) = + 1 ∀t results in a sequence of state switches from the first to the last potential well. More complicated sequences of state switches can be achieved by setting Γ(t) (with Γ ( t ) = 1 ∀ t ) appropriately. For example, choosing the sign of a set of random numbers drawn from some distribution centered around zero results in a random sequence of state switches. The sign of Γ(t) may also be chosen depending on the current state of the system. In the upper part of Figure 2, we show for a 4-state system exemplary phase space representations of the averaged dynamical variables for each state. In the lower part, we show for two choices of Γ(t) exemplary time series of the dynamical variable z(t) of the potential landscape and of the averaged dynamical variable x ̄ ( 1 ) ( t ) of sub-network 1.

There are alternative ways to model a potential landscape (Hänggi et al., 1990), for example using an n-th order polynomial. Advantages of using a succession of inverted Gaussian functions include the simple and intuitive way of adding further potential wells, thereby retaining the order of the wells. A disadvantage is the smooth barrier between potential wells, which may result in a rapid escape from a well once the above mentioned escape condition is fulfilled, and thus to short transition times (cf. Figure 2). This can be avoided by using, for example, fragmented barriers that can be constructed using the classic Cantor fractal construction process (see, e.g., Mandelbrot (1982); Omelchenko et al. (2015)). Such fragmented barriers may also mimic riddled basin of attractions (Alexander et al., 1992). Another way to achieve a non-smooth barrier would be adding e.g., colored noise to the potential landscape.

An example of a potential landscape with fragmented barriers is shown in Figure 3 along with time series of the averaged dynamical variable x ̄ ( 1 ) ( t ) of sub-network 1 of the FHN-NoN and of the dynamical variable z(t) of the potential landscape. The inclusion of a fragmented barrier can be regarded as adding “intermediate states” that temporarily trap the system. Note that the dynamics within these intermediate states differs from the ones observed in the potential wells. Figure 4 provides a synopsis of the accumulated waiting times of the FHN-NoN dynamics within each state and demonstrates how the steepness of fragmented barriers impacts on the transition time between states.

(A): same potential landscape as in Figure 2 but with fragmented barriers, generated by using the classic Cantor fractal construction process. We bridge “gaps” in the barriers that result from the fractal construction process by inserting segments with adjustable steepness α. The fragmented barrier begins and ends outside of the local minimum of a potential well such that an amplitude-based transition can be induced from the potential well into the barrier. (B,C): Exemplary time series of the averaged dynamical variable x ̄ ( 1 ) ( t ) of the FHN-NoN and of the dynamical variable z(t) of the potential landscape.

Distributions of accumulated waiting times of the FHN-NoN dynamics within the potential wells (states S1–S4) of the landscape with fragmented barriers shown in Figure 3 and for Γ(t) chosen state-dependent such that the system evolves alternating from the first to the last state and back. Distributions (kernel density estimates) were derived from 10 realizations with random initial conditions of oscillators and of the dynamical variable z(t) of the landscape. Other control parameters as in Figure 3. The larger waiting times within states S2 and S3 result from these states being visited, on average, more often than the other states as they each can be reached via two transitions. The inset shows the dependence of the average transition time between potential wells on the steepness α of the fragmented barrier. The average transition time was estimated over 20 realizations with different initial conditions. Lines are for eye-guidance only.

 2.3 Modifying the System’s Resilience: An Example

As already mentioned above, several control parameters allow one to modify the resilience of our multistable system which is briefly illustrated in the following. We consider a system as in Figure 3 with three desired states (S1, S2, S3) representing its normal functioning and one undesired state (S4) representing an aberrant functioning. The parameters controlling the distance between wells S1, S2 and S3 as well as the height of the barriers between these wells are identical but differ from those of well S4. We allow for a state-dependent switching between states, and by gradually moving S4 closer to S3 (i.e., decreasing the distance Δz between these states), we mimic a progressive loss of the system’s resilience. Since the height of the barrier between S3 and S4 is enlarged, the system is trapped in S4 once it enters this state. A schematic of this modification along with exemplary excerpts of time series of observables from some oscillators are presented in Figure 5.

Schematic of a modification of resilience of a multistable system with three desired states (S1, S2, S3) and one undesired state (S4) (A) and excerpts of exemplary time series of system observables from various states (B). Setting of control parameters of the potential landscape: a G = [36.5, 36.5, 36.5, 36.5]; μ G = [ − 18.3, − 6.2, 6.2, 22.3]; ( σ G ) 2 = [ 12,12,12,5 ] . Coupling constant for all states: C w ( 1 ) = 0.115 , for state S1: C w ( 2 ) = 0.116 and C b = 1.045 ⋅ 10−5; for state S2: C w ( 2 ) = 0.116 and C b = 1.425 ⋅ 10−5; for state S3: C w ( 2 ) = 0.125 and C b = 1.045 ⋅ 10−5; for state S4: C w ( 2 ) = 1.25 and C b = 1.5 ⋅ 10−5. Initial configuration at time interval t ∈ 0,1 0 5 : Self-induced transitions are possible between states S1, S2, and S3 but a transition into the undesired state S4 is not possible. Starting at time t = 105, the position μ 4 G is linearly decreased up to μ 4 G = 19.3 at time t = 2.5 ⋅ 10-5, which decreases the distance Δz S3,S4 between S3 and S4. The lower part of the figure exemplifies excerpts of time series of the x-components of oscillators 1, 3, 5, 10. and 25 (increasing value of control parameter b i ) of sub-network 2 at time intervals indicated in the upper part. The excerpts are labeled according to the system’s state, and T1 and T2 indicate excerpts from transitory phases.

Before closing this section, we briefly summarize the main aspects of our testbed that allow us to modify resilience of a multistable networked dynamical system in a controlled manner. Our testbed provides the means to simulate the dynamics of a multistable system with the help of a network of networks of FitzHugh-Nagumo oscillators coupled to a potential landscape that consists of a succession of a number of potential wells with smooth or fragmented barriers. Various control parameters allow one to generate distinguishable dynamics for each (desired or undesired) state, to adjust the waiting time of the system within each state, as well as the transition time between states. Our testbed also allows for a generation of time series of system observables, and these time series may serve as input to data-driven indicators of resilience.

 3 An Exemplary Evaluation of a Data-Driven Indicator of Resilience

In the following, we utilize time series generated by our testbed for an exemplification of a data-driven indicator of resilience. Rings et al. (2019) proposed a time-series-based and non-perturbative approach to efficiently monitor dynamical resistance, an indicator of resilience of a networked dynamical system. The approach is fully data-driven since it takes into account the units’ individual signals only and consists of the following three central steps of analysis:

1 Probe with high temporal resolution the dynamical coupling structure between interacting system units;

2 identify dynamical regimes (here: states) from similar time-dependent coupling structures;

3 define dynamical resistance R as the minimum distance between all accessible dynamical regimes.

Step of analysis # 1: The dynamical coupling structure (the second term on the r.h.s. of Eq. 1; coupling strength, coupling structure, and coupling function) can be probed with bivariate time-series-analysis techniques developed in statistics, nonlinear dynamics, information and synchronization theory as well as in statistical physics (Pikovsky et al., 2001; Kantz and Schreiber, 2003; Pereda et al., 2005; Hlaváčková-Schindler et al., 2007; Marwan et al., 2007; Stankovski et al., 2017; Tabar, 2019). Here we use three widely-used techniques, namely the zero- and maximum-lag cross-correlation (Brillinger, 1981; Rosenblum et al., 1997) as well as the (normalized) mutual information (Kraskov et al., 2004). These techniques allow one to estimate the similarity/interdependence ρ uv between pairs of time series u and v each of length T (with T much smaller than the total observation time). If appropriately normalized, ρ uv assumes values between 0 and 1, indicating either complete independence or complete dependence. We use a sliding window approach to calculate ρ uv between all pairs of units in a time-resolved manner which results in a temporal sequence of interaction matrices ρ.

Step of analysis # 2: In order to identify dynamical regimes, one can define similarity between two interaction matrices ρ (t l ) and ρ (t m ) at times t l and t m as ξ (t l , t m ) ≡‖ρ (t l ) − ρ (t m )‖, where ‖… ‖ denotes the Euclidean norm (Münnix et al., 2012). The similarity matrix ξ—estimated for all times t m and t l —then contains pertinent information about the system’s dynamics, and recurrent patterns in the similarity matrix indicate dynamical regimes (Marwan et al., 2007). In order to identify these regimes and their number, Rings et al. (2019) proposed to use a time-resolved hierarchical clustering analysis of coupling structures in an abstract space spanned by all pairwise interactions. For our investigations, we use a k-means algorithm (MacQueen, 1967) given that the number of different dynamical regimes (clusters k = N s) is known a priori.

Step of analysis # 3: The minimum Euclidean separation between cluster centroids is taken as the minimum distance between dynamical regimes and is an estimate for dynamical resistance R: the larger this distance between regimes the higher is the capacity of a system to absorb disturbances and to reorganize while undergoing dynamical changes so as to still retain essentially the same functionality.

For our investigations, we consider the example from Section 2.3 and simulate a gradual loss of resilience of the system, which we assume to result from a “perturbation” mediated by the undesired state. We again incrementally decrease the distance Δz S3,S4 between states S3 and S4, and for each increment we record time series of observables of each FHN oscillator for 106 time steps thereby starting from identical initial conditions of each oscillator. For our analyses, we use the oscillators’ x-components from sub-network 2 that we observe using the identity as measurement function.

In Figure 6, we show how the shortening of the distance Δz S3,S4 impacts on the waiting times within each state. As expected, the median waiting time within state S4 increased upon decreasing Δz S3,S4, while within state S3 the median waiting time gradually decreased. Waiting times within states S1 and S2 remained largely unaffected by the perturbation. In the upper part of Figure 7, we summarize our findings for dynamical resistance R. Depending on the bivariate time-series-analysis technique employed to estimate R, we observe the initial resilience of the system to be diminished by about 10% as S4 gets closer to S3. Our interpretation of this loss of resilience due to a “perturbation” mediated by an undesired state is further corroborated by the distinct increased area between cluster centroids reflecting a deformation of the initial configuration of the system’s dynamical regimes (lower part of Figure 7).

Impact of decreasing the distance Δz S3,S4 between states S3 and S4 on the waiting times in each of the four states (S1–S4). Medians and variances (lengths of error bars) obtained from 10 realizations of the simulation setup. A zero value of the waiting time of S4 indicates that this state is never reached because it is too far away from S3. The vertical dashed line indicates onset of perturbation. Other lines are for eye-guidance only.

Relative change of dynamical resistance [δR; (A)] and of area between cluster centroids [δV; (B)] upon decreasing the distance Δz S3,S4 between states S3 and S4. Dynamical resistance R estimated with mutual information (filled circles), zero-lag cross-correlation (filled triangles), and maximum-lag cross-correlation (filled squares). Area V between cluster centroids calculated from the distances between the three cluster centroids representing states S1–S3 (symbols as above). Data normalized to the respective values for large Δz S3,S4. Medians and variances (lengths of error bars) obtained from 10 realizations of the simulation setup. The vertical dashed line indicates onset of perturbation. Other lines are for eye-guidance only.

 4 Conclusion

We develop a testbed that allows one to modify resilience of a multistable networked dynamical system in a controlled manner and to generate time series of observables that may be used to evaluate the suitability of data-driven indicators of resilience. The testbed presented here was designed in such a way that it provides a means—with the help of an adjustable potential landscape, sufficiently many and, more importantly, contextual control parameters—to simulate a multistable system with a number of system states and with self-induced switching between them as well as to simulate distinguishable dynamics for each state. Waiting times within states are sufficiently long to allow data-driven indicators of resilience to reliably characterize the different states. With the inclusion of fragmented barriers between potential wells, transition times between states can be preset. We here considered a potential landscape—consisting of a succession of adjustable wells—driven by the rich dynamics of a network of networks of diffusively coupled FitzHugh-Nagumo (FHN) oscillators (Rydin Gorjão et al., 2018). The network’s dynamics is chaotic and can exhibit different dynamical patterns such as low-amplitude oscillations, nonlinear waves, and rare recurring high-amplitude phenomena (Ansmann et al., 2013). The network can also exhibit self-induced switchings between these patterns (states) without a change of control parameters (Ansmann et al., 2016). Our testbed allows for short computation times and can modify resilience during run time. As an example, the generation of time series of observables of a fully connected network of two networks of 50 FHN oscillators with 106 data points each requires about 3 min on a 64-bit architecture with a single CPU at 2.2 GHz.

Using time series data generated by our testbed for a multistable system gradually perturbed by an undesired state, we performed an exemplary evaluation of a data-driven indicator of resilience of a networked dynamical system (Rings et al., 2019). Our findings indicate that this indicator—dynamical resistance R—appears to be capable of tracking changes in resilience, at least to some extent and for the scenario considered here. Nevertheless, findings also indicate that its performance appears to depend on the bivariate time-series-analysis technique employed to characterize couplings between system units. Future studies would need to address the question as to which extent the influence of a unit’s self-dynamics can be neglected when estimating resilience of a networked dynamical system. In addition, future studies would need to tackle the largely unsolved issue of how to reliably interpret findings obtained with data-driven indicators particularly with respect to Holling’s definition of resilience.

We foresee various extensions to our testbed, also in view of evaluating other data-driven indicators. For example, one may consider other configurations of the potential landscape, other models of networked dynamics, coupling and measurement functions that are of relevance for a given research field. Time-dependent control parameters [see, e.g., Nicolis and Nicolis (2014); O’Regan and Burton (2018)] for both the network dynamics and for the potential landscape will introduce various non-stationarities, thus bringing our testbed closer to natural systems. At a similar token, the introduction of stochasticity (Freidlin and Wentzell, 1984; Arnold, 1998) into our testbed may allow for various noise-related phenomena such as noise-induced transitions (Horsthemke and Lefever, 1984), stochastic resonance (Gammaitoni et al., 1998), or noise-induced tipping (Ritchie and Sieber, 2017; Wunderling et al., 2021). Their time-series-analysis-based investigation in networked dynamical system may, however, require more refined and better adapted analysis techniques (Rydin Gorjão et al., 2019, 2021; Tabar, 2019; Aslim et al., 2021).

 Data Availability Statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation. The code of this work has been provided as the open source software package resiland written in the programming language Python. It is freely available on github under the: https://github.com/tobfischer/resiland.

 Author Contributions

All authors conceived the research project and wrote the paper. All authors contributed to the article and approved the submitted version.

 Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with one of the authors (KL).

 Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

We are grateful to Timo Bröhl for interesting discussions and for critical comments on earlier versions of the manuscript.

 References Alexander J. C. Yorke J. A. You Z. Kan I. (1992). Riddled Basins. Int. J. Bifurcation Chaos 02, 795–813. 10.1142/s0218127492000446 Ansmann G. Karnatak R. Lehnertz K. Feudel U. (2013). Extreme Events in Excitable Systems and Mechanisms of Their Generation. Phys. Rev. E 88, 052911. 10.1103/PhysRevE.88.052911 Ansmann G. Lehnertz K. Feudel U. (2016). Self-induced Switchings between Multiple Space-Time Patterns on Complex Networks of Excitable Units. Phys. Rev. X 6, 011030. 10.1103/PhysRevX.6.011030 Ansmann G. (2018). Efficiently and Easily Integrating Differential Equations with JiTCODE, JiTCDDE, and JiTCSDE. Chaos 28, 043116. 10.1063/1.5019320 Arnold L. (1998). Random Dynamical Systems. Berlin: Springer. Ashwin P. Wieczorek S. Vitolo R. Cox P. (2012). Tipping Points in Open Systems: Bifurcation, Noise-Induced and Rate-dependent Examples in the Climate System. Phil. Trans. R. Soc. A. 370, 1166–1184. 10.1098/rsta.2011.0306 Aslim E. Rings T. Zabawa L. Lehnertz K. (2021). Enhancing the Accuracy of a Data-Driven Reconstruction of Bivariate Jump-Diffusion Models with Corrections for Higher Orders of the Sampling Interval. J. Stat. Mech. 2021, 033406. 10.1088/1742-5468/abe59a Barabási A. Posfai M. (2016). Network Science. 1st edn. Cambridge, UK: Cambridge University Press. Benettin G. Galgani L. Giorgilli A. Strelcyn J.-M. (1980a). Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems; a Method for Computing All of Them. Part 1: Theory. Meccanica 15, 9–20. 10.1007/BF02128236 Benettin G. Galgani L. Giorgilli A. Strelcyn J.-M. (1980b). Lyapunov Characteristic Exponents for Smooth Dynamical Systems and for Hamiltonian Systems; a Method for Computing All of Them. Part 2: Numerical Application. Meccanica 15, 21–30. 10.1007/BF02128237 Boettiger C. Hastings A. (2013). No Early Warning Signals for Stochastic Transitions: Insights from Large Deviation Theory. Proc. R. Soc. B. 280, 20131372. 10.1098/rspb.2013.1372 Bonhoeffer K. F. (1948). Activation of Passive Iron as a Model for the Excitation of Nerve. J. Gen. Physiol. 32, 69–91. 10.1085/jgp.32.1.69 Brillinger D. (1981). Time Series: Data Analysis and Theory. San Francisco, USA: Holden Day. Chernihovskyi A. Lehnertz K. (2007). Measuring Synchronization with Nonlinear Excitable Media. Int. J. Bifurcation Chaos 17, 3425–3429. 10.1142/s0218127407019159 Clements C. F. McCarthy M. A. Blanchard J. L. (2019). Early Warning Signals of Recovery in Complex Systems. Nat. Commun. 10, 1681. 10.1038/s41467-019-09684-y Dai L. Vorselen D. Korolev K. S. Gore J. (2012). Generic Indicators for Loss of Resilience before a Tipping point Leading to Population Collapse. Science 336, 1175–1177. 10.1126/science.1219805 Dakos V. Carpenter S. R. van Nes E. H. Scheffer M. (2015). Resilience Indicators: Prospects and Limitations for Early Warnings of Regime Shifts. Phil. Trans. R. Soc. B 370, 20130263. 10.1098/rstb.2013.0263 Diks C. Hommes C. Wang J. (2019). Critical Slowing Down as an Early Warning Signal for Financial Crises? Empir. Econ. 57, 1201–1228. 10.1007/s00181-018-1527-3 Ditlevsen P. D. Johnsen S. J. (2010). Tipping Points: Early Warning and Wishful Thinking. Geophys. Res. Lett. 37, L19703. 10.1029/2010gl044486 FitzHugh R. (1961). Impulses and Physiological States in Theoretical Models of Nerve Membrane. Biophysical J. 1, 445–466. 10.1016/S0006-3495(61)86902-6 Freidlin M. I. Wentzell A. D. (1984). Random Perturbations of Dynamical Systems. New York: Springer Gammaitoni L. Hänggi P. Jung P. Marchesoni F. (1998). Stochastic Resonance. Rev. Mod. Phys. 70, 223–287. 10.1103/RevModPhys.70.223 Gerster M. Berner R. Sawicki J. Zakharova A. Škoch A. Hlinka J. (2020). FitzHugh–Nagumo Oscillators on Complex Networks Mimic Epileptic-Seizure-Related Synchronization Phenomena. Chaos: Interdiscip. J. Nonlinear Sci. 30, 123130. 10.1063/5.0021420 Glass L. Hunter P. McCulloch A. (1991). Theory of Heart. New York: Springer. 10.1007/978-1-4612-3118-9 Hänggi P. Talkner P. Borkovec M. (1990). Reaction-rate Theory: Fifty Years after Kramers. Rev. Mod. Phys. 62, 251. 10.1103/revmodphys.62.251 Hagemann A. Wilting J. Samimizad B. Mormann F. Priesemann V. (2021). Assessing Criticality in Pre-seizure Single-Neuron Activity of Human Epileptic Cortex. Plos Computat. Biol. 17, e1008773. 10.1371/journal.pcbi.1008773 Hernández-Navarro L. Faci-Lázaro S. Orlandi J. G. Feudel U. Gómez-Gardeñes J. Soriano J. (2021). Noise-driven Amplification Mechanisms Governing the Emergence of Coherent Extreme Events in Excitable Systems. Phys. Rev. Res. 3, 023133. 10.1103/physrevresearch.3.023133 Hlaváčková-Schindler K. Paluš M. Vejmelka M. Bhattacharya J. (2007). Causality Detection Based on Information-Theoretic Approaches in Time Series Analysis. Phys. Rep. 441, 1–46. 10.1016/j.physrep.2006.12.004 Holling C. S. Goldberg M. A. (1971). Ecology and Planning. J. Am. Planners 37, 221–230. 10.1080/01944367108977962 Horsthemke W. Lefever R. (1984). Noise-Induced Transitions. Theory and Applications in Physics, Chemistry and Biology. Berlin: Springer. Kantz H. Schreiber T. (2003). Nonlinear Time Series Analysis. 2nd edn. Cambridge, UK: Cambridge University Press. 10.1017/CBO9780511755798 Karnatak R. Ansmann G. Feudel U. Lehnertz K. (2014). Route to Extreme Events in Excitable Systems. Phys. Rev. E 90, 022917. 10.1103/PhysRevE.90.022917 Koch C. (1999). Biophysics of Computation: Information Processing in Single Neurons. New York: Oxford University Press Computational Neuroscience. Kraskov A. Stögbauer H. Grassberger P. (2004). Estimating Mutual Information. Phys. Rev. E 69, 066138. 10.1103/PhysRevE.69.066138 Kuehn C. Zschaler G. Gross T. (2015). Early Warning Signs for Saddle-Escape Transitions in Complex Networks. Sci. Rep. 5, 13190. 10.1038/srep13190 Kuehn C. (2011). A Mathematical Framework for Critical Transitions: Bifurcations, Fast-Slow Systems and Stochastic Dynamics. Physica D 240, 1020–1035. 10.1016/j.physd.2011.02.012 Lehnertz K. Bröhl T. Rings T. (2020). The Human Organism as an Integrated Interaction Network: Recent Conceptual and Methodological Challenges. Front. Physiol. 11, 1694. 10.3389/fphys.2020.598694 Lenton T. Livina V. Dakos V. Van Nes E. Scheffer M. (2012). Early Warning of Climate Tipping Points from Critical Slowing Down: Comparing Methods to Improve Robustness. Phil. Trans. Roy. Soc. A: Math. Phys. Eng. Sci. 370, 1185–1204. 10.1098/rsta.2011.0304 Lyapunov A. M. (1892). The General Problem Of Stability Of Motion (In Russian) Doctoral dissertation Kharkov: University of Kharkov, Kharkov Mathematical Society. MacQueen J. B. (1967). “Some Methods for Classification and Analysis of Multivariate Observations,” in Fifth Berkeley Symposium on Mathematical Statistics and Probability, June 21–July 18, 1965 and December 27, 1965–January 7, 1966. Editors Cam M. L. Neyman J. (Berkeley, USA: Statistical Laboratory of the University of California), 281–297. Mandelbrot B. B. (1982). The Fractal Geometry of Nature. San Francisco: Freeman.: Marconi M. Métayer C. Acquaviva A. Boyer J. Gomel A. Quiniou T. (2020). Testing Critical Slowing Down as a Bifurcation Indicator in a Low-Dissipation Dynamical System. Phys. Rev. Lett. 125, 134102. 10.1103/physrevlett.125.134102 Marwan N. Romano M. C. Thiel M. Kurths J. (2007). Recurrence Plots for the Analysis of Complex Systems. Phys. Rep. 438, 237–329. 10.1016/j.physrep.2006.11.001 Masoliver M. Malik N. Schöll E. Zakharova A. (2017). Coherence Resonance in a Network of FitzHugh-Nagumo Systems: Interplay of Noise, Time-Delay, and Topology. Chaos: Interdiscip. J. Nonlinear Sci. 27, 101102. 10.1063/1.5003237 Meyer K. (2016). A Mathematical Review of Resilience in Ecology. Nat. Resour. Model. 29, 339–352. 10.1111/nrm.12097 Mitra C. Kurths J. Donner R. V. (2015). An Integrative Quantifier of Multistability in Complex Systems Based on Ecological Resilience. Sci. Rep. 5, 16196. 10.1038/srep16196 Münnix M. C. Shimada T. Schäfer R. Leyvraz F. Seligman T. H. Guhr T. (2012). Identifying States of a Financial Market. Sci. Rep. 2, 644. 10.1038/srep00644 Nagumo J.-I. S. Arimoto S. Yoshizawa S. (1962). An Active Pulse Transmission Line Simulating Nerve Axon. Proc. IRE 50, 2061–2070. 10.1109/jrproc.1962.288235 Nicolis C. Nicolis G. (2014). Dynamical Responses to Time-dependent Control Parameters in the Presence of Noise: a normal Form Approach. Phys. Rev. E 89, 022903. 10.1103/PhysRevE.89.022903 O’Keeffe P. E. Wieczorek S. (2020). Tipping Phenomena and Points of No Return in Ecosystems: beyond Classical Bifurcations. SIAM J. Appl. Dyn. Syst. 19, 2371–2402. 10.1137/19m1242884 O’Regan S. M. Burton D. L. (2018). How Stochasticity Influences Leading Indicators of Critical Transitions. Bull. Math. Biol. 80, 1630–1654. 10.1007/s11538-018-0429-z Omelchenko I. Provata A. Hizanidis J. Schöll E. Hövel P. (2015). Robustness of Chimera States for Coupled FitzHugh-Nagumo Oscillators. Phys. Rev. E 91, 022917. 10.1103/PhysRevE.91.022917 Pereda E. Quian Quiroga R. Bhattacharya J. (2005). Nonlinear Multivariate Analysis of Neurophysiological Signals. Prog. Neurobiol. 77, 1–37. 10.1016/j.pneurobio.2005.10.003 Pikovsky A. S. Rosenblum M. G. Kurths J. (2001). Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge, UK: Cambridge University Press. Plotnikov S. Lehnert J. Fradkov A. Schöll E. (2016). Synchronization in Heterogeneous FitzHugh-Nagumo Networks with Hierarchical Architecture. Phys. Rev. E 94, 012203. 10.1103/PhysRevE.94.012203 Quinlan A. E. Berbés-Blázquez M. Haider L. J. Peterson G. D. (2016). Measuring and Assessing Resilience: Broadening Understanding through Multiple Disciplinary Perspectives. J. Appl. Ecol. 53, 677–687. 10.1111/1365-2664.12550 Ramlow L. Sawicki J. Zakharova A. Hlinka J. Claussen J. C. Schöll E. (2019). Partial Synchronization in Empirical Brain Networks as a Model for Unihemispheric Sleep. EPL (Europhys. Lett.) 126, 50007. 10.1209/0295-5075/126/50007 Rings T. Ansmann G. Lehnertz K. (2017). How Important Are Hubs for the Generation of Extreme Events in Networks of Excitable Units? Eur. Phys. J.-Spec. Top. 226, 1963–1970. 10.1140/epjst/e2017-70021-3 Rings T. Mazarei M. Akhshi A. Geier C. Tabar M. R. R. Lehnertz K. (2019). Traceability and Dynamical Resistance of Precursor of Extreme Events. Sci. Rep. 9, 1744. 10.1038/s41598-018-38372-y Ritchie P. Sieber J. (2017). Probability of Noise-And Rate-Induced Tipping. Phys. Rev. E 95, 052209. 10.1103/PhysRevE.95.052209 Rocşoreanu C. Georgescu A. Giurgiţeanu N. (2000). The FitzHugh–Nagumo Model: Bifurcation and Dynamics. Dordrecht: Kluwer Academic Publishers. Romero-Ortuño R. Martínez-Velilla N. Sutton R. Ungar A. Fedorowski A. Galvin R. (2021). Network Physiology in Aging and Frailty: The Grand Challenge of Physiological Reserve in Older Adults. Front. Netw. Physiol. 1, 712430. 10.3389/fnetp.2021.712430 Rosenblum M. G. Pikovsky A. S. Kurths J. (1997). From Phase to Lag Synchronization in Coupled Chaotic Oscillators. Phys. Rev. Lett. 78, 4193–4196. 10.1103/physrevlett.78.4193 Rydin Gorjão L. Saha A. Ansmann G. Feudal U. Lehnertz K. (2018). Complexity and Irreducibility of Dynamics on Networks of Networks. Chaos: Interdiscipilinary J. Nonlinear Sci. 28, 106306. 10.1063/1.5039483 Rydin Gorjão L. Heysel J. Lehnertz K. Tabar M. R. R. (2019). Analysis and Data-Driven Reconstruction of Bivariate Jump-Diffusion Processes. Phys. Rev. E 100, 062127. 10.1103/PhysRevE.100.062127 Rydin Gorjão L. Witthaut D. Lehnertz K. Lind P. G. (2021). Arbitrary-order Finite-Time Corrections for the Kramers–Moyal Operator. Entropy 23, 517. 10.3390/e23050517 Saha A. Feudel U. (2017). Extreme Events in FitzHugh-Nagumo Oscillators Coupled with Two Time Delays. Phys. Rev. E 95, 062219. 10.1103/PhysRevE.95.062219 Scheffer M. Bascompte J. Brock W. A. Brovkin V. Carpenter S. R. Dakos V. (2009). Early-warning Signals for Critical Transitions. Nature 461, 53–59. 10.1038/nature08227 Scheffer M. Bolhuis J. E. Borsboom D. Buchman T. G. Gijzel S. M. Goulson D. (2018). Quantifying Resilience of Humans and Other Animals. Proc. Natl. Acad. Sci. (U.S.A.) 115, 11883–11890. 10.1073/pnas.1810630115 Schoenmakers S. Feudel U. (2021). A Resilience Concept Based on System Functioning: A Dynamical Systems Perspective. Chaos: Interdiscip. J. Nonlinear Sci. 31, 053126. 10.1063/5.0042755 Stankovski T. Pereira T. McClintock P. V. E. Stefanovska A. (2017). Coupling Functions: Universal Insights into Dynamical Interaction Mechanisms. Rev. Mod. Phys. 89, 045001. 10.1103/RevModPhys.89.045001 Tabar M. R. R. (2019). Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems: Using the Methods of Stochastic Processes. Cham-Switzerland: Springer. 10.1007/978-3-030-18472-8: van der Bolt B. van Nes E. H. Scheffer M. (2021). No Warning for Slow Transitions. J. Roy. Soc. Interf. 18, 20200935. 10.1098/rsif.2020.0935 van der Pol B. van der Mark J. (1928). The Heartbeat Considered as a Relaxation Oscillation, and an Electrical Model of the Heart. Phil. Mag. 7 (6), 763. 10.1080/14786441108564652 Vanselow A. Wieczorek S. Feudel U. (2019). When Very Slow Is Too Fast-Collapse of a Predator-Prey System. J. Theor. Biol. 479, 64–72. 10.1016/j.jtbi.2019.07.008 Walker B. Holling C. S. Carpenter S. Kinzig A. (2004). Resilience, Adaptability and Transformability in Social–Ecological Systems. Ecol. Soc. 9, 5. 10.5751/es-00650-090205 Weinans E. Quax R. van Nes E. H. van de Leemput I. A. (2021). Evaluating the Performance of Multivariate Indicators of Resilience Loss. Sci. Rep. 11, 9148. 10.1038/s41598-021-87839-y Wilkat T. Rings T. Lehnertz K. (2019). No Evidence for Critical Slowing Down Prior to Human Epileptic Seizures. Chaos: Interdiscip. J. Nonlinear Sci. 29, 091104. 10.1063/1.5122759 Wunderling N. Krönke J. Wohlfarth V. Kohler J. Heitzig J. Staal A. (2021). Modelling Nonlinear Dynamics of Interacting Tipping Elements on Complex Networks: the PyCascades Package. Eur. Phys. J. Spec. Top. 230, 3163–3176. 10.1140/epjs/s11734-021-00155-4