Inter-layer competition in adaptive multiplex network

We consider competition between layers in adaptive multiplex networks of phase oscillators, where adaptation principles (which cause intra-layer topology evolution) are inspired by real world homophily and homeostasis phenomena. Our model yields the emergence of both scale-free topologies and meso-scale structures in the layers, for an appropriate choice of the control parameters. We further report that the growth of the number of interacting layers leads to a decrease of the global order, due to inter-layer structural competition. However, the increase of the system’s scale can effect local synchronization between neighboring (or strongly coupled) nodes. Such unforeseen phenomena is connected with the nature of the competitive mechanism, which implies the rivalry for optimal structure within the whole system, a situation occurring in a variety of natural systems.


Introduction
Since the beginning of complex networks theory the main motivation of the researchers was the aim at describing processes that take place in the real world systems [1]. As always in science, new ideas and approaches emerged as a consequence of the need of a more accurate correspondence between results from models and real data. This is why more and more attention is currently paid to multi-layer networks [2], as they very elegantly furnish a representation of systems where multiple relations exist between networking elements.
In particular, multiplex networks (the sub-class of multi-layer networks where nodes are the same on all layers) found applications in technological [3][4][5], biological [6], and social systems [7,8], as well as in human brain, artificial intelligence and robotics [9,10]. Multiplex networks are furthermore useful tools for the analysis of climatological data, e.g. to represent and investigate the topology of statistical relationships between the fields of distinct climatological variables [11]. In [12], four kinds of immunization strategies for multiplex networks are proposed to describe the spread of diseases in the real world, and the effectiveness of those strategies on different types of topologies is investigated. Using multiplex networks is also appropriate for describing various transportation networks: [13] considers the influence of the topology of each network layer on the emergence of structural features in the whole network using the European air transportation network, where each layer represents a commercial airline. Kurant et al [14] used a layered model to study the load distribution in three transportation systems, where the lower layer is the physical infrastructure and the upper layer represents the traffic flows. Urban bus-transport and railway systems can be mapped in two spaces, where the nodes represents a bus stop or a railway station [15]. At the same time, different levels of country organization can be presented in a multiplex framework: [16] reveals the features of (and interactions between) spatially independent and spatially distributed networks by multiplexing the administrative structure and geo-population relationships of a given country. Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
On the other hand, research on collective dynamics in multiplex networks concentrated so far on synchronization processes between the network's layers, as emergence of synchronized states plays a significant role in many natural phenomena [2,17,18]. Furthermore, competition between layers is ubiquitous in real world networks, e.g. in urban mobile networks, and in neuroscience when describing the processes occurring both in local neuron networks and during the interaction between different parts of the brain [19][20][21], in particular interactions between dynamical timescales in the neural network [22]. In this paper we study how the competition between network layers can influence the structural and dynamical features of multiplex networks.
Despite some obvious differences among them, real world networks feature a certain number of similar properties [23], e.g. a power-law distribution of the connection weights and a strong modularity, i.e. the presence of structures at a meso-scopic scale. The emergence of these structural and dynamical features can be explained as a consequence of a self-organization processes involving structure-dynamics adaptation of two fundamental principles [24]: (a) the connections between nodes which have a similar (synchronized) behavior, have a trend to be enhanced; (b) the available resources for each node to sustain its connections with other elements of the network are limited. The first principle is well-known under the terms of homophily [7,25] or Hebbian learning [26,27], and the second one models what is called homeostasis [28]. We will show that the competition, balance, and trade-off between these two mechanisms are the key ingredients to describe intralayer pattern formation in adaptive multiplex networks.
To construct a consistent presentation, the paper is divided into several sections. Section 2 contains a specification of the mathematical model under study, including the description of the adaptive processes that we consider for the evolution of the multiplex network. Section 3 presents the obtained results, specified in the case of a two-layer network (section 3.1), where we give a good sight on how the network topology changes in response to such competition process, and describes how the size of the system effects its structural and dynamical properties (section 3.2). Section 4 concentrates on some analytical results, and the over all obtained scenarios are discussed in the conclusive section.

Model under study
The main goal of the current study is understanding how the competition between network's layers can reshape their topology in the presence of homophily and homeostasis. We use the Kuramoto model, which is a well proved paradigm for describing various forms of collective dynamics in a variety of physical, biological, chemical and social systems. Real-life systems have been represented as networks of coupled phase oscillators in studies of synchronization phenomena [29][30][31][32][33][34][35]. We here aim at showing that two key features of real world networks (a scale-free distribution of the weights, and the formation of meso-scale structures) are the results of selforganization under the presence of the above-mentioned mechanisms, homophily, and homeostasis.
We consider M layers, and N oscillators (nodes) on each layer. Each oscillator i (i = 1... N) on the layer l (l = 1... M) is characterized by a natural frequency ω i (equal for all layers, and in this sense one speaks of a multiplex network) and by a phase i l j (which may, instead, vary from layer to layer), and interacts with all other nodes in the lth layer. The phase dynamics is ruled by where w t ij l ( ) stands for the weight of the connection between nodes i and j on layer l, and λ is the coupling strength. At the same time, the weights are here considered as co-evolving quantities, which follow (for each pair i, j of nodes and in each layer) the equation: where N i l is the neighborhood of node i on layer l, i.e. the set of indices which denote the nodes connected to node i on layer l.
The time-dependent parameters p t ij l ( ) are the measures of phase coherence between the nodes i and j on layer l averaged over the time interval [t−T, t], which are described by the following equation: with T being a characteristic adaptation time. To implement the homeostatic property in our model we set the total weight of all incoming connections of each node in each layer to be a constant value over all times. Initially, the weights are randomly assigned and normalized as whereas all phases, i l j , and natural frequencies, ω i , are chosen randomly in the interval [−π; π]. It should be noted that the second and the third terms in the right side of equation (2) reflect adaptive interlayer interactions and competitive interactions between the layers for optimal topology, respectively. The homophily property is accounted for by the first term in the right side of equation (2): the high degree of phase coherence between the nodes i and j on layer l enhances the strength of their inter-layer connection. At the same time, the second term in the equation reflects homeostasis in that it subtracts the interaction between all over neighboring nodes in the layer l. Finally, the third term is responsible for the competitive interaction between the layers of the network for optimal topology: the increase of synchrony between nodes i and j in the layer l leads to a decrease of the weight between the corresponding nodes in other layers. A schematic representation of these adaptive mechanisms is reported in figure 1.

Numerical results
In this section we present numerical results from simulations of the model described above, with the aim of describing how structural and dynamical features emerge as functions of our two key parameters: the coupling strength λ, and the adaptation time T. The section is organized in two subsections: section3.1 depicts the case of a two layers network, while section3.2 examines how such emergent properties depend on the size of the network, i.e. on the number of interacting layers.
In order to monitor the system's features, it is useful to rely on the calculation of three order parameters. The first is the time-dependent measure of the degree of global synchronization: The second parameter quantifies the average value of the order parameters inside the networks layers: Finally, the average synchronization between connected nodes in the graph can be defined as the local order parameter, and is measured as:

Two layers network
Let us first concentrate on the case M=2, and N=500 oscillators on each of the two layers. Initially (i.e. for t ä [0, 500]) we evolve the phases of oscillators by means of equation (2) with weights w ij fixed on their initial conditions, i.e. we let the system unfold without adaptation. Then (i.e. at t=500 time units), the homophily and homeostasis terms are activated, and we let the system progress for the next 2,000 time units with the co-evolving weights of equation (3). To warrant attainment of the network topology to its asymptotic state, we also monitor the total change in the adjacency matrix over one unit step of integration: All order parameters are calculated starting from the time at which γ<10 −4 , and their values are accumulated (and averaged) over the next 500 time units. The three panels of figure 2 report the global, intra-layer and local order parameters in the (T, λ) parameter space. In the first panel, one notices that for low coupling strengths (below the critical value for global synchronization) the adaptive mechanism enhances global order in the network, in that the increase of the adaptation time causes the growth of r g . Further, the second panel indicates that the layer order parameter r grows with λ faster than r g . The areas where the calculated values of r g and r become significantly different correspond to an asynchronous behavior between nodes with same index located on different layers. Finally, one notices that r link is strongly dependent on T, and the growth of adaptation time generically enhances synchronization of neighboring nodes. Such a feature is the hallmark of the emergence of modular synchronization, i.e. the network comes out to be segregated into distinct clusters.
In order to grasp more details on the network's adaptive dynamics, we concentrate the attention on three specific configurations, marked by the points A=(0.2, 20), B=(1.20, 45) and C=(2.80, 55) in figure 2. As weights are evolving, the inequality in the two layers' topology can be quantified by the final (i.e. at the last calculation time) total difference between their adjacency matrices:  Points labeled with A, B and C correspond to three particular system's configurations that will receive specific attention in our paper.
However, at moderate values of λ and T the intra-layer order parameter exceeds the global order, and therefore there is evidence of formation of different intra-layer topologies as a result of strengthening the connections inside the layer. Figure 4 reports the visualizations of each layer's topology at the points A, B and C. In the configuration A, the system selects a inhomogeneous topology, with weakly pronounced structural clusters. This topology has a resemblance with a power-law weight distribution, and can be considered as the closest to real systems. When considering configuration B, a strongly modular (clustered) structure begins to appear with growing of λ and T. Figure 4 allows to visualize the inequality in the layers' topologies: the structure of the second layer is much more modular. Finally, when control parameters approach a critical value (like close to configuration C), clusters begin to disintegrate, and the network returns to a state similar to that of configuration A, but this time the weight distribution strongly deviates from a power-law, and becomes much more homogeneous.

Multilayered networks
It is typical for real systems to have more than two types of relations between their constituent units. Therefore, it is useful to examine the case of more than two layers. To this purpose, we let M vary from M=1 to M=20, and examine how the in the number of layers influences the behavior of the system for the three configurations of control parameters which were discussed in the previous section. Figure 5 presents the order parameters versus the number of layers in the configurations A, B and C. In all configurations, the global order ( figure 5(a)) has a tendency to decrease when increasing the size of the system, indicating a weakening of synchronization with the growth of M. However, its dynamics strongly varies in different configurations: while at A one observes an almost linear dependence, at B and C a drastic decrease is observed already at the transition from a one layer to a four-layers multiplex. These results suggest that the effects of competitive mechanisms between layers become more and more pronounced with both the increase of intra-layer synchronization and the size of the system. At the same time, the intra-layer order parameter (see figure 5(b)) practically does not depend on the number of layers in the network, and only displays different values at points A, B and C.
The most unexpected result concerns the local order parameter, r link , when averaged over layers (see figure 5(c)). Despite the stability of intra-layer order, indeed, one sees that increasing the number of layers results in a considerable decrease of the synchronization between neighboring nodes, i.e. increasing the size of network depletes synchronization at the local scale, and makes the competitive mechanisms more expressed.
It is important also to elucidate the impact of the intra-layer connections on the formation of different topologies among the layers. Figure 6(a) presents the adjacency matrices of a five-layer network with nodes sorted by their ranking value in the second smallest eigenvector of the 1st layer's Laplacian matrix. Such an eigenvector, also known as Fiedler's vector, is connected indeed with network's modularity. Therefore, by calculating the second smallest eigenvector, one can encapsulate the partition of each of the five-layers of the network (presented on figure 6(b)), and the adjacency matrices give a representation on the structure of connections within the clusters that are formed.
The plots shown in figure 6 illustrate the network structure for the configuration B, which features prominent modular effects. One can notice the presence of communities on each layer, which is demonstrated

Analytical treatment
Let us prove our numerical results in the framework of analytical treatment of the network behavior. For this purpose we expand the approach proposed in [23] for analytical estimation of the network weight distribution in the uncoupled limit ( = + D = D ( )in the uncoupled limit, where Δ ij =ω i −ω j . Let us analyze the entire system by taking as reference frequency that of the ith, ω i node on the strength of Since initially all other natural frequencies ω j are uniformly distributed from −π to π, then Δ j =ω j −ω i are also independent identically distributed (i.i.d.) random variables, whose conditional distribution p j i w D D ( | ) is uniform in the range [−π−ω i ; π−ω i ]. Thus, all L j D ( )(or simply L j ) are also i.i.d., and their probability distribution function (PDF) conditioned on ω i is In accordance with equation (2), the weights are asymptotically stable values. Taking into account both global intra-layer interaction and multiplex inter-layer influence between the elements of the network, one can derive the equation for the weights in the following form According to the properties of the PDF for L j [23], one can apply the central limit theorem to G, and finds that it converges to a Gaussian random variable, whose expectation value and variance are (in the limit of large N) G N L j á ñ = á ñand var[G]=N var[L j ], respectively. Hence, following [23], one can obtain the analytical expression for the distribution of w ij j at 0 l  as:  Figure 7 shows the comparison of the weight distribution from numerical simulations at negligibly small coupling strength λ=10 −4 (T=100 and M=2) with the analytical curve from equation (10). One can see that equation (10) gives indeed a rather good approximation of the weight distribution also for very small values of λ. In case of large λ clustering takes place in the network under study. Thus, the appearance of wellpronounced peak in network weights distribution indicates that network behavior is is no longer subject to the analytical model.

Conclusion
In summary, we have considered inter-layer competitive phenomena in an adaptive multiplex network of phase oscillators with topology of connections between the units established in accordance with mechanisms of homeostasis and homophily. We used the model based on Kuramoto paradigm, and examined the effects of coupling strength, adaptation time and size of the network at different stages and configurations.
In the process of adaptation, the model acquires topological properties which are common for a rather vast class of real world networks, such as scale-free distribution of the weights, and presence of cluster structures at a meso-scopic scale. These properties emerge only for moderate values of the control parameters (the coupling strength λ, and the adaptation time T), in connection with partial synchronization and yet rather high values of local order parameter. We also revealed that the increase in the number of interacting layers leads to a decrease of global order due to inter-layer structural competition, while the average intra-layer order remains the same.
Our results may pave the way to understand some typical phenomena of natural systems, like the rivalry of different layers of interactions for optimal structure within a whole system.