A Inter-layer correlations

# Non-linear growth and condensation in multiplex networks

## Abstract

Different types of interactions coexist and co-evolve to shape the structure and function of a multiplex network. We propose here a general class of growth models in which the various layers of a multiplex network co-evolve through a set of non-linear preferential attachment rules. We show, both numerically and analytically, that by tuning the level of non-linearity these models allow to reproduce either homogeneous or heterogeneous degree distributions, together with positive or negative degree correlations across layers. In particular, we derive the condition for the appearance of a condensed state in which one node in each layer attracts an extensive fraction of all the edges.

###### pacs:
89.75-k, 89.75.Fb, 89.75.Hc

Various complex systems are well described by multiplex networks of nodes connected through links of distinct types, which constitute separate yet co-evolving and interdependent layers (1); (2); (3); (4). Examples of multiplex structures can be found in social, technological, transportation and communication systems, and in general wherever a certain set of elementary units is bound by different kinds of relationships (7); (6); (5). In these systems, links of different types are intertwined in non-trivial ways, so that it is not possible to study each layer separately. In particular, a node can have different degrees at the various layers, so that a hub at one layer might not be a hub in another layer (8) or, conversely, the hubs might tend to be the same across different layers (9); (10). Also, it has been shown that the presence of an edge at a certain layer of a multiplex network is often correlated with the presence of the same edge on another layer, which corresponds to a significant overlap of links (4); (1); (5); (11). Some recent studies have focused on dynamical processes on multiplexes, including percolation (12); (13); (14); (15); (16), diffusion (17); (18); (19); (20), spreading (21); (22); (23); (24), traffic (25), cascades (26), and cooperation (27); (28); (29), and a few recent works have suggested that degree correlations (32); (33); (34) as well as overlap of links (13); (30); (31) may have a substantial impact on the emergence and stability of collective behaviors in multiplex systems.

It is therefore interesting to investigate the mechanisms responsible for the appearance of inter-layer correlations in multiplexes. A few different approaches for the modelling of multiplex networks have been recently proposed. Some of them aim at defining appropriate static null-models for multiplexes (7); (6); (5), while some other focus on capturing the non-equilibrium nature of multiplexes and on providing possible physical explanations for their formation (9). However, until now, all the existing models for growing multiplexes with homogeneous and heterogeneous degree distributions allow for positive inter-layer degree correlations only.

In this Article, we propose and study a general growth model of multiplex networks based on a non-linear preferential attachment mechanism. Using both analytical and numerical arguments, we show that this model generates different regimes and displays a transition towards a condensed state where only a few hubs dominate the degree distribution of each layer. Moreover, in the non-condensed regime the model can generate multiplexes with homogeneous or heterogeneous degree distributions, having either positive or negative inter-layer degree correlations. Finally we notice that in the multi-layer version of non-linear preferential attachment the structure of the network dramatically depends on fluctuations, and that the mean-field approach, which was fundamental to understand network growth in single-layer networks, actually fails to a large extent in predicting the dynamics of the growth process.

The paper is organized as follows. In Sec. I we define a general class of non-linear preferential attachment growth models for multi-layer networks focusing, as an example, on the case of a 2-layer multiplex. In this case, the growth is completely determined by the relative values of two attaching exponents, called and . In Sec. II we investigate the role of the exponent when one of the two terms of the attaching kernel is linear, i.e. . This is a first generalization of the classical linear preferential attachment model (36). In Sec. III we derive a mean-field solution for the proposed class of models, and we show that the mean-field approximation fails to account for most of the observed structural properties of the multiplex, in particular regarding the possibility to obtain negative inter-layer degree correlations. In Sec. IV we present the master equation of the model and we solve it to derive the conditions for the appearance of a condensed state. In Sec. V we show and discuss the full phase diagram of the model based on numerical simulations, which is in perfect agreement with the analytical predictions obtained by solving the master equation. In Sec. VI and in Sec. VII we discuss, respectively, the effect of the parameters on the role played by hubs, by means of the recently introduced multiplex cartography (3), and the appearance of mixed degree correlation patterns. In Sec. VIII we focus on the values of characteristic path length and multiplex interdependence obtained as a function of and , while in Sec. IX we show how the model can be calibrated in order to reproduce some of the structural properties of two real-world multiplex networks. In Sec. X we discuss three possible generalizations of the model to the case of -layer multiplex networks, providing also the analytical solution for the boundary of the condensed phase. Finally, in Sec XI we draw our conclusions and we discuss possible future directions of research in the field of multiplex network modelling.

## I Model

Let us consider a multiplex network consisting of layers, one for each type of relationship among nodes, defined by the vector of adjacency matrices , where and if and only if node and node are connected by an edge on layer . A node of the network is characterized by the vector of the degrees of its replicas at each layer, where . We are interested in the mechanisms which might be responsible for the growth of the multiplex. We start from a connected graph with nodes and we assume that, at each time , a new node arrives in the graph, carrying new links in each layer, and that the probability for node to attach on layer to an existing node is a function of the degrees of at all layers:

 Π[ℓ]i→j∝f[ℓ](k[1]j,k[2]j,…,k[M]j) (1)

For the sake of clarity, and without loss of generality, we focus on a multiplex network with two layers, where we denote by the degree of node in layer , and by the degree of in layer , and we assume that

 Π[1]i→j∝f(kj,qj),    and    Π[2]i→j∝f(qj,kj). (2)

In the context of single-layer networks, non-linear attachment kernels of the form , with , have been introduced in Ref. (35), as a generalization of linear preferential attachment models (36). We extend this idea to networks with multiple layers, also allowing for negative exponents to mimic the case in which new nodes prefer to avoid linking to high-degree nodes. We adopt the general expression

 f(kj,qj)=kαjqβj (3)

where, by tuning the two exponents , we can model different attachment strategies. If the exponents and in Eq. (3) are both positive (negative), then new nodes will preferentially link to nodes which are well-connected (poorly connected) on both layers. Conversely, if and (resp. and ), a new node will be preferentially linked, in layer 1, with nodes which are well-connected (resp. poorly-connected) in layer 1 and poorly connected (resp. well-connected) in layer 2. A specular interpretation holds for the attachment probability on layer 2. As we will show in the following, the attachment probabilities in Eqs. (2) and (3) are general enough to produce multiplex networks with different degree distributions, and with positive and negative correlations between the degrees of a node at the two layers. There are several possible ways to generalize this model to the case of more than two layers, and some of them are discussed in Sec. X.

## Ii Semi-nonlinear attachment

Let us first consider the case and , i.e. when the probability to attach to node at layer (resp., at layer ) is proportional to (resp., to ). In particular, when and , we recover the uncorrelated linear preferential attachment kernel, which has been extensively studied in Ref. (9). In this case, the degree distribution in each layer is a power law with , and the multiplex exhibits positive inter-layer degree correlations, the degree of a node being essentially determined by its age.

When the growth process can produce multiplex networks with homogeneous, heterogeneous or condensed degree distribution on each layer, characterized by either assortative or disassortative inter-layer degree correlation patterns, depending on the sign of . In Fig. 1 we report the results obtained by simulating the growth of a multiplex for and in the range . In order to characterize the degree distributions of the two layers we plot, as a function of , the variance of the degree distribution, the maximum degree , the number of different degree classes present in each layer, and the participation ratio . Given a degree sequence , is defined as (38)

 Y−12=⎡⎣∑i(ki∑jkj)2⎤⎦−1. (4)

It is easy to show that when for all , i.e. for homogeneous degree distributions, while if most values of are equal, except for a few nodes for which we have , i.e. in the presence of a condensate state where a few nodes connect to nearly all the other nodes of a layer. When and is positive, we observe a transition to a condensed state, characterized by small , large , , and , signalling the existence of a few dominant nodes. Conversely, for negative values of we obtain heterogeneous degree distributions (large values of , relatively large values of , and ). In particular, these distribution are power-laws. In fact, as shown in. Fig. 2(a)-(b), for the degree of the largest hub of the graph at time scales as and the fluctuations of the degree distributions scale as .

In Fig. 2(c) we report the exponent of the power law distribution of the two layers, whose value clearly depends on . In particular, for we recover , as in the standard single-layer linear preferential attachment. When then converges to . When is negative and close to zero, we observe a strange phenomenon, which is also responsible for the peaks in , , and shown in Fig. 1. In this region, as we increase the distribution becomes first more homogeneous (with a peak of for ) and then again more heterogeneous, up until for .

This apparently strange behavior can be explained by considering that for the two layers are competing, i.e. a node having high degree on one layer will tend to have small degree on the other layer. In this case, a small negative value of actually reduces the heterogeneity of the attachment probability distribution that we have for , allowing small-degree nodes (for which the effect of layer competition is mitigated by the fact that is negative and close to zero) to acquire more edges. Conversely, when the value of becomes smaller then local fluctuations start to play a fundamental role, and the distribution becomes more heterogeneous again.

Three typical examples of degree distributions obtained for different values of are reported in Fig.3(a), in particular for . When we observe a homogeneous distribution for small values of , and one node acquires a finite fraction of the edges, i.e. the network is condensed. For the degree distribution is a power-law with exponent  (9). Finally, for the degree distribution can be fitted by a power-law with exponent .

Interestingly, the value of also determines the sign and value of inter-layer degree correlations, as confirmed by the plot of the average degree , at layer , of nodes having degree at layer , shown in Fig. 3(b)-(d) (9). It is clear that, by tuning , one can obtain either positive (, and ) or negative () inter-layer degree correlations. Finally, in Fig. 3(e) we plot as a function of the value of the Kendall’s rank correlation coefficient computed on the degree sequences of the two layers (see Appendix). For we have disassortative inter-layer degree correlations (), meaning that a hub on one layer is a poorly-connected node on the other layer, while for the degrees of the two replicas of the same node are positively correlated ().

## Iii Mean-field approximation

The mean-field approach has been proven to be extremely good in making qualitative predictions on the degree distribution of growing network models. Here we report convincing evidence that this approach is not able to capture essential properties of the proposed non-linear growth model. In fact in this model stochastic effects are fundamental to describe the evolution of the system. As we will see in a moment, the conclusion of the mean-field theory is that the expected degrees of the same node in the two layers are equal, irrespective of the values of the two exponents and . However, such a conclusion is in clear disagreement with the results obtained by numerical simulation and reported, for instance, in Fig. 3, which indeed confirm that for some combinations of the exponents and the degrees of a node at the two layers can be negatively correlated.

In the mean-field approximation the degree of a node at time acquires a deterministic value equal to its average degree in the stochastic model. If we indicate by and by the average degree of node on layer 1 and on layer 2 respectively, the mean-field approximation assumes that the degree of a node in layer 1 is equal to , i.e. and similarly that the degree in layer 2 of node at time is given by . Since, in this approximation, the average number of links that at time a node acquires in layer 1 is given by

 mf(κi,ϕi)∑jf(κj,ϕj), (5)

while the average number of links that a node acquires at time in layer 2 is given by

 mf(ϕi,κi)∑jf(ϕj,κj), (6)

when and as in Eq.(3), the mean-field equations for and at large times read

 dκidt=καiϕβiCt,             dϕidt=ϕαiκβiCt, (7)

with the constant to be self-consistently determined as

 C=limt→∞∑ti=1καiϕβimt. (8)

Assuming that is a constant, the Eqs. can be rewritten as

 κβ−αidκidlnt=1C(κiϕi)β=ϕβ−αidϕidlnt. (9)

Therefore we find for

 d[κβ−α+1i−ϕβ−α+1i]dlnt=0, (10)

while we have for

 dlnκi−lnϕidlnt=0. (11)

Therefore, if we consider the initial conditions the mean-field approach implies always . Inserting this solutions in the Eqs. we get

 dκidt=κα+βiCt (12)

yielding the solution

 κi(t)=m(tti)1/C (13)

for , and the solution

 κi(t)=[m1−(α+β)+1−α−βCln(tti)]1/(1−(α+β)) (14)

for .

For we observe a singularity in the solution for indicating the fact that the self-consistent equation for cannot be satisfied. By studying the master equation we will show that for we observe a condensation phase transition. Starting form the solution given by Eq. (13) and Eq.(14) the predicted degree distribution is scale free with power-law exponent for and a Weibull distribution for .

Overall we can say that the mean-field approach provides a solution that reflects the symmetry of the model in the two layers. Nevertheless this approach mostly fails in characterizing the correlations between the degrees of the same node in different layers. As we said before, the behavior of the model and its predictions are not supported by the simulations because the dynamics of the model is strongly affected by stochasticity and noise. In particular we argue that the strong deviations from the mean-field behavior that we observe in the simulations are due to the fundamental role played by stochastic effects on the degrees of the nodes recently arrived in the network. In fact these nodes will have a small degree in both layers and the fluctuations on these quantities will strongly affect the linking probability distribution.

## Iv Master equation

More theoretical insights about the model defined by Eq. (2) and Eq. (3) come from the solution of the master equation of the system, which accounts for the expected number of nodes with links in layer 1 and links in layer 2. Let us consider for simplicity the case . The master equation needs to take into account that at any time one of the following events can occur:

1. The number of nodes with degree in layer 1 and degree in layer 2 increases by one if the new node links in layer 1 but not in layer 2 to a node of degree in layer 1 and degree in layer 2.

2. The number of nodes with degree in layer 1 and degree in layer 2 increases by one if the new node links in layer 2 but not in layer 1 to a node of degree in layer 1 and degree in layer 2.

3. The number of nodes with degree in layer 1 and degree in layer 2 increases by one if the new node links in layer 1 and also in layer 2 to a node of degree in layer 1 and degree in layer 2.

4. The number of nodes with degree in layer 1 and degree in layer 2 decreases by one if the new node links in layer 1 to a node of degree in layer 1 and degree in layer 2.

5. The number of nodes with degree in layer 1 and degree in layer 2 decreases by one if the new node links in layer 2 to a node of degree in layer 1 and degree in layer 2.

Moreover, for and the average number of nodes with degree in layer 1 and degree in layer 2 increases by one at each time step, since the newly arrived node has degrees and . Taking into account all these possibilities, we can write the master equation as

 Nk,q(t+1) = Nk,q(t)+δk,mδq,m+f(k−1,q)M(t)(1−f(q,k−1)M(t))Nk−1,q(t)(1−δk,m) (15) +f(q−1,k)M(t)(1−f(k,q−1)M(t))Nk,q−1(t)(1−δq,m)−[f(q,k)+f(k,q)M(t)−f(k,q)f(q,k)[M(t)]2]Nk,q(t) +f(q−1,k−1)f(k−1,q−1)[M(t)]2Nk−1,q−1(t)(1−δq,m)(1−δk,m)

where , , is the Kronecker delta and is given by:

 M(t)=∑k,qf(k,q)Nk,q(t)=∑k,qf(q,k)Nk,q(t) (16)

Eq. (15) can be solved by using techniques similar to those adopted for single-layer networks or for multiplex networks with linear or semi-linear attachment kernels (35); (37); (9). In particular, by solving the master equation we obtain an analytical explanation for the appearance of a condensed phase. In fact, the master equation depends on the quantity which satisfies, in the thermodynamic limit , the relation

 M(t)=∑k,qf(k,q)Nk,q(t)=∑k,qf(q,k)Nk,q(t). (17)

Assuming that the normalization sum scales like , i.e. with constant, we can derive a recursive expression for . However, the hypothesis depends on the value of the exponents and in general is not satisfied. A deviation from this scaling indicates that in each layer we have a node that is grabbing an extensive number of links , i.e. we are in a condensed network phase.

### iv.1 Conditions for condensation

In order to show in which region of the phase space condensation occurs we first find a sufficient condition for condensation and then we will show that this condition is also a necessary one. We make use of the master equation to estimate , respectively for and , by considering (without loss of generality) the case . We observe that for at each time there are no vertices that in layer 1 have degree greater than , therefore the master equation given by Eq. (15) becomes

 Nk,1(k)=(k−1)αM(k−1)(1−(k−1)βM(k))Nk−1,1(k−1) (18)

But for large times . Moreover the fractions since only the first node of the network can have degree equal to the time . Consequently

 M(t)≥tα. (19)

Instead, if then at time there are no nodes that have at the same time degree in layer 1 greater than and degree in layer 2 greater then . In this case the master equation given by Eq. (15) becomes

 Nk,k(k)=(k−1)2(α+β)[M(k−1)]2Nk−1,k−1(k−1). (20)

Notice that the fractions since only the first node of the network can have degrees equal to , where is the time. Therefore we get that

 M(t)≥{tαifβ≤0tα+βifβ>0

This means that for and or for and

 M(t)≥tξ,∀ξ>1. (21)

This is a sufficient condition to have condensation, since in this case the expected number of nodes that at time have degrees scales with , i.e. . In fact, starting from the master equation, satisfies the following relation

 dN1,1(t)dt=−21M(t)N1,1(t)+1, (22)

where in writing this equation we have neglected higher order terms in . If Eq. (21) is satisfied, then the first term in the right-hand side of Eq. (22) is negligible and we have for large . This implies that the number of nodes with degrees different from is negligible, so that in this region we have a condensation phenomenon with few nodes grabbing an extensive number of connections.

Let us now show that the condition , and , is also necessary for condensation. Let us assume that we have a condensation of the links. In this scenario, we will have for one node with degree on layer 1, say node , and another node with degree on layer 2, say node ; conversely, for we will expect to have exactly one node, say node , having degrees . Since we have condensation then we can write an upper bound to , by taking into account only the contribution of the condensed nodes:

 M(t)≤{tαforβ≤0tα+βforβ>0 (23)

Putting Eq. (23) together with the lower bound given by Eq. (21) we find that satisfies the scaling

 M(t)≃{tαforβ≤0tα+βforβ<0

But we know that with , therefore we confirm that if the condensation transition occurs then either and or and . Therefore the condensation transition occurs only in the region or in the region , . In particular, for the same node will be the condensate node in both layers, while for the condensate node in one layer will not be the condensate node in the other layer. When the condensate nodes in the two layers might be either the same node or different nodes in different realizations.

### iv.2 Solution of the master equation in the non condensed phase

We consider now the master equation in the non condensed phase where with independent on , for . As we have seen above, this implies that the parameters satisfy the conditions: and or and . In this region of the phase space, we have always and therefore we can neglect the terms proportional to in the rate equation, finding the master equation for evolving multiplex in the non condensed phase, i.e.

 dNk,q(t)dt =Ak−1,qtNk−1,q(t)+Bk,q−1tNk,q−1(t)+ (24) −[Ak,q+Bk,qt]Nk,q(t)+δk,mδq,m

where we have put

 Ak,q=kαqβC,               Bk,q=qαkβC, (25)

and is a constant that can be determined self-consistently as

 C=limt→∞1t∑k,qkαqβNk,q(t). (26)

Assuming valid in the large time limit, we can solve for and we get

 Pm,q=(∏qj=mBk,j−11+Ak,j+Bk,j)Pm,mPk,q=∑qr=1(∏qj=r+1Bk,j−11+Ak,j+Bk,j)Ak−1,r1+Ak,r+Bk,rPk−1,r

These recursive equations can be used to solve numerically for the joint degree distribution of the degrees in the two layers, but unfortunately for there is no closed form analytical solution to these equations.

## V Numerical Results

The predictions obtained by solving the master equation of the model are in very good agreement with the phase diagram of the system obtained through simulations, reported in Fig. 4(a)-(b) (, , ). In these figures we show, for each value of the two parameters and , the corresponding values of (a) and (b), which allow us to visualize the two separate regions of the phase space. In region I the degree distribution is not condensed, while in region II we observe condensation as indicated by both a small value of and of . The shape of the boundary between the two regions agrees very well with the analytical prediction provided by the solution of the master equation in the thermodynamic limit (indicated by the solid lines in panel (a) and panel (b)). We notice that region I can be further divided into two separate sub-regions, according to the fact that the resulting degree distribution at each layer is homogeneous (region ) or heterogeneous (region ).

It is interesting to analyze the transition to condensation as a function of at fixed , i.e. region . In particular, we are interested in checking whether the degree distribution becomes a power–law before we reach the condensation transition (we already know that at the boundary of the condensation transition the degree distribution is a power-law, with an exponent which depends on , as discussed in Sec. II). Therefore, we analyzed the scaling of the degree of the largest hub and of the fluctuations of the degree distribution as a for increasing values of . The results corresponding to are reported in Fig.5(a) and Fig. 5(b). Notice that for we observe homogeneous degree distributions, i.e. no scaling of fluctuations with and a logarithmic scaling of , while for we have , which corresponds to . However, we observe that scales as already for , and in particular in the region . Also, in this region scales as , indicating that in region the degree distribution of each layer is a power-law.

Concerning the sign of inter-layer correlations, in Fig. 4(c) we report the Kendall’s correlation coefficient of the degree sequences at the two layers, where the two regions where inter-layer degree correlations are respectively positive (region ) and negative (region ) are separated by a solid black line. It is interesting to note that a multiplex can exhibit either positive or negative inter-layer correlations independently of the fact that its layers have homogeneous or heterogeneous distributions. While from the linking probabilities given by Eq. (2)-(3) we expect when , when the degrees of a node in the two layers tend to be negatively correlated. However, the interpretation of the phase diagram of for negative is less trivial, and the shape of the boundary between the regions and needs some explanation. In fact, when Eq. (3) implies that if a node has high degree in one layer, it will have low probability to acquire new links in the other layer, so that the degrees of the old nodes of the network will be negatively correlated. This is clear by looking at Fig. 4(d), which confirms that for the inter-layer degree correlations of older nodes are always negative.

However, for some values of the value of computed on the whole network could be positive, due to the presence of a large majority of younger nodes having small degrees on both layers (i.e. fickle nodes), whose values are mostly determined by stochastic fluctuations. In general, for large negative values of the fraction of fickle nodes is reduced, until it becomes zero for (the dashed line in Fig. 4(c) corresponds to the values of ), and in this case all the nodes have negative correlated degrees, resulting in a negative value of . We notice that the existence of two sub-regions in the phase diagram of for is not a finite-size effect, as confirmed by the results shown in Appendix Fig. A-1.

## Vi Multiplex cartography

The authors of Ref. (3) have recently introduced the concept of multiplex cartography, which is in the same spirit of the network cartography proposed by Guimerá and Amaral in Ref. (39); (40). Multiplex cartography is based on two measures, namely the Z-score of the overlapping degree of a node:

 z(oi)=oi−⟨o⟩σo (27)

where while and are the average and standard deviation of over all the nodes, and the multiplex participation coefficient:

 Pi=MM−1⎡⎢⎣1−M∑α=1⎛⎝k[α]ioi⎞⎠2⎤⎥⎦. (28)

The multiplex participation coefficient of a node characterizes its involvement in the layers of the multiplex. In fact, tends to if node has exactly the same degree on all the layers, while if node is isolated on all the layers but one. With respect to the Z-score of their overlapping degree, we distinguish hubs, for which , from regular nodes, for which . With respect to the multiplex participation coefficient, we call focused those nodes for which , mixed the nodes having and truly multiplex (or even simply multiplex) the nodes for which . The scatter-plot of and provides information about the patterns of participation across nodes of different degree classes, and gives insight about the different roles played by nodes.

In Fig. 6(a)-(c) we report the multiplex cartography diagrams for different values of (). It is interesting to notice that layer competition (i.e., ) enhances the variability of the multiplex cartography but produces multiplexes in which hubs are predominantly focused (top-left corner of the plots) while poorly-connected nodes are predominantly multiplex. Conversely, strong layer concordance (i.e., ) tends to produce multiplexes in which nodes belong to just a few different classes, i.e. either multiplex hubs or multiplex nodes.

## Vii Mixed correlations

Since the combination of and allows to produce multiplex graphs having either assortative or disassortative intra-layer degree-degree correlations and positive, null or negative inter-layer degree correlations, it is interesting to look at the combination of intra-layer and inter-layer correlations. In particular, we might ask whether a node being a hub on layer is preferentially connected on layer with other hubs or instead with leaves. So in general we can be interested in assessing whether:

1. a hub tends to be connected with other hubs or to poorly-connected nodes (intra-layer correlations);

2. a hub on one layer tends to be either a hub or a poorly-connected node in the other layer (inter-layer correlations);

3. a hub in one layer has neighbors in the other layer who are connected either to other hubs or poorly-connected nodes(type-1 mixed correlations).

4. the neighbors of a hub in one layer are either hubs or poorly-connected nodes in the other layer (type-2 mixed correlations).

We measure type-1 mixed correlations using the quantity:

 knn(q)=∑kP(k=ki|q=qi)1ki∑ja[1]ijkj (29)

which is the average degree of first neighbors on layer of a node having degree at layer . Similarly, we can define the dual quantity:

 qnn(k)=∑qP(q=qi|k=ki)1qi∑ja[2]ijqj (30)

If the plot of is an increasing (decreasing) function of , then we say that the mixed correlations of layer with respect to layer are positive (negative), or assortative (disassortative).

Type-2 mixed correlations can be quantified through the following expression:

 ˜qnn(k)=1Nk∑iδ(ki,k)1ki∑ja[1]ijqj (31)

which corresponds to the average degree at layer of the neighbors on layer of a node having degree on layer , and by the dual expression:

 ˜knn(q)=1Nq∑iδ(qi,q)1qi∑ja[2]ijkj (32)

Here (resp. ) indicates the number of nodes having degree equal to (resp. ) on layer (resp. on layer ). In Fig. 6(d)-(f) we show the intra- inter- and mixed correlation patterns obtained for several values of . Interestingly, for different values of the parameters one obtains different intra- and inter-layer correlation patterns, but also assortative or disassortative mixed correlations.

## Viii Distance and interdependence

Despite the main focus of the present work is on the properties of degree distribution and inter-layer degree correlations, we have also explored the distribution of shortest path length and the actual organization of shortest paths in the multiplex as a function of the two parameters and . It is important to notice that in a multiplex network the shortests paths between any pair of nodes are not limited to just one layer but can instead span both layers. Therefore, aside with the classical measure of characteristic path length:

 ⟨l⟩=1N(N−1)∑i∑j

which is just the average over all possible pairs of nodes of the distance between node and node , we also computed the multiplex interdependence(25); (9):

 λ=1N∑i∑j≠iψijσij. (34)

The quantity is the average ratio between the number of shortest paths between node and node which use edges lying on both layers and the total number of shortest paths between and in the multiplex. When then almost all shortest paths run in just one layer, while at the other extreme all shortest paths use edges in both layers.

In Fig. 7(a)-(b) we report the value of and for a synthetic multiplex of nodes as a function of and . Notice that the behaviour of closely mirrors that of the participation ratio reported in Fig. 4. As expected, the characteristic path length is smaller in the condensed phase, due to the presence of condensed nodes which are connected to virtually all the other nodes, and is larger in the non-condensed phase. The behavior of is more interesting. In fact, the non-condesed phase is characterised by a relatively high interdependence, and its value (which is almost always confined in the interval ) does not heavily depend on the actual value of and . In the condensed phase, instead, we spot two different sub-regions. For we have , which is expected since in this regime the same node is the condensed one on both layers, and half of the shortest paths can indeed run on both layers. When the condensed nodes on the two layers are distinct, so that all the shortest path run on just one layer, i.e. through the condensed node of that layer, and consequently .

## Ix Model calibration

Here we discuss the possibility of calibrating the non-linear preferential attachment model defined in Eq. (3), i.e. of choosing appropriate values of and , in order to reproduce some of the structural properties of a real-world multiplex network. As an example, we consider two 2-layer multiplex networks constructed from the IMDB costarring multi-layer network data set described in Ref. (8). In this data set, each layer corresponds to a different movie genre. The first multiplex consists of the nodes (actors) who have acted both in Adult and in Western movies, while the second one includes actors who have starred both in Adult and Mystery movies. Both systems are characterized by negative inter-layer degree correlations ( and , respectively). Let us now imagine that we want to set the values of and in order to construct a synthetic network having the same inter-layer degree correlation pattern of each of the two multiplexes.

In Fig. 8 we report the curves in the - plane corresponding to the values of measured in the Adult-Western and in the Adult-Mystery multiplexes. All the points of the solid red curve are pairs of values which produce a 2-layer multiplex with (Adult-Mystery), while the pairs indicated by the dashed blue line correspond to (Adult-Western). We notice that each of these pairs of parameters produces a synthetic multiplex network having exaclty the same value of , but in general different values of characteristic path length and interdependence . In the same plot we show, for each network, the range of values which guarantee, respectively, a value of or a value of compatible with those observed in the real multiplexes. This example shows that, despite being elegant and analytically solvable, the non-linear attachment model cannot reproduce, at the same time, all the structural properties of real-world networks. This fact suggests that the preferential attachment mechanism is just one among several ingredients responsible for the formation of multiplex networks.

## X General models for M layers

The model defined in Eq. (2) and Eq. (3) can be generalized to the case of multiplexes with layers in at least three different ways. We review them in the following, and for the first two generalization we also give a sketch of the the derivation of the conditions for condensation.

### x.1 One vs. All

A simple extension would be to consider an attaching function

 fa(→k)=(k[a])α∏b≠a(k[b])β (35)

which says that the probability for a new node to attach on layer to a node having degree depends on the power of and on the product of the powers of the degrees of the same node at the other layers . In this case, each layer can either compete with all the others () or cooperate with all of them (), and the behavior of any two layers will be exactly the same of that studied in the previous Sections.

Following a similar approach used to determine the condensation phase diagram for the model of two layers, it is easy to show that the condensation occurs in a multiplex of layers satisfying the attachment rule given by Eq. (35) under the following conditions

 α>1, β<0      or        α+β(M−1)>1, β>0. (36)

In the case and there are nodes in which the condensation occurs, exactly one in each of the layers. Each of these condensed nodes has degree in exactly one layer (say layer ), while its degree on all the other layers is equal to . Instead for and the condensation occurs on a single node that has degree in all the layers of the multiplex.

### x.2 Two groups of layers

Another possible extension of Eq. (3) to the case of layers considers layers divided into two groups, say and . We denote by the group of layers to which layer belongs, and by their cardinality and . We define the attaching function:

 f[a](→k)=∏b:Γ(b)=Γ(a)(k[b])α∏b:Γ(b)≠Γ(a)(k[b])β (37)

meaning that the probability for a new node to connect on layer with a node of degree depends on the product of the -power of the degrees of the destination node at all layers belonging to the same group of layer multiplied by the product of the powers of the degrees of the destination node at all layers belonging to the other group. Also in this case the dynamics of pairwise relationships between layers belonging to different groups is similar to that observed in the 2-layer case discussed in the previous Sections. Though, the phase diagram is not exactly the same. In fact, the condensation could occur either only on the layers belonging to , or only on the layers belonging to or on all the layers at the same time.

Following a similar approach used to determine the condensation phase diagram for the model of two layers, it is possible to show that the condensation occurs in a multiplex of layers satisfying the attachment rule given by Eq. (37) under the following conditions

 αM1>1, β<0or αM2>1, β<0or ξ=α(M21+M22)+2βM1M2>2, β>0. (38)

In the case and there is a node in which the condensation occurs. This node has all the degrees in layers given by . Similarly for one node becomes the condensate in layers if . If both and these two nodes where the condensation occurs coexist in the multiplex and are distinct. Instead for and the condensation occurs on a single node that have all the degrees in every layer .

### x.3 More complex layer interconnections

Finally, we consider the case in which the degree of a node at each single layer might interact with the degree of the same node at any other layer by means of a power or . We define a interaction matrix , such that if layer interacts with layer through the exponent , while if interacts with through the exponent . Notice that in general , i.e. is not necessarily symmetric. In this case the attaching function reads:

 fa(→k)=∏b:ca,b=1(k[b])α∏b:ca,b=−1(k[b])β (39)

This model is very general and allows a pretty rich interplay between the degree distributions of the layers. In this case the conditions for condensation depend on the structure of the interconnection matrix , and the derivation is left as a future work.

## Xi Conclusions

In this Article we have introduced a general class of non-linear models to grow multiplexes which display a rich variety of behaviors, including the appearance of positive, null and negative inter-layer degree correlations and the transition to a condensed phase. We have shown that the model is highly sensitive to stochasticity, so that the mean-field approach, which has been fundamental to study growth processes on single-layer networks, fails to give account for some of its most interesting properties. Conversely, the solution of the master equation of the system gives some general theoretical insights which will certainly prove to be a useful guide in the exploration of real-world multiplex networks.

We would like to stress the fact that the class of growth models proposed in this work includes only some of the ingredients which might be responsible for the formation and evolution of multi-layer networks. As a matter of fact, real networked systems rarely evolve only by the addition of new nodes and edges at discrete time-steps. Depending on the structure and function of the multiplex system under study, nodes can also disappear and re-join the network again, with a different number of edges on each layer, and edges might be rewired, severed and re-created, sometimes according to the state of some dynamical processes occurring on the network. Also, the arrival and departure of nodes and the creation and rewiring of edges might be affected by different levels of topological and temporal correlations. All these ingredients should be taken into account for a more accurate modelling of real-world multiplex networks, and this will certainly be the subject of future research in this novel field of network science. Nevertheless we believe that, despite the few simplifying assumptions introduced to make the model analytically tractable, the present work clearly points out that multiplex networks are indeed characterized by new, additional and somehow unexpected levels of complexity, and that the multiplex perspective might reveal interesting aspects of real-world complex systems which have remained unnoticed until now.

###### Acknowledgements.
V.N. and V.L. acknowledge support from the Project LASAGNE, Contract No.318132 (STREP), funded by the European Commission. M.B. is supported by the FET-Proactive project PLEXMATH (FP7-ICT-2011-8; grant number 317614) funded by the European Commission. This research utilised Queen Mary’s MidPlus computational facilities, supported by QMUL Research-IT and funded by EPSRC grant EP/K000128/1.

## Appendix A Inter-layer correlations

### a.1 Coefficients to quantify inter-layer degree correlations

To detect and quantify the presence of inter–layer degree correlations we have evaluated the Pearson’s linear correlation coefficient , the Spearman rank correlation coefficient and the Kendall’s rank correlation coefficient of the degree distributions at the two layers. If we denote as and the degrees of node respectively at layer and layer , the Pearson’s correlation coefficient of the two degree sequences is defined as:

 r=⟨kq⟩−⟨k⟩⟨q⟩σkσq (A-40)

where the averages are taken over all the nodes in each layer, and the are the corresponding standard deviations. Similarly, if we denote by the rank of the degree of node on the first layer, and by the rank of the degree of node on the second layer, the Spearman’s correlation coefficient is defined as:

 ρ=∑i(r(ki)−¯¯¯¯¯¯¯¯¯¯r(k))(r(qi)−¯¯¯¯¯¯¯¯¯r(q))√∑i(r(ki)−¯¯¯¯¯¯¯¯¯¯r(k))2∑i(r(qi)−¯¯¯¯¯¯¯¯¯r(q))2 (A-41)

where and are the averages respectively at layer and layer .

If we consider node and and we call the ranking induced at each layer by the degree sequence, we say that is a concordant pair with respect to if the ranks of the two nodes agree, i.e. if both and or both and . If a pair of nodes is not concordant, then it is said discordant. The Kendall’s coefficient measures the correlation between two rankings by looking at concordant and discordant pairs:

 τ=nc−nd√(n0−n1)(n0−n2) (A-42)

where is the number of concordant pairs, is the number of discordant pairs, and is the total possible number of pairs in a set of elements. The terms and account for the presence of rank degeneracies. In particular, let us suppose that the first ranking has tied groups, i.e. sets of elements such as all the elements in one of this set have the same rank. If we call the number of nodes in the tied group, then is defined as:

 n1=m∑i=112ui(ui−1).

Similarly, is defined as follows:

 n2=n∑j=112vj(vj−1)

where we have made the assumption that the second ranking has tied groups, and that the tied group has elements.

The Kendall’s coefficient is equal to when the rankings induced by the degree sequence at each layer are perfectly concordant, while if one of the two rankings is exactly the opposite of the other.

### a.2 Pearson’s coefficient in the non-condensed phase

Using the master equation we can derive several relations between the moment of the degree distribution at long times. In particular it can be shown that for