Our method is based on the assumption that the information structure of a system captures every effective interaction among its agents and thereby reflects the condition of the system. The abstract nature of information theory allows us to analyse systems independently of their specific instantiation. We only rely on the availability of longitudinal data reflecting the strength of the system’s individual components in a very broad sense. Hence, in general, our method can be applied to any complex system. The only condition is that for a given period of time and for every component of the system, a time series of quantitative data reflecting the outcome of interactions exists. Such time series could exemplarily be biomass of a plant species, number of individuals of an animal species, or sales of a company. The data type can differ among the components of the system, i.e. be heterogeneous.
In a first step, networks of information transfer are inferred via pairwise estimation of transfer entropy9 among all agents. Considering these networks and, in particular, their development over time, offers insights into functional interactions.
In the second step, potential, connectedness and resilience are computed solely using the networks of information transfer (see Supplementary A for a review of the adaptive cycle and its defining variables). Here, we utilize capacity and ascendency as being defined by Ulanowicz in the context of ascendency theory11. Note that Ulanowicz also used information theory to define capacity as an entropy of flows and ascendency as mutual information between inflow and outflow. While the first one is a measure of the average indeterminacy in the fluxes of the network, the latter quantifies the efficiency the system has in making use of its capacity12. However, being rooted in systems ecology, Ulanowicz always considered flows of physical quantities, such as energy or resource fluxes. In contrast, we will derive the quantities from networks of information transfer, abstracting from the physical representation of the interaction. Thus, potential is the capacity of the network of information transfer, and connectedness the corresponding ascendency.
The challenging part of our approach is to find an appropriate measure of resilience. There exist various conceptions and following definitions of resilience13. For our purposes, Holling’s definition of resilience fits best, namely to define resilience as ”the magnitude of disturbance that can be absorbed before the system changes the variables and processes that control behavior” (Ref.1, p. 28). There have been various approaches to make this notion measurable, however, all of them either depending on the specific system under observation14,15,16,17 or requiring deep knowledge of the system dynamics18. Resilience has also been studied from a network perspective (see e.g.19,20). Since we are modeling complex systems as networks of information transfer, our definition is inspired by a common concept in spectral graph theory. We use the so-called graph Laplace operator, which captures vulnerability of a network with respect to perturbations of its topological structure.
Taken together, the development of these three variables displays the system’s course through the adaptive cycle, helping to better understand system maturation and to evaluate its current condition. We will now provide a detailed description of our method, its implementation, and its application in the three case studies presented in this paper.
Step 1: estimation of networks of information transfer
Let (mathscr {V}) be a collection of variables, quantifying the state of agents defining a system. Let (I = (i_1, dots , i_N)) and (J = (j_1, dots , j_N)) be two sets of samples of states for the components I and J, say. For example, I and J can be identified with abundances of two interacting species at time points (1, dots , N). We consider the time series I and J as realisations of two approximately stationary discrete Markov processes. This allows us to compute Schreiber’s transfer entropy9, serving as a measure of their effective interaction. Transfer entropy from J to I is defined as
$$begin{aligned} T_{J rightarrow I} = sum _{n = 1}^{N-1} p left( i_{n+1},i_{n}, j_{n} right) cdot log left( frac{p left( i_{n+1}|i_{n}, j_{n} right) }{p left( i_{n+1} | i_{n} right) } right) . end{aligned}$$
(T_{J rightarrow I}) quantifies the average reduction in uncertainty about the future of I given the past of J. In other words, how much additional information do we gain about the next state of I, knowing not only the past of I itself, but the past of J as well. In the literature, a multitude of studies on the interpretation of transfer entropy in general and in specific contexts can be found21,22,23.
As the probabilities occurring in the definition of transfer entropy are in general not known, we have to estimate transfer entropy on the basis of given realizations of the random variables, i.e. the data given as samples of the time series. Typically, we do not use all available samples to estimate transfer entropy at time t but samples falling within a certain window of time preceding time t. The size (w_t) of this windows can either be fixed, or depend on the time t, e.g. (w_t = t). In the first case, the window is “shifted” going along with t to guarantee transfer entropy always being estimated on the same number of samples. In the second case, the window starts at the beginning of the time series and is extended with increasing t. In this case, the full history of the time series is considered for estimating transfer entropy. The choice of the window size depends on the system under consideration. In any case, it should be at least as large as the assumed order of the underlying Markov process. We then compute the information transfer from J to I at time t estimating transfer entropy over the period (t-w_t+1,dots ,t). To be precise,
$$begin{aligned} T_{J rightarrow I}^t = sum _{n = t-w_t+1}^{t} p left( i_{n+1},i_{n}, j_{n} right) cdot log left( frac{p left( i_{n+1}|i_{n}, j_{n} right) }{p left( i_{n+1} | i_{n} right) } right) . end{aligned}$$
Depending on the size of (w_t) and the data being available, it can be useful to increase the number of data points falling within every window by interpolation. For our calculations, we used the Matlab function pchip. Interpolation stabilizes the estimation in case of small window sizes. At the same time, interpolating too many points reduces stochasticity in the time series due to the deterministic component being introduced by the interpolation model. Thus, there is a trade-off between stochasticity and stability which has to be taken into account.
We estimated (T_{J rightarrow I}^t) using the Kraskov-Stögbauer-Grassberger (KSG) estimator TransferEntropyCalculatorKraskov as being provided with the JIDT toolkit24. For all our calculations, the function has been called using the data ((j_{t-w_t+1},dots ,j_t)) and ((i_{t-w_t+1},dots ,i_t)) in the mode computeAverageLocalOfObservations. For all other parameters of the estimation procedure, we used the default values (k=k_{tau}=l=l_{tau}=delay=1). Other choices of these parameters can be reasonable depending on the specific system to be analysed. To distinguish actual interactions from random noise, we tested all estimates via hundred-fold bootstrapping using the function computeSignificance(100) incorporated in the toolkit. Tests passing below a given significance level have been accepted and thus lead to an edge between the corresponding components with the estimated transfer entropy defining the corresponding weight. Estimating and testing for all pairs of components at fixed time t, we finally obtained a weighted, directed graph
$$begin{aligned} G^t = left( mathscr {V},{T_{J rightarrow I}^t|(J,I) in mathscr {V} times mathscr {V} } right) end{aligned}$$
as being our inferred model of interaction at time t. Given time series of abundances of length N for each component, this results in a sequence of interaction networks for time points (w_1,dots ,N).
Summarizing, the first step infers models of interaction among the given variables in form of a series of networks capturing the interaction patterns and their strengths. These network models can then be used in the second step to actually determine the position of the system within the adaptive cycle.
Step 2: determining potential, connectedness, and resilience
As mentioned before, our definitions of potential and connectedness are based on Ulanowicz’s notions of capacity and ascendency11. Ulanowicz provides further information on the theoretical background of these measures. Let
$$begin{aligned} T^t = sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t end{aligned}$$
be the total transfer of the system at time t. We further introduce the following shorthand notation
$$begin{aligned} T_{J}^{text {out},t} = sum _{I in mathscr {V}} T_{J rightarrow I}^t qquad text{ and } qquad T_{I}^{text {in},t} = sum _{J in mathscr {V}} T_{J rightarrow I}^t. end{aligned}$$
Define
$$begin{aligned} P^t = – sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t cdot log left( frac{T_{J rightarrow I}^t}{T^t} right) qquad hbox {as the system’s} ,potential ,hbox {at time} ,t end{aligned}$$
and
$$begin{aligned} C^t = sum _{(J,I) in mathscr {V} times mathscr {V}} T_{J rightarrow I}^t cdot log left( frac{T_{J rightarrow I}^tT^t}{T_{J}^{text {out,t}}T_{I}^{text {in,t}}} right) qquad ,hbox {as its} ,connectedness ,hbox {at time}, t. end{aligned}$$
Being essentially a sum over the indeterminacy in each transfer within the system, potential can be interpreted as a measure of the system’s power for evolution and its ability to develop. Recall that development of the system as a whole necessarily relies on communication, i.e., transfer of information among its components. In contrast, connectedness measures the degree of internal coherence of the system by contrasting information leaving one component with information arriving at another component.
In order to define resilience, we need to capture vulnerability of the system with respect to unforeseen perturbation. In terms of graph theory, this can be achieved by studying the eigenvalues of a certain matrix, being associated with the graph. Indeed, the smallest non-trivial eigenvalue of the so-called graph Laplacian of an undirected graph quantifies the vulnerability of the graph with respect to disturbance of the topology of the graph25,26. In our case, we need to transfer this idea to the case of the directed graphs (G^t).
Thus, given (G^t=left( mathscr {V},T^t_{J rightarrow I}|(J,I) in mathscr {V} times mathscr {V}right)) be a non-empty, weighted, directed graph with vertex set (mathscr {V}). Let further (c > 0) be a constant. Let (D_{out}) and (D_{in}) be the diagonal matrix of out-degrees and in-degrees, respectively, and A the weighted adjacency matrix. We then define the following Laplace type operators of (G^t):
$$begin{aligned} L_{out} = c cdot D^{-frac{1}{2}}_{out} left( D_{out}- A right) D^{-frac{1}{2}}_{out}, quad hbox { and }quad L_{in} = c cdot D^{-frac{1}{2}}_{in} left( D_{in} – A right) D^{-frac{1}{2}}_{in}, end{aligned}$$
following the convention that (D^{-frac{1}{2}}_*(u,u) = 0) for (D_*(u,u) = 0). Note that, for the sake of readability, we omitted the superscript t in these definitions. For all case studies presented in this paper, we used
$$begin{aligned} c = frac{1}{max { T^t_{J rightarrow I}|(J,I) in mathscr {V} times mathscr {V} }} end{aligned}$$
as standardization constant.
Since A is no longer symmetric, the spectrum of (L_{out}) and (L_{in}) is complex in general. Nevertheless, the distance of the spectrum to the imaginary axis in the complex plane still determines the stability of the graph. Therefore, we define resilience of the graph G as the smallest, non-trivial absolute value of the real parts of all eigenvalues of its two Laplacian matrices, i.e.
$$begin{aligned} R^t = min left{ |mathfrak {R}sigma | :sigma in {{,mathrm{Spec},}}(L_{out}) cup {{,mathrm{Spec},}}(L_{in}), sigma ne 0right} . end{aligned}$$
See Supplementary E for a more detailed explanation motivating this definition as well as for an alternative definition of the involved Laplacian matrices.
Our definitions of the three systemic variables are summarized in Table 1. In addition, Table 2 displays basic information concerning the data sets and parameters of our case studies. Note that, for visualization purposes, Figs. 3, 4, and 5 show a smoothed version of the estimated variables as being obtained by applying the R functions smooth.spline and splinefun.
Schematic representation of our quantification method. In the first step, time series of abundance data (a) are transferred into networks of information transfer (b). In the second step, the three systemic variables (c) are computed on basis of the networks. The figure depicts the window shifting method.
We created the R package QtAC in order to enable a straightforward application of our method27. The package comprises all functions required to compute a system’s course through the adaptive cycle and to visualize the results.
Figure 2 illustrates the key idea of our approach. Figure 2a shows randomly generated abundances of five components (A,B,C,D,E). To estimate the position of this small system within the adaptive cycle at time t and (t+1), we estimate transfer entropy for all pairs of components based on the samples within the window ((t-w+1, dots , t)) and ((t-w+2, dots , t+1)), respectively, and test for significance. This results in two inferred interaction networks shown in Fig. 2b. Using these networks, we can compute potential, connectedness and resilience at these two points in time. Figure 2c depicts the shift the system has made in the coordinate system spanned by the three characteristic variables.
The decrease in resilience from t to (t+1) mainly follows from loosing the edge (Drightarrow C) at (t+1). With component D being connected with the rest of the system by one edge, only, the system becomes more vulnerable, since perturbation of the edge (Drightarrow E) would fully decouple D from the rest of the system. Similarly, the loss of this edge also leads to a decrease of potential. Heuristically, the more edges a system has, the more potential there is to change from one state to another. Note that the even distribution of weights also added to the system’s potential, as for example edges (Arightarrow E) and (Arightarrow C) both loose weight. The moderate decline of connectedness follows from the loss of the edge (Drightarrow C) as well as from the smaller capacity of the edges (Arightarrow C) and (Arightarrow E), decreasing the overall total edge weight.
Source: Ecology - nature.com
