Dynamic interdependence and competition in multilayer networks

TL;DR

This paper introduces a general framework for modeling interdependence and competition in dynamic multilayer networks, revealing complex behaviors like explosive transitions, hysteresis, and chaos in interconnected systems.

Contribution

It develops a novel, broadly applicable dynamic dependency framework and applies it to oscillator networks and epidemic models, uncovering new phenomena.

Findings

Discovered explosive synchronization transitions.

Identified hysteresis and multi-stability in interdependent systems.

Revealed chaotic behavior in complex network interactions.

Abstract

From critical infrastructure, to physiology and the human brain, complex systems rarely occur in isolation. Instead, the functioning of nodes in one system often promotes or suppresses the functioning of nodes in another. Despite advances in structural interdependence, modeling interdependence and other interactions between dynamic systems has proven elusive. Here we define a broadly applicable dynamic dependency link and develop a general framework for interdependent and competitive interactions between general dynamic systems. We apply our framework to studying interdependent and competitive synchronization in multi-layer oscillator networks and cooperative/competitive contagions in an epidemic model. Using a mean-field theory which we verify numerically, we find explosive transitions and rich behavior which is absent in percolation models including hysteresis, multi-stability and chaos. The framework presented here provides a powerful new way to model and understand many of the interacting complex systems which surround us.

Figures

Equations

$${\dot{x}}_i=g(x_i)+\lambda \sum _{j=1}^N{A}_{ij}h(x_i,x_j)$$

$${\lambda }_B\to {\lambda }_B{F}_i^{A\to B}(t)$$

$$F_i^{A\to B}(t)=\{\begin{array}{ll}\left|z_i^A\right|(t)& \text{mbox}interdependent,\\ 1-\left|z_i^A\right|(t)& \text{mbox}competitive,\\ 1& \text{mbox}decoupled\end{array}$$

$$\dot{x}_{i}^{\sigma}=g(x_{i}^{\sigma})+\lambda_{\sigma}\prod_{\mu=1}^{M}F_{i}^{\mu\to\sigma}\sum_{j=1}^{N}A_{ij}^{(\sigma)}h\big{(}x_{i}^{\sigma},x_{j}^{\sigma}\big{)},$$

$$F_i(t)=\left(\begin{array}{cc}1& F_i^{B\to A}(t)\\ F_i^{A\to B}(t)& 1\end{array}\right)$$

$$\text{mathpzc}{Z}_{\sigma }=\int \frac{k{P}_{\sigma }(k)}{⟨k{⟩}_{\sigma }}{G}_{\sigma }({\lambda }_{\sigma }\text{mathpzc}{Z}_{\sigma }{F}^{\mu \to \sigma },k)\mathrm{d}k,$$

$${\dot{\theta }}_i^{\sigma }={\omega }_i^{\sigma }+{\lambda }_{\sigma }{F}_i^{\mu \to \sigma }\sum _{j=1}^N{A}_{ij}^{(\sigma )}\sin ({\theta }_j^{\sigma }-{\theta }_i^{\sigma }),$$

$$z_i(t)=r_i(t)e^{i{\psi }_i}=\frac{1}{k_i}\sum _j{A}_{ij}{e}^{i({\theta }_j)}$$

$${\dot{\theta }}_i^{\sigma }={\omega }_i^{\sigma }-{\lambda }_{\sigma }{r}_i^{\mu }{r}_i^{\sigma }{k}_i^{\sigma }\sin ({\theta }_i^{\sigma }-{\psi }_i^{\sigma }),$$

$${\dot{\Delta \theta }}_i^{\sigma }={\omega }_i-{\lambda }_{\sigma }{k}_i^{\sigma }{R}^{\mu }{R}^{\sigma }\sin \Delta {\theta }_i^{\sigma }.$$

$${\dot{\Delta \theta }}_i^{\sigma }={\omega }_i-{\lambda }_{\sigma }{k}_i^{\sigma }(1-R_{\mu })R_{\sigma }\sin \Delta {\theta }_i^{\sigma }.$$

$${\dot{x}}_i=-\gamma x_i+\beta (1-x_i)\sum _{j=1}^N{A}_{ij}{x}_j,$$

$${\Theta }_i=\frac{1}{{\bar{k}}_i}\sum _j{A}_{ij}{x}_j.$$

$$\dot{x}_{i}^{\sigma}=-x_{i}^{\sigma}+\beta_{\sigma}F_{i}^{\mu\to\sigma}k_{i}^{\sigma}\Theta_{i}^{\sigma}\big{(}1-x_{i}^{\sigma}\big{)},$$

$${\Theta }_{\sigma }=\frac{{\beta }_{\sigma }{\Theta }_{\sigma }{F}^{\mu \to \sigma }}{⟨k{⟩}_{\sigma }}\underset{1}{\overset{+\infty }{\int }}{k}^2{P}_{\sigma }(k)\frac{1}{1+{\beta }_{\sigma }k{\Theta }_{\sigma }{F}^{\mu \to \sigma }}dk.$$

$${\dot{\theta }}_i^{\sigma }={\omega }_i^{\sigma }+{\lambda }_{\sigma }{F}_i^{\mu \to \sigma }\sum _{j=1}^N{A}_{ij}^{(\sigma )}\sin ({\theta }_j^{\sigma }-{\theta }_i^{\sigma }),$$

$$\text{mathpzc}{z}_i^{\sigma }(t)=R_{\sigma }(t)e^{ı{\Psi }_{\sigma }(t)},\forall i\in {\mathfrak{G}}_{\sigma },$$

$${\dot{\phi }}_i^{\sigma }={\omega }_i-{\lambda }_{\sigma }{k}_i^{\sigma }{F}^{\mu \to \sigma }{R}_{\sigma }\sin {\phi }_i^{\sigma },$$

$$R_{\sigma}=\lambda_{\sigma}R_{\sigma}\int_{\mathbb{R}}\mathpzc{x}\,\mathcal{I}_{\sigma}\big{(}\lambda_{\sigma}R_{\sigma}\mathpzc{x}\big{)}\Gamma^{\mu\to\sigma}(\mathpzc{x})\operatorname{\mathrm{d}}\mathpzc{x},$$

$${\Theta }_i:=\frac{1}{k_i}\sum _{j=1}^N{A}_{ij}{x}_j,$$

$${\dot{x}}_i=-\gamma x_i+\beta (1-x_i)k_i{\Theta }_i.$$

$$x_i=\frac{\beta k_i{\Theta }_i}{\gamma +\beta k_i{\Theta }_i},$$

$$x_i=\frac{\beta k_i\Theta }{\gamma +\beta k_i\Theta },$$

$$\Theta =\beta \Theta \underset{1}{\overset{+\infty }{\int }}\frac{k^{2}P(k)}{⟨k⟩\left(\gamma +\beta k\Theta \right)}dk,$$

$$\begin{cases}\begin{aligned} \dot{x}_{i}^{\mathrm{A}}&\,=-\gamma_{\mathrm{A}}x_{i}^{\mathrm{A}}+\beta_{\mathrm{A}}\Theta_{i}^{\mathrm{B}}\big{(}1-x_{i}^{\mathrm{A}}\big{)}k_{i}^{\mathrm{A}}\Theta_{i}^{\mathrm{A}},\\ \dot{x}_{i}^{\mathrm{B}}&\,=-\gamma_{\mathrm{B}}x_{i}^{\mathrm{B}}+\beta_{\mathrm{B}}\Theta_{i}^{\mathrm{A}}\big{(}1-x_{i}^{\mathrm{B}}\big{)}k_{i}^{\mathrm{B}}\Theta_{i}^{\mathrm{B}},\end{aligned}\end{cases}$$

$${\Theta }_{\sigma }=\frac{{\beta }_{\sigma }}{⟨k{⟩}_{\sigma }}\underset{1}{\overset{+\infty }{\int }}{k}^2{P}_{\sigma }(k)\frac{{\Theta }_{\sigma }{\Theta }_{\mu }}{{\gamma }_{\sigma }+{\beta }_{\sigma }k{\Theta }_{\sigma }{\Theta }_{\mu }}dk,$$

Full text

Many real-world complex systems include macroscopic subsystems which influence one another. This feature arises, for example, in competing or mutually reinforcing neural populations in the brain Fox et al. (2005); Stumpf et al. (2010); Mišić et al. (2015), opinion dynamics among social groups Halu et al. (2013), and elsewhere Angeli et al. (2004); Helbing (2013). It is therefore important to understand the possible consequences that different types of inter-system interactions might have. In 2010, substantial progress was made when the theory of percolation on interdependent networks was introduced Buldyrev et al. (2010). This model showed that when nodes in one network depend on nodes in another to function, catastrophic cascades of failures and abrupt structural transitions arise, as observed in real-world systems Buldyrev et al. (2010); Leicht and D’Souza (2009); Kivelä et al. (2014); Boccaletti et al. (2014); Danziger et al. (2016a).

However interdependent percolation is limited to systems where functionality is determined exclusively by connectivity: either to the largest connected component Buldyrev et al. (2010), backbone Danziger et al. (2015) or a set of source nodes Min and Goh (2014). It thus provides only a partial understanding of real-world systems, where the network serves as the base upon which a dynamic process occurs Boccaletti et al. (2002); Barzel and Barabási (2013).

Here, we propose a general framework for modelling interactions between dynamical systems. Two fundamental and ubiquitous ways in which nodes in one system can influence nodes in another one are interdependency or cooperation, as in critical infrastructures Rosato et al. (2008); Shekhtman et al. (2016); Danziger et al. (2016a) or among financial networks Kenett and Havlin (2015); Brummitt and Kobayashi (2015), and competition or antagonism, which is common in ecological systems Hibbing et al. (2009); Coyte et al. (2015), social networks Halu et al. (2013), or in the human brain Fox et al. (2005); Nicosia et al. (2017). It is not uncommon to find interdependent and competitive interactions simultaneously, in predator-prey relationships in ecological systems Blasius et al. (1999), in binocular rivalry in the brain Fries et al. (1997), or even in phenomena like “frenemies” and “coopetition” in social systems Lee et al. (2012). Special cases of competitive interactions between networks have been studied, but without a general framework or ability to uncover universal patterns between systems Masuda and Konno (2006); Gómez-Gardeñes et al. (2015); Valdez et al. (2016).

We model the cross-system interaction by multiplying the coupling strength of a node to its neighbors in one network by a function of the local order of a node in another network. If the function is increasing, the potential for local order of the two nodes is positively correlated, reflecting an interdependent interaction; if it is decreasing, then the potential for local order is anti-correlated, reflecting a competitive coupling. Because local order often reflects the instantaneous local functionality and can be meaningfully defined for a wide range of complex systems, this framework can capture an unprecedented variety of inter-system interactions.

After presenting the general equations for dynamical dependence, we examine in detail a system of two networks of Kuramoto oscillators and a system of two susceptible-infected-susceptible (SIS) epidemic processes, with different combinations of competitive and interdependent interactions between them. We find that under an interdependent interaction, the system exhibits behavior familiar from interdependent networks: abrupt transitions from order to disorder and cascade plateaus at criticality. Furthermore, because of the added richness of the dynamical models, new features such as a forward transition (and hysteresis), and a metastable region are observed. Similarly under a competitive interaction, we find coexistence, hysteresis and multi-stability. When the two types of interactions are asymmetrically implemented, we observe novel oscillatory states and chaotic attractors. Because the new interactions are expressed via local order terms, we are able to perform a mean-field approximation of the exact equations, which are solved numerically and verified against extensive GPU-accelerated simulations on large synthetic networks.

I Model

We begin by considering isolated dynamical systems composed of $N$ nodes which evolve according to

[TABLE]

where $x_{i}$ is the dynamic state of node $i$, $g$ and $h$ are scalar-valued functions of self-dynamics and pairwise interactions, respectively, $\lambda$ is the coupling strength and the network of interactions is represented by $A_{ij}$ which is equal to one if nodes are connected and zero otherwise. An extremely wide range of physical models can be described by equations of this form, including oscillators, epidemic processes, opinion models and many others. The question that we raise here is, how can a dependency relationship be defined between systems such as these? Specifically, we want to include the effects of local functionality in one system promoting or suppressing the onset of functionality in another system. To address this lack, we propose to define a link between systems via multiplication of the coupling strength $\lambda$ of a node in one network by a function of the local order $z_{i}(t)$ around a node in another network. Local order is defined as the weighted average of the ordered state of the neighbors of $i$, with the precise definition dependent on the model in question. In general we replace:

[TABLE]

where $F_{i}$ is a function of the local order parameter $z_{i}(t)$ of node $i$ in layer $A$. Because the numbering of the nodes is arbitrary, we assume that node $i$ in layer $A$ affects node $i$ in layer $B$. As seen in Fig. 1, a suitable choice of $F_{i}$ can represent an interdependent or competitive interaction. We thus consider

[TABLE]

and leave more exotic interactions for future study.

In this manner, we can describe the evolution of the dynamic state in layer $\sigma$ of an $M$-layer ensemble of interdependent, competitive or mixed interacting dynamical systems with the equations:

[TABLE]

where we assume, for simplicity, that the $g$ and $h$ functions are the same in each layer, i.e., the same process is taking place in each layer. We further assume that the time scale of the processes is the same. The product in Eq. (4) reflects the assumption that if the interaction term $F_{i}^{\mu\to\sigma}$ goes to zero in any of the layers, it suppresses the coupling of node $i$ in every layer, reflecting mutual interdependence (or competition).

We can thus consider the $F_{i}^{\mu\to\sigma}$ terms as elements of a supra-adjacency matrix $F_{i}(t)$ describing the interactions between layers at node $i$, which for the case of two interacting networks $A$ and $B$ would be represented as:

[TABLE]

where we have assumed that there are no self-interactions.

Because $F_{i}$ is determined entirely by the local order parameter, it is straightforward to analyze the system with a mean-field approach by simply replacing $z_{i}$ with $\mathpzc{Z}$, the corresponding global order parameter. For many systems described by Eq. (1), the sum on neighbors can be rewritten as a function of the local order parameter, and this is used as the basis for a mean-field theory. Because we also define the interaction in terms of local order, the mean-field approach simultaneously solves the intra-layer and inter-layer dynamics. For two networks, we can then determine a general self-consistent equation for the global order parameter:

[TABLE]

where $F^{\mu\to\sigma}$ is equal to the mean-field interaction term: $\mathpzc{Z}_{\mu}$ for interdependent and $1-\mathpzc{Z}_{\mu}$ for competitive and the $G$ function is a dynamics-dependent function based on the mean-field solution of the single layer case. One feature of this function is that it always has a multiplicative factor of $\mathpzc{Z}_{\sigma}$ and therefore a zero solution, which may or may not be stable. This equation holds when all of the nodes are coupled. In the methods section, we provide the generalized equation for the case where only a fraction $f<1$ are coupled and the remaining nodes are decoupled. Thus for any mean-field solvable dynamical system, our framework enables study of the evolution of the coherent state under any combination of competitive and interdependent links.

In the following, we will demonstrate how this new type of cross-layer dynamic link impacts the onset of macroscopic order in two archetypal systems: Kuramoto oscillators and reversible (SIS) epidemics. We examine these systems in three configurations of cross-layer links: (1) interdependent both ways, (2) competitive both ways and (3) one way interdependent and one way competitive (henceforth referred to as “hybrid”). All three cases have real-world motivations and are studied here with full coupling ($f=1$) and partial coupling ($f<1$).

II Interdependent and competitive synchronization

Synchronization is a common feature of diverse physical systems. It has been observed that competitive and cooperative interactions between neural populations play key roles in vision Fries et al. (1997) and elsewhere in the brain Fox et al. (2005).

By convention, we refer to the dynamic phase of oscillator $i$ as $\theta_{i}$, take $g(\theta_{i})=\omega_{i}$ as the constant function mapping each node to its natural frequency $\omega_{i}$, and $h(\theta_{i},\theta_{j})=\sin(\theta_{j}-\theta_{i})$, we obtain from Eq. (4), for the case of two networks:

[TABLE]

where $F_{i}$ is defined as in Eq. (3) and the local order parameter $z_{i}$ is defined as:

[TABLE]

where $\psi_{i}$ is the instantaneous average phase of the neighbors of $i$.

For the fully interdependent case, we have:

[TABLE]

where we have rewritten the interaction term using the local order (see Methods for details). This can be reduced to $N$ decoupled equations using a mean-field approximation $r_{i}^{\mu}=R^{\mu}=\lvert\sum\exp(i[\theta_{i}-\Psi])\rvert$ for global average phase $\Psi$ and moving to the rotating frame by letting $\Delta\theta=\theta_{i}-\Psi$:

[TABLE]

Using a continuum approximation we then obtain a self-consistent integral equation for the global order in both networks which defines the fixed points, stability and phase flow for this system, as described in the methods section. This leads to a self-consistent equation of the form of Eq. (5) with $G(x,k)=x\int_{|\omega|<x}d\omega g_{\rm freq}(\omega)\sqrt{1-\frac{\omega}{xk}}$, the familiar equation from the self-consistent solution following Kuramoto Kuramoto (1975). We note that the special case of two adaptively coupled Kuramoto networks with $\lambda^{\sigma}=\lambda^{\mu}$ has been considered previously in the context of explosive synchronization Filatrella et al. (2007); Zhang et al. (2015); Danziger et al. (2016b).

We find that, similar to the mutual giant component of interdependent percolation, the global synchronization level undergoes discontinuous backward transition as the coupling strength is decreased. An interesting and significant departure from the percolation models is the existence of a novel forward transition from desynchronized to synchronized state. When $f=1$ the zero solution is always nominally stable but becomes unstable only due to the existence of a nearby saddle-point (Fig. 2a and Supp. Fig. 1-2) and fluctuations of the size of the basin of stability Zou et al. (2014). However, for $f<1$, as the coupling strength is increased, there is a point at which the desynchronized solution ceases to be stable and the system spontaneously jumps to the synchronized branch (Fig. 2c and Supp. Fig. 3-4). This forward transition does not exist in percolation models and its existence allows us to delineate a clear metastable region which we can roughly predict by assuming a characteristic fluctuation size (Fig. 2b,d). Metastable synchronization is absent in the standard Kuramoto model but has been observed, for example, in Josephson junctions Barbara et al. (1999); Filatrella et al. (2007).

Turning to the case of competitive synchronization, which is relevant for competing neural populations Fox et al. (2005), in particular the visual processing of optical illusions Eagleman (2001), we can follow the same steps as above but with $F_{i}=1-z_{i}$ to obtain:

[TABLE]

This system shows pronounced metastability without coexistence for $f=1$ (Supp. Fig. 5-7) and particularly interesting phenomena when $f<1$. In addition to the simple states of desynchronization or total domination by one network or the other, we find chimera-like states in which the decoupled nodes synchronize while the competitively coupled nodes remain desynchronized leading to a global synchronization level bounded by $f$, the fraction of nodes with the competitive interaction. We find that the system can transition continuously between the state where only one network synchronizes to the state where one dominates but the other has chimera-like partial synchronization and then discontinuously as the fully synchronized network drops to partial synchronization and the partially synchronized network becomes fully synchronized (Fig. 3). We further find multistability of up to four solutions for a small region of the phase space (see Supp. Fig. 8-9). These features demonstrate the rich synchronization patterns that can be caused by competing networks.

Finally, we consider the hybrid case, in which the synchronization of nodes in network $A$ promote the onset of synchronization of nodes in network $B$, even as that onset of synchronization in network $B$ suppresses the ability of nodes in $A$ to synchronize. This type of behavior is observed in binocular rivalry, in which neurons associated with the dominant eye synchronize more strongly when the weak eye is stimulated, but the weak eye synchronizes less strongly when the dominant eye is stimulated Fries et al. (1997). In this case network $A$ will be described by Eq. (10) while network $B$ will be described by Eq. (9). In such a case, we find convergent oscillatory solutions when $f<1$ (Fig. 4a-c) and chaotic attractors when $f=1$ (Fig. 4d,e). The source of the chaotic behavior is related to the fact that when the order parameter oscillates, it causes the effective coupling strength to also oscillate and oscillatory coupling strengths have been shown to lead to the onset of chaos So and Barreto (2011). The figures here have been derived assuming uniformly distributed natural frequencies. Results for Lorentzian frequency distribution are shown in Supp. Fig. 10.

A unique advantage of our model is that, because the coupling between the networks is on the level of order and not of phase, we are now able to model the cooperative onset of synchronization at different frequency bands. This is significant in light of the complex interactions between neural populations of different frequencies Gans et al. (2009); Lowet et al. (2015) which cannot be captured by existing multilayer models and the potential for networks of synchronizing oscillators to represent learning tasks Reichert and Serre (2013), modeling fundamentally multi-frequency phenomena such as hearing Wang et al. (2016) and physiological synchronization between organ systems in the body which may be very different from one another Bashan et al. (2012); Ivanov et al. (2016). These phenomena may be more naturally modelled using the idea of order affecting order rather than as oscillators representing different entities summing their phases directly. Mathematical treatment of multimodal interdependent synchronization is presented in Supp. Sec. S3.

III Epidemics

Because the cross-system dynamic relationships defined in Eq. (4) are general, it is simple to adapt them to other scenarios in which the dynamic state of a process in one network affects its activity in another layer. In epidemics, multiple diseases or strains can spread cooperatively Newman and Ferrario (2013); Chen et al. (2013); Azimi-Tafreshi (2016), with the onset of one disease increasing the susceptibility to other diseases. At the same time, information about vaccines can spread on social networks similar to those upon which the disease spreads, leading to competitive contagions Newman (2005); Sanz et al. (2014); Sahneh and Scoglio (2014). If exposure to the disease in one layer induces an individual to vaccinate in another layer, we may find that the disease layer enhances the vaccine layer, even as the vaccine layer suppresses the disease layer (the assymetric “hybrid” case of our model) Ahn et al. (2006). Our model of dynamic dependence offers a unified framework to study a broad spectrum of possible interactions between epidemics on multilayer networks. For simplicity, we treat a probabilistic reversible (SIS) model which, in the single layer case evolves according to:

[TABLE]

where $\gamma$ is the recovery rate and $\beta$ is the transmission rate. For simplicity, we renormalize the time such that $\gamma=1$ and define the local order parameter as

[TABLE]

This immediately leads to the epidemics analog of Eq.(4):

[TABLE]

where $F_{i}$ is defined as in Eq. (3) for interdependent or competitive interactions. Using the standard mean-field theory Pastor-Satorras and Vespignani (2001) we let $\Theta_{i}\rightarrow\Theta$ and obtain the general self-consistent equations (special case of Eq. (5) with $G(x,k)=\frac{x}{1+kx}$):

[TABLE]

Here, $F^{\mu\to\sigma}$ equals $\Theta_{\sigma}$ for the interdependent case, and $1-\Theta_{\mu}$ for the competitive case. Letting $F^{\mu\to\sigma}=1$ we recover the known decoupled single-layer case Pastor-Satorras and Vespignani (2001). Analagous equations for $f<1$ are provided in the Methods section.

Looking at the results of these interactions, we find a number of novel and realistic phenomena. In the interdependent case, we find that there is a first-order transition in the forward direction and an abrupt transition in the backward direction (Fig. 5a,b) when $f<1$ and no forward transition at all when $f=0$. In contrast to the case of synchronization, where a finite level of fluctuations (due to the quenched disorder of the natural frequency distribution) made the zero-solution spontaneously jump to the synchronized state, the system leaves the zero solution only when it becomes a proper unstable fixed point. In practical terms, we see that for interdependent epidemics it is much harder (i.e., requires a much larger reduction in $\beta$) to stop an epidemic once it has become endemic than it is to keep it from breaking out. Another practical consequence of the bistability is that, in contrast to non-cooperative epidemics, the zero solution has a finite basin of stability even after the endemic phase has emerged. This means that in cooperative/interdependent epidemics, small outbreaks are expected to die out even for comparatively high transmission rates, but that a large outbreak can become endemic.

In the competitive case we also see new behaviors that differ from synchronization. Instead of metastability, we find a broad regime of coexistence (Fig. 5c,d). This indicates that it is possible for the disease and its cure to coexist in the same system, presenting a challenge for any attempt to completely eradicate it. The hybrid case (shown in Supp. Fig. 13) does not display the chaotic behavior of the analagous synchronization system but instead has oscillating but convergent attractors (Supp. Fig. 14).

IV Connection with interdependent percolation

We argue here that the new framework developed above represents a generalization of interdependence from percolation to general dynamic systems (those which follow Eq. (1)). The reason we describe the interaction in Eq. (4) as representing a dynamic dependency link is that it has the same impact on dynamic order that the percolation dependency link has on the connectivity-based order of interdependent percolation. This is because, when we implement the dependency interaction from Eq. (3), the result is that as long as node $i$ is locally disordered in network ${B}$, the coupling strength of node $i$ in network $A$ is multiplied by $|\mathpzc{z}_{i}^{\mathrm{B}}|\approx 0$, effectively cutting it off from the other nodes in $A$ and suppressing the onset of order in $A$. This reflects a dynamical generalization of interdependent percolation Buldyrev et al. (2010): in percolation the local order of a node is defined as whether any of its neighbors leads to the largest connected component; if it has no local order, then it suppresses the ability of the node that depends on it from becoming ordered. As with interdependent percolation, dynamic interdependence leads to increased vulnerability as reflected in a higher threshold and the abruptness of the transition. There have been several studies of competitive percolation as well, but the absence of a natural temporal evolution has led to difficulties in describing the dynamics in each system Zhao and Bianconi (2013); Kotnis and Kuri (2015); Watanabe and Kabashima (2016).

One of the most intriguing findings that we report here is that all of the systems: percolation, synchronization and epidemics, undergo a remarkably similar transition in the backward direction (from order to disorder). In all of the cases, the transition is characterized by a slow cascading process with a “plateau” (Fig. 6). The system quickly falls from a highly ordered state to an intermediate state and then spends the bulk of its relaxation time with the order parameter maintaining an approximately constant value, before dropping to zero at an exponential rate. While this cascade has been studied in the context of a branching process in percolation Zhou et al. (2014); Lee et al. (2016), here we observe it through a fixed point that transitions from stable to marginally unstable (see phase flows in Supp. Fig. 2 and 4). The macroscopic signature of this process is remarkably similar across these highly diverse processes, reflecting a universality to interdependent transitions across diverse dynamics.

V Conclusion

Here we have presented a universal theoretical framework for interdependent and competitive links between general dynamic systems. Whereas previous attempts to include new effects in multi-layer systems have required ad-hoc definitions and assumptions, we now present a simple and realistic approach that can be easily adapted to a wide variety of systems. Furthermore, because the connections are based on the fundamental statistical physical properties of local-order, standard mean-field methods can be applied straightforwardly to obtain analytic results.

This generalization of interdependence to dynamical systems allows for the development of new models for neural, social and technological systems that better capture the subtle ways in which different systems can affect one another. Though we focus here on the fundamental interactions of interdependence and competition, more exotic interactions can be described by this framework by replacing the simple linear functions of $z_{i}$ in Eq. 3 with more complex functions. We find that the phenomenology recovers key features of interdependent percolation but also uncovers new phenomena like oscillatory and chaotic states, metastability and coexistence, which were not known to be caused by inter-system interactions.

Materials & Methods

Mean-field approximation for interacting synchronization.

We consider therefore a network of Kuramoto models, each describing a population of oscillatory elements, and model the cross-system interactions according to our framework. For simplicity, let us focus on a system of $M=2$ dynamically dependent Kuramoto models on arbitrary networks of equal size. Assuming no self-dynamical dependencies, the dynamics of the overall system is described by the set of $2N$ coupled equations

[TABLE]

where $i=1,\dots,N$, and $\sigma=\mathrm{A},\,\mathrm{B}$ with $\mu\neq\sigma$.

In order to identify all the possible asymptotic coherent states of the model (15), we can follow the traditional self-consistent method developed by Kuramoto Kuramoto (1975, 2012). In particular, let us consider the case of two uncorrelated random graphs with prescribed degree sequences Bollobás (2001); Newman (2010); Cohen and Havlin (2010). Adopting the so-called annealed network approximation Bianconi (2002); Dorogovtsev (2010); Boccaletti et al. (2014), we replace the entries $A_{ij}^{(\sigma)}$ with their ensemble averages, i.e. the probabilities $p_{ij}^{\sigma}=k_{i}^{\sigma}k_{j}^{\sigma}/(\langle k\rangle_{\sigma}N)$ that vertices $i$ and $j$ with degree $k_{i}^{\sigma}$ and $k_{j}^{\sigma}$, respectively, are connected in layer $\sigma=\mathrm{A},\,\mathrm{B}$, being $\langle k\rangle_{\sigma}$ its corresponding mean degree. Within this approach, each local order parameter (7) becomes independent of the node indices, so that

[TABLE]

meaning that each node feels the global mean-field measured over the entire network it belongs to. Let us further assume that in the thermodynamic limit both populations reach asymptotically steady states with constant mean-field amplitudes and constant frequencies, so that $\dot{R}_{\sigma}=0$ and $\ddot{\Psi}_{\sigma}=\dot{\Omega}_{\sigma}=0$ for $\sigma=\mathrm{A},\mathrm{B}$, and choose equal, symmetric and unimodal frequency distributions on layers $\mathfrak{G}_{\mathrm{A}}$ and $\mathfrak{G}_{\mathrm{B}}$. In this case, the equations governing the phase dynamics (15) can be rewritten as

[TABLE]

where $\varphi_{i}^{\sigma}\equiv\theta_{i}^{\sigma}-\Omega_{\sigma}t$ are local phase differences Basnarkov and Urumov (2008); Petkoski et al. (2013); Iatsenko et al. (2013) defined after moving to rotating reference frames on each unit circle comoving with frequencies $\Omega_{\sigma}$ (see Supp. Sec. 1 for more details).

Solutions of the system (17), shows that the population of oscillators in each layer splits into two groups, namely drifting and frequency locked oscillators Strogatz (2000). The phases of the latter ones are entrained by the mean-field and correspond to fixed point solutions of the system (15), i.e. $\dot{\boldsymbol{\varphi}}^{\sigma}=\boldsymbol{0}$. Drifting oscillators, on the other hand, are not entrained and never reach steady states. Nonetheless, based on the model’s assumptions, drifting oscillators have vanishing contributions to the complex order parameters in the thermodynamic limit (see Supp. Sec. 1), so that the overall collective behaviour of the system is dominated entirely by the locked oscillators. We find hence that all the possible coherent stationary states of the system (15) can be obtained by solving self-consistently a system of two integro-transcendental equations

[TABLE]

where $\mathcal{I}_{\sigma}$ is an integral function accounting for the degree distributions of the network $\mathfrak{G}_{\sigma}$ and the natural frequency distribution of the Kuramoto oscillators sitting on its nodes (see Supp. Sec. 1), whilst $\Gamma^{\mu\to\sigma}(\mathpzc{x}):=f\delta(\mathpzc{x}-F^{\mu\to\sigma})+(1-f)\delta(\mathpzc{x}-1)$ is the distribution of the dynamical dependencies from $\mathfrak{G}_{\mu}$ to $\mathfrak{G}_{\sigma}$, with $\sigma,\,\mu=\mathrm{A},\,\mathrm{B}$ and $\mu\neq\sigma$.

Notice that the null solution describing the mutual incoherent phase always satisfies Eqs. (18), though it not always stable. Non-vanishing solutions of (18) can be found numerically for general graph topologies and frequency distributions Lafuerza et al. (2010); Omel’chenko and Wolfrum (2012), once the strategy for the dynamical interactions between the layers has been chosen.

Mean-field approximation for interacting epidemics.

Consistently with the general framework (7)–(4), let us introduce a global order parameter for SIS epidemics, here defined as $\Theta:=\sum_{i}k_{i}\Theta_{i}/\sum_{i}k_{i}$, where

[TABLE]

are local order parameters measuring the local spread of the disease in the neighbourhoods of the $i^{\mathrm{th}}$ node Barrat et al. (2008); Newman (2010); Pastor-Satorras et al. (2015). Inserting the latter into the above equations, one finds

[TABLE]

Stationary solutions can then be found by setting $\dot{x}_{i}=0$ in the latter expression, which yields the local self-consistent equations

[TABLE]

for every $i=1,\dots,N$. The above identifies what is in the literature known as individual-based mean field approach. In what follows, we will consider instead the degree-based approximation introduced by by Pastor-Satorras et al. Pastor-Satorras and Vespignani (2001); Pastor-Satorras et al. (2015), which allows us to solve just one self-consistent equation for the global order parameter. Moreover, it is straightforward to prove that this approximation is actually equivalent to the annealed network approximation for uncorrelated random graphs Newman (2010); Pastor-Satorras et al. (2015), in which case the local order parameters $\Theta_{i}$ become node-independent so that, $\Theta\equiv\Theta_{i}=\frac{1}{\langle k\rangle N}\sum_{i=1}^{N}k_{j}x_{j}$. Stationarity implies then

[TABLE]

so that one can find all the possibly stationary states by solving the self-consistent equation

[TABLE]

where we replaced the sums with definite integrals over the degree distribution after moving to the thermodynamic limit. Moving to dynamically dependent epidemic processes, let us consider the the moment the case of fully interdependency links between the replica nodes. In such case, node $i$ in layer $\mathfrak{G}_{\mathrm{A}}$ feels the local spread of the disease of its replica node in layer $\mathfrak{G}_{\mathrm{B}}$, which adaptively changes its local rate of transmission, i.e. $\beta_{\mathrm{A}}\mapsto\beta_{\mathrm{A}}\Theta_{i}^{\mathrm{A}}$ and viceversa. Hence

[TABLE]

and one arrives eventually to the following self-consistent system of equations

[TABLE]

for $\sigma=\mathrm{A},\,\mathrm{B}$ and $\mu\neq\sigma$.