A Novel Measure Inspired by Lyapunov Exponents for the Characterization of Dynamics in State-Transition Networks

The combination of network sciences, nonlinear dynamics and time series analysis provides novel insights and analogies between the different approaches to complex systems. By combining the considerations behind the Lyapunov exponent of dynamical systems and the average entropy of transition probabilities for Markov chains, we introduce a network measure for characterizing the dynamics on state-transition networks with special focus on differentiating between chaotic and cyclic modes. One important property of this Lyapunov measure consists of its non-monotonous dependence on the cylicity of the dynamics. Motivated by providing proper use cases for studying the new measure, we also lay out a method for mapping time series to state transition networks by phase space coarse graining. Using both discrete time and continuous time dynamical systems the Lyapunov measure extracted from the corresponding state-transition networks exhibits similar behavior to that of the Lyapunov exponent. In addition, it demonstrates a strong sensitivity to boundary crisis suggesting applicability in predicting the collapse of chaos.


Introduction
Complex network theory has had many interdisciplinary applications in different domains of social sciences, epidemiology, economy, neuroscience, biology etc. [1]. In recent years different network approaches have been also developed for nonlinear time series analysis. For a detailed review see [2]. Proper mapping between a discrete time series and a complex network in order to apply the tools of network theory in an efficient manner is not a trivial question. In case of continuous-time dynamical systems it can be even more complicated. There are several approaches to this problem, here we mention three large categories [2]: (1) Proximity networks are created based on the statistical or metric proximity of two time series segments. The most studied variant of proximity networks are recurrence networks [3]. These have found many applications, in the characterization of discrete [4] and continuous dynamical systems [5,6], in the classification of medical signals [7], and in the analysis of two-phase flows [8]. A special version are joint recurrence networks, which were developed for the detection of synchronization phenomena in coupled systems [9].
(2) Visibility graphs capture the convexity of subsequent observations [10]. The methods of natural and horizontal visibility graphs belong to this class. Visibility graphs were used for the analysis of geophysical time series [11], for the characterization of seismic activity [12], two-phase fluid flows [13], and for algorithmic detection of autism [14]. (3) State-transition networks (STN) represent the transition probabilities between discretized states of the dynamics. These can be threshold-based networks [15] or ordinal partition networks [16]. These also have found many applications in the domain of biological regulatory networks [17], study of signals of chaotic circuits [18], electrocardiography [19], economic models [20] and climate time series [21]. STNs are in fact an equivalent representation of discrete-time finite-state Markov chains with a time-homogeneous transition matrix [22].
The adjacency matrix of the STN is the transition matrix of the Markov chain, and it is a right stochastic matrix (a real square matrix, with each row summing up to one [23]), having many known properties [22], see Section 2.1 for details. One of the frequently used entropy-type measures for characterizing Markov chains is the Kolmogorov-Sinai entropy [24].
While these approaches have been extremely useful -as shown by the wide range of applications, people have always concentrated on applying graph theory tools to obtain a new understanding of the dynamics, employing several traditional network measures [1,2]. Here we would like to take an inverse step and generalize a prominent measure from nonlinear dynamics theory on STNs. The Lyapunov exponent is one of the most used quantities for analyzing dynamical systems [25]. However, its estimation for time series, when the dynamical equations are not available, is seldom trivial [26]. Here we introduce an analogous measure on STNs that can become an effective measure also in time series analysis. Inspired by the theory of dynamical systems we look at the trajectory length of a system evolving over discrete time according to the transition probabilities defining its STN. We show that by associating the appropriate length measure to the transitions the combined ensemble and time average of the length of a single step yields the Kolmogorov-Sinai entropy. The new quantity, we name it Lyapunov measure, is defined in analogy to the Lyapunov exponent by estimating the variance of trajectory lengths during a random walk over the network. We find that the Lyapunov measure is able to distinguish between periodic and chaotic time series, and detect furthermore crisis-type bifurcations by presenting pronounced peaks in the vicinity of these parameters.
After a short description of the STNs, we present the new Lyapunov measure and its properties. We test and compare its behaviour with the Kolmogorov-Sinai entropy on a theoretical network model with cyclic properties, the discrete-time Henon map [27], and the continuous-time Lorenz system [28].

State-Transition Networks
Mapping a dynamical system into an STN requires us to assign the different states of the dynamics to certain nodes of the network [15,16]. Directed edges of the network correspond to transitions between the discrete (or discretized) states of the system, characterized by the transition probability p ij from state i to j, where the probability of leaving the node equals to one (we have a right stochastic matrix), Trajectories consisting of several consecutive timesteps of the dynamics determine a path between distant nodes i and k of the STN, with the probability of selecting a particular path, e.g., i → j → k, given that the trajectory starts from node i is hence being given by the product of the respective transition probabilities. The transition probabilities, p ij introduced in Equation (1) are conditional by their definition assuming that the starting node of the transition is i. Therefore, the non-conditional probability of visiting edge i → j can be expressed by the Bayes-formula, where x i is the probability of residing in node i. Analogously to geometric distance, one may assign a length l ij as weights to the respective edges [29], which takes low values for high probabilities and vice-versa: The total length L ijk of such a i → j → k path is then given by the sum of lengths of the links along the path. A time series can hence be encoded by an STN as described above. Nodes of the network represent the spatial structure, while the time-like behavior is encoded by weighted and directed edges. Mathematically, the weighted adjacency matrix of transition probabilities, (P) ij = p ij , determines the time evolution of an ensemble of trajectories on the STN: where , and x i (t) denotes the probability of finding the system in state (node) i at time t. This process may also be seen as a time-homogeneous discrete-time Markov chain with a finite state space [24]. For STNs the stationary distribution is given by the renormalized eigenvector corresponding to the unit eigenvalue, Note that since the evolution operator P of the STNs considered here is a stochastic irreducible matrix [23], its largest eigenvalue is always 1 and the existence of the corresponding positive eigenvector is guaranteed by the Perron-Frobenius theorem [30]. For aperiodic transition matrices the long-term distribution is independent of the initial conditions, the system evolving over time to the stationary state x * . In case of periodic solutions however the long-time averaged distribution is also given by the stationary solution,

Lyapunov Measure
In the context of dynamical systems theory, Lyapunov exponents are the best known quantities used to characterize the system's behavior [25,26]. Here, we introduce an analogous quantity for STNs. Given a STN, one may define trajectories similarly to random walks on graphs: for any given initial state i the next state j is chosen randomly, using the transition probabilities p ij . As described above, the weighted length of such a trajectory segment is given by Equation (3). By concatenating t subsequent segments one may construct an ensemble of different trajectories each of length where l k is the length associated to the kth segment in the trajectory. For a large enough ensemble, in the limit of long times, both the average of the total path length L and its variance ∆L scales linearly with the total number of steps, t (see the Erdős-Rényi-type weighted random graphs in Figure 1A and the respective average lengths L and variance ∆L as a function of time in Figure 1B): As we shall see in the following the first observation allows for the measurement of the Kolomogorov-Sinai entropy [24] for STNs, while the second one enables us to define a novel network measure, similar in spirit to the Lyapunov exponents. While the above scaling behavior offers an intuitive picture of the properties of average path length and its variance, an alternative approach allows a straightforward mathematical treatment and facilitates computation to a great extent. In the limit of infinitely large ensemble of infinitely long walks the problem reduces to a basic Markov process. Using the definition of the total trajectory length (7), one can easily show that the asymptotic behavior, that is when t → ∞, of the ensemble averaged path length, grows linearly in time. The scaling factor, given by the average edge length of the STN, l , weighted by the non-conditional occurrence probability, q * ij = x * i p ij , of visiting the edge i → j in the asymptotic limit (compare Equation (2)), turns out to be equivalent with the Kolmogorov-Sinai entropy for Markov chains [24], measuring the average entropy of transition probabilities, with respect to the stationary distribution x * i given by Equation (5).

Figure 1.
Dependence of the Lyapunov measure defined in Equation (14) on the cyclicity parameter, c, for complete graphs of different sizes, N, edges associated with probabilities defined in Equation (15) and lengths expressed by Equation (3). (A) Examples of graphs with N = 10 nodes using different cyclicity parameters, c. The probabilities, p ij , associated with edges are represented by their width. (B) Dependence of the mean, L , and variance, ∆L, of trajectory length (see Equation (8)) on the number of steps, t, during simulated random walk on graphs similar to those in panel (A). The size of the networks is N = 50. Data points correspond to ensemble averages of 1000 independent trajectories. (C) The Lyapunov measure as a function of the cyclicity parameter, c, for graphs similar to those in panel (A). Noisy lines with markers represent simulation results of random walk on these graphs. The measure is calculated according to Equation (14) using ensemble averages of 1000 trajectories. Each trajectory includes 1000 steps. Starting nodes are chosen randomly with uniform probability and the first 500 "thermalization" steps are ignored when estimating trajectory lengths. Dashed lines represent simulation results on similar graphs by estimating the variance of the edge length according to Equation (13) with C = 0. Continuous lines represent the analytical model (see Equation (20)).
Similarly, due to the intrinsic diffusive nature of the envisioned process the variance of the trajectory length, ∆L, also grows linearly with the number of steps, t. The average squared path length, can be given in terms of the average lengths l and squared lengths l 2 , and of the asymptotic average of the correlation function C(k, k ) for timesteps k, k → ∞, The variance of the total path lengths is hence proportional to the variance of edge lengths for the whole network, σ 2 l = l 2 − l 2 , weighted according to transition frequencies during the stationary Markov process and the time averaged correlation function C .
Traditional Lyapunov exponents measure the exponential divergence rate of initially close-by trajectories during the time evolution of a dynamical system. For STNs, we can define a network measure bearing similarities with the Lyapunov exponents by considering the scaling of the length difference of pairs of random paths started from the same initial node of the network. As a further simplification, we define the Lyapunov measure for STNs using a proportional quantity, the variance of the total length of single trajectories, For uncorrelated processes, viz when C = 0, we can estimate the mean L and the variance ∆L of the trajectory lengths by the overall mean l and variance σ 2 l in the length "covered" in a single step, a methodology which will be first tested for random networks.

Properties of the Lyapunov Measure
The STN analog of limit cycles and strange attractors in dynamical systems would be circular networks and high degree networks, respectively. Let us investigate the dependence of a STN's Lyapunov measure, Λ, on the degree of cyclicity of the network. We shall build our description on a simple theoretical model consisting in a complete directed graph of N nodes identified by their label i ∈ 0, 1, . . . , N − 1. In order to study the effect of cyclicity the transition probabilities are associated to the edges in three steps: first, to each of the N(N − 1) links we assign a random u ij value distributed uniformly between zero and one; second, we rescale these probabilities such that each node has a privileged target and source node with transition probabilities: where "mod" denotes the modulo operation and c ∈ [0, 1] is the parameter quantifying the cyclicity of the graph (see Figure 1A). At c = 1 we have a purely cyclic graph where each node is only linked to the next node in the list while for c = 0 we retain the original u ij weights. As a last step, we normalize the outgoing weights according to Equation (1).
Computing the average and the variance for an ensemble of trajectory lengths, started from random initial nodes, we obviously recover the theoretically calculated linear scaling behavior (compare Equations (9) and (13)).
The results of the simulation implementing Equation (14) are represented by the noisy lines with markers in Figure 1C. There is an apparent dependence on the size of the network. However, a complete graph with all links equivalent (c = 0) with weights distributed uniformly yields Λ = 1/4 irrespective of the size of the network. A simple cycle with all nodes having a single incoming and outgoing link (c = 1) eliminates all randomness and hence the variance based Lyapunov measure goes to zero. For intermediate values of the cyclicity parameter, c, the Lyapunov-measure evolves smoothly exhibiting a maximum somewhere in the upper half of the unit interval.

Analytical Study of the Lyapunov Measure
For a conceptual grasp of the behavior shown in Figure 1C we further simplify our model. In a fully connected graph of size N all l i -s in Equation (7) can be regarded as identical random variables, hence where σ 2 l ≡ l 2 − l 2 is the variance of the variables. We can work with a single vertex and extrapolate the results to all vertices. One neighbor of this "generic" vertex is privileged as it is the next in the cycle (see Figure 1A). Let us assume that the rest of the N − 2 neighbors are equivalent (their links have equal weights): where n = N − 2, and with the cyclicity parameter c ∈ [0, 1] the probability p = 1 n+1 + n n+1 c will change between 1/(n + 1) (all links equivalent) and one (fully cyclic). The length of link i is set to − ln p i , the probability of choosing that link is p i : The 1/4 offset at c = 0 is a consequence of the fixed transition probability values in Equation (17). Appendices A and C account for this offset. A detailed investigation of the maximum is under Appendix B.

Time Series Analysis with STNs
The random network model in the previous section has the benefit of offering a basic understanding of the Lyapunov measure. In order to test the scope of its applicability we apply it to the well-known Henon map and the Lorenz system and compute Λ for the STNs constructed for time series with different control parameters. We show that the network measure is able not only to distinguish between periodic and chaotic regimes but also presents pronounced peaks before crisis-type bifurcations. To that end we need to map "real-world" time series to STNs. In this section we introduce a methodology which is generic enough to be applied both on discrete-and continuous-time dynamical systems.

Construction of STNs
To construct the STN corresponding to the dynamics of the time series generated by a dynamical system one needs to discretize space and time adequately. In case of discrete-time systems time is inherently defined in terms of integer-valued timesteps. For continuous-time systems this step is however not as straightforward. A time sampling with a constant sampling rate may not be beneficial for all cases since the trajectory may evolve with various speeds in the phase space, having segments where it slows down while at other parts advancing abruptly orders of magnitude further within the same time interval. A too high sampling frequency would lead to a very high number of self-loops in the STN, since for any finite spatial mesh there exist several consecutive time points which fall within the same grid cell. On the other hand, a very low sampling rate leads to the loss of time correlations. Here we propose to use the well-defined method of Poincare sections, which naturally generates the equivalent map, i.e., the discrete-time version of the dynamical system. The Poincare map can be constructed by tracking the consecutive unidirectional intersection points of the trajectory and the Poincare surface, a 2D surface in case of 3D systems.
On the other hand, for an STN one needs to assign the nodes of the network to different states during the dynamics. One may hence discretize the phase space by dividing it into equal sized bins, corresponding to the nodes of the network. Here we consider the effective phase space of the system, viz the largest rectangle in the phase plane of the map/Poincare map which fully covers the attractor, and construct a mesh, each of the newly created bins defining a coarse-grained state of the original dynamical system. Directed edges of the network correspond then to transitions between the discrete states of the phase space, characterized by the transition probabilities p ij .
To illustrate the methodology we present examples of different attractors with the corresponding STNs for the discrete-time Henon map [27] and the continuous-time Lorenz system [28] (see Figures 2 and 3). The two systems are briefly surveyed in Section 3.1. For the Henon-map we use n max = 10 4 steps long trajectories started from randomly chosen initial conditions while discarding the first n trans = 10 3 transient steps to let the dynamics settle onto the attractor of the system. As a function of the a parameter the attractor may have either a fractal structure of folded filaments or it may consist of discrete points separated in the phase plane (see the top panels of Figure 2). For the Lorenz system we used the method of Poincare sections, with the Poincare surface being a predefined 2D plane. As an example, in Figure 3 we illustrate the projection of the x = 15 Poincare plane by a dashed line in the (z, x) phase plane (top row). The corresponding Poincare map is constructed using t max = 5000 time unit long trajectories, discarding the initial t trans = 300 time units. The choice of the Poincare section does not effect the results qualitatively. For comparison, a second, perpendicular Poincare section, defined by y = 0, have also been tested (not shown here). We found that the results are robust with respect to the choice of the Poincare plane. The topology of the resulting Poincare sections also reflect faithfully the structure of the attractors (middle row).
For both the Henon map and Lorenz system we selected parameters corresponding to qualitatively different dynamical regimes: traditional chaotic attractors with extended fractal filaments (orange and red colors), partially predictable chaotic attractors with small patches of Cantor-sets [31] (green color), and individual points for periodic attractors (blue).
The STNs are constructed using a 20 × 20 mesh in the phase planes of the respective maps. For computing the transition probabilities p ij we generate long discrete-time trajectories as described above and count the number of i → j transitions along these trajectories. The p ij probabilities are finally given by normalizing the number of jumps between consecutive states according to Equation (1). Chaotic dynamics generates densely wired graphs (left columns of Figures 2 and 3) and widely distributed global transition frequencies ∝ q * ij (as illustrated by the width of the network connections), in contrast to periodic motion for which the network collapses to a simple cyclic chain of edges and vertices (right column of Figures 2 and 3). Though, seemingly two of the selected chaotic attractors are identical (see the orange and red attractors in Figure 2 and 3), looking more closely one realises that for the red colored STNs the global transition probabilities (the widths of the edges) are more heterogeneous than for the orange colored STNs, allowing for a more probable cyclic path in the network. Note that the random network model introduced in Section 2.2.1 resembles this network structure for large (but smaller than 1) cyclicity parameters, e.g., when c ∼ 0.75.

Network Measures
To characterize the network with a single scalar parameter we compare the two STN measures discussed in Section 2.2. First, we adopt here the well-known Kolmogorov-Sinai (KS) entropy for Markov chains, as defined by Equation (9). Second, we compare the above introduced novel measure (14) similar in idea to the well-known Lyapunov exponents. As a basis of comparison we also provide the bifurcation diagrams of the Henon and Lorenz systems, together with the classical largest Lyapunov exponent λ m (see Figure 4).
As changing the control parameters a and ρ, for the Henon map (21) and for the Lorenz system (22) respectively, both exhibit a whole series of complex bifurcations (see the top panels of Figure 4). These bifurcation scenarios are, as expected, faithfully reflected by the largest Lyapunov exponent values: λ m < 0 and λ m = 0 for periodic behavior of the Henon map and Lorenz equations, respectively, while λ m > 0 corresponding to chaotic dynamics in both cases.
The Kolmogorov-Sinai entropy, S KS , is defined for the STNs obtained from the dynamical systems for the actual control parameter values. Here we implement both the original definition based on the global transition probability matrix q * ij and the algorithm based on generating random trajectories in the STN and measuring their linear scaling factor with time (compare Equations (9) and (10)). As expected the exact matching of the two definitions is reflected by the numerical results as well (see Figure 4): S KS = 0 denoting exactly cyclic networks, low positive values of S KS correspond to partially predictable chaotic motion, while chaotic dynamics generates high entropy networks.
The Lyapunov measure, Λ, is computed here based on the estimation of the variance in trajectory lengths for large ensembles of random paths, according to Equation (14) (see Section 3.3 for details on ensemble statistics). In the examined parameter region, while bearing many similarities with Kolmogorov-Sinai entropy, the Lyapunov measure exhibits interesting behavior in the vicinity of crisis-type bifurcations points (abrupt disappearance or reduction in size of the attractor, viz in the size of the black shaded regions in the bifurcation diagram) [32]. Approaching the crisis from the direction of the chaotic region Λ presents an abrupt peak, forecasting in some sense the collapse of the chaotic attractor (see the red colored sharp peaks in the bottom panels of Figure 4).
As it is demonstrated using the random network model in Section 2.2.1, the Lyapunov measure has a maximum point with respect to the cyclicity parameter (see also Appendix B). The appearance of the peaks is closely related to this special cyclic topology in real systems as well (as shown in Appendix D). The height of the peaks is furthermore boosted by the correlations of edge lengths along the paths of random walks (see Equation (14)).
Interestingly these precursor peaks seem to be more pronounced for boundary crisis than for the case of interior crisis (compare for example the red-colored peak of the Lorenz system with the peak around the parameter value denoted by the red dot).
Our results, summarized in Figure 4, demonstrate that STNs can encode all the relevant information about the dynamics of the system. For chaotic dynamics with positive maximal Lyapunov exponent λ m > 0 the network measure Λ is also positive. For periodic motion both quantities are zero, λ m = Λ = 0. Interestingly, while the Lyapunov exponent decreases when approaching the bifurcation point where chaos disappears, the here introduced Lyapunov measure shows a diverging tendency.  Figures 2 and 3. To construct the STNs the Henon-map is iterated for 3 × 10 4 timesteps omitting the initial 1000 transient steps, respectively for the Lorenz system 5000 timeunits were used ignoring the first 500 units long transients. The Kolmogorov-Sinai entropy is computed both using the adjacency matrix of global transition probabilities q * ij (compare Equations (9) and (10)) and by generating and ensemble of random paths using the same number of steps as for the Lyapunov measure. The Lyapunov network measure, defined by Equation (14), is calculated over an ensemble of 100 random trajectories of 10 4 steps, neglecting the initial 5000 steps (see Section 3.3). As an example of boundary-crisis precursor, a single peak in the Lyapunov measure is colored in red for both systems.

Discrete-and Continuous-Time Dynamical Systems
The Henon-map is considered as a prototype system for studying bifurcations and chaos in discrete-time dynamical systems [32]. The mapping from state (x n , y n ) to (x n+1 , y n+1 ) is defined by The effect of the nonlinear term x 2 n in the dynamics may be tuned by changing parameter a, hence in most of the studies it is considered the control parameter, while choosing b = 0.3 standard value for the other parameter. As a function of the a parameter the attractor may have either a fractal structure of folded filaments or it may consist of discrete points separated in the phase plane (see the top panels of Figure 2).
The Lorenz system is probably the most known continuous-time dynamical systems exhibiting chaos on a butterfly-shaped attractor. Being a 3D system defined aṡ it is one of the simplest examples which allows for the presence of chaotic behavior. While we consider the standard σ = 10 and β = 8/3 parameters, we select with ρ ∈ [180, 182] a parameter interval for which one finds not only chaotic and periodic behaviors, but also partially predictable chaos (PPC), as introduced in [31]. The respective attractors are illustrated in Figure 3.

Phase-Space Discretization and Poincare Sections
The binning of the phase space is a cornerstone in the construction of STN networks. For a low number of bins one loses information regarding the dynamical states of the system. On the other hand for a too high resolution the STN collapses to chain-like networks even for an otherwise complex dynamics, requiring exponentially long time series to overcome the statistical unreliability of the data. For chaotic time series, the number of nodes of the obtained STN increases with the binning resolution. Yet the measure proposed here turns out to be reliable since it is a monotonous function of the number of nodes (compare Equation (20) and Figure 1). The choice of the Poincare section also does not affect the results qualitatively. For comparison, two perpendicular Poincare sections have been tested, defined by x = 15 and y = 0, respectively. We found that the results are robust with respect to the choice of the Poincare plane used for the construction of the STN.

Ensemble Averages and Asymptotic Behavior
Given a STN, the Lyapunov measure (14) is computed numerically by constructing random trajectories of total length L i (t) at time t, and measuring their variance over an ensemble of n paths, i = 1, . . . , n. For a large enough ensemble the ∆L/t converges to a steady state value, termed here as the Lyapunov network measure and denoted by Λ. For the STNs obtained from the Lorenz-system generated time series for ρ = 180.7 the fluctuations become relatively small for an ensemble average computed over n > 10 3 random trajectories and for t > 5000 time units (see the left panel of Figure 5). In order to reduce the fluctuations one may further average over the stationary ∆L/t values in time. Different dynamical behaviors correspond to separated plateaus in the time series of ∆L/t, leading to well-defined and different Lyapunov-measure values for four parameter settings presented in Figure 3 (compare the right panel of Figure 5).

Discussion
Complex dynamics in continuous phase space and time bears similarities and differences compared to those occurring on networks. Inspired by the utility of the Lyapunov exponent in dynamical systems theory we introduced a somewhat analogous network measure for STNs based on the Kolmogorov-Sinai entropy for graphs. The random STN model shows that using the variance instead of the mean length, equivalent with the Kolmogorov-Sinai entropy, contributes to a surge in the measure as we approach cyclicity (see Figure 1C). Our analytical study explains the important properties of the new measure as they manifest in the case of random networks. In order to assess the connection between the new measure and the classical Lyapunov exponent we also introduced a novel procedure for converting time series to STNs. Our examples include both discrete-time (Henon-map) and continuous-time (Lorenz) systems. The Kolmogorov-Sinai entropy of the so obtained Henon/Lorenz STNs reproduces the control-parameter dependence of the Lyapunov exponent of the corresponding dynamical systems. The newly introduced Lyapunov measure, however, exhibits an additional and unexpected property. It appears to be sensitive to boundary crisis where a chaotic attractor is suddenly created or destroyed (see bottom panels of Figure 4). This extreme jump in the measure during crisis is due to the cumulation of the basic increase experienced also with random networks (see Figure 1C) and the effect of correlations due to the heterogeneity of the Henon/Lorenz STNs. In present paper the Lyapunov measure for these networks was estimated based on trajectory simulations. However, accounting for the correlations as introduced in Equation (12) should be possible based solely on the transition probabilities, p ij . In this sense, similarly to the case of the Lyapunov exponent, it may allow two slightly different but -for the case of ergodic systems-equivalent interpretations: one may either follow trajectories and compute the quantities of interest along them or determine the natural/stationary distributions underlying the dynamics and calculate spatial averages according to these probabilities [32]. The pronounced peaks during boundary crisis suggests possible applications for forecasting transitions from chaos to periodicity.
The proposed method can be applied to any stationary time series and any dynamical system with attractors. It might not always provide as relevant output as for the two examples given in the paper. We expect, however, similar results for most of the commonly studied systems. In case of multistable dynamical systems the Lyapunov measure has to be computed separately for each attractor, similarly to the traditional Lyapunov exponent. For non-stationary time series the proposed methods and quantities may need further developments.

Appendix D. Insight into the Emergence of Precursor Peaks
The surge in the Lyapunov measure as we approach the collapse of chaos (see the region in red in the bottom panels of Figure 4 is a cumulative result of three factors: (i) the random network model in Section 2.2.1 shows that as cyclicity grows the measure attains a maximum (see Figure 1C); (ii) below we provide evidence for a superlinear growth in the cyclicity itself close to the bifurcation point; (iii) nontrivial correlations (see Equation (12)) present in real STNs can further boost the effect. Within the presented conceptual framework cyclicity can be quantitatively defined only through the random network model of Section 2.2.1. Its scope can be extended to general networks through Equations (9) and (18). Following Equation (9) we can compute the Kolmogorov-Sinai entropy, ergo, the mean length l of an edge for any STN. By inverting the relationship l (c) from Equation (18) we can extract information on the cyclicity, c, of these networks. Subsequently, the Lyapunov measure can be estimated from Equation (20). There is an inherent mismatch between the two approaches to the mean lengths formalized in Equations (10) and (18). As opposed to a general STN the analytical model assumes a complete graph making the degree of the nodes and the size of the network coincide. This discrepancy can be to some extent bridged by using the average node degree as n when inverting the relationship in Equation (18). The results depicted in Figure A2 are the counterparts of those shown in the bottom right panel of Figure 4.  (9). The cyclicity parameter is estimated by inverting numerically the l (c) relationship in Equation (18). As a last step, the Lyapunov measure was approximated using the analytical expression from Equation (20). The random walk simulation counterpart of the image is presented in the bottom right panel of Figure 4. For more details, refer to Appendix D.