Economic networks: Heterogeneity-induced vulnerability and loss of synchronization

To frame our study, we use the well-established definitions of resilience proposed in ecology—where the concept originated¹¹ and adapted to the analysis of dynamical systems in general. Resilience can be evaluated using one or more of three metrics: the size of the basin of attraction of a specific dynamical regime;^12,13 the speed at which a system, after a given perturbation, bounces back to its attractor;¹⁴ and the intensity of the perturbation that a system can sustain without switching to another, less desirable attractor.¹⁵ The word “metric” here is not used in its strict mathematical acceptance as a distance between two points in an arbitrary space but in the broader sense used in measuring performance (e.g., Ref. 16).

We study these three metrics in a toy model for production networks. Dynamical models have been proposed to examine the propagation of supply disruptions in interfirm or intersector input–output networks. The emergence and the amplitude of cascading failures^17,18 and their potential consequences for production location¹⁹ have been investigated. The relationship between the network structure and vulnerability has recently been explored,^20,21 as well as the role of adaptive agent behavior²² and that of delays.²³ We note that the questions addressed by these models do have conceptual connections with epidemic modeling, which focuses on the spreading of human diseases or of computer viruses across a network of agents.²⁴ 

In this paper, we focus on the role of heterogeneities that are present at the network, node, and link levels. First, production networks are reported to have a high heterogeneity in connectivity, with few highly connected firms interacting with more peripheric ones. Fujiwara and Aoyama²⁵ have analyzed the power-law distributions of the in- and out-degrees of empirical production networks. The heterogeneity of these distributions influences how localized shocks aggregate.²⁶ The vulnerability of various network structures has been studied in other types of models: in the general network theory (e.g., Ref. 27), in load network dynamics,^28–31 in epidemic-spreading models^32,33 and, most recently, in networks of networks.³⁴ 

Furthermore, many node-level characteristics are heterogeneous in economic networks, such as their size.³⁵ Gabaix³⁶ found that a power-law distribution of firm size can affect the aggregate volatility of economies. Few models have incorporated firm-level heterogeneities in economic dynamics, with the exception of Hallegatte,³⁷ who found that heterogeneity in capacity or in inventory, can amplify supply-side shocks.

Finally, all links are not equivalent. In particular, they are associated with delays, which primarily originate from decision-making, production and transportation. They represent a key concern in operations research,³⁸ since they can lead to mismanagement of inventories, synchronization issues, and oscillations. A key example, which has been empirically and experimentally verified, is the bullwhip effect.³⁹ Still, with the exception of the stock market model developed by Miśkiewicz,⁴⁰ no economic network model incorporates link delays. In epidemic models, delays have been shown to foster the outbreak of epidemics,^41,42 but the effect of delay heterogeneity has not been tackled so far. Buzna, Peters, and Helbing⁴³ have developed a model of load dynamics on networks with uniformly distributed delays, but the role of inhomogeneity in the delays was not included.

Understanding the role of delay heterogeneity in networks is therefore of interest in a variety of fields. For economic networks Coluzzi et al.²³ took a first step in this direction, and their contribution is the basis of the present paper. We investigate here the influence of heterogeneity in connectivity, inventory, and delays on the resilience of production networks. To do so, we formulate a model governed by a coupled system of Boolean delays equations (BDEs), whose Boolean variables evolve in continuous time.^44,45

The BDE framework is useful in building simple models of complex nonlinear systems. They allow one to simplify many components while fully retaining some core dynamical ingredients: threshold behavior, multiple feedbacks, and distinct time delays.⁴⁶ In fact, the behavior of Boolean networks with homogeneous delays that use exclusively the “AND” operator is well documented.⁴⁷ More general BDEs—with arbitrary Boolean operators and with heterogeneous delays—have been successfully applied to the study of climate dynamics (e.g., Refs. 48 and 49), earthquake physics,^50,51 and genomics.⁵² 

Our BDE model is based on Coluzzi et al.,²³ extended to incorporate inventories. It encompasses four elements: (i) production based on complementary inputs; (ii) heterogeneous delays in production and delivery; (iii) heterogeneous inventories; and (iv) supplier–buyer networks with heterogeneous connectivity. The addition of a simple economic process allows us to go one step beyond the static effect of network structures on resilience.

Production, which is usually modeled as a real or integer-valued variable, is represented here by a Boolean variable. An agent with production 1 is said to be fully active: it meets the needs of its customers and maintains its inventories. A production of 0 does not imply that the agent is completely inactive, only that it is affected by some sort of disruption, such as quality issues, delayed or lower-than-expected production. Thanks to this simplification, we can focus on the dynamics of propagation that follows initial perturbations. We thereby choose to not monitor the numerical value of economic damages, in order to better track the spreading of localized shocks, and analyze the role of the key nonlinear components of this economic dynamics: production interdependencies, heterogeneous delays, and inventories.

To analyze the effects of heterogeneous connectivity, we developed a class of scale-free networks that exhibits some empirical features of real economic networks. Specifically, Fujiwara and Aoyama²⁵ studied a large dataset of over 1 million Japanese firms and their supply relationships, provided by the company Tokyo Shoko Research Ltd. They found that the distributions of in- and out-degrees have a power-law distribution for degrees between 10 and 1000, followed by an exponential decay for higher degree, and an average connectivity of 4. The exponents of the cumulative distributions $P > ( k )$ of in- and out-degree are 1.35 ± 0.02 and 1.26 ± 0.02, respectively, where the error corresponds to 1.96 times the estimated standard deviation. We generate a class of networks that reproduce these features, which we call SF–FA networks—SF for scale-free and FA for Fujiwara and Aoyama—and compare their behavior with that of Erdős–Rényi (ER) networks. This type of nationwide datasets has recently been used to assess the propagation of disruptions sparked by natural disasters.⁵³ 

Specific linkages between two or more sectors may severely affect recovery. For instance, Hallegatte⁸ reports that, following hurricane Katrina, the lack of housing and health care services hindered the return of workers and the reconstruction, which were crucial, however, for the recovery of these services. The timing of such cross-sectorial interactions is likely to play a major role in the joint recovery of both sectors. This type of problem motivates our applying BDE models to networks of networks, whose usefulness in representing interdependence between specific economic sectors has been already demonstrated.³⁴ The setting of nested networks allows us to study how synchronization and the loss thereof may affect network resilience.

We provide results on network resilience to local perturbations in the following setting. Initially, all agents fully produce except one that gets temporarily disrupted. This initial perturbation may spark a wave of disruptions that either vanish or induce permanent disruptions. This phenomenon is well illustrated by the so-called Albuquerque accident: a localized and short fire sparked by a storm generated major disruptions in the mobile phone industry.⁵ Through extensive simulations, we characterize the network behavior and evaluate the influence of heterogeneity in inventories and in delays. We use the three previously mentioned metrics to gauge resilience: size of the attractor basin of the full activity state, transient speed, and perturbation intensity. We then propose a resilience indicator and discuss its implications for resilience building. Our results also include intriguing synchronization behavior, such as metastable relaxation oscillations, recalling the importance of delays in network synchronization.⁵⁴ 

In Sec. II, we first formulate the BDE model and its delays. Next, the simulation routines and the algorithm generating the classes of networks investigated are specified. In Sec. III, we start our investigation by studying how a single perturbation propagates in fundamental network motifs. We build a probabilistic model to analyze the propagation in simple acyclic structures in Sec. III A, and in Secs. III B–III D, we identify the key mechanisms that can prevent network recovery.

Section IV provides results for large interfirm networks. Guided by the results of Sec. III, we study in Sec. IV A, the global and local structures of SF–FA networks that influence resilience, and contrast this analysis with ER networks. Section IV B presents the effects of heterogeneity in inventories and in delays.

In Sec. V, we turn to the systems of interdependent networks, and show how differences between the networks can affect vulnerability via loss of synchronization. We measure in Sec. VI, the critical role played by individual agents and derive an aggregate resilience indicator. Finally, we discuss in Sec. VII, potential implications for resilience building, suggest several extensions, and comment on the role of synchronization in network resilience.

The economy is modeled here as a network, in which each agent is represented as a node and produces using inputs from other agents. These dependence relationships are represented by directed links. The topology of the network is captured by an n × n connectivity matrix $M = ( m i j )$⁠, where n is the number of agents, and $m i j = 1$ if agent i is a supplier of agent j, and $m i j = 0$ otherwise. No agent can supply itself, i.e., $m i i = 0$⁠. Some agents are primary producers, and produce without supplier.

Each of the n agents can only be in one state, modeled as a Boolean variable x_i with $i = 1 , .. , n$⁠: either it is fully active, x_i = 1, or it is disrupted and x_i = 0. A fully active agent is able to deliver its outputs to its clients.

The Boolean variable x_i is continuously updated over time via delay equations. Each supplier–buyer relationship is characterized by a delay, denoted by $τ i j$⁠; this delay represents the time it takes agent j to produce and deliver inputs to agent i. With the exception of primary producers, all agents need inputs from each one of their suppliers, and they produce upon the simultaneous reception of all inputs.

In the absence of inventories, full activity of an agent i depends exclusively on the set $S ( i )$ of its suppliers. When agent i is a primary producer, the set $S ( i )$ is empty and the agent is always active. Otherwise, activity of agent i is maintained at t only when all suppliers j were active at $t − τ j i$⁠. The n variables that represent activity are updated as follows:

Here, the operator $∏$ is the repeated Boolean product “AND,” and we use the usual notation AND $≡ ∧$ and OR $≡ ∨$⁠.

We assume that agents hold inventories that allow them to stay active when their supply is disrupted. A second Boolean variable, y_i, is attached to each agent: it represents the availability, $y i ( t ) = 1$⁠, or unavailability, $y i ( t ) = 0$⁠, of agent i's inventory at time t. Disruption, $x i ( t ) = 0$⁠, only occurs if agent i does not receive some supply while its inventory is depleted; see Eq. (2a).

Inventory is created through production: when agent i receives all its supply simultaneously, i.e., when $∏ j ∈ S ( i ) x j ( t − τ j i ) = 1$⁠, then i produces and it both delivers outputs to clients and replenishes its own inventory; the latter will be available after λ_i time units; see Eq. (2b) below. This real-valued quantity represents the size of the inventory: agents with large λ_i can stand longer external perturbations, because their production can be stored for a larger period of time.

For non-primary producers, the dynamical system for $x i ( t )$ and $y i ( t )$ is

Equation (2b) can be substituted into (2a), yielding a single dynamical equation for activity

Note that the presence of the Boolean product operator AND, as well as of the addition operator OR, generates time irreversibility. Two variables $x i ( t )$ and $x j ( t )$ that are coupled by AND will tend to 0 after finite time, while coupling by OR will push them towards 1. In either case, it is not possible to unambiguously infer past states from future states. The system of BDEs in Eq. (2) is, therefore, said to be dissipative, in the precise sense defined by Ghil and Mullhaupt⁴⁵ and reviewed by Ghil, Zaliapin, and Coluzzi.⁴⁶ We thus expect the solutions for all the network topologies and combinations of delays and inventories that will be considered herein to contain initial transients and lead to asymptotic behavior—whether stationary, periodic or aperiodic—in finite time. Some of the transients, though, may turn out to be very long, as found already by Coluzzi et al.²³ 

The n continuous equations (2) are simulated using an ad-hoc C++ program. They are evaluated every Δ units of time, which sets the time resolution of the simulations. We typically used a unit of time $Δ t = 0.01$⁠, which implies that all time-valued parameters and variables have two decimals. Each agent is an object, or node, identified by an integer i, in which the activity $x i ( t )$⁠, inventory duration λ_i, the identification numbers of the suppliers $S ( i )$ and the associated time delays $τ j i , j ∈ S ( i )$ are stored. The graph structure is given by the list of $S ( i )$⁠. Note that, to calculate the state x_i at time t, as in Eq. (2), we need the activity x_j of its suppliers at $t − τ j i$ and $t − τ j i − λ i$⁠. We therefore need to store within each agent the trajectory of its activity x_i between $t − τ j i − λ i$ and $t − Δ t$⁠.

This dynamical system has a trivial fixed-point attractor: if all agents are active, i.e., if $x i ( t ) = 1$ for $i = 1 , … , n$⁠, the system remains stationary. We study the resilience of this attractor, called full activity, to small perturbations: how quickly and how far do external perturbations propagate? Does the network always get back to full activity, or can perturbations push the system into other alternative states? How do the network topology and the heterogeneity of delays and inventories influence such behavior?

An external perturbation is defined as the shutdown of a single agent's production over a certain time interval, denoted by δ₀. We primarily focus on the consequences of such perturbations that occur in isolation. This case represents for instance an accident occurring at a single production facility, which can then propagate to other locations, as in the Albuquerque accident mentioned in Sec. I.⁵ For an initial perturbation δ₀ affecting agent r, the initial state is formulated as follows: for all $i = 1 , … , n$ and for $t ≤ 0 , x i ( t )$ equals one. From t = 0 onward, all the $x i ( t )$'s are updated according to Eq. (3), except $x r ( t )$⁠, which is maintained at 0 for $0 ≤ t < δ 0$⁠. Unless otherwise specified, we set $δ 0 = 1$⁠.

We will monitor how such initial perturbations propagate throughout the network. In particular, they may create propagating “waves” of disrupted agents that are either damped out or continue to propagate indefinitely. The key observable is the relative number of fully active agents over time, denoted by $ρ ( t )$

The asymptotic regimes can reach alternative attractors, which are characterized by their dimension in the space $ℤ 2 n$⁠, where $ℤ 2 = { 0 , 1 }$—the attractors can be fixed points, as well as periodic or aperiodic attractors. To examine the extent of economic damages in the asymptotic regime, we measure the average asymptotic value of fully active agents, denoted by $ρ ∞$⁠. If the asymptotic regime is a fixed point, $ρ ∞$ is the value of $ρ ( t )$ at this point. If it is periodic, $ρ ∞$ is the average value of $ρ ( t )$ over one period. When $ρ ∞$ is below 0.5, most agents are not fully active anymore, and the network is said to have collapsed.

Transients are characterized by the speed and amplitude of the propagating wave. We observe, for a given agent k, the ensemble $A ( k )$ of agents that get perturbed at least once following the initial perturbation, i.e., $A ( k ) = { i : ∃ t ≥ 0 , x i ( t ) = 0 }$⁠. The cardinality of $A ( k )$⁠, denoted by $α ( k )$⁠, measures the propagation amplitude. We also record the time at which each agent in $A ( k )$ gets perturbed for the first time. The time of first perturbation, denoted by t_i, is such that $x i ( t ) = 1$ for $t < t i$ and $x i ( t i ) = 0$⁠. The average time of first perturbation, denoted by θ, is a proxy for propagation speed: for comparable amplitudes α, a wave with a small θ propagates faster than the one with a high θ.

We start in Sec. III by a detailed investigation on how perturbations propagate in fundamental network motifs, such as lines, trees, and loops. We particularly focus on interacting multiples loops. To create each of these types of networks, we specify for each agent i the set of suppliers $S ( i )$ that leads to the desired topology. For instance, for a line of 3 agents, we set $S ( 1 ) = 0$ and $S ( 2 ) = 1$⁠. The studies of these simple motifs allow us to identify the key propagation mechanisms that can lead, in more complex networks, to alternative attractors. Next, for each motif, we build a probabilistic model to analyze how particular distributions of delays and inventories influence propagation. We use in particular, lognormal, exponential, and uniform probability distributions.

In Sec. IV, we analyze how the distributions of delays and inventories affect the resilience of SF–FA and ER networks. To create these more complex networks, we use or adapt existing algorithms that generate C++ network objects from the igraph library.⁵⁵ We then transform these objects into a list of sets of suppliers $S ( i )$⁠.

To generate SF–FA networks, we use the method first proposed by Goh, Kahng, and Kim,⁵⁶ in the form used by Chung and Lu.⁵⁷ The number of incoming and outgoing links attributed to each node is proportional to a weight, or load factor; this method is implemented in the so-called fitness game algorithm of the igraph library. For a network with n nodes, the factor of the i^th node is $( i + i 0 − 1 ) − 1 / ( γ − 1 )$⁠, where γ is the targeted power-law exponent and i₀ an integer used to introduce an exponential decay for high degrees. To match the results of Fujiwara and Aoyama²⁵ for their huge data set of Japanese firms, we use γ-values of 1.35 for in-degrees and of 1.26 for out-degrees, while i₀ is set to 12. The connectivity c is defined as the average number of in- and out-degrees per agent, and it is set to 4, i.e., there are 4 times more links than nodes.

Directed ER networks are generated through the so-called $G ( n , m )$ algorithm, which is directly accessible in the igraph library; here n denotes the number of nodes and m the number of links.⁵⁸ To compare with SF–FA networks, we choose $m = 4 n$⁠. Unless otherwise specified, all networks have $n = 10 3$ agents. This number is large enough to be in qualitative agreement with the behavior of larger networks, and low enough to be computationally feasible when exploring large regions of the parameter space.

We design a simulation routine to comprehensively assess the behavior of the network from all possible initial perturbations. For networks with n agents, this routine consists of n simulations. Each of these simulations corresponds to the initial perturbation of a different agent. For instance, in the k^th simulation, agent k is initially perturbed, and we measure the propagation amplitude $a ( k )$⁠, the average time of first disruption $θ ( k )$⁠, and the asymptotic average density of fully active agents $ρ ∞ ( k )$⁠. The statistics of the n-member ensembles ${ α ( k ) } , { θ ( k ) }$ and ${ ρ ∞ ( k ) }$ captures the system's response to small perturbations. We denote the mean of these three ensembles by A, Θ, and $R ∞$⁠, respectively.

We first study the homogeneous case, where all delays and inventories are equal, i.e., $τ i j ≡ τ$ and $λ i ≡ λ$⁠. We track how the three ensembles ${ α ( k ) } , { θ ( k ) }$ and ${ ρ ∞ ( k ) }$ change in the two-dimensional (2-D) parameter space $( τ , λ )$⁠. Then, to assess the effect of heterogeneities in delays and inventories, we draw each delay $τ i j$ and each inventory λ_i from two lognormal distributions with means $( μ τ , μ λ )$ and standard deviations $( σ τ , σ λ )$⁠, respectively. We then track how the dynamical response changes along an increase in $σ τ$ and $σ λ$⁠. We use as a baseline the homogeneous case where τ  =  1 and $λ = 0.2$⁠.

In Sec. V, we study how heterogeneities in delays and inventories affect propagation in acyclic networks, where disruptions never persist indefinitely, i.e., where $R ∞ = 1$⁠. Weisbuch and Battiston,¹⁹ for instance (see Fig. 1 there), used homogeneous delays and a regular lattice—which can be seen as a particularly simple case of acyclic topology—to model production networks.

We generate random acyclic networks (RAs) of connectivity c, in two steps. First, we build an ER network of connectivity 2c, and extract its connectivity matrix; then we use the upper triangle of this matrix to generate the RA network. In RAs, there is a natural direction of the production flow from primary producers, i.e., agents without a supplier, to final producers, i.e., agents without a customer.

We characterize the sensibility of RAs to a small perturbation by evaluating the extent to which final producers are impacted by initial perturbations of primary producers. We measure, for each primary producer, the average time final producers are disrupted. The average of this quantity, divided by the duration of the initial perturbation δ₀, defines the transfer coefficient, denoted by TC. To give a precise definition to TC, we denote the ensemble of primary producers by I and its cardinality by n_I, the ensemble of final producers by J and its cardinality by n_J, and the time series of the activity of agent i following the initial perturbation δ₀ of agent j by $x i j ( t )$⁠. The transfer coefficient TC is then defined by

Finally, we examine the resilience of a network of interdependent RA networks. We focus on the case of two RA networks, where the primary producers of one network are supplied by the final producers of the other. To generate these interlocked networks, we produce two RA networks, and each primary producer of the first RA network randomly picks a supplier among the final producers of the second one; reciprocally, primary producers of the second RA network pick suppliers among the final producers of the first one. Figure 1 illustrates all classes of networks studied in the paper. Table I summarizes the parameters, variables, Boolean operators and observables used in the model.

Consider a linear network, illustrated in Fig. 1(a), with $n ≫ 1$ agents. It represents an isolated supply chain in its simplest form. Agent 0 is the primary producer, it supplies agent 1, which itself supplies agent 2, and so on up to agent n – 1. With this structure, the dynamical system of Eq. (3) becomes, for i > 0

Let us first consider a simple case without inventory, i.e., when $λ i = 0$⁠, and with homogeneous and unitary delays, i.e., $τ i − 1 , i ≡ 1$⁠. If agent 0 is externally perturbed at t = 0 during δ₀ units of time, agent 1 will be disrupted at t = 1 during δ₀ units of time as well, so will agent 2 at t = 2, and so on. There is a wave of disruptions that progresses steadily along the line. If the delays are heterogeneous, the rhythm of this progression depends on the particular sequence of delays.

Inventories can reduce the wave of disruptions. For instance, if agent 1 holds inventories, say $λ 1 = 0.5$ unit of time, and the initial perturbation δ₀ is lower than λ₁, agent 1 remains active and the wave stops. If, however, δ₀ is higher than λ₁, then at $t = 1 + λ 1$ the inventory is exhausted and agent 1 is disrupted during $δ 0 − λ 1$ units of time.

In the general case, by applying Eq. (3) to agents 0 and 1, we can demonstrate that, if agent 0 is externally perturbed during δ₀ units of time, then agent 1 will be temporarily disrupted provided the initial perturbation is longer than its inventory, i.e., if $δ 0 > λ 1$⁠. In this case, agent 1 will be perturbed at $t = τ 01 + λ 1$ for $δ 0 − λ 1$ units of time. One can iteratively apply this result to agents 1 and 2, then 2 and 3, and so on, and derive the two following general results.

First, agent i will be perturbed if the duration of the external perturbation exceeds the sum of the inventories of all agents between agent 1 and itself, i.e., if

Next, if inequality (7) is satisfied, then the duration of the perturbation affecting agent i is given by the difference between the duration of the initial perturbation and the sum of the inventories $δ 0 − ∑ k = 1 i λ k$⁠. Both the delays and the inventories determine the time of first perturbation

These results correspond to most real-world observations of supply chain disruptions, in which a perturbation propagates and is gradually damped by inventories. Formulas (7) and (8) above allow us to easily characterize the resilience of such linear networks in the homogeneous case. When all delays and inventories are equal to τ and to λ, respectively, the initial perturbation propagates up to agent $z ≡ ⌊ δ 0 / λ ⌋$⁠, where $⌊ x ⌋$ denotes the integer part of x. In linear networks, this ratio also corresponds to the number of disrupted agents α, which is thus inversely proportional to the inventories λ, and is not affected by delays. The time t_i of first disruption of agents within the set $A$ is $t i = ( τ + λ ) i$⁠. The average time θ of first disruption is given by $θ = ( τ + λ ) ( a + 1 ) / 2$⁠: it varies linearly with the sum of the delays and the inventories, $τ + λ .$ In such a linear network, there is no permanent disruption: perturbations eventually decay.

In reality, supply chains are more reticulated, because agents have multiple suppliers and customers. So-called trees, which are illustrated in Fig. 1(b), are the other building blocks of such networks. The dynamics of disruptions is qualitatively similar, with decaying waves of disruptions propagating across the tree. Quantitatively, we can readily extend the results obtained for linear networks to trees, as follows: Layer 0 is occupied by one root agent, which supplies d agents in layer 1, d² in layer 2, and so on. Ternary trees, with d = 3, were used, for instance, by Zaliapin, Keilis-Borok, and Ghil^50,51 in studying colliding cascades of rupture and healing in earthquake dynamics. In trees, the initial perturbation of the root agent can lead to significant propagating waves. The formula $z ≡ ⌊ δ 0 / λ ⌋$⁠, derived for linear networks, can be applied, where z now represents the level to which the perturbation propagates. Given that layer k is occupied by d^k agents, we have

For the heterogenous case, we build a probabilistic model to assess the effect of a change in key features of the delay and inventory distribution on propagation. Letting each delay and inventory be an independent, identically distributed random variable, our probabilistic model is based on a recursive algorithm that yields the expected values of the amplitude and average time of first disruption, $E ( α )$ and $E ( θ )$⁠. The algorithm is presented in  Appendix A and its results are compared with those of the numerical experiments. We find that, for both a uniform and a lognormal distribution, the more heterogeneous inventories are, the further the perturbation tends to propagate. In these simple network structures, in which pathways never interesect, delay heterogeneity does not affect the propagation amplitude.

These results on linear networks and trees will help us understand the propagation of disruptions in large RA networks in Sec. V A. Such networks can represent the supply chain of large industries taken in isolation, such as electronic or automotive industries.

In fact, the supply chains of different industries interact, leading to large and complex production networks, as documented by Fujiwara and Aoyama.²⁵ In particular, they exhibit circularity. For instance, transportation firms depend strongly on the fuel supply, which itself depends on transportation services.

While propagating disturbances eventually vanish in acyclic networks, they can be sustained or even amplify in purely cyclic networks. In the latter, asymptotic behavior involves recovery, periodic disruptions or collapse; such a simple, isolated loop is illustrated in Fig. 1(c). Based on the findings of Subsection III A on linear networks, the following results are easy to demonstrate.

In the absence of inventory, an initial perturbation propagates indefinitely around the loop: it periodically affects each agent, and the period is equal to the sum of the delays. If other agents depend directly or indirectly on this loop, they are also affected by the periodic wave of disruptions.

When some agents hold inventories, propagating waves are gradually damped and eventually die out; their behavior can be derived using the analytic approach taken for the linear network in Sec. III A. This damping can however be canceled by short time delays: if the initial perturbation propagates around the loop so rapidly that the first agent has not yet recovered from its initial perturbation, then all agents remain permanently disrupted. In an isolated loop with agents labeled $k = 0 , 1 , … , n − 1$⁠, the condition for such fast collapse is

More complex types of behavior occur in networks in which multiple loops interact. In this case, a new type of phenomenon arises, which we call “fine detuning”. Essentially, it involves two or more periodicities associated with the distinct loops being close to commensurable i.e., to having an exact rational relationship. The precise definition is given in Eq. (11) below. The present subsection describes this mechanism and its origin. A rigorous proof is given in  Appendix B.

When two loops interact through one agent, as in Figs. 1(d) and 1(e), any perturbation reaching this pivotal agent duplicates, inducing one propagating wave in each one of the two loops. In the absence of inventory, these two waves propagate freely, duplicate again at the pivotal agent, leading therewith to new propagating waves in each loop, and so on.

If the propagation times through either loop are equal, then the two waves return to the pivotal agent at the same time. This exact overlapping leads to a periodic regime, as it does in an isolated loop. More generally, if one propagation time is an integer multiple of the other then, in the absence of inventory, a periodic regime sets in. In the absence of such an integer-multiple relation between propagation times, any initial perturbation gradually amplifies through duplication and eventually leads to permanent disruption of all agents.

Contrary to isolated loops, explored in Sec. III B, inventories do not guarantee recovery in interacting loops, even for moderate initial perturbations. We refer to the case in which the sum of the delays in one loop and that in the other differ by a small amount as fine detuning. In this case, a single perturbation can permanently disrupt all agents through a process of successive duplication and recombination of propagating waves.

When the difference in delays between both loops—i.e., the detuning—is large and incommensurable, the wave that propagates in one loop behaves independent of the one in the other loop; both waves gradually decay due to inventories, as seen in the isolated loops in Sec. III B. When both loops are tuned, i.e., when the sums of delays in each loop are equal, the waves are perfectly superimposed at the pivotal agent, and they synchronously decay. Only when the difference in delays is small, i.e., when detuning is fine, do the waves superimpose with merely a small phase shift and thus amplify.

More precisely, for an initial perturbation δ₀ that affects the pivotal agent, fine detuning occurs when the difference in total delays $Δ T$ between the two loops falls within the following range:

where L_a is the sum of inventories in the faster loop, labeled a, L_b is the sum of inventories in the slower loop, and λ₀ is the inventory of the pivotal agent. For this condition to hold, it is furthermore necessary that $δ 0 > max { L a , L b } − λ 0$⁠, in other words that the inventories be low enough to allow the initial perturbation to propagate back to agent 0 in both branches. Under these conditions, the network collapses.  Appendix B provides the details of the proof, along with relevant examples of time series.

For two loops as in Figs. 1(d) and 1(e), Fig. 2 shows how the asymptotic density $ρ ∞$ of fully active agents varies with the total delay T_b of one of the two loops. The phenomenon of fine detuning occurs in the blue areas, on either side of T_a, i.e., for T_b > T_a, as in the inequalities (11), as well as for T_b < T_a, which also corresponds to these inequalities when the loops labeled a and b are interchanged.

We notice the appearance of periodic solutions on the less detuned side of a collapsed interval, i.e., closer to the “resonance” between the two “periods,” T_a and T_b. On the farther side, there is an abrupt change from collapse to recovery, cf. Figure 13 in  Appendix B.

This detuning also occurs, to a lesser extent, when one of the total delays, T_a or T_b, is close to an integer multiple of the other, as shown by the narrower green areas of Fig. 2: when $( T a = 3 , T b = 6 )$ and when $( T a = 3 , T b = 1.5 )$⁠. In the latter case, given the parameter set chosen, fine detuning only occurs for $T b > T a / 2$⁠. The red area of Fig. 2 corresponds to the phenomenon of rapid collapse described in Sec. III B: delays and inventories are so short that the propagating waves return to agent 0, while it has not yet recovered from the initial perturbation.

When more than two loops interact, the process of wave duplication and recombination can lead to even more complex dynamical behavior. In particular, given three interacting loops, as illustrated in Figs. 1(f) and 1(g), we observe that periodic regimes occur for continuous parameter intervals, while for two interacting loops such regimes only set in at isolated parameter values, as indicated by the vertical dotted lines in Fig. 2.

The dynamical response to initial perturbations of a specific network is shown in Fig. 3 for six sets of parameters. This network has n = 4 agents, structured into three loops connected through agent 0, as in Fig. 1(f). As illustrated in panels (a) and (b), permuting delays within a loop modifies the finite-length transient but does not alter the asymptotic periodic solution. A minor change in one delay can lead to a discontinuous jump in period length: from 1 for $τ 02 = 2$ in panels (a) and (b) to 2 for $τ 02 = 1.9$ in panel (c). Such period doubling is well known in differentiable dynamical systems and has been discussed in the BDE context by Saunders and Ghil.⁴⁹ Changing inventories does not modify the period length but it does affect the actual periodic pattern; see panels (a), (d), (e) and (f).

In this network, additional numerical simulations (not shown) suggest that, in addition to full activity or collapse, periodic solutions occur in large regions of the 10-dimensional parameter space $( τ 01 , τ 10 , τ 02 , τ 20 , τ 03 , τ 30 , λ 0 , λ 1 , λ 2 , λ 3 )$⁠. These regions have a finite volume in $ℝ 10$⁠, in other words they are not points as in the case of two interacting loops, but are a union of several disconnected sets. Within such a set, the period length can remain stable or change discontinuously with a change in one delay. Additional results are shown in  Appendix C and require further study.

Throughout Sec. III, we have identified the qualitative properties of the dynamics in simple network motifs, which are the building blocks of larger, complex and more realistic structures of economic production. These properties are summarized in Table II. In Sec. IV, we will show how the dynamics of these motifs helps us understand damage propagation in larger networks.

As shown in Sec. III, the presence of loops in a network, and in particular, that of interacting loops, enable permanent disruption. In complex networks, loops are the constitutive motif of their so-called strongly connected component (SCC).^23,60 An SCC is a maximal subset of agents within which each agent can reach and be reached by any other agent through directed links. SCCs are embedded in a weakly connected component (WCC), where the WCC is the maximal subset of agents within which each agent can reach and be reached be any other agent, disregarding the direction of the links. Fujiwara and Aoyama²⁵ found that 99% of the 1019 854 firms in their dataset were present in the WCC and 45.3% in the SCC.

If a permanent disruption reaches an SCC, all agents in the SCC itself, along with those that source their inputs in the SCC, get permanently disrupted. The whole disrupted region is called the out-component of the SCC, labeled OUT. Identifying all agents in OUT allows one to predict the potential extent of permanent disruptions. Conversely, an initial perturbation can only trigger permanent disruptions if it is located in the IN of an SCC, i.e., among the agents producing inputs used by any agents within the SCC. Investigating all IN agents allows one, therefore, to estimate the maximum probability that an initial perturbation triggers a permanent disruption. The major WCC, SCC, IN, and OUT components of a network can be usefully mapped into a so-called bow-tie diagram; see Fig. 13 in Coluzzi et al.²³ and discussion therein.

Although ER networks fail to capture many features of real production networks, they are the most studied types of complex networks, and are therefore a useful basis of comparison for our analysis of the more realistic SF–FA networks. In ER networks, a giant SCC arises when the average connectivity is larger than 1; see Fig. 4(a). The expected sizes of the IN and OUT components of the giant SCC are equal, and reach half of the network for c = 1.39 when $n = 10 5$⁠. At this specific connectivity, and in the absence of inventory, there is one chance in two that a random perturbation triggers a permanent disruption. Such a disruption would then affect half of the agents. In a scale-free network, the high concentration of links ensures that a giant SCC emerges for very low connectivities; this is the case also in our SF–FA networks, as seen in Fig. 4(b). The asymmetry between in and out-degree distributions in SF–FA networks leads, moreover, to different sizes of the IN and OUT components.

According to the results of Sec. III, the asymptotic response is determined by the presence of interactions between multiple short loops. We investigate next the density of agents that are engaged in at least two loops of five agents or less, as a function of network structure; this quantity is denoted by $D SIL$⁠, where the subscript SIL stands for “short interacting loops.” Longer loops are expected to play only a minor role when the average inventory exceeds 0.2.

For both ER and SF–FA networks, $D SIL$ increases with the average connectivity (not shown), but tends to decreases with the network size, as shown in Fig. 4(c). This decrease is, however, much sharper for ER networks; hence, for large n, almost no agent is engaged in multiple short loops. This result is consistent with the findings of Bianconi, Gulbahce, and Motter,⁶¹ who showed that the number of short loops becomes negligible in ER networks in the limit of large n.

For SF–FA networks, though, the density of agents that are engaged in at least two short loops, $D SIL$⁠, stabilizes for large n and remains above 50%. Following Bianconi, Gulbahce, and Motter,⁶¹ this difference is likely to be related to the positive correlation in SF–FA networks between in-degree and out-degree.²⁵ 

The high occurrence of short interacting loops in SF–FA networks suggests that they are more prone to permanent disruptions than ER networks. Moreover, for connectivity below 1, SF–FA networks still exhibit SCCs, and can therefore permanently collapse, while ER networks cannot. For larger connectivities, SF-FA networks have a smaller giant OUT component than ER networks, suggesting that the extent of a possible collapse, in terms of the number of inactive nodes, is expected to be lower.

In both ER and SF–FA networks, three asymptotic regimes are possible: full activity, periodic disruption or permanent disruption. In the homogeneous case with τ  =  1 and $λ = 0.2$⁠, both ER and SF–FA networks fully recover from any initial perturbation. However, even a moderate level of heterogeneity in delays reduces resilience fairly abruptly; see Fig. 5(a). This effect is due to the occurrence of fine detuning, cf. Sec. III C. Slight variability in the delays of multiple interacting loops can induce the collapse of these loops, which then extends to the whole giant SCC and OUT components.

For ER networks, the transition between full recovery and collapse occurs for a standard deviation of delays $σ τ$ between 0.06 and 0.12; see Fig. 5(a). The rapid drop of $R ∞$ is due, in fact, to a change in the structure of the distribution of the individual values of ${ ρ ∞ ( k ) }$⁠. Before the transition, the network always recovers and all $ρ ∞ ( k )$ equal 1. During the transition, a second mode—in the statistical sense of the term—appears in the distribution. It corresponds to the collapse of the giant OUT component, and it is equal to 1 minus the relative size of the giant OUT component; for ER networks with n = 1000 and c = 4, this mode equals $1 − ( 0.98 ± 0.008 )$⁠. The weight of this mode in the statistical distribution of ${ ρ ∞ ( k ) }$ starts from 0% and reaches $82.5 % ± 3.0 %$ after the transition.

In SF–FA networks, the transition is likewise shown in Fig. 5(a) and it occurs for even lower heterogeneity levels, because of the higher density of short interacting loops; see the values of $D SIL$ in Fig. 4(c) at $n = 10 3$⁠. These networks have smaller giant OUT components than ER networks, namely 0.90 ± 0.02, but permanent disruption occurs more often, following $88.2 % ± 0.9 %$ of the initial perturbations. Furthermore, the $σ τ$-value at which transition occurs is more variable among different SF–FA than for ER networks; see the blue vs. the black interquartile ranges in Fig. 5(a).

As $σ τ$ increases, no persistent periodic orbits arise, in either ER or SF–FA networks. Transients can, however, exhibit oscillations for moderate delay heterogeneity, $σ τ ≲ 0.15$⁠, as illustrated by the time series of Figs. 5(b)–5(d). The period of these small oscillations corresponds to the mean delay $μ τ$⁠, and they quickly vanish for higher $σ τ$-values. Delay heterogeneity tends to accelerate the disruption waves and shorten the transients. After the transition from full recovery to collapse, a further increase in $σ τ$ reduces the average time of first disruption Θ. In SF-FA networks, it is reduced by about 11% between $σ τ = 0.1$ and $σ τ = 0.5$⁠, from 4.4 ± 0.09 to 3.9 ± 0.09.

Inventory heterogeneity also deteriorates the ability of the system to recover, i.e., its resilience. As seen in Fig. 5(e), a similar transition occurs as $σ λ$ increases; this transition, however, occurs for much larger $σ λ$-values. For instance, for ER networks, the decrease in $R ∞$ starts when the standard deviation of inventories $σ λ$ reaches 0.15, i.e., for a coefficient of variation of 75%, as opposed to a coefficient of variation of 6% for delay heterogeneity.

In this case, though, the transition between full recovery and collapse is accompanied by the emergence of asymptotically stable periodic orbits: these are represented by the dashed lines in Fig. 5(e), and illustrated by the time series of Fig. 5(g). In addition to full recovery and collapse of the giant OUT component, the distribution of $ρ ∞ ( k )$ has intermediate values, which correspond to the time average of $ρ ( t )$ during one period of these oscillations. The presence of these periodic orbits tends to soften the transition from high to low $〈 ρ ∞ 〉$-values, which is less abrupt than in Fig. 5(a).

As for the delays, the transition in SF-FA networks occurs for lower $σ λ$⁠. The emergence of asymptotically periodic orbits is related to the interaction between multiple short loops, as illustrated in Sec. III D, and it is largely fostered by the homogeneity of the delays. As soon as delays become heterogeneous, fine detuning occurs and most oscillations become transient and lead to collapse. In this case, the counterpart (not shown) of the diagram in Fig. 5(e) does not exhibit periodic orbits.

The detrimental role of heterogeneity on resilience can be observed by comparing panels (a) and (c), in which delays and inventories are homogeneous, with panels (b) and (d) of Fig. 6, in which they are heterogeneous with a coefficient of variation of 50%. For both ER and SF–FA networks, the portion of the parameter plane $( μ τ , μ λ )$ corresponding to full recovery (dark blue) is significantly reduced. As previously illustrated, this detrimental effect of heterogeneity on resilience mostly comes from the heterogeneity of delays. When delays are heterogeneous, the exact mean value of the distribution has almost no impact on resilience; see Figs. 6(b) and 6(d). This result provides additional evidence for permanent disruption occurring through fine detuning, in which only the difference in delays matters and not their absolute value, in agreement with inequalities (11).

In the homogeneous case, collapse occurs for short delays only. Since all delays are equal, fine detuning can only occur between loops of distinct length, e.g., between a 2-loop and a 3-loop, or between a 2-loop and a 4-loop. The smallest difference in delays between such interacting loops—i.e., the quantity $Δ T$ in Eq. (11)—is therefore τ: that is why small τ-values are necessary for fine detuning. For instance, using Eq. (11), if the initial perturbation hit an agent connecting a 2-loop and a 3-loop, fine detuning occurs for $λ ≤ τ < 1 − λ$⁠.

Small delays also enable fast collapse when the sum of delays and inventories in one loop falls below the initial perturbation, see Sec. III B. In loops of two agents, for instance, fast collapse occurs when $τ < δ 0 / 2 − λ$⁠. Further mechanisms are also likely to induce collapse, such as those underlying the behavior of three interacting loops, as illustrated in Sec. III D.

Comparing Figs. 6(b) and 6(d) shows that SF–FA networks are less resilient than ER networks. The transition between full recovery—the dark region—and permanent disruption—the light—occurs for larger inventory levels. According to the results of Sec. IV A, SF-FA networks have a larger number of short interacting loops, and thus a higher probability for fine detuning.

In the homogeneous case, Figs. 6(a) and 6(c) are very similar. However, because of the same topological feature, the parameter region of full activity is slightly reduced for SF–FA networks in the general area $( 0 ≤ τ ≤ 0.5 , 0.25 ≤ λ ≤ 0.45 )$ of low delays and relatively high inventories. On the other hand, because SF–FA has a smaller giant OUT component, the extent of collapse is lower than in ER networks. This can be seen in Figs. 6(b) and 6(d), where the light region is slightly darker than in Figs. 6(a) and 6(c), indicating larger $R ∞$-values.

Instead of considering a whole production network at once, it may be useful to first investigate the supply chain of a specific industry taken in isolation, such as the ones of the transportation or energy sector, and then investigate their interactions. A sectorial supply chain can be represented by acyclic networks. In such a topology, no permanent disruption occurs. Still, heterogeneities in delays and inventories do have an impact on certain features of the transient regime. An RA network can be seen as a sequence of linear chains that intersect without forming a loop. Recalling the study of linear networks in Sec. III A, the amplitudes of the propagating waves induced by the perturbation of primary producers are likely to increase as inventory heterogeneity increases.

This prediction is confirmed in Fig. 7(a), which shows that the transfer coefficient TC increases with $σ λ$⁠. This indicator measures the average disruption time of final producers following the external perturbation of a primary producer. In other words, a broader distribution of inventories makes an RA network more permeable, i.e., final producers are more likely to be disrupted due to the failure of primary producers. Since more agents get perturbed, full recovery takes longer.

Figure 7(b) shows a positive, yet moderate, impact of delay heterogeneity on TC. This phenomenon is due to a process similar to fine detuning, which occurs between parallel pathways. The initial perturbation generates several propagating waves, some of which get recombined further down the network. This recombination can potentially amplify the disruption, as shown in Sec. III C for fine detuning between loops. This explanation is confirmed by the fact that TC increases only for small $σ τ$-values; this amplifying effect then vanishes, as in the transition for cyclic networks illustrated in Fig. 5(a), where fine detuning only requires a small degree of delay heterogeneity.

The interdependencies between two sectorial supply chains, or two types of infrastructures, can be represented by linking RA networks. In a system of two interdependent RA networks—where each primary producer of one network is supplied by one final producer of the—heterogeneity in the delays significantly reduces resilience. When the delays in both networks are drawn from the same lognormal distribution, the behavior (not shown) is qualitatively similar to that of the ER and SF–FA networks in Fig. 5(a), with a quick transition to collapse for moderate $σ τ$⁠.

Two time series are displayed in Fig. 8, where the red curves correspond to the network where the initial perturbation occurs. When the delays are homogeneous or almost homogeneous, as in Fig. 8(a), the propagating waves fade out, while in Fig. 8(b) they amplify through fine detuning due to delay heterogeneity. Inventory heterogeneity can also lead to collapse, but for even larger $σ λ$-values than in ER and SF–FA networks, while periodic orbits arise to a smaller extent during the transition (not shown).

In all propagation phenomena in these paired networks, whether leading to recovery or to collapse, we observe a back-and-forth pendulum movement: disruptions are alternately stronger in one network, then in the other, and so on. In Figs. 8(a) and 8(b), the red curve starts below the blue curve, then crosses above it, and so on, until the asymptotically stationary regime is reached. Additional simulations show that the time between two crossings lasts between 3.5 and 5.5 time steps, which is slightly higher the total delays along an average pathway from the primary to the final producers. We also observe a slight acceleration of these crossing throughout the transients leading to collapse.

When the delay distributions are distinct from one network to the other, collapse can occur even for distributions that are both homogeneous. In the so-called dichotomous case, where all delays are equal to τ_a in network a and to τ_b in network b, a phenomenon of resonance occurs due to fine detuning, as shown in Fig. 9. When τ_a or τ_b is close to being an integer of the other, both networks recover; otherwise, they collapse. This process involves longer loops formed by the interconnection between the two networks, and thus occurs for lower inventories— $λ = 0.1$ in Fig. 9. Rebounds are observed for intermediate multiplicity, for $τ b ≈ 1.5 or 2.6$⁠. This observation suggests finer mechanisms of interactions between propagation waves that could be examined in further work.

An important aspect of resilience is the maximum intensity of perturbations that the system can handle without leaving the attractor basin of full activity; this question is related to the original inquiry of Albert, Jeong, and Barabási²⁷ into the “tolerance” of complex networks. In our Boolean setting, the intensity of an external perturbation is apprehended by its duration.

For SF–FA and ER networks, we measure for each agent the critical duration of external perturbation beyond which the giant OUT component collapses. For all agents $i = 1 , .. , n$⁠, we denote by $δ c ( i )$ this duration, also called the critical perturbation. When the relative size of the OUT component equals R_O agents, we define the critical perturbation as follows:

The values of $δ c ( i )$ were determined here by numerical network simulations. External perturbations lasting less than $δ c ( i )$ are called subcritical perturbations: they do not threaten the stability of the network.

Agents that are not in the giant IN component cannot induce permanent disruptions: they have no critical perturbation, i.e., $δ c ( i ) = ∞$ for all agents $i ∉ IN$⁠. All other agents in the giant IN component have a finite $δ c ( i )$-value, which depends on the topology, the delays, and on the inventories of the other agents. The lower $δ c ( i )$⁠, the more important agent i is for global resilience. Figure 10(a) shows an example of how critical durations are spread over an SF–FA network with n = 50 agents. The spread of $δ c ( i )$-values is relatively limited.

This indicator is now compared with the canonical centrality measures developed in the network theory.⁶⁰ We determine whether central nodes, according to closeness, betweenness, and eigenvector centrality, have also low critical perturbations $δ c ( i )$⁠. Additional simulations are performed for five ER and five. SF–FA networks with lognormally distributed delays with parameters $( μ τ = 1 , σ τ = 0.5 )$ and inventories with $( μ λ = 0.2 , σ λ = 0.1 )$⁠.

We observe that, on average, the ensemble ${ δ c ( i ) : i ∈ IN }$ is satisfactorily correlated with closeness centrality: $− 0.72 ± 0.066$ for ER networks and $− 0.68 ± 0.031$ for SF–FA networks. Closeness centrality is the inverse of the average geodesic distance linking one node to all other nodes, and therefore captures the potential damages that the failure of this node could cause. The correlation coefficient is not higher, in absolute terms, when the geodesic distance is replaced by the propagation time, i.e., by the sum of delays and inventories.

The remaining variability that is not captured by closeness centrality is due to the particular propagation phenomena studied in Sec. III, especially fine detuning. These mechanisms, enabled by the interaction between short loops, occur more often in SF–FA networks. This behavior is a plausible explanation for the correlation coefficient being lower for SF–FA networks than for ER ones.

Correlations are weaker for betweenness centrality: $− 0.36 ± 0.06$ for ER networks and $− 0.37 ± 0.04 %$ for SF–FA networks; they are also weaker for eigenvector centrality: $− 0.01 % ± 0.06 %$ for ER networks and $− 0.36 % ± 0.04 %$ for SF–FA networks. Betweeness centrality, which measures how important a node is in connecting each pair of nodes, captures wider information on the structure of the incoming links, which are not relevant for propagation. Eigenvector centrality, based on in and out-degree concentration, has a particularly low correlation coefficient for ER networks, where degrees are more evenly spread.

The ensemble ${ δ c ( i ) }$ condenses information on one of the three aspects of resilience. We study how this ensemble varies with a change in the topology and the distribution of delays and inventories. In particular, the sum $Δ c$ of the $δ c ( i )$'s over the complement of the IN set can be used as a scalar indicator for resilience: the higher this indicator, the larger perturbations a network can sustain without collapsing. The following results present quantities that are averaged over five repetitions of the same class of networks, and confirm the findings in Sec. IV concerning the attractor basins and propagation speed.

First, the $Δ c$-values of SF–FA networks are lower by almost 50% than those of ER networks with similar delay and inventory distributions. In networks of 1000 agents with heterogeneous delay and inventory distribution, $Δ c = 817 ± 34$ in the ER networks and $Δ c = 440 ± 37$ in the SF–FA networks.

To illustrate this result, two distributions of $δ c ( i )$'s are shown in Fig. 10(b), one for SF–FA networks, the other for ER networks. First, the shape of the two histograms is similar and their widths are approximately equal. Second, as expected from the results of Sec. IV, delay heterogeneity tends to reduce $Δ c$⁠. For ER networks with $σ τ = 0.05$ (not shown), i.e., before the sharp drop in $R ∞$ apparent in Fig. 5(a), $Δ c = 910 ± 51$⁠, while it drops to 817 ± 34 in Fig. 10(b), where $σ τ = 0.5$⁠.

The networks under consideration here recover from any initial subcritical perturbation. However, the critical perturbations $δ c ( i )$ were determined for initial perturbations occurring one by one only. They apply to situations where, after each perturbation, the system has enough time to recover. What happens when additional subcritical perturbations occur during the recovery process?

Additional simulations were carried out, where, instead of imposing a unique perturbation at the beginning of the run, random perturbations occur. Specifically, agents have a certain probability to get perturbed per unit of time; this failure rate is denoted by f. The duration of theses subcritical perturbations is set to $0.8 δ c ( i )$ for agents that are in the giant IN component and to 1 for agents that are not in it. The following results are based on a limited set of numerical simulations, and should therefore be considered as merely indicative.

SF–FA networks with 1000 agents and lognormally distributed delays and inventories—where $( μ τ = 1 , σ τ = 0.5 )$ and $( μ λ = 0.2 , σ λ = 0.1 )$—are robust when f is lower than 0.3%, i.e., when 3 subcritical perturbations occur on average per unit of time. For higher rates, collapse occurs within 500 units of time; this latency time interval decreases dramatically as f increases. For smaller networks, collapse occurs for even lower f-values. For instance, in SF–FA networks with 100 agents, collapse often occurs within 500 units of time for f as low as 0.03%. The two situations are illustrated in Figs. 11(a) and 11(b).

We observe that, at such low rates, collapse occurs following two perturbations that are very close to each other in time, within one or two units of time. It suggests that the combined perturbations, even subcritical, of interacting agents could lead to collapse. The different behavior observed in large and small networks could originate from the higher density of short interacting loops in smaller networks, as shown in Sec. IV A. The probability that two perturbations concomitantly affect the same loops or two loops in interaction is therefore higher.

These analyses show that large networks are robust to random subcritical perturbations occurring at a moderate rate—less than 3 per time steps for n = 1000. However, if perturbations are not random and exhibit some spatial correlations, even lower failure rates can induce collapse. A combination of localized perturbations occurring in specific regions of the network, whether ER or SF–FA, can have significant economy-wide impacts.

The point in time at which dynamical processes occur on a network's nodes does matter. In this paper, we have demonstrated that heterogeneity in delay values significantly influences how economic networks react to small, localized perturbations. Several aspects of a network's resilience are affected by the heterogeneity of its delays: the size of the economically desirable basin of attraction is reduced, the intensities of the perturbations from which the network can recover are smaller, and the speeds at which perturbations propagate are higher. We argue that the distribution of delays is crucial in influencing the dynamics of input–output networks, a factor well worth considering when building or analyzing such economic models.

In particular, the destabilizing effect of short delays in production and delivery is rather counterintuitive. Unlike feedback control systems, in which short delays prevent oscillations and instabilities, quick interactions between interdependent processes speed up the propagation of a disruption, resulting in a longer recovery or a collapse.

We found that inventory heterogeneity can also play a major role in acyclic supply chains, by increasing the permeability of the network, i.e., the probability that final producers get affected by the failures among primary producers. This finding confirms the conclusions of Hallegatte,³⁷ who found an adverse effect of heterogeneity in either capacity or inventories upon a network's recovery from natural disasters. Imbalances in inventories, though, have a weaker effect on a network's vulnerability than the one induced by delay heterogeneity.

We identified the key mechanisms underpinning these results. Of particular significance is the process of fine detuning, whereby a single perturbation induces multiple disturbance waves in a network with interacting loops; these waves subsequently recombine and form new, longer perturbations. The latter originate from interactions between loops with short and slightly distinct time delays. The aggregate behavior of large, complex networks was thus found to be related to their local topological features. This finding allowed us to explain why real production-and-supply networks, represented here as SF–FA networks, are consistently less resilient than comparable ER networks. The higher vulnerability of scale-free networks is consistent with the conclusions drawn from epidemic-spreading models,^32,33 a consistency that suggests conceptual linkages between contagion dynamics and disruption propagation. Our findings on SF–FA networks could be refined by incorporating other features of real economic networks identified by Fujiwara and Aoyama,²⁵ in particular disassortativity, clustering, and community structures.

The process of fine detuning was also observed in networks of networks, and it is induced by heterogeneity in delays between networks. An acyclic network with homogeneous delays can safely interact with other acyclic structures, as long as their respective delays are integer multiples of each other. But the absence of integer multiplicity between delays leads to a loss of synchronization between the disruption waves, and may thus prevent the joint recovery of the networks Buldyrev et al.⁶² have shown, in effect, that interacting networks can exhibit very different failure behavior than single networks, and our findings here contribute therewith to the study of networks of networks.

Finally, while aiming mainly at gaining insights into their resilience, we explored many aspects of a network's behavior, such as the size of its attractor basins, as well as the transient speed of damage propagation and the effects of the intensity of exogenous perturbations. Studying these various aspects separately, as well as in combination with each other, allowed us to increase the robustness of our conclusions.

We argue that, although sometimes it may be presented as ill-defined, the concept of resilience can usefully guide system analyses and stimulate original investigations.

In Sec. VI, we evaluated how critical individual agents are for overall resilience. We built an agent-level indicator $δ c ( i )$ based on the maximum perturbation intensity—measured here by its duration—that an agent i can handle without causing permanent disruptions to a large fraction of the network. The importance of individual agents can be partly predicted by closeness centrality, but not by two additional, commonly used centrality measures that we considered—betweenness and eigenvector centrality.

It appears, therewith, that no one-size-fits-all indicator can measure the importance of nodes, independent of the context. Hence, vulnerability indicators need to be specifically tailored to capture the relevant dynamical processes at stake. In a real-world context, to evaluate how critical a specific firm is using the concept of a critical perturbation, one would need to quantify the levels of inventory of its direct partners, and the extent to which they depend on each other, as well as on the travel times of the damages.

How powerful can such vulnerability indicators be in mitigating risks? Section VI C shows that our indicator $δ c ( i )$ can be used to map the vulnerability of the production networks to localized perturbations. From this information, critical nodes could be targeted with measures that speed up their post-shock recovery. Such a strategy could help alleviate the system-wide consequences of events that are uncorrelated, both in time and in space.

Our study suggests, however, that vulnerability assessments based on single-failure scenarios do not capture the whole risk landscape: they have serious limitations, for instance, when it comes to large-scale events, such as natural disasters, or to coordinated attacks, which can lead to concomitant failures of specific combinations of agents. Comprehensively assessing the individual role of firms for such low-probability high-impact risks is particularly challenging.

We argue that progress can be made by identifying critical network structures, as done in Sec. III. We found that it is the interactions between basic structures, more than the structures themselves taken in isolation, that give rise to vulnerabilities. For large, complex networks, it is likely that no single indicator can cover all types of damage mechanisms. From a regulatory perspective, then, several vulnerability or resilience indicators should be combined in applications, without underestimating the remaining uncertainties.

Besides our results on network resilience, this work has revealed some intriguing aspects of Boolean processes with delays. Rosin et al.^63,64 observed somewhat similar characteristics of such processes in experimental realizations of BDEs using electronic logic circuits.

In large regions of parameter space, we found relaxation oscillations emerging in the density $ρ ( t )$ of undamaged agents. This was the case, in particular, for networks consisting of three interacting loops, as apparent in the time series of Figs. 3 and 14, as well as for complex networks with homogeneous delays and heterogeneous inventories, cf. Figs. 5(f)–5(h) and Fig. 15 in  Appendix C below. These oscillations asymptote, after a transient, to a periodic limit cycle, over whose period the initial perturbation propagates and synchronously affects a large number of nodes.

Within a period, inventories first attenuate the beginning of the disruption, but are not large enough to stop it completely. Next, they are gradually emptied out, thus leading to the reappearance of disruptions. Finally, inventories get replenished, the disruptions temporarily disappear again entirely, and the cycle can start all over.

Since delays are homogeneous, the replenishment of the inventories occurs simultaneously. This particular timing leads to the periodic but sudden disappearance of almost all disruptions, leading in turn to a very steep upward slope of $ρ ( t )$ within each period; see, for instance, Fig. 5(g). The pattern of the relaxation dynamics, i.e., the downward slope of $ρ ( t )$⁠, is directly related to the fact that inventories become used gradually, one agent after the other.

We observe in Fig. 14 that modifying one inventory changes the relaxation pattern. The higher the average level of inventories, the longer the relaxation phase within each period.

Delays determine the return of the inventories, and therefore the timing of the steep upward peak in $ρ ( t )$⁠. They therefore determine the periodicity. As seen in Sec. III, these oscillations are due to interactions between multiple loops. A change of delay in one loop can desynchronise different loops, and lead to a different period. A small change in delay can have a strong discontinuous effect on the length of the period; see Figs. 14(a)–14(c). However, delays that are too heterogeneous tend to desynchronize the loops, which can lead either to recovery or to complete collapse. These phenomena are further discussed in  Appendix C.

The patterns seen in Fig. 15 below suggest that metastable oscillations can stabilize and persist indefinitely. The fact that, within a given network, the same propagation phenomena can lead at first to a relaxation pattern, that lasts for several periods, and then to a different pattern, with sometimes a distinct periodicity, suggests that relaxation oscillations may temporarily set in over different regions of the network. The transition from one oscillatory regime to another is marked by a temporary desynchronization, and can span multiple periods; see Figs. 15(b)–15(f). This behavior gives rise to increasingly long transients. D'Huys et al.⁶⁵ found experimentally a different phenomenon that also generates very long transients in Boolean loops with “NOT” operators.

The fact that oscillations can move throughout the network is of particular interest for networks with several centers, such as the system of interdependent networks studied in Sec. V. This movement is likely to be at the origin of the pendulum alternation between two acyclic networks observed in Fig. 8.

In this paper, we have formulated a simple model of an economic network, and analyzed how heterogeneity in delays, inventories and connectivity negatively affect resilience to small perturbations. We identified the key underpinning mechanisms leading to different dynamical responses, and established robust correspondences between the local structure of the networks and their aggregate behavior. Empirically calibrated networks were used to represent nationwide production systems. Our analysis was then extended to networks of networks, which may be used to study the interactions between sectorial supply chains or infrastructures. We investigated resilience using three major network characteristics: size of the attractor basins, duration of the transients, and intensity of the perturbations. These three characteristics allowed us to compare and combine our partial results throughout this paper and provide more robust final findings.

The use of BDEs allowed us to capture the highly nonlinear processes at work without modeling the exact production levels. Thanks to this key simplification, we were able to perform extensive numerical simulations on complex networks and study large portions of the models' parameter space. Even if the quantitative results herein cannot be tested directly on real data, the qualitative understanding gained allows us to identify new directions for empirical studies.

We argue that variability in production and delivery delays and inventory unbalances could amplify the economy-wide impacts of localized supply disruptions. This hypothesis is empirically verifiable using econometric methods. We also suggest that, beyond the simple message that interdependencies between key industries increase vulnerability, it is the time scale over which the connections actually happen that matters, as well as the existence and nature of a given connection. Case studies could reveal whether, indeed, short interactions between interdependent industries can lead to sizeable disruption risks.

Albeit restricted to Boolean variables, it was possible to satisfactorily model inventories. The BDE framework seems therefore flexible enough to integrate more detailed representations of economic processes. Our model could thus incorporate, for instance, a certain degree of substitutability between inputs, known to alleviate propagation phenomena. Such a modification could be achieved by replacing the operator AND by OR, or by introducing other Boolean operators. Agents could also be allowed, when disruptions become severe, to switch to alternative suppliers within the same network, or, as implemented by Coluzzi et al.,²³ to receive external provisions of inputs.

The authors thank F. D'Andrea and G. Weisbuch for the numerous suggestions and G. Weisbuch for reading attentively and providing constructive comments on the earlier version of the paper. This work was supported by a Monge Fellowship of the Ecole Polytechnique (CC), as well as by Grant Nos. N00014–12-1–0911 and N00014–16-1–2073 from the Multidisciplinary University Research Initiative (MURI) of the Office of Naval Research, and by the Energy and Prosperity Chair of the Institut Louis Bachelier (MG).

We consider a linear network of $n ≫ 1$ agents, where inventories and delays are random variables, respectively, denoted by $Λ 0 , … , Λ n − 1$ and $T 01 , … , T 1 − 2 n − 1$⁠. Each of these variables is independent, and follows a certain probability distribution. The duration of the disruption affecting each agent is also a random variable, denoted by $Δ 0 , … , Δ n − 1$⁠. We further denote by f_X the probability distribution function of a random variable X, and by g_X its cumulative distribution function. The duration of the initial perturbation, $Δ 0$⁠, is exogenous, and its probability distribution function is given, along with those of each inventory and delay. We denote by Ω_k with $1 ≤ k < n$ the occurrence of a disruption at agent k, following the initial perturbation of agent 0. We add to this set the event $Ω 0$⁠, which represents the occurrence of the exogenous perturbation on agent 0. The probability of this event, $Pr ( Ω 0 )$⁠, is 1.

The probability that the perturbation affecting agent 0 propagates to agent 1 corresponds to the probability that the inventory of agent 1 is lower than the initial perturbation. This probability, denoted by $Pr ( Ω 1 | Ω 0 )$⁠, and called the transmission probability, is obtained through the following formula:

We assume that the perturbation reaches agent 1. To determine whether it propagates further and reaches agent 2, we need to gauge the relative value of $Δ 1$ and $Λ 2$⁠. To that end, we need to determine the probability distribution of $Δ 1 = Δ 0 − Λ 1$ conditional to the occurrence of $Ω 1$⁠. It is the probability distribution for $Δ 1$ for positive values of $Δ 1$⁠. We can write

We can then use the convolution product only to determine the probability distribution of $Δ 1 | Ω 1$

We can use those results to initiate an iterative algorithm to compute the duration of all perturbations Δ_k conditional to their occurrence Ω_k. The following set of equations allows one to compute the probability distribution of $( Δ k | Ω k )$ given the probability distribution of $( Δ k − 1 | Ω k − 1 )$ and the a priori probability distribution chosen for inventories

This algorithm yields as a by-product all transmission probabilities $Pr ( Ω k | Ω k − 1 )$ for k > 0. We compound these probabilities to obtain the probability that the k^th agent gets disrupted

The expected amplitude of the propagation is simply the sum of all probabilities

The discrete probability distribution of A can also be obtained. The probability that A equals an integer a is the probability that the initial perturbation propagates up to agent a without propagating further

For trees, we make the simplifying assumption that each pathway connecting the root agent to the final agents is independent. In a tree composed of L layers, with each agent supplying to d agents, the expected amplitude is $E ( A ) = ∑ k = 0 K d k Pr ( Ω k )$⁠. Figure 12 shows how the average amplitude changes with the mean and standard deviation of inventories for different probability distribution, using numerical simulations and the algorithm described in this  Appendix.

Suppose that the duration of the initial perturbation is a random variable which follows an exponential distribution probability with rate parameter r, and that all inventories follow the same exponential distribution with rate parameter q. The probability distribution function and cumulative distribution function of inventories are

Using Eqs. (A6) and (A7), we have

In other words, if agent 1 gets disrupted, the probability distribution of the duration of this disruption is the same as the probability distribution of the initial perturbation. This result applies to all agents. Consequently, for all $k = 1 , .. , n$⁠, the transmission probability is constant and equal to $q / ( r + q )$⁠.

Fine detuning can occur in a network composed of two interacting loops, whose topology is illustrated in Figs. 1(d) and 1(e). Due to this phenomenon, an initial perturbation can amplify and lead to the full and permanent disruption of the network. In this  Appendix, we determine the conditions that lead to such small detuning. We denote by T_a and T_b as the sum of the delays in each loop and assume, without loss of generality, that $T a ≥ T b$⁠. Similarly, we denote by L_a and L_b the sum of inventories in each loop, including the inventory of the pivotal agent, called agent 0, whose inventory equals λ₀. We assume that L_a and L_b are strictly positive.

When agent 0 gets perturbed at t = 0, two waves of disruption emerge, one in each loop. If we consider each loop as independent, then Eq. (7) states that agent 0 gets disrupted through loop a. if the duration of the initial perturbation exceeds L_a. We call this event a direct disruption (a). Conversely, a direct disruption (b) occurs if the duration of the initial perturbation exceeds L_b. If either one of these two conditions is not met, then the study of the propagation in this system boils down to the one in isolated loops, where the propagating wave eventually vanishes.

If both propagating waves, though, reach agent 0 again, they can interact. We denote by $a ′$ the supplier of agent 0 belonging to loop a, and by $b ′$ its supplier from loop b. The dynamical equation governing the behavior $x 0 ( t )$ of agent 0 is

where $τ a ′ 0$ and $τ b ′ 0$ are the delays between agent $a ′$ and agent 0, and between agent $b ′$ and agent 0, respectively.

Direct disruption (a) corresponds to the case in which the inventory is emptied due the disruption of supplier $a ′$⁠, and the production is stopped due to the same disruption, and likewise for direct disruption (b). But another disruption can also occur when the inventory of agent 0 is emptied by the disruption of one supplier, and the production is then stopped by the perturbation of another one. Specifically, agent 0 is disrupted at t if agent $a ′$ is disrupted at $t − τ a ′ 0 − λ 0$ and agent b at $t − τ b ′ 0$⁠. We call this phenomenon a (ab) cross-disruption. Conversely, a (ba) cross-disruption occurs when agent 0 is disrupted at t because agent $b ′$ is disrupted at $t − τ b ′ 0 − λ 0$ and agent $a ′$ at $t − τ a ′ 0$⁠.

To study these cross-disruptions, we need to know when the propagation waves hit suppliers $a ′$ and $b ′$⁠. From Eq. (8), supplier $a ′$ is perturbed at $( T a − τ a ′ 0 ) + ( L a − λ 0 )$ until $( T a − τ a ′ 0 ) + δ 0$⁠, and supplier b at $( T b − τ b ′ 0 ) + ( L b − λ 0 )$ until $( T b − τ b ′ 0 ) + δ 0$⁠. The conditions for a (ab) cross-disruption to occur are therefore

They can be rewritten as

Both conditions can be realized at the same time if

Since T_a < T_b and $L a < δ 0$⁠, the first condition is always met. The second condition can be rewritten as: $Δ T < δ 0 − L b + 2 λ 0$⁠.

The same analysis for a (ba) cross-disruption, yields that the first necessary condition for such a disruption to occur is $T b + L b ≤ T a + δ 0$⁠; under the present assumptions, the latter condition is not met. Hence, if T_a < T_b, only a (ab) cross-disruption is possible.

The starting and ending times of the cross-disruption (ab) are, respectively, $max { T a + L a , T b + L b − λ 0 }$ and $min { T a + λ 0 + δ 0 , T b + δ 0 }$⁠. Since the direct disruption (a) occurs between $T a + L a$ and $T a + δ 0$⁠, the cross-disruption (ab) always occurs after the direct disruption (a). Both disruptions overlap if disruption (ab) starts before the end of disruption (a), i.e., if $Δ T < L a − L b + λ 0$ or if $Δ T < δ 0 − L b + λ 0$⁠. Since $δ 0 > L a$⁠, the second condition $Δ T < δ 0 − L b + λ 0$ is sufficient for overlapping to occur. Similarly, given that the direct disruption (b) occurs between $T b + L b$ and $T b + δ 0$⁠, the cross-disruption (ab) always occurs before. Both disruptions overlap if disruption (ab) ends after the start of disruption (b), i.e., if $Δ T < λ 0$ or if $Δ T < δ 0 − L b + λ 0$⁠. Since $δ 0 > L b$⁠, the second condition $Δ T < δ 0 − L b + λ 0$ is sufficient for overlapping to occur.