Method for Measuring Resilience of Complex Systems Using Network Theory and Graph Energy

Because of societies’ dependence on systems that have increasing interconnectedness, governments and industries have an increasing interest in managing the resilience of these systems and the risks associated with their disruption or failure. The identification and localization of tipping points in complex systems is essential in predicting system collapse but exceedingly difficult to estimate. At critical tipping-point thresholds, systems may transition from stable to unstable and potentially collapse. One of the approaches to measuring a complex system's resilience to collapse has been to model the system as a network, reduce the network behavior to a simpler model, and measure the resulting model's stability. In particular, Gao et al. introduced a method in 2016 that includes a resilience index that measures precariousness, the distance to tipping points. However, those mathematical reductions can cause the model to lose information on the topological complexity of the system. Using computational experimentation, a new method has been formulated that more-accurately predicts the resilience and location of tipping points in networked systems by integrating Gao et al.’s method with a measurement of a system's topological complexity using graph energy, which arose from molecular orbital theory. Herein, this method for measuring and managing system resilience is outlined with case studies involving ecosystem collapse, supply-chain sustainability, and disruptive technology. The precariousness of these example systems is found using the dimension reduction, and a graph-energy correction is quantified to fine-tune the measurement. Lastly, the integration of this method into systems engineering processes is explored to provide a measurement of precariousness and give insight into how a complex system's topology affects the location of its tipping points.


I. INTRODUCTION
So much depends on the nature of tipping points in critical systems that must not fail.From power grids to food supply chains, from socio-economic to ecological systems, modern societies rely on the continuous steady operation of the systems that support their needs.When a tipping point is reached, the system becomes unstable and spirals into a different behavior pattern.Severe tipping points can cause collapse.For example, in ecologies, a tipping point can be triggered by or cause extinctions that lead to the complete collapse of the ecosystem.In supply chains, a tipping point can cause a sustainable supply chain to become overloaded and decimate supplies.
To actively monitor the health of these systems, their resilience can be measured.The challenge of accurately modeling the resilience of networked dynamical microscale to macroscale systems is common to economics, biology, ecology, chemistry, physics, and engineering [1], [2], [3] and is essential in predicting risk [4].Resilience can be measured in several ways, including the time for the system to recover, the  decrease in performance during a disturbance, and the recovered performance after a disturbance [5].When it comes to resilience with respect to tipping points, one resilience metric is the precariousness of the system, which is its distance to the nearest tipping point [6].Fig. 1 illustrates this concept.When a system is perturbed away from its nominal operating zone and moves closer to a tipping point, its resilience decreases because of its precariousness.If precariousness could be measured, then it could provide system managers, designers, and stakeholders with information on how to avoid approaching undesired tipping points to sustain critical systems.
Until recently, estimating resilience relied on lowdimensional models [7], [8] that did not adequately capture the complexity of the system [9].In 2016, a revolutionary resilience metric was introduced by Gao et al. through a paper published in Nature [9].This first-of-a-kind tool uses a new mathematical dimension, named the effective plane, to show how close complex systems are to potential collapse.Fig. 2, redrawn from Gao et al. [9], shows a tipping point model for ecosystems using their method.A healthy ecosystem begins at a point along line x H , with a high number of species captured in a state parameter x eff .As the modeled ecosystem is weakened, Gao et al.'s resilience index, β eff , drops and eventually reaches a critical threshold, β c eff .Below this threshold, the ecosystem can become unstable and drop to the x L line, which means that the ecosystem completely collapsed to zero population.Because of this relationship, the difference between β eff and β c eff can be used as a measurement of precariousness.When the resilience index β eff reaches β c eff , the system can trigger a tipping point.
One fascinating aspect of the effective plane that was introduced by Gao et al. is that the behavior of large complex networks is reduced to that of a simpler system that has only one or several tipping points.The resulting model enables a measurement of precariousness for large networks, but it also loses some of the structure of the system and the complex dynamics in that topology.We theorize that there is an impact from this loss of topological complexity in this model reduction, and if we can measure and adjust for that topological complexity, then we could improve the prediction of tipping points in large complex systems.
To begin, we theorize that when large-scale systems are modeled as networks, their structure can be similar to that of molecular systems, and their levels of structural complexity can be measured using adjacent methods from molecular science [10].For example, researchers in complexity theory and systems engineering [11], [12] have correlated the topological complexity of systems with graph energy (E), which was originally a topological measurement for molecules [13], [14], [15], [16].The origin of graph energy was a measurement of binding energy in Hückel molecular orbital (HMO) theory, which quantifies total π -electron energy [17] and was shown to correlate with thermodynamic properties and resonance stabilization of organic molecules [18], [19].Furthermore, graph energy has been shown to correlate with tipping-point behaviors in chemistry and biology [20], including quantitative structure-activity relationship (QSAR) models [21], [22], entropy [23], [24], and properties of proteins [25], [26].Hypothesizing that graph energy could provide similar insight into the impact of system complexity on tipping points in large-scale systems, we adapt molecular orbital theory in the context of graph energy to estimate a system's topological complexity.
Ultimately, we seek insight into the larger fundamental question: does the complexity of a system affect its resilience to transition or collapse?This question is explored in this research by first discussing the previous work that this research builds on, in measuring and managing system complexity, tipping points, and resilience.Then the equations for graph energy are adapted to become a measurement for system complexity and interface with the resilience index from Gao et al. that measures a system's precariousness to tipping points.The combination of these methods is illustrated in Fig. 3.These methods are then tested in three case studies of systems undergoing tipping points to assess the impact of the system complexity on the system's resilience to triggering those tipping points.

II. LITERATURE REVIEW A. MEASURING AND MANAGING SYSTEM COMPLEXITY
Complex nonlinear systems make up many of the systems that our society depends on today.One definition is that complexity is the degree of difficulty in accurately predicting future behavior of a system [27], [28].The MIT Engineering Systems Division defined a complex system as a system "with components and interconnections, interactions, or interdependencies that are difficult to describe, understand, predict, manage, design, or change.It includes dynamic, static, and structural aspects" [29].Other attributes that make a system complex include nonlinearity, nonequilibrium dynamics, and interrelated dynamics [27], [30], [31], leading to tipping points that are difficult to predict.
For measuring complexity, many sources of quantitative complexity measures from the literature focus on structural, dynamic, organizational, and information/computational complexity [32], [33].One complexity metric from the physical sciences is graph energy, E .The origin of graph energy was an equation derived by Hückel in 1931 [17] to measure the electron energy for a certain class of organic molecules.The theory was expanded by Graovac, Gutman, and Trinajstić using graph theory [13].The graph-theory form of HMO theory measures the absolute energy levels (ε i ) of the ith molecular orbit were related to the eigenvalues (λ i ) of the adjacency matrix (A) of the molecular graph (G) using the following: The energy of a graph, or the graph energy (E ), was initially represented as E = N i=1 |λ i |, so that (1) became ε π = nα + nβE.Nikiforov extended this theory to any matrix [34] by replacing the eigenvalues with singular values.The resulting graph energy is represented by the following, where n i=1 σ i is the sum of the singular values of a graph's adjacency matrix A: ( The graph-energy concept was adapted by Sinha et al. [11] to measure the structural complexity, C, of a system using (3).In this equation, the complexity contributions are summed for each component i, with α i as the individual complexity of each component, β i j as the complexity of the interface between components i and j, A as the adjacency matrix of the graph, γ as a scaling factor, and E (A) as the graph energy using ( 2) ( Sinha et al.'s research showed that the structure of the system was more centralized with low graph energy, E , and became more distributed with higher E as shown in Fig. 4, redrawn from [11].
Building from the previous work of Sinha et al., which showed that graph energy can be used as a measure of topological complexity for systems, the research in this article implements graph energy to determine how topological complexity affects resilience to tipping points.

B. MEASURING AND MANAGING TIPPING POINTS
One remarkable kind of emergent behavior of complex systems is phase transitions, more colloquially named "tipping points" [35].Phase transitions occur when a system "tips" from one state of behavior into another state of behavior.The measurement of tipping points has been advanced by Solé [36] and van Voorn [37] in the categorization of types of tipping points, illustrated in Fig. 5 redrawn from [36].A bifurcation, or first-order phase transition, occurs when the stable steady state transitions to an unstable state, shown as a dashed line.A separatrix, or second-order phase transition, is when the state of a system crosses a boundary that causes there to be two possible states.
The case studies that were selected for this research are bifurcation phase transitions, where the systems begin with a stable initial state, but then can cross a threshold where they can either continue to operate in a stable state or become unstable.

C. MEASURING AND MANAGING SYSTEM RESILIENCE
In the engineering domain, one definition of resilience, per IEEE Technical Report PES-TR65 and FERC Docket No. AD18-7-000, is "the ability to withstand and reduce the magnitude and/or duration of disruptive events, which includes the capability to anticipate, absorb, adapt to, and/or rapidly recover from such an event" [38], [39].One metric that can be incorporated into resilience management is the magnitude of disturbance the system can absorb before its steady-state shifts to a new equilibrium [40].This type of resilience measure is used frequently in ecology research, leading to some articles naming it an "ecological resilience" [41], [42].Or this approach can be seen as a measurement of precariousness, or how close a system is to its tipping points [6], as shown in Fig. 1.
One method for measuring precariousness is β eff , the resilience index by Gao et al. [9] discussed in Section I.The research that has built on that work mostly includes investigating the practical application of Gao et al.'s methods, e.g., [43], [44], and [45], or comparing other methods or models for measuring large network stability or resilience, e.g., [46] and [47].There are fewer works that directly use and expand on Gao et al.'s resilience index method.They include a model that expands on the mutualistic ecological model used by Gao et al. to include predator-prey dynamics [48] and an investigation of the assumptions and limitations of Gao et al.'s method [49].The following research is the first to investigate the impact of graph energy on this resilience index.

III. METHOD FOR ASSESSING IMPACT OF SYSTEM ENERGY ON RESILIENCE TO TIPPING POINTS
As introduced in Section I, the question for this research is: "does a system's complexity impact its resilience to triggering tipping points?"The approach selected for this research was to investigate how system topological complexity, as measured by graph energy, impacts a system's resilience to tipping points, as measured by the resilience index from Gao et al.The method used to explore this problem is to repeat the method used by Gao et al. in the development and validation of their resilience index, but adding an assessment of how the system's graph energy, used as a measure of topological complexity, impacts the location or nature of tipping points, and applying the new combined method to three case studies.The approach has the following steps illustrated in Fig. 6.
The method developed by Gao et al. [9] starts with a network-based dynamics model of a system.Each node, i, has a differential equation for dx i dt that describes the behavior of the system.In this equation, F (x i ) models the self-dynamics of the node, and G(x i , x j ) models the effect of node i on node j, which is summed over N nodes This set of dynamics equations is translated into the effective plane using the following conversion, in which s in = (s in 1 , . . ., s in u ) T is the vector of incoming weighted degrees in adjacency matrix A β eff becomes a "resilience index" that is dependent on the topology, or structure, or the network.And x eff portrays a combined description of the dynamics of the system in the effective plane.Through this conversion, the set of N number of differential equations is reduced to one differential equation The tipping points are found by solving for the stability points of (7).The critical resilience index located at those stability points is β c eff .To adapt graph theory to interface with this model as a measurement of topological complexity, this new approach begins with a form of (1) adapted into (8), where the total energy contribution of each component (ε i ) includes the energy of each component, α i , the energy of each interface, β i j , and adjustment for the topological energy of the system, E (A binary ).In graph theory, the components are nodes, the interfaces are edges The graph energy, E (A binary ), is the sum of the singular values, σ i , of the binary form (01) of the adjacency matrix A binary Of note, the adjacency matrix used in Gao's method is weighted, A. The methods can begin to be integrated by setting A i j = A i j,binary β i j .A new term, an effective system energy, ε eff , can be calculated by translating (8) into the effective plane using , and ε eff = 1 T Aε 1 T A1 .Then the weighted adjacency matrix element reduces to the resilience index, β eff , as shown in the following: This resulting equation provides insight into how graph energy, E (A binary ), can be adjacently applied to systems in other fields.Conceptually, the right side of this effective energy, ε eff , equation can be understood as a measurement of the energy stored in the structure of the interfaces between the components of a system.For example, when it is applied to an ecosystem model, it represents the topological structure in the interactions between species.One of the effects of the reduction developed by Gao et al. [9] is that the system's model loses information about the system's topological structure, and the nuances in the dynamics caused by that structure.The research in this article predicts that the measurement of graph energy may be used to assess the impact of reduced topological complexity in the system.

IV. CASE STUDIES AND RESULTS
Three use cases were then simulated to assess the relationship between energy and the resilience index calculated with these methods.The first case study repeats the ecological analysis performed by Gao et al., but adds an assessment of how graph energy impacts their method.The second case study analyzes a supply-chain model, to show that these combined methods are extensible to managing tipping points in other kinds of systems that include engineered systems.The third case study assesses a tipping point when a disruptive technology takes over a market.The translation of these case studies into the effective plane and the calculation of their tipping points are in references [9], [50], [51], [52].
The ecological case studies in Gao et al. [9] were repeated, but with calculation of the graph energy, E (A binary ), of each ecosystem.The dynamics were simulated using the following: (11) where the first term, B i , is the migration rate.The second term describes the growth of the species, x i , due to Allee effect and migration [9], with C as the Allee constant, and K as the environment carrying capacity.The third term models the interaction between species i and j, with the parameters D, E, and H used to saturate that part of the function when populations are large.For the ecosystem case studies, these parameters were set equal to the same values as Gao et al., with E i = E = 0.9, and H j = H = 0.1 [9].Following Gao et al., the dynamics equation is recast in terms of the effective plane where β eff is the resilience index and x eff shows population health.The stability points for this equation are found by setting (12) equal to zero and solving for β eff following the steps described by Gao et al.This gives a predicted bifurcation at β c eff = 6.97 [9].Adjacency matrices were formed to describe species interactions in real ecosystems, using the data from [9] [see Fig. 7(f)].The adjacency matrices capture the relationships between species in a mutualistic bipartite ecosystem [see Fig. 7(a) and (b) redrawn from Gao et al. [9]].The simulation from Gao et al. [9], which uses a fourth-order Runge-Kutta stepper (MATLAB function ode45), was adapted to also calculate graph energy using (9).While monitoring graph energy scaled by network size, E/N, where N is the number of nodes, the ecosystems were perturbed by removing network nodes and edges, as illustrated in Fig. 7(c) and (d).To simulate species extinction, a row or column of the adjacency matrix was deleted.In this case, both the graph energy and β eff decreases [see Fig. 7(c)].To simulate ecological (e.g., geographic) separation, an edge between species was removed without necessarily causing species extinction, and the β eff decreases, but the graph energy initially increases [see Fig. 7(d)].After many ecological separations, the ecosystem is so weakened that it has few interspecies connections left.At this point, every subsequent separation causes cascading extinctions among species as they lose their mutualistic support.Both the graph energy and β eff values decrease and the ecosystem collapses.For each iteration on the adjacency matrix, the species interaction was simulated to determine if the ecosystem sustained itself or collapsed to zero population [labeled with black stars, Fig. 7(e)].Color dots show results of 50-100 simulation runs.
Results show that for systems with higher graph energy, the ecosystems did not start collapsing at β c eff = 6.97, as predicted by Gao et al., but instead weakened further into a lower β eff region.Ecosystems with higher graph energy exhibit increased resilience compared to those with lower graph energy.Using Fig. 4 as a reference, the ecosystems with a higher E are more distributed in their topology, such that in the mutualistic ecosystem simulation they would likely maintain more supportive connections between species than a more-centralized lower E ecosystem, which would have fewer supportive connections between species.Those supportive connections would contribute to that increased region of resilience for the higher graph energies.For the supply-chain case study, the calculation for graph energy was incorporated into the simulations of the randomly generated supply-chains from references [50], [51].Using the linearized supply chain dynamics model from Helbing et al. [53] with the following variables.1) n i , x i , and y i are the linearized deviation from the steady state inventory, production, and external flow of goods of kind i with respect to the steady state.2) c i j is an element of the supply matrix C which is a ratio of the goods of kind i needed to produce goods of kind j. 3) A is a scalar parameter representing the negative of the partial derivative of the desired production rate with respect to the current inventory.A is assumed constant across all goods i. 4) B is a scalar parameter representing the negative of the partial derivative of the desired production rate with respect to the current inventory change rate.B is assumed constant across all goods i.The translation of these equations into the effective plane is discussed in Edwards et al. [51] and uses the following conversions: T is the vector of incoming weighted degrees in matrix C. The translation into the effective plane results in 16) With these equations, the dynamics of the entire network are modeled with two variables that describe the deviation from the steady state inventory and production in the effective plane (n eff and x eff ), one driving function (y eff ), and two adjustable constants to model the response of all nodes to changing demand and inventory stockpiles (A and B).The stability points of these equations are at β c eff = 1, above which the supply chain becomes overloaded with maximized production rates and zeroed inventories [50].
The simulations for the case study were run with adjacency matrices from the production processes for polystyrene egg trays and recycled egg trays [51], [54].The percentage of reuse for recycled egg cartons was randomly varied as an integer between 0% and 100%, and simulated to a distant time of t = 1000 weeks with a MATLAB function ode45 solver [51].
Results from these simulations reveal a similar impact of graph energy on the location of the tipping point as the results from the ecological modeling.Fig. 8 plots the graph energy, E, of the binary form of each supply matrix, normalized by the matrix size N, against its β eff .Systems that have unstable supply chains are plotted with a mauve star.Results show that when the resilience index is above the predicted tipping point, when β eff > 1, the supply chain could not sustain itself because it had an internal dependency on reusing more material than it was creating.In this respect, β eff could be a measure of material conservation in a supply chain.If this model were expanded to include the big-data of all internal cycles and external limits for a large complex supply chain, then it could be a powerful tool to detect changes in the data that would lead to resource depletion.Similar to the ecosystem analysis, the region of system instability shifts with the energy level of the system, and could benefit from a correction like the illustrated β c eff,corrected .
For the disruptive technology case study, the calculation for graph energy was incorporated into the simulations of the technology markets from reference [52].The model is adapted from the following predator-prey model from Pielou [55]: Ünver [56] performed a study of technology industry dynamics using this Lotka-Volterra model, in which the disruption of digital cameras in the film-camera industry was simulated.Building on that previous work, the problem can be decomposed into a set of equations, with x 1 representing the units produced with incumbent technology and x 2 representing the units produced with disruptive technology.The adjacency matrix A i j captures the sales that each technology is taking away from the other This overall equation in the effective plane can be further decomposed to enable the study of the interactions of different parameters or subsystems.In this case, the behavior of the incumbent (x 1 ) and disruptive (x 2 ) technologies is separated into x 1eff and x 2eff to distinguish the interaction between these two parameters.The system is decomposed using the following conversion equations: (21) resulting in the following equations in the effective plane: with the following stability points: The case study simulated was the same film-camera disruptive case as Ünver [56], but with the simplification of a = a 1 = a 2 = 0.The results of 100 000 simulations are shown in Fig. 9.When the technology is successful and has significant sales at the end of the simulation, then its resilience index and graph energy is shown as a blue dot.However, when the incumbent technology is unsuccessful and has zero sales at the end of the simulation, or when the disruptive technology has less than 8e 5 sales, then the results are sown as a black star.The actual tipping point of these systems is located on the edge of the black stars, a boundary across which the technology can go from always succeeding to potentially failing.This tipping point boundary clearly shifts away from the predicted tipping point, β c eff , as the graph energy of the system increases, with a similar pattern to the previous case studies.
Combining the results of each of the case studies, Fig. 10 plots E against β eff − β c eff , so that all the graphs have the same tipping-point location at 0. The graph energy, E, was not normalized by N, which was theorized by Sinha et al. [10] and included in the first two case studies in this research, because in the derivation of ( 8) and ( 9) the authors determined that there was no such factor.The data are labeled with a star ( * ) when the ecosystem collapsed, the supply chain was unsustainable, or the technology failed to gain or maintain a market share.
Based on these results, the estimate for tipping-points (β c eff ) for all systems using this method can be adjusted to slope inversely with graph energy, E. Assuming an inverse relationship between β eff and E [based on the effective system energy (8)], a corrected resilience index β c eff,corrected is found using the enveloping edges between stable and unstable regions.A MATLAB script was used to "grab" the data points on the upper edge of the unstable ecosystems and failed disruptive technology, and the lower edge of the unstable supply chains and failed incumbent technology.These data points, circled in blue, define the tipping-point boundary.Then MATLAB was used to calculate an exponential fit of these points on the boundary between stable and unstable systems.The results show that increasing the graph energy leads to a greater offset in the system tipping-point threshold, following this exponential curve.This correction defined by ( 25) is shown in Fig. 10 as a blue line Ultimately, systems engineers who are studying the tippingpoints of systems want to either prevent the triggering of an undesired tipping-point, or to control a system through a desired tipping point.To enable this kind of systems management, a measurement of precariousness, or how far a system is from a tipping point, is desirable.The adjusted resilience index, δ eff in the following, provides a direct estimate of precariousness, the distance of a system to its tipping points, and it includes the correction derived in (25): When the method from Gao et al. is applied to a system modeled as a network, (25) can be applied to incorporate a correction to the location of the predicted tipping point, and (26) can be implemented to calculate how close the system is operating to that tipping point.Returning to the original research question, "does the complexity of a system affect its resilience to transition or collapse?"These results provide a measurement for how the topological complexity, as measured by graph energy, affects a system's precariousness to tipping points, as measured by the resilience index β eff .

V. INTEGRATION INTO SYSTEMS ENGINEERING PROCESSES
This method can be integrated into systems engineering processes for measurement and management of system resilience.
During system development, the dynamics of the system can be modeled for scenario analysis, as shown in Fig. 11.The dynamics of the system under development can be modeled, and the disruption potential can be measured using the adjusted resilience index in (26).If a system under development shows a particularly low adjusted resilience index, changes to the design can be assessed for enabling a higher adjusted resilience index.The results can then be factored into system development requirements.Through this consideration, the model insights shape the development and ultimately influence the resulting system.Similarly, this approach can be applied throughout the development and after deployment, potentially incorporating it into risk and change management processes.Fig. 12 depicts how this method would integrate as part of a cyclical process, which can be repeated as necessary to update the model with new information and, as a result, reduce uncertainty and risk.As depicted, the method in a loop enables continuous prediction and adjustment that can be aligned with other cyclical systems engineering processes.

VI. CONCLUSION
In summary, this research expanded on a theory for a universal resilience index measurement from Gao et al. [9] by introducing a complementary measurement of system topological complexity using graph energy that, when combined, shows how topological complexity impacts the location of tippingpoints in large-scale systems when modeled as networks.Simulations from case studies in different fields provided data to investigate this relationship.The results were the introduction of a correction to Gao et al. [9] that adjusts for the topological complexity initially lost in their method's modeling reductions and reveals the location of tipping points with greater accuracy for the measurement and management of system resilience.Also, an adjusted resilience index was introduced that more directly estimates precariousness, the distance of a system to its tipping points.These adjustments were fit to the data from the three case studies, and the resulting equations can be used as a correction for the method when applied to these kinds of systems models.Future studies can investigate the generalizability of this method with a wider set of experiments, potentially resulting in a refined fit or modified form of equations.
Overall, this approach can bring new breakthrough insights to domains spanning from mass extinctions to critical infrastructure systems, shedding light on the status of their strength and vulnerability to potential collapse.It provides a means to better understand the tipping points of systems that are so critical to sustain, with a measurement that can be used in system management and risk mitigation to actively avoid or potentially control a system through a tipping point.

FIGURE 1 .
FIGURE 1. Precariousness of a system to a tipping point.

FIGURE 3 .
FIGURE 3. Combination of system resilience and complexity theories to answer research question.

FIGURE 4 .
FIGURE 4. When graph energy is calculated from a graph model of a system, it can be used to measure system complexity [11].

FIGURE 7 .
FIGURE 7. Bipartite ecological networks (a) are formed into projection networks, which are the basis for the weighted adjacency matrix (b).Graph energy patterns of real ecosystems as they are perturbed using node loss (c) and edge loss (d).The actual region where ecosystems collapse (labeled with black * ) is offset from the predicted tipping point using the universal resilience index, β c eff (labeled with a red line) indicating a correction (β c eff,corrected , labeled with a blue line) (e).Each color represents a different ecosystem (f).

Fig. 7 (
e) shows how the offset between where the Gao et al. prediction of the tipping point, β c eff , could be corrected with a new prediction, β c eff,corrected .

FIGURE 8 .
FIGURE 8. Tipping point decreases as graph energy of the supply matrix increases in supply-chain model.

FIGURE 10 .FIGURE 11 .
FIGURE 10.Boundary between stable and unstable systems shifts with E, suggesting a correction to the tipping point estimate, β eff,corrected c .

FIGURE 12 .
FIGURE 12. Integration of adjusted resilience index with cyclical systems engineering processes.