aa r X i v : . [ phy s i c s . s o c - ph ] J u l Growing multiplex networks
V. Nicosia, G. Bianconi, V. Latora,
1, 2 and M. Barthelemy School of Mathematical Sciences, Queen Mary University of London, Mile End Road, E1 4NS, London (UK) Dipartimento di Fisica e Astronomia, Universit`a di Catania and INFN, 95123 Catania, Italy Institut de Physique Th´eorique, CEA, CNRS-URA 2306, F-91191, Gif-sur-Yvette, France
We propose a modelling framework for growing multiplexes where a node can belong to differentnetworks. We define new measures for multiplexes and we identify a number of relevant ingredientsfor modeling their evolution such as the coupling between the different layers and the distribution ofnode arrival times. The topology of the multiplex changes significantly in the different cases underconsideration, with effects of the arrival time of nodes on the degree distribution, average shortestpath length and interdependence.
PACS numbers: 89.75.Fb, 89.75.Hc and 89.75,-k
Many different physical, biological and social systemsare structured as networks, and their properties are now,after a decade of efforts, well understood [1–4]. How-ever, a complex network is rarely isolated, and some ofits nodes could be part of many graphs, at the sametime. Examples include multimodal transportation net-works [5, 6], climatic systems [7], economic markets [8],energy-supply networks [9] and the human brain [10].In these cases, each network is part of a larger systemin which a set of interdependent networks with differentstructure and function coexist, interact and coevolve. Sofar network scientists have investigated these systems bylooking at one type of relationship at a time, e.g., by ana-lyzing collaboration networks and email communicationsas separate graphs. However, the structural properties ofeach of these networks and their evolution can depend ina non-trivial way on that of other graphs to which theyare interconnected. Consequently, these systems are bet-ter represented as multiplexes , i.e. graphs composed by M different layers in which the same set of N nodes canbe connected to each other by means of links belongingto M different classes or types. Despite some early at-tempts in the field social network analysis [11], the char-acterization of multiplexes is still in its infancy, mainlydue to the lack of multiplex data. However, some re-cent works have already proposed suitable extensions tomulti-layer graphs of classic network metrics and mod-els [12–14]. Preliminary results show that multiplexicityhas important consequences for the dynamics of processesoccurring in real systems, including routing [12, 15], dif-fusion [16], cooperation [17], election models [18], andepidemic spreading [19]. Nowadays, an increasing num-ber of new data sets of multiplex systems, e.g. com-ing from large online social networks [20, 21], tradingnetworks [22] and human neuroimaging techniques [23],are rapidly becoming available and demand for adequatemodels to understand their structure and evolution.In this Letter we propose and study a generic modelof multiplex growth, inspired by classical models basedon preferential attachment, in which the probability fora newly arrived node to establish connections to existing nodes in each of the layers of a multiplex is a function ofthe degree of other nodes at all layers. We define two newmetrics to characterize the structure of multiplexes andwe study the effect of different attachment rules and theimpact of delays in the arrival of nodes at different layerson the structure of the resulting network. We provideclosed forms for both the degree distributions at eachlayer and the inter-layer degree-degree correlations, andwe show how different attachment kernels can change thedistributions of distances and interdependence.More precisely, a multiplex is a set of N nodes whichare connected to each other by means of edges belong-ing to M different classes or types. We represent eachclass of edges as a separate layer , and we assume thata node i of the multiplex consists of M replicas, one foreach layer. We denote by V [ α ] the set of nodes in layer α and by E [ α ] the set of all the edges of a given type α .An M -layer multiplex is therefore fully specified by thevector A = [ A [1] , A [2] , . . . , A [ M ] ], whose elements are theadjacency matrices A [ α ] = { a [ α ] ij } , where a [ α ] ij = 1 if node i and node j are connected by an edge of type α , whereas a [ α ] ij = 0 otherwise. We denote by k [ α ] i = P j a [ α ] ij the de-gree of node i at layer α , i.e. the number of edges of type α of which i is an endpoint, and by k i the M -dimensionalvector of the degrees of the replicas of i . In general, thedegrees of the replicas of i are distinct, and some replicascan also be isolated (i.e. k [ α ] i = 0 for some value of α ).In the following we consider all the edges at all layers tobe undirected and unweighted. As in the case of classi-cal ‘singlex’ graphs, we can characterize each layer α ofa multiplex by studying the degree distribution P ( k [ α ] ),and the joint-degree distribution P ( k [ α ] , k ′ [ α ] ). However,we are interested here in the structural properties of themultiplex as a whole, so we propose to quantify the corre-lations between the degrees of replicas of the same nodeat two different layers α and α ′ , by constructing the inter-layer joint-degree distributions P ( k [ α ] , k [ α ′ ] ), or the con-ditional degree distributions P ( k [ α ′ ] | k [ α ] ). In particular,we can look at the projection of the conditional distribu-tion obtained by considering the average degree ¯ k [ α ′ ] atlayer α ′ of nodes having degree k [ α ] at layer α :¯ k [ α ′ ] ( k [ α ] ) = X k [ α ′ ] k [ α ′ ] P ( k [ α ′ ] | k [ α ] ) (1)By plotting this quantity as a function of k [ α ] we can de-tect the presence and the sign of degree correlations be-tween the two layers. For a multiplex with no correlationsbetween layers α and α ′ we expect ¯ k [ α ′ ] ( k [ α ] ) = h k [ α ′ ] i and ¯ k [ α ] ( k [ α ′ ] ) = h k [ α ] i . If ¯ k [ α ′ ] ( k [ α ] ) increases with k [ α ] we say that the degrees of the two layers have positive(assortative) correlations, while if ¯ k [ α ′ ] ( k [ α ] ) is a decreas-ing function of k [ α ] we say that the degrees on layer α and α ′ are anticorrelated (or disassortatively correlated).We notice that a similar concept of inter-network assor-tativity was already defined in Ref. [24] for the case ofinterdependent graphs, while the authors of Ref. [13] pro-posed to measure inter-layer assortativity by means of thePearson’s linear correlation coefficient of degrees [25]. -6 -4 -2 ~k -3.0 k ~ m(k+2)/(m+1)
10 100 100010 -6 -4 -2 ~ k -3.0 k ~ k b)a)kP(k) k(k)k FIG. 1. (color online)
Synchronous linear attachment.
Panel a-b: the degree distribution P ( k ) (left) and the projec-tion ¯ k ( k ) of the inter-layer degree correlations (right) closelyfollow the theoretical curves (solid black lines) and are rela-tively insensitive to the coupling matrix. In addition to the assortativity, we can also charac-terize the ‘multiplex reachability’ of a node i , e.g., bycomputing the average distance L i from i to any othernode of the multiplex, and comparing this average dis-tance with that measured on each layer separately. Thepresence of more than one layer in a multiplex producesan increase in the number of available paths, so that thedistance between two nodes of a multiplex will be, in gen-eral, smaller than or at most equal to that measured oneach layer separately. A better measure to quantify thevalue added by the multiplexicity to the reachability ofnodes is the interdependence [12] which for a node i isdefined by λ i = X j ∈ N j = i ψ ij σ ij (2) where ψ ij is the number of shortest paths between node i and node j which use edges lying on more than one layer,while σ ij is the total number of shortest paths between i and j in the multiplex. The interdependence of a mul-tiplex is computed as the average node interdependence λ = 1 /N P i λ i with λ ∈ [0 , λ is close to zero, thenmost of the shortest paths among nodes lie on just onelayer, while if λ is close to 1 the majority of the shortestpaths exploit more than one layer.The few models of multiplexes proposed so far arebased on the juxtaposition of random graphs [13]. How-ever, networks usually result from a growing process con-sisting in the addition of nodes and edges over time. Forthis reason, we introduce here a model of growing mul-tiplex networks. Most of the classical growing modelsfor single-layer networks start from an initial connectedgraph with m nodes and assume that new nodes arrivein the graph one by one, carrying m edge stubs, and con-nect with other existing nodes according to a prescribedattachment rule. In that case, each node i has a unique arrival time t i , but in multiplexes, instead, each layercan exhibit a different edge-formation dynamics, and ingeneral the edges of the M replicas of a new node arenot created at the same time. For instance, a face-to-face interaction relationship is usually established beforetwo individuals become friend, while two locations areusually connected by a road before a direct railway linebetween them is constructed. Consequently, we assumethat a newly arrived node has exactly m stubs on eachlayer of the multiplex (in Appendix we briefly discuss thecase where m is a random variable), but the replica of anode i on layer α can connect its m stubs at a differenttime t [ α ] i . We denote by t i the vector of arrival timesof the replicas of node i . In order to make the modelanalytically tractable, we make two simplifying assump-tions. The first is that there exists a layer ¯ α so that t [¯ α ] i ≤ t [ α ] i ∀ i, ∀ α = ¯ α . This is equivalent to saying thata newly arrived node must first create its connections onlayer ¯ α before any of its replicas can create connectionson any other layer α = ¯ α . We call ¯ α the master layer (in Appendix we briefly discuss the case in which thisassumption does not hold, and each node can arrive firston any of the M layers of the multiplex). The secondassumption is that nodes arrive one by one on the mas-ter layer, at equal discrete time intervals t = { , , . . . , } .We label the nodes of a growing multiplex according tothe ordering induced by their arrival on the master layer.Without loss of generality, in the following we assumethat the the master layer is the first one, i.e. ¯ α = 1, andthat the arrival times of the replicas of node i have theform t [ α ] i = T ( t [1] i , ξ [ α ] ( τ )) (3)where T is a certain function of t [1] i and of the randomvariable ξ [ α ] ( τ ). By appropriately choosing T and ξ [ α ] ( τ ) k(t) L i P(L i ) λ i P( λ i )t δ δ
40% 45% 50% 55% 60% H i P(H i ) β =1.1 β =1.5 β =2.0 β =2.510 -6 -3 P(k) ~ k -3 -6 -3 a)b) c) d) e) FIG. 2. (color online)
Delayed linear attachment. (a)Degree distribution on the first layer. When β is close to 1,super-hubs appear. (b) The degree k ( t ) of the largest hub ofthe first layer as a function of time scales as ( t/t ) δ , where δ approaches 0 . β increases (insets). The value of β tunes the shape of the distribution of average shortest pathlengths (c) and node interdependence (d). (e) The percentageof times H i that the maximal-degree node in a shortest pathfrom node i belongs to the first layer. When β is small, super-hubs in the first layer are more abundant in the shortest paths. we can model different arrival behaviors, including a) si-multaneous arrival ( T = t [1] i ); b) power-law delayed ar-rival ( T = t [1] i + ξ ( τ ) and P ( ξ = τ ) = ( β − τ − β for τ ≥ β > i connects to m existing nodes in the master layer, accord-ing to a certain attachment rule. As in the preferentialattachment models [27], we assume that the attachmentprobability depends on the degree of a node. However,in a multiplex the probability for node i to connect tonode j on each layer α can depend not only on k [ α ] j butalso on the degrees of j ’s replicas on the other layersΠ [ α ] i → j = F [ α ] j ( k j ) P l F [ α ] l ( k l ) (4)For the sake of clarity and without loss of generality,we focus in the following on 2-layer multiplexes with α = 1 ,
2. We begin with the simplest case of linear at-tachment which is the natural extension of the Barab´asi-Albert model [27]. In this case, we consider that theprobability for a newborn node i to connect to an existingnode j on layer α is proportional to a linear combinationof the degrees of j at all layers. The attachment kernelscan then be expressed as (cid:20) F [1] [ k, q ] F [2] [ k, q ] (cid:21) = C (cid:20) kq (cid:21) = (cid:20) c [1 , c [1 , c [2 , c [2 , (cid:21)(cid:20) kq (cid:21) (5)where we use here and in the following the notations k [1] = k and k [2] = q [28]. The coefficients c [ r,s ] tune the dependence of the attachment probability atlayer r on the degrees of nodes at layer s . In the case of 2-layer multiplexes we can represent the set ofcoefficients C = { c [ r,s ] } using the compact notation { c [1 , , c [1 , , c [2 , , c [2 , } . The dynamics can be easilysolved in mean-field (see Appendix for details) and insome specific cases we can fully characterize the degreecorrelations within the two different layers by analyticallysolving the master-equation. If we denote by N k,q ( t ) thenumber of nodes having, at time t , degree k on the firstlayer and degree q on the second layer, and by Π [ α ] k,q theprobability that one of these N k,q ( t ) nodes acquires oneof the m new links on layer α at time t + 1, the masterequation can be written as [29] N k,q ( t + 1) = N k,q ( t ) + G − L (6)where G = m h Π [1] k − ,q N k − ,q ( t ) + Π [2] k,q − N k,q − ( t ) i + δ k,m δ q,m L = m h Π [1] k,q + Π [2] k,q i N k,q ( t )represent, respectively, the expected increase ( G ) and theexpected decrease ( L ) of N k,q at time ( t + 1). Assumingthat N k,q = tP ( k, q ) for large t , the solution of Eq. (6)is obtained by solving the corresponding recursive ex-pression (see Appendix for details). In the following wesummarize the master-equation solution in some partic-ularly interesting cases. First of all let us consider simul-taneous arrival of the nodes in the two layers. If we set C = { , , , } then the attachment probability at eachlayer will depend only on the degree of the nodes in thesame layer. In this case the degree distribution in thefirst layer reads [30, 31] P ( k ) = 2 m ( m + 1) k ( k + 1)( k + 2) , k > m (7)and the degree distribution in the second layer is identi-cally equal. This distribution goes as P ( k ) ∼ k − γ with γ = 3. If we solve the master-equation for the multi-plex evolution we obtain the analytical expression for theinter-layer joint degree probability P ( k, q ) P ( k, q ) = 2Γ(2 + 2 m )Γ( k )Γ( q )Γ( k + q − m + 1)Γ( m )Γ( m )Γ( k + q + 3)Γ( k − m + 1)Γ( q − m + 1)(8)The average degree ¯ k ( q ) at layer 1 of nodes having degree q at layer 2 reads: ¯ k ( q ) = m ( q + 2)1 + m (9)Notice that even if the two layers grow independently,the simultaneous arrival introduces non-trivial inter-layerdegree correlations. In fact, in the mean-field approach,the degree of a node on each layer increases over timeas k [ α ] i ( t ) = m (cid:16) t/t [ α ] i (cid:17) / (see Appendix for details), so -4 -2 P(k)
10 10010 k(k)
10 100 10 1003 4 L i P(L i ) t λ( t) λ i P( λ i ) {0,1,0,1}{0,1,1,0}{0,1,1,1}{1,0,0,1}{1,0,1,0}{1,0,1,1}{1,1,0,1}{1,1,1,0}{1,1,1,1} a)b)c) d) e) kk k FIG. 3. (color online)
Semi-linear attachment.
Degreedistributions (a) and inter-layer degree correlations (b) forsemi-linear attachment. (c) The distribution of the averageshortest path length from one node to all the other nodesheavily depends on the coupling pattern. Similarly, the in-terdependence of a node λ ( t ) is always a sublinear functionof the arrival time t but its shape depends on the couplingpattern at work (d). In general, older nodes have smallerinterdependence. (e) The coupling pattern also affects thedistribution of node interdependence. The smallest averageinterdependence is observed when the two layers are indepen-dent (yellow curve). that the degrees of the two replicas of a node i depend, forlarge t , only on their arrival time. If both replicas havethe same arrival time, i.e. t [1] i = t [2] i then the degree ofthe two replicas will be positively correlated. In Fig. 1 wereport the degree distribution and the values of ¯ k ( q ) fortwo coupling patterns, which are in good agreement withthe theoretical curves [32]. It is clear from the figure thatin the synchronous arrival case the shape of the couplingmatrix is actually not very relevant and that the value ofthe degree distribution exponent and strong assortativityare robust features of these multiplexes.If we consider a power-law delayed arrival time on thesecond layer, the results are significantly different. InFig. 2 we illustrate how the exponent of the delay distri-bution β affects the structure of the obtained multiplex.The bulk of the degree distributions are still power laws P ( k ) ∼ k − γ with γ = 3, but the shape of the far taildepends now on β : for small β , a few nodes are predom-inant and become super-hubs (as also shown in Fig. S-1in Appendix). The average shortest path and the inter-dependence are also significantly affected, as shown inFig. 2c-d. In particular, when β is closer to one the pres-ence of more predominant ‘old’ hubs lowers the averageshortest path and the inter-layer assortativity. Moreover,broader delays cause a lower participation of hubs of thesecond layer in shortest paths, as shown in Fig. 2e.So far we have considered the case of two scale-freegrowing networks, but it would be interesting to con-struct multiplexes in which a scale-free network is cou-pled to a network with a peaked degree distribution. Inthis respect, we introduce a semi-linear attachment ker- nel which allows to grow multiplexes in which the twolayers have different topological structures. The model isdefined as follows (cid:20) F [1] [ k, q ] F [2] [ k, q ] (cid:21) = C (cid:20) k (cid:21) (10)where C is still a 2 × C = { , , , } we can analytically solve the masterequation and the degree distributions of the two layersread: P [1] ( k ) ∼ k − , P [2] ( q ) ∼ e − q (11)while the inter-layer joint degree distribution is equal to P ( k, q ) = a ( k ) k − m X n =0 (cid:16) k − mn (cid:17) (cid:18) m m + k − n (cid:19) q − m +1 ( − k − m + n (12) where a ( k ) = Γ( k )Γ( m +1)Γ( k − m +1) . The function ¯ k ( q ) isgiven by ¯ k ( q ) = m (cid:18) m + 1)1 + 2 m (cid:19) q − m +1 (13)Similar relations can be derived for the other couplingpatterns. In Panel a) and b) of Fig. 3 we report thedegree distribution and the value of ¯ k ( q ) for three dif-ferent coupling patterns, which are in good agreementwith the theoretical curves. In the semi-linear model thecoupling pattern has a dramatic impact on other struc-tural properties of the multiplex such as the distributionof the average shortest path length from each node andthe distribution of node interdependence. In particular,the interdependence is smaller for older nodes, and growssublinearly with time. This implies that navigation forold nodes is easier within a single layer while youngernodes will have to resort to the different layers to reach atarget. In addition, a sublinear growth implies that thesystem performance increases very slowly.V.N. and V.L. acknowledge support from the ProjectLASAGNE, Contract No.318132 (STREP), funded bythe European Commission. M.B. thanks Queen MaryUniversity of London for its warm welcome at the startof this project, and is supported by the FET-Proactiveproject PLEXMATH (FP7-ICT-2011-8; grant number317614) funded by the European Commission. [1] R. Albert and A.-L. Barabasi, Rev. Mod. Phys. , 47(2002).[2] M. E. J. Newman, SIAM Review , 167–256 (2003).[3] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez and D.-U. Hwang, Phys. Rep. , 175–308 (2006). [4] A. Barrat, M. Barthelemy and A. Vespignani,
DynamicalProcesses on Complex Networks (Cambridge UniversityPress, Cambridge, England, 2008).[5] M. Kurant and P. Thiran,
Phys. Rev. Lett. , 138701(2006).[6] S.-R. Zou, T. Zhou, A.-F. Liu, X.-L. Xu and D.-R. He, Phys. Lett. A , 4406–4410 (2010).[7] J. Donges, H. Schultz, N. Marwan, Y. Zou and J. Kurths,
Eur. Phys. J. B , 635–651 (2011).[8] J. Yang, W. Wang and G. Chen, Physica , 2435–2449 (2009).[9] S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley andS. Havlin,
Nature (London) , 1025–1028 (2010).[10] E. Bullmore and O. Sporns,
Nat. Rev. Neurosci. , 186–198 (2009).[11] S. Wasserman and K. Faust Social Network Analysis (Cambridge University Press, Cambridge, 1994)[12] R. G. Morris and M. Barthelemy,
Phys. Rev. Lett. ,128703 (2012).[13] K.-M. Lee, J. Y. Kim, W. kuk Cho, K.-I. Goh and I.-M.Kim,
New J. Phys. , 033027 (2012).[14] G. Bianconi, Phys. Rev. E , 062806 (2013).[15] Y. Zhuo, Y. Peng, C. Liu, Y. Liu and K. Long, PhysicaA , 2401–2407 (2011).[16] S. G´omez, A. D´ıaz-Guilera, J. G´omez-Garde˜nes, C. J.P´erez-Vicente, Y. Moreno and A. Arenas,
Phys. Rev.Lett. , 028701 (2013).[17] J. G´omez-Garde˜nes, I. Reinares, A. Arenas and L. M.Flor´ıa,
Sci. Rep. , 620 (2012).[18] A. Halu, K. Zhao, A. Baronchelli and G. Bianconi, EPL–Europhys. Lett. , 16002 (2013).[19] A. Saumell-Mendiola, M. ´A. Serrano and M. Bogu˜n´a,
Phys. Rev. E , 026106 (2012).[20] M. Szell, R. Lambiotte and S. Thurner, Proc. Natl. Acad.Sci. USA , 13636 (2010).[21] M. Magnani, B. Micenkova and L. Rossi, preprint,arXiv:1303.4986 (2013).[22] M. Barigozzi, G. Fagiolo and D. Garlaschelli,
Phys. Rev.E , 046104 (2010).[23] V. Nicosia, M. Valencia, M. Chavez, A. D´ıaz-Guilera, V.Latora, Phys. Rev. Lett. , 174102 (2013).[24] R. Parshani, C. Rozenblat, D. Ietri, C. Ducruet and S.Havlin,
EPL–Europhys. Lett. , 68002 (2010).[25] Notice that, as shown in [26], the Pearson’s coefficientbecomes infinitesimal in N when degree fluctuations areunbounded (e.g., in the case of power-law distributions P ( k ) ∼ k − γ with γ < Phys. Rev. E ,022801 (2013).[27] R. Albert, H. Jeong and A.-L. Barabasi, Nature ,130–131 (1999).[28] It is easy to show that the time required to simulate thegrowth of a multiplex with N nodes and M layers in-creases as O ( N M ). See Appendix for additional details[29] Notice that for t ≫ [ α ] k,q ≪ ∀ k, q ; conse-quently, aymptotically in time we can neglect the prob-ability that the same node acquires, at time t , a link inlayer 1 and in layer 2.[30] S. N. Dorogovtsev, J. F. F. Mendes and A. N. Samukhin, Phys. Rev. Lett. , 4633 (2000).[31] P. L. Krapivsky and S. Redner, Phys. Rev. E , 066123 (2001).[32] All the results shown in the figures correspond to 2-layermultiplexes with N = 10000 nodes, m = 3, m = 3, andare averaged over 50 realizations. Mean-field theory
In this section, using a mean-field approach, we discuss the time evolution of the degree of the nodes of a multiplexin the different layers and we derive long-time expressions for the degree distribution at each layer and for the inter-layer degree-degree correlations. We first consider the linear attachment case (A) and then proceed to the semi-linearcase (B). We also provide a concise discussion of the case in which new nodes bring a random number of new edges(C).
Linear attachment kernel on both layers
According to the growth model discussed in the main text, the probability that a newly arrived node i in themultiplex creates a link to node j on layer α can be written asΠ [ α ] i → j = F [ α ] j ( k j ) P l F [ α ] l ( k l ) (S-1)where F [ α ] j ( k j ) is a certain function of the degrees of the replicas of node j . If F [ α ] j ( k j ) is a linear function of k j ∀ α ,and all the replicas of the new node arrive at the same time, the temporal evolution of the degree of a node on eachlayer is governed by the equations " dk [1] dtdk [2] dt = 12 t C (cid:20) k [1] k [2] (cid:21) = 12 t (cid:20) c [1 , c [1 , c [2 , c [2 , (cid:21)(cid:20) k [1] k [2] (cid:21) (S-2)with the constraints c [1 , + c [1 , = 1 and c [2 , + c [2 , = 1. Since the matrix elements are real and non-zeros, themaximal eigenvalue is real. Moreover since we have just two layers then both eigenvalues λ and λ are real. Wenotice that if we impose that each row of matrix C must sum to 1, and that the coefficients c [ r,s ] are non-negative (toensure that Π [ α ] i → j is a probability distribution ∀ α ), then we can write C = (cid:20) a − ab − b (cid:21) , ≤ a ≤ , ≤ b ≤ b − a + 1 = 0 the matrix C has eigenvalues λ = 1 and λ = ( a − b ) = 1, with eigenvectors u = [1 ,
1] and u = [1 , − b/ (1 − a )] (the degenerate case b − a + 1 = 0 is considered below). Since the eigenvalues aredistinct, then C is diagonalizable, i.e. it is similar to the diagonal matrixΛ = (cid:20) a − b (cid:21) , ≤ a ≤ , ≤ b ≤ C . The system in Eq. (S-2) can be also written in the form˙ k ( t ) = A ( t ) k ( t ) = α ( t ) C k ( t ) (S-3)where α ( t ) = t is a scalar function, and C is a constant matrix. This is a homogeneous time-varying linear dynamicalsystem, whose temporal evolution is fully determined by the initial state k s = m u and by the state transition matrixΦ( t, s ), where s is the time at which a node is added to the graph k ( t ) = Φ( t, s ) k s = m Φ( t, s ) u . (S-4)Since C is diagonalizable, then the transition matrix can be written asΦ( t, s ) = e R ts A ( τ ) d τ = e C R ts d τ τ = e Cσ . (S-5)To compute the transition matrix Φ( t, s ) we first compute the exponential matrix e Cσ and then we substitute σ with R ts d τ τ . Since the eigenvalues of C are distinct, the corresponding eigenvectors form a base of R , so we have e Cσ = V e Λ V − (S-6)where V is the matrix whose columns are the eigenvectors of C and Λ is the diagonal matrix of the eigenvalues of C .After some simple algebra we obtain e Cσ = 11 − ( a − b ) " be σ + (1 − a ) e ( a − b ) σ (1 − a ) e σ − (1 − a ) e ( a − b ) σ be σ + b (1 − a ) a − b e ( a − b ) σ (1 − a ) e σ − b (1 − a )( a − b ) e ( a − b ) σ . (S-7)By substituting σ with R ts d τ τ = (log t − log s ), we get:Φ( t, s ) = 11 − ( a − b ) " b ( ts ) + (1 − a )( ts ) a − b (1 − a )( ts ) − (1 − a )( ts ) a − b b ( ts ) + b (1 − a ) a − b ( ts ) a − b (1 − a )( ts ) − b (1 − a ) a − b ( ts ) a − b (S-8)and the solution of the system reads k s ( t ) = m Φ( t, s ) u = m (cid:18) ts (cid:19) u . (S-9)In this case we have h k [1] | k [2] i ∝ k [2] , (S-10) h k [2] | k [1] i ∝ k [1] . (S-11)Let us now consider the degenerate case b − a + 1 = 0. Since we imposed that each row of the matrix C has to sumto 1, then b − a + 1 = 0 only if a = 1 and b = 0. In this case the two layers evolve independently, the matrix C isdiagonal and the time evolution of the degree on each layer reads k [ α ] s ( t ) = m (cid:18) ts (cid:19) , (S-12)with α = 1 ,
2. This means that in the case of linear attachment kernel on both layers without delay the degreedistribution of each layer is a power-law P ( k ) ∼ k − γ with exponent γ = 3 and we have h k [1] | k [2] i ∝ k [2] , (S-13) h k [2] | k [1] i ∝ k [1] . (S-14) Semi-linear attachment kernel
For the semi-linear attachment kernel we have " dk [1] dtdk [2] dt = 1 t (cid:20) am am +1 − a bm bm +1 − b (cid:21) (cid:20) k [1] k [2] (cid:21) + 1 t " m (1 − a )2 am +1 − a m (1 − b )2 bm +1 − b (cid:21) (S-15)which is in the form ˙ k ( t ) = A ( t ) k ( t ) + B ( t ) w ( t ) (S-16)where A ( t ) = " amt (2 am +1 − a ) bmt (2 am +1 − a ) , B ( t ) = " m (1 − a ) t (2 am +1 − a ) m (1 − b ) t (2 bm +1 − b ) , w ( t ) = (cid:20) (cid:21) . The system in Eq. (S-16) is a non-homogeneous time-varying linear dynamical system where w ( t ) represents anexternal forcing function. If we call Φ( t, s ) the state transition matrix of the corresponding homogeneous system˙ k ( t ) = A ( t ) k ( t ), it is possible to show that the unique solution of Eq. (S-16) is given by k s ( t ) = Φ( t, s ) k s + Z ts d σ Φ( t, σ ) B ( σ ) w ( σ ) . (S-17)The form of the transition matrix Φ( t, s ) associated to the homogeneous system depends on the value of a . When a = 0 then Φ( t, s ) readsΦ( t, s ) = " (cid:0) ts (cid:1) β (2 abm +(1 − a ) b )(2 abm − ab + a ) (cid:0) ts (cid:1) β + (( a − b − abm )(2 abm − ab + a ) , where β = am am − a + 1 . (S-18)By plugging Eq. (S-18) into Eq. (S-17) one obtains the mean-field temporal evolution of k [1] s and k [2] s : k [1] s ( t ) = am − a + 1 a (cid:18) ts (cid:19) β + a − a ,k [2] s ( t ) = δ (cid:18) ts (cid:19) β + η (log( t ) − log( s )) + ǫ, (S-19)where δ = 2 a bm + (3 a − a ) bm + ( a + 1) b a bm − a b + a ,η = ( a − ab ) m a bm − a b + a ,ǫ = ((2 a − a ) b + a ) m − ( a − b a bm − a b + a . In general, if b = 0 then for t → ∞ we have (cid:0) ts (cid:1) β ≫ (log( t ) − log( s )). Consequently, Eqs. (S-19) can be written as k [1] s ( t ) ≃ (cid:18) m + 1 − aa (cid:19) (cid:18) ts (cid:19) β , (S-20) k [2] s ( t ) ≃ δ (cid:18) ts (cid:19) β (S-21)so that the degree distribution on both layers reads P ( k [ ℓ ] ) ∼ k − ( β +1) = k − (3+ − aam ) (S-22)and we have h k [1] | k [2] i ∝ k [2] , h k [2] | k [1] i ∝ k [1] . (S-23)Instead, if b = 0 the solution for the degree of nodes on the second layer reads k [2] s ( t ) = η (log( t ) − log( s )) + ǫ, (S-24)while k [1] s ( t ) is expressed by Eq. (S-20). In this case, the degree distribution on the first layer is the same as inEq. (S-22), while for the second layer we have P ( k [2] ) ∼ e − kη (S-25)and in the limit of large k [1] ( t ) , k [2] ( t ) we obtain h k [1] | k [2] i ∝ e βk [2] η , h k [2] | k [1] i ∝ log( k [1] ) . (S-26)Eqs. (S-18—S-24) are valid when a = 0. When a = 0 the state transition matrix readsΦ( t, s ) = (cid:20) bm (log t − log s )2 bm − b +1 (cid:21) (S-27)and the generic solutions for k [1] s ( t ) and k [2] s ( t ) are k [1] s ( t ) = m (log( t ) − log( s )) + m,k [2] s ( t ) = bm (log( t ) − log( s )) + (cid:0) (2 − b ) m + 2 bm (cid:1) (log( t ) − log( s )) + 4 bm + (2 − b ) m bm − b + 2 . (S-28)In this case, the degree distribution on the first layer is exponential P ( k [1] ) ∼ e − km . On the second layer, the functionalform of the degree distribution depends on the value of b . It is easy to verify that when b = 0 then P ( k [2] ) ∼ e − km ,and we have in the limit of large k [1] ( t ) , k [2] ( t ) h k [1] | k [2] i ∝ k [2] , h k [2] | k [1] i ∝ k [1] . (S-29)Conversely, when b > P ( k [2] ) ∼ e −√ µk + ν √ µk + ν (S-30)where µ = m (cid:2) b m + (cid:0) b − b (cid:1)(cid:3) ,ν = m (cid:2) b m + (cid:0) b − b (cid:1) m + (3 b − b + 1) (cid:3) . In this case, in the limit of large k [1] ( t ) , k [2] ( t ), we have h k [1] | k [2] i ∝ p k [2] , h k [2] | k [1] i ∝ (cid:16) k [1] (cid:17) . (S-31) Fluctuations in the number of edges
In principle, the mean-field approach could be also applied to the case in which the number of edges brought onlayer α by each new-born node is not fixed but is a random variable ξ [ α ] drawn from a given distribution P ( ξ [ α ] ). Inthis case we should solve the system of stochastic differential equations: " dk [1] dtdk [2] dt = 12 (cid:0) κ [1] ( t ) + κ [2] ( t ) (cid:1) (cid:20) ξ [1] ( t ) 00 ξ [2] ( t ) (cid:21) (cid:20) a − ab − b (cid:21) (cid:20) k [1] k [2] (cid:21) (S-32)where: κ [ α ] ( t ) = Z t d τ ξ [ α ] ( τ )The random variable ξ [ α ] is a positive integer with average h ξ [ α ] i and the dominant term at large times of κ [ α ] is thengiven by κ [ α ] ( t ) ≃ t h ξ [ α ] i + o ( t )This implies in particular that at large times, the effect of randomness in the number of edges is negligible, and thebehavior of the system is governed by the average number of edges h ξ [ α ] i added in each layer. Master Equation approach for the model without delay
We provide here the derivation of exact expressions of P ( k ) and P ( k, q ) starting from the master equation of thesystem. We denote by N k,q ( t ) the average number of nodes that at time t have degree k in layer 1 and degree q in0layer 2. We start from a small connected network and at each time we add a node which brings, at same time, m new edges in layer 1 and m new edges in layer 2. We assume that, when we add the new node i to the network, theexpected number of new links in layer 1 attached to a node j of degree k in layer 1 and degree q in layer 2 is givenby m Π [1] i → j = A k,q t . Similarly, the expected number of new links in layer 2 attached to a node j of degree k in layer 1and degree q in layer 2 is given by m Π [2] i → j = B k,q t . In addition to that, we work in the hypothesis that in the large t limit, t ≫
1, we have A k,q /t ≪ B k,q /t ≪ N k,q ( t + 1) = N k,q ( t ) + A k − ,q t N k − ,q ( t ) + B k,q − t N k,q − ( t ) − (cid:20) A k,q t + B k,q t (cid:21) N k,q ( t ) + δ k,m δ q,m (S-33)for k ≥ m and q ≥ m , as long as N m − ,q ( t ) = N k,m − ( t ) = 0. Assuming that N k,q = tP ( k, q ) is valid in the largetime limit t ≫
1, we can solve for the combined degree distribution P ( k, q ) indicating the probability that a node hasat the same time degree k in layer 1 and degree q in layer 2. We get the master equations P ( m, q ) = q Y j = m +1 B m,j − A m,j + B m,j P ( m, m ) ,P ( k, q ) = q X r = m q Y j = r +1 B k,j − A k,j + B k,j A k − ,r A k,r + B k,r P ( k − , r ) (S-34) Solution of the master equation in three simple cases (i) Linear attachment kernel –
Let us first consider a linear preferential attachment kernel, in which c [1 , = c [2 , = 1and c [1 , = c [2 , = 0. In this case we have A k,q = k ,B k,q = q . (S-35)The recursive Eqs. (S-34) read P ( m, q ) = Γ( q )Γ(3 + 2 m )Γ( m )Γ(3 + q + m ) P ( m, m ) ,P ( k, q ) = q X r = m (cid:18) Γ( q )Γ(3 + r + k )Γ( r )Γ(3 + q + k ) (cid:19) k −
12 + k + q P ( k − , r ) (S-36)where P ( m, m ) is fixed by the normalization condition P ∞ k = m P ∞ q = m P ( k, q ) = 1. Using the relation q X r = m Γ( k + r − m )Γ( r − m + 1)Γ( k − m ) = Γ( k + q − m + 1)Γ( k − m + 1)Γ( q − m + 1) (S-37)it can be proved recursively that P ( k, q ) takes the following expression P ( k, q ) = 2Γ(2 + 2 m )Γ( m )Γ( m ) Γ( k + q − m + 1)Γ( k + q + 3) Γ( q )Γ( q − m + 1) Γ( k )Γ( k − m + 1) (S-38)Summing over the degree in layer 2 we can find the degree distribution P ( k ) in layer 1, i.e. P ( k ) = P ∞ q = m P ( k, q )obtaining the known result for a single layer, P ( k ) = 2 m (1 + m ) k ( k + 1)( k + 2) (S-39)The function h k ( q ) i is given by h k ( q ) i = P ∞ k = m kP ( k, q ) P ∞ k = m P ( k, q ) = m m ( q + 2) (S-40)1Similar expressions are obtained for P ( q ) and h q ( k ) i , by summing Eq. (S-38) over k . (ii) Uniform attachment kernel – Let us now consider a uniform attachment kernel, in which every target node j ischosen with probability Π [1] i → j = t in layer 1 and with probability Π [2] i → j = t in layer 2, so that A k,q = mB k,q = m. (S-41)In this case the recursive Eqs. (S-34) read P ( m, q ) = (cid:18) m m (cid:19) q − m P ( m, m ) ,P ( k, q ) = q X r = m (cid:18) m m (cid:19) q − r m m P ( k − , r ) (S-42)where P ( m, m ) is again fixed by the normalization condition P ∞ k = m P ∞ q = m P ( k, q ) = 1. Using again the relationprovided in Eq. (S-37) it is easy to prove recursively that P ( k, q ) is given in this case by P ( k, q ) = 1 m (cid:18) m m (cid:19) k + q − m +1 Γ( k + q − m + 1)Γ( k − m + 1)Γ( q − m + 1) (S-43)Moreover the degree distribution P ( k ) = P ∞ q = m P ( k, q ) and P ( q ) = P ∞ k = m P ( k, q ) of a single network are given by P ( k ) = 11 + m (cid:18) m m (cid:19) k − m and P ( q ) = 11 + m (cid:18) m m (cid:19) q − m (S-44)while the function h k ( q ) i = P ∞ k = m kP ( k, q ) /P ( q ) is given by h k ( q ) i = ( q + 2) m m (S-45) (iii) Semi-linear attachment kernel – Finally we analyze the case of semi-linear attachment with c [1 , = c [2 , = 1and c [1 , = c [2 , = 0. We have: A k,q = k B k,q = m. (S-46)In this case the Eqs. (S-34) read as P ( m, q ) = (cid:18) m m (cid:19) q − m P ( m, m ) ,P ( k, q ) = q X r = m (cid:18) m k + 2 m (cid:19) q − r k −
12 + k + 2 m P ( k − , r ) (S-47)where P ( m, m ) is fixed by the normalization condition P ∞ q = m P ∞ k = m P ( k, q ) = 1. It can be shown recursively thatthese equations have the following solution, P ( k, q ) = 1Γ( m + 1) Γ( k )Γ( k − m + 1) k − m X n =0 (cid:18) k − mn (cid:19) (cid:18) m m + k − n (cid:19) q − m +1 ( − k − m − n , with the associated degree distributions P ( k ) = P ∞ q = m P ( k, q ) and P ( q ) = P ∞ k = m P ( k, q ) given by P ( k ) = 2 m (1 + m ) k ( k + 1)( k + 2) ,P ( q ) = 1Γ(1 + m ) ∞ X k = m k − m X n =0 (cid:18) k − mn (cid:19) (cid:18) m m + k − n (cid:19) q − m +1 ( − k − m + n = (cid:18) m m (cid:19) q − m
11 + m h k ( q ) i = P ∞ k = m kP ( k, q ) /P ( q ) is given by h k ( q ) i = (cid:18) m + 1)1 + 2 m (cid:19) q − m +1 m (S-48) Role of β in the delayed arrival In Fig. S-1 we show the time evolution of the maximum degree k M ( t ) on the first layer, for different values of β .Notice that k M ( t ) ∼ ( t/s ) δ . The effect of the exponent β tuning the width of the delay distribution is evident: thelarger the value of β , the closer δ is to 0 .
5, the value observed in the case of synchronous arrival. Consequently, therightmost part of the degree distribution is broader when β is close to 1 and becomes more similar to P ( k ) ∼ k − when β increases. Finite size effects
It is interesting to investigate how the properties of the multiplexes generated using the model we propose dependon the number of nodes N . For instance, most of the mean-field predictions for the degree distributions and inter-layerdegree correlations are valid in the limit of large N . However, as shown in Fig. S-2 and in Fig. S-3, the properties ofthe degree distributions and of inter-layer degree-degree correlations are similar to those predicted for large N evenfor relatively small multiplexes, e.g. with N = 1000. Time complexity
The most efficient algorithm for the construction of a simplex networks based on preferential attachment takesadvantage of random sampling with rejection and runs in O ( N m ). However, in the case of a multiplex the procedureto sample a candidate neighbour j of a newly-arrived node i is a bit more complicated. Let us first consider a two-layermultiplex described by Eq. (S-2). When we sample the candidate neighbours of node i at layer 1 at time t , each node j should be sampled with a probability proportional to ak [1] j + (1 − a ) k [2] j . The simplest way to implement suchsampling is to construct a vector S whose n -th entry is equal to P nj =1 ak [1] j + (1 − a ) k [2] j (the first element S [0] of thearray is set equal to zero); then we sample a real number ζ in the interval (0 , S [ t ]] and we choose the node j suchthat ζ ≤ S [ j ] and ζ > S [ j − S at each time t requires O ( t ) operations while thesampling of a single node j can be efficiently implemented by binary search, requiring at most O (log( t )) operationsper edge, so that the sampling of m edges requires at most O ( m log( t )) steps. Thus, the total number of operationsneeded to sample a layer of a multiplex is: N X j = m +1 O ( t ) + O ( m log( t )) = O ( N ) + O ( mN log( N )) = O ( N + mN log( N )) (S-49)It is easy to verify that the construction of a M -layer multiplex requires a number of steps O ( M ( N + mN log( N )) (S-50)which, for m fixed, is dominated by O ( M N ). Therefore, the time complexity of this algorithm is linear in thenumber of layers and quadratic in the number of nodes. We notice that in principle it is possible to construct betteralgorithms to sample growing multiplexes by implementing a smart policy to update the array S . Randomly-chosen master layer
In the main text we made the simplifying assumption that each node arrives first on the master layer and thenon the other layers, after a certain delay. We call this assumption “Equal master layer” (EML). In this Section webriefly comment on the case in which this assumption does not hold, i.e. when a node first arrives either on the first3 t k M (t) β =1.1 - δ =0.56 β =1.5 - δ =0.53 β =2.0 - δ =0.51 β =2.5 - δ =0.50 FIG. S-1.
Effect of β on the maximum degree of a layer. When the exponent β of the delay distribution is close to 1,the maximum degree of a layer scales as ( t/s ) δ , with δ > /
2. Consequently, the rightmost side of the degree distribution isbroader when β is close to 1 and converges to P ( k ) ∼ k − as β increases. k -4 N(k)
N=1000N=5000N=10000N=50000N=100000~ k -3 FIG. S-2.
Degree distribution on the first layer for different values of N . Even for relatively small multiplexes, i.e. N = 1000, the exponent of the degree distribution is almost equal to γ = 3. or on the second layer, and then arrives on the other layer after a power-law distributed delay. We call this case“Randomly-chosen master layer” (RML), to stress the fact that the master layer of each node is chosen at randomamong the M layers of the multiplex. In particular, we are interested in the case in which a newly arrived node4 k k(k) N=1000N=5000N=10000N=50000N=100000k(k) ~ k
FIG. S-3.
Inter-layer degree correlation as a function of N . The shape of the inter-layer degree-degree correlations isalready evident even for small values of N . The curves have been vertically displaced to facilitate visual comparison. selects one of the M layers of the multiples as its master layer with uniform probability p = 1 /M . In Fig. S-4, S-5and S-6 we report, respectively, the degree distributions, the temporal scaling of the degree of the largest hub andthe distribution of shortest path lengths and node interdependence for RML with M = 2. The plots suggest that therandom choice of the master layer produces a more balanced distribution of super hubs between the two layers, whichhas a relevant impact on the distribution of shortest path lengths and node interdependence (Fig. S-6). Conversely,the degree distributions and the temporal scaling of the degree of the largest hub are practically indistinguishablefrom those observed in EML (Fig. S-4 and S-5).5 -6 -3 P(k) ~ k -3 -6 -3 -6 -3 P(k) ~ k -3 k -6 -3 a)b) FIG. S-4.
Degree distribution of mixing master layers . The degree distribution of the first layer for the EML (panel a)and the RML (panel b) for β = 1 . β = 2 .
100 1000 10000 t k M (t) β=1.1 − δ=0.57β=1.5 − δ=0.53β=2.0 − δ=0.52β=2.5 − δ=0.50 FIG. S-5.
Temporal scaling of degree for RML . The degree k M of the largest hub in the first layer scales as ( t/s ) δ , where δ depends on the exponent β of the delay distribution and approaches 0 . β increases. Notice that the values of δ are prettysimilar to those observed in the case of EML, reportes in Fig. S-1. L i P(L i ) λ i P( λ i ) L i P(L i ) λ i P( λ i ) a)b) FIG. S-6.