Nonlinear growth and condensation in multiplex networks
Abstract
Different types of interactions coexist and coevolve 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 coevolve through a set of nonlinear preferential attachment rules. We show, both numerically and analytically, that by tuning the level of nonlinearity 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.75k, 89.75.Fb, 89.75.HcVarious complex systems are well described by multiplex networks of nodes connected through links of distinct types, which constitute separate yet coevolving and interdependent layers Cardillo2013_Emergence (); Cardillo2013_Airports (); battiston (); Thurner (). 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 arenas (); Aleva:2013 (); Bianconi2013 (). In these systems, links of different types are intertwined in nontrivial 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 Nicosia2014_corr () or, conversely, the hubs might tend to be the same across different layers Nicosia2013 (); kim (). 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 Thurner (); Cardillo2013_Emergence (); Bianconi2013 (); Halu2013 (). Some recent studies have focused on dynamical processes on multiplexes, including percolation Buldyrev:2010 (); Cellai2013 (); Baxter2013 (); Min2013_robustness (); Bianconi2014_percolation (), diffusion Gomez2013 (); Sole2013 (); DeDomenico2013_RW (); Halu2013_Pagerank (), spreading Saumell2012 (); Granell2013 (); Min2013_epidemic (); Cozzo2013_epidemic (), traffic Morris:2012 (), cascades Morris:2013 (), and cooperation GomezGardenes2012 (); Jiang2013 (); Wang2013 (), and a few recent works have suggested that degree correlations Parshani2010 (); Lee2012 (); Valdez2013 () as well as overlap of links Cellai2013 (); Li (); Hu () 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 interlayer 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 nullmodels for multiplexes arenas (); Aleva:2013 (); Bianconi2013 (), while some other focus on capturing the nonequilibrium nature of multiplexes and on providing possible physical explanations for their formation Nicosia2013 (). However, until now, all the existing models for growing multiplexes with homogeneous and heterogeneous degree distributions allow for positive interlayer degree correlations only.
In this Article, we propose and study a general growth model of multiplex networks based on a nonlinear 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 noncondensed regime the model can generate multiplexes with homogeneous or heterogeneous degree distributions, having either positive or negative interlayer degree correlations. Finally we notice that in the multilayer version of nonlinear preferential attachment the structure of the network dramatically depends on fluctuations, and that the meanfield approach, which was fundamental to understand network growth in singlelayer 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 nonlinear preferential attachment growth models for multilayer networks focusing, as an example, on the case of a 2layer 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 ba (). In Sec. III we derive a meanfield solution for the proposed class of models, and we show that the meanfield approximation fails to account for most of the observed structural properties of the multiplex, in particular regarding the possibility to obtain negative interlayer 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 battiston (), 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 realworld 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:
(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
(2) 
In the context of singlelayer networks, nonlinear attachment kernels of the form , with , have been introduced in Ref. krapivsky (), as a generalization of linear preferential attachment models ba (). 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 highdegree nodes. We adopt the general expression
(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 wellconnected (poorly connected) on both layers. Conversely, if and (resp. and ), a new node will be preferentially linked, in layer 1, with nodes which are wellconnected (resp. poorlyconnected) in layer 1 and poorly connected (resp. wellconnected) 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 Seminonlinear 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. Nicosia2013 (). In this case, the degree distribution in each layer is a power law with , and the multiplex exhibits positive interlayer 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 interlayer 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 Derrida:1987 ()
(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 powerlaws. 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 singlelayer 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 smalldegree 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 powerlaw with exponent Nicosia2013 (). Finally, for the degree distribution can be fitted by a powerlaw with exponent .
Interestingly, the value of also determines the sign and value of interlayer 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) Nicosia2013 (). It is clear that, by tuning , one can obtain either positive (, and ) or negative () interlayer 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 interlayer degree correlations (), meaning that a hub on one layer is a poorlyconnected node on the other layer, while for the degrees of the two replicas of the same node are positively correlated ().
Iii Meanfield approximation
The meanfield 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 nonlinear 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 meanfield 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 meanfield 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 meanfield 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
(5) 
while the average number of links that a node acquires at time in layer 2 is given by
(6) 
when and as in Eq.(3), the meanfield equations for and at large times read
(7) 
with the constant to be selfconsistently determined as
(8) 
Assuming that is a constant, the Eqs. can be rewritten as
(9) 
Therefore we find for
(10) 
while we have for
(11) 
Therefore, if we consider the initial conditions the meanfield approach implies always . Inserting this solutions in the Eqs. we get
(12) 
yielding the solution
(13) 
for , and the solution
(14) 
for .
For we observe a singularity in the solution for indicating the fact that the selfconsistent 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 powerlaw exponent for and a Weibull distribution for .
Overall we can say that the meanfield 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 meanfield 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:

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.

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.

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.

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.

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
(15)  
where , , is the Kronecker delta and is given by:
(16) 
Eq. (15) can be solved by using techniques similar to those adopted for singlelayer networks or for multiplex networks with linear or semilinear attachment kernels krapivsky (); Dorogovtsev (); Nicosia2013 (). 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
(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
(18) 
But for large times . Moreover the fractions since only the first node of the network can have degree equal to the time . Consequently
(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
(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
This means that for and or for and
(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
(22) 
where in writing this equation we have neglected higher order terms in . If Eq. (21) is satisfied, then the first term in the righthand 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:
(23) 
Putting Eq. (23) together with the lower bound given by Eq. (21) we find that satisfies the scaling
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.
(24)  
where we have put
(25) 
and is a constant that can be determined selfconsistently as
(26) 
Assuming valid in the large time limit, we can solve for and we get
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 subregions, 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 powerlaw, 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 powerlaw.
Concerning the sign of interlayer 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 interlayer 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 interlayer 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 interlayer 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 subregions in the phase diagram of for is not a finitesize effect, as confirmed by the results shown in Appendix Fig. A1.
Vi Multiplex cartography
The authors of Ref. battiston () 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. Guimera05 (); GuimeraNature05 (). Multiplex cartography is based on two measures, namely the Zscore of the overlapping degree of a node:
(27) 
where while and are the average and standard deviation of over all the nodes, and the multiplex participation coefficient:
(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 Zscore 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 scatterplot 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 (topleft corner of the plots) while poorlyconnected 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 intralayer degreedegree correlations and positive, null or negative interlayer degree correlations, it is interesting to look at the combination of intralayer and interlayer 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:

a hub tends to be connected with other hubs or to poorlyconnected nodes (intralayer correlations);

a hub on one layer tends to be either a hub or a poorlyconnected node in the other layer (interlayer correlations);

a hub in one layer has neighbors in the other layer who are connected either to other hubs or poorlyconnected nodes(type1 mixed correlations).

the neighbors of a hub in one layer are either hubs or poorlyconnected nodes in the other layer (type2 mixed correlations).
We measure type1 mixed correlations using the quantity:
(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:
(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).
Type2 mixed correlations can be quantified through the following expression:
(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:
(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 interlayer 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 interlayer 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:
(33) 
which is just the average over all possible pairs of nodes of the distance between node and node , we also computed the multiplex interdependenceMorris:2012 (); Nicosia2013 ():
(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 noncondensed phase. The behavior of is more interesting. In fact, the noncondesed 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 subregions. 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 nonlinear 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 realworld multiplex network. As an example, we consider two 2layer multiplex networks constructed from the IMDB costarring multilayer network data set described in Ref. Nicosia2014_corr (). 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 interlayer 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 interlayer 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 AdultWestern and in the AdultMystery multiplexes. All the points of the solid red curve are pairs of values which produce a 2layer multiplex with (AdultMystery), while the pairs indicated by the dashed blue line correspond to (AdultWestern). 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 nonlinear attachment model cannot reproduce, at the same time, all the structural properties of realworld 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 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
(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
(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:
(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 2layer 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
(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:
(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 nonlinear models to grow multiplexes which display a rich variety of behaviors, including the appearance of positive, null and negative interlayer degree correlations and the transition to a condensed phase. We have shown that the model is highly sensitive to stochasticity, so that the meanfield approach, which has been fundamental to study growth processes on singlelayer 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 realworld 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 multilayer networks. As a matter of fact, real networked systems rarely evolve only by the addition of new nodes and edges at discrete timesteps. Depending on the structure and function of the multiplex system under study, nodes can also disappear and rejoin the network again, with a different number of edges on each layer, and edges might be rewired, severed and recreated, 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 realworld 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 realworld 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 FETProactive project PLEXMATH (FP7ICT20118; grant number 317614) funded by the European Commission. This research utilised Queen Mary’s MidPlus computational facilities, supported by QMUL ResearchIT and funded by EPSRC grant EP/K000128/1.Appendix A Interlayer correlations
a.1 Coefficients to quantify interlayer 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:
(A40) 
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:
(A41) 
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:
(A42) 
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:
Similarly, is defined as follows:
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 noncondensed 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 we have
(A43)  
In particular we have,
(A44) 
Therefore the Pearson’s linear coefficient , defined as
(A45) 
with can be also written as
(A46) 
Appendix B Stability of negative interlayer correlations
It is interesting to investigate whether the existence of two subregions in the phase diagram of for (see Fig. 4(c) in the main text) is indeed due to finitesize effects or not. To this aim, we computed for networks whose size varied across three orders of magnitude. The results are reported in Fig. A1, for different values of and . As made clear by the figures, the values of measured for a certain pair do not depend on the size of the multiplex, and therefore the shape of the region is not an artifact due to finite size effects.
References
 (1) A. Cardillo et al., Sci. Rep. 3, 1344 (2013).
 (2) A. Cardillo et al., Eur. Phys. J. Spec. Top. 215, 2333 (2013).
 (3) F. Battiston, V. Nicosia and V. Latora, Phys. Rev. E 89, 032804 (2014).
 (4) M. Szell, R. Lambiotte, S. Thurner, Proc. Natl. Acad. Sci., USA, 107, 13636 (2010).
 (5) G. Bianconi, Phys. Rev. E 87, 062806 (2013).
 (6) M. Kivela, et al, arXiv:1309.7233 (2013).
 (7) M. De Domenico et al., Phys. Rev. X 3, 041022 (2013).
 (8) V. Nicosia, V. Latora, arxiv:1403.1546 (2014).
 (9) V. Nicosia, G. Bianconi, V. Latora and M. Barthelemy, Phys. Rev. Lett. 111, 058701 (2013).
 (10) J, Y. Kim and K.I. Goh, Phys. Rev. Lett. 111, 058702 (2013).
 (11) A. Halu, S. Mukherjee, G. Bianconi, Phys. Rev. E 89, 012806 (2014).