Towards a Data-Driven Estimation of Resilience in Networked Dynamical Systems: Designing a Versatile Testbed

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.


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 resiliencedynamical resistance  takes into account that the dynamics of some elementary unit of a Ndimensional networked dynamical system can be described by its self-dynamics as well as by interactions with other units: 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 selfdynamics 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.

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.

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: 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 selfinduced 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: where C (k) w and C (k,l) b 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).

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): ( 4 ) Here a G n , μ G n , and σ G n 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: where 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 (k) w 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.

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 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). 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. 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 datadriven indicators of resilience.

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 nonperturbative approach to efficiently monitor dynamical FIGURE 3 | (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.
FIGURE 4 | 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.
Frontiers in Network Physiology | www.frontiersin.org 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 timedependent 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 timeseries-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;. Here we use three widelyused techniques, namely the zero-and maximum-lag crosscorrelation (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 timeresolved 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.
Frontiers in Network Physiology | www.frontiersin.org March 2022 | Volume 2 | Article 838142 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 10 6 time steps thereby starting from identical initial conditions of each oscillator. For our analyses, we use the oscillators' xcomponents 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). FIGURE 6 | 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.
FIGURE 7 | 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. 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 datadriven 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 10 6 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 . 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 timeseries-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 selfdynamics 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 noiserelated 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-seriesanalysis-based investigation in networked dynamical system may, however, require more refined and better adapted analysis techniques (Rydin Gorjão et al., 2019Aslim 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.