Controllability of multiplex, multi-timescale networks
Márton Pósfai, Jianxi Gao, Sean P. Cornelius, Albert-László Barabási, Raissa M. D'Souza
CControllability of multiplex, multi-timescale networks
M´arton P´osfai ∗ Complexity Science Center and Department of Computer Science,University of California, Davis, CA 95616, USA andDepartment of Physics of Complex Systems,E¨otv¨os University, Budapest, H-1117, Hungary
Jianxi Gao and Sean P. Cornelius
Center for Complex Network Research, Department of Physics,Northeastern University, Boston, MA 02115, USA
Albert-L´aszl´o Barab´asi
Center for Complex Network Research, Department of Physics,Northeastern University, Boston, MA 02115, USACenter for Cancer Systems Biology, Dana-Farber Cancer Institute,Harvard University, Boston, MA 02215, USADepartment of Medicine, Brigham and Women’s Hospital,Harvard Medical School, Boston, MA 02115, USA andCenter for Network Science, Central European University, Budapest, H-1051, Hungary
Raissa M. D’Souza
Complexity Science Center, Department of Computer Scienceand Department of Mechanical and Aerospace Engineering,University of California, Davis, CA 95616, USA andSanta Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA (Dated: September 30, 2016) a r X i v : . [ phy s i c s . s o c - ph ] A ug bstract The paradigm of layered networks is used to describe many real-world systems – from biologicalnetworks, to social organizations and transportation systems. While recently there has been muchprogress in understanding the general properties of multilayer networks, our understanding of how tocontrol such systems remains limited. One fundamental aspect that makes this endeavor challengingis that each layer can operate at a different timescale, thus we cannot directly apply standard ideasfrom structural control theory of individual networks. Here we address the problem of controllingmultilayer and multi-timescale networks focusing on two-layer multiplex networks with one-to-oneinterlayer coupling. We investigate the practically relevant case when the control signal is appliedto the nodes of one layer. We develop a theory based on disjoint path covers to determine theminimum number of inputs ( N i ) necessary for full control. We show that if both layers operate onthe same timescale then the network structure of both layers equally affect controllability. In thepresence of timescale separation, controllability is enhanced if the controller interacts with the fasterlayer: N i decreases as the timescale difference increases up to a critical timescale difference, abovewhich N i remains constant and is completely determined by the faster layer. We show that thecritical timescale difference is large if Layer I is easy and Layer II is hard to control in isolation. Incontrast, control becomes increasingly difficult if the controller interacts with the layer operating onthe slower timescale and increasing timescale separation leads to increased N i , again up to a criticalvalue, above which N i still depends on the structure of both layers. This critical value is largelydetermined by the longest path in the faster layer that does not involve cycles. By identifying theunderlying mechanisms that connect timescale difference and controllability for a simplified model,we provide crucial insight into disentangling how our ability to control real interacting complexsystems is affected by a variety of sources of complexity. ∗ [email protected] . INTRODUCTION Over the past two decades, the theory of networks proved to be a powerful tool forunderstanding individual complex systems [1, 2]. However, it is now increasingly appreciatedthat complex systems do not exist in isolation, but interact with each other [3, 4]. Indeed,an array of phenomena – from cascading failures [5, 6] to diffusion [7] – can only be fullyunderstood if these interactions are taken into account. Traditional network theory is notsufficient to describe the structure of such systems, so in response to this challenge, theparadigm of multilayer networks is being actively developed. Here we study a fundamental,yet overlooked aspect of multilayer networks: each individual layer can operate at a differenttimescale. Particularly, we address the problem of controlling multilayer, multi-timescalesystems focusing on two-layer multiplex networks. Recently significant efforts have been madeto uncover how the underlying network structure of a system affects our ability to influenceits behavior [8–16]. However, despite the appearance of coupled systems from infrastructureto biology, the existing literature – with a few notable exceptions [17–20] – has focused oncontrol of networks in isolation, and the role of timescales remains unexplored.Control of multilayer networks is important for many applications. For example, considera CEO aiming to lead a company consisting of employees and management. Studying thenetwork of managers or the network of employees in isolation does not take into account im-portant interactions between the different levels of hierarchy of the company. On the otherhand, treating the system as one large network ignores important differences between thedynamics of the different levels, e.g. management may meet weekly, while employees are indaily interaction. In general, the interaction of timescales plays an important role in orga-nization theory [21]. Or consider gene regulation in a living cell. External stimuli activatesignaling pathways which through a web of protein-protein interactions affect transcriptionfactors responsible for gene expression. The activation of a signaling pathway happens on thetimescale of seconds, while gene expression typically takes hours [22]. As a third example,consider an operator of an online social network who wants to enhance the spread of certaininformation by interacting with its users. However, a user may subscribe to multiple socialnetworking services and may opt to share news encountered in one network through a differ-3nt one – out of reach of the operator. The dynamics of user interaction on different websitescan be very different depending on user habits and the services offered [23–25]. For examplethe URL shortening service Bit.ly reports that the half-life of shared links depends on thesocial networking platform used: half the clicks on a link happened within 2.8 hours afterposting on Twitter, within 3.2 hours on Facebook and within 7.4 hours on Youtube [26].Common features of these examples are that (i) each interacting subsystem is describedby a separate complex network; (ii) the dynamics of each subsystem operate on a different,but often comparable timescale and (iii) the external controller directly interacts with onlyone of the subsystems. Here we study the control properties of a model that incorporatesthese common features, yet remains tractable. More specifically, we study discrete-time lin-ear dynamics on two-layer multiplex networks, meaning that we assume one-to-one couplingbetween the nodes of the two layers. This choice ensures both analytical tractability and theisolation of the role of timescales from the effect of more complex multilayer network struc-ture. Identifying the underlying mechanisms that govern the controllability of this simplemodel provides crucial insight into disentangling how our ability to control real interactingcomplex systems is affected by a variety of sources of complexity.So far only limited work investigated controllability of multilayer networks. Menichettiet al. investigated the controllability of two-layer multiplex networks governed by lineardynamics such that the dynamics of the two layers are not coupled, but the input signals inthe two layers are applied to the same set of nodes [18]. Yuan et al. identified the minimumnumber of inputs necessary for full control of diffusion dynamics, allowing the controller tointeract with any layer [19]. Zhang et al. investigated the controllable subspace of multilayernetworks with linear dynamics without timescale separation if the controller is limited tointeract with only one layer; showing that it is more efficient to directly control peripheralnodes than central ones [20]. Here we also limit the controller to one layer, yet by exploringthe minimum input problem, we offer a direct metric which allows us to compare our findingsto previous results for single-layer networks [10]. More so, the key innovation of our work isthat we take into account the timescale of the dynamics of each layer, a mostly overlookedaspect of multilayer networks.It is worth mentioning the recent work investigating the related, but distinct problem of4ontrollability of networks with time-delayed linear dynamics [27]. The key difference betweentime-delay and timescale difference is that for time-delayed dynamics the state of a node willdepend on some previous state of its neighbors; however, the typical time to change the stateof a node remains the same throughout the system. While in case of timescale difference, thetypical time needed for changes to happen can be different in different parts of the system.In the next section, we introduce a simple model that captures some common propertiesof multilayer networks and we describe the problem setup. In Sec. III, we develop a theory todetermine the minimum number of inputs required for controlling multiplex, multi-timescalenetworks with discrete-time linear dynamics relying on graph combinatorial methods. InSec. IV, we use networks with tunable degree distribution to systematically uncover therole of network structure and timescale separation. We study three scenarions: no timescaleseparation, Layer I operates faster and Layer II operates faster. Finally, in Sec. V we providea discussion of our results and we outline open questions.
II. MODEL DEFINITION
We aim to study the controllability of coupled complex dynamical systems with the follow-ing properties: (i) each subsystem (layer) is described by a complex network; (ii) the operationof each layer is characterized by a different timescale and (iii) the controller only interactsdirectly with one of the layers. We propose a model that satisfies these requirements andyet is simple enough to remain tractable. We focus on two-layer multiplex systems, meaningthat there is a one-to-one correspondence between the nodes of the two layers.The model is defined by a weighted directed two-layer multiplex network M which consistsof two networks L I and L II called layers and a set of links E I , II connecting the nodes ofthe different layers. Each layer L α (where α ∈ { I , II } ) consists of a set of nodes V α = { v α , v α , . . . , v αN } and a set of links E α , where a directed link ( v αi , v αj , w αij ) ∈ E α is an orderednode pair and a weight representing that node v αi influences node v αj with strength w αij . Thetwo layers are connected by link set E I , II = { ( v I i , v II i , w I,II i ) | i = 1 , . . . , N } , in other words,there is directed one-to-one coupling from Layer I to Layer II (Fig. 1a). Although the linksare weighted, the exact values of the weights do not have to be known for our purposes.5ur goal is to control the system by only interacting directly with Layer I. We studylinear discrete-time dynamics x I ( t ) = A I x I ( t − τ I ) + Bu ( t − τ I ) x II ( t ) = A II x II ( t − τ II ) + ∆ τ I ( t ) Dx I ( t − τ I ) if ( t mod τ I ) = 0 , if ( t mod τ II ) = 0 , (1)where x I ( t ) and x II ( t ) ∈ R N represent the state of nodes in Layer I and II; the matrices A I and A II ∈ R N × N are the transposed weighted adjacency matrices of Layer I and II,capturing their internal dynamics. The weighted diagonal matrix D ∈ R N × N captures howLayer I affects Layer II.Vector u ( t ) ∈ R M provides the set of independent inputs and the matrix B ∈ R N × M defines how the inputs are coupled to the system. To differentiate between the function u ( t )and an instance of the function at a given time step, we refer to a component u i ( t ) of vector u ( t ) as an independent input, and we call its value at time step t (cid:48) , u i ( t = t (cid:48) ), a signal.Finally, τ I , τ II ∈ { , , . . . } are the timescale parameters of each subsystem, meaning thatthe state of Layer I is updated according to Eq. (1) every τ I th time step; and Layer II isupdated every τ II th time step. And∆ τ I ( k ) = k mod τ I ) = 0 , k mod τ I ) (cid:54) = 0 , (2)is the Kronecker comb, meaning that Layer I directly impacts the dynamics of Layer II if thetwo layers simultaneously update. We investigate three scenarios: (i) the subsystems operateon the same timescale, i.e. τ I = τ II = 1; (ii) Layer I updates faster, i.e. τ I = 1 and τ II > τ I > τ II = 1.We seek full control of the system as defined by Kalman [28], meaning that with theproper choice of u ( t ), we can steer the system from any initial state to any final state infinite time. To characterize controllability, we aim to design a matrix B such that the systemis controllable and the number of independent control inputs, M , is minimized. The minimumnumber of inputs, N i , serves as our measure of how difficult it is to control the system.To find a robust and efficient algorithm to determine N i , we rely on the framework ofstructural controllability [29]. We say that a matrix A ∗ has the same structure as A , if6he zero/nonzero elements of A and A ∗ are in the same position, and only the value ofthe nonzero entries can be different, in other words, in the corresponding network the linksconnect the same nodes, only the link weights can differ. A linear system of Eq. (1) de-fined by matrices ( A I , A II , D , B ) is structurally controllable if there exists matrices with thesame structure ( A ∗ I , A ∗ II , D ∗ , B ∗ ) such that the dynamics defined by ( A ∗ I , A ∗ II , D ∗ , B ∗ ) arecontrollable according to the definition of Kalman. Note that ultimately we are interestedin controllability and not structural controllability. Yet, structural controllability is a usefultool because (i) if a linear system is structurally controllable, it is controllable for almost alllink weight combinations [29] and (ii) determining structural controllability can be mappedto a graph combinatorial problem allowing for efficient and numerically robust algorithms. III. MINIMUM INPUT PROBLEM
Before addressing the minimum input problem of multiplex networks, we revisit the case ofsingle-layer networks by providing an alternative explanation of the Minimum Input Theoremof Liu et al. [10]. This new approach readily lends itself to be extended to multiplex, multi-timescale networks. Thus providing the basis for Sec. III B, in which we develop an algorithmto determine N i for two-layer multiplex networks. A. Single-layer networks
The linear discrete-time dynamics associated to a single-layer weighted directed network L are formulated as x ( t + 1) = Ax ( t ) + Bu ( t ) , (3)where x ( t ), A , B and u ( t ) are defined similarly as in Eq. (1) (Fig. 2a). To obtain a graphcombinatorial condition for structural controllability we rely on the dynamic graph D T ,which represents the time evolution of a system from t = 0 to t = T [30–32]. Each node v i in L is split into T + 1 copies { v i, , v i, , . . . , v i,T } , each copy v i,t represents the state of node v i at time step t . We add a directed link ( v i,t → v j,t +1 ) for t = 0 , , . . . , T − v i → v j ) in the original network, representing that the state of7ode v j at time t + 1 depends on the state of its in-neighbours at the previous time step. Toaccount for the controller, for each independent input we create T nodes u i,t ( i = 1 , , . . . , M ; t = 0 , , . . . , T −
1) each representing a control signal (i.e. the value of the i th input at timestep t ). We draw a directed link ( u i,t → v j,t +1 ) for t = 0 , , . . . , T − b ji (cid:54) = 0, where b ji isan element of matrix B .According to Theorem 15.1 of Ref. [30], a linear system ( A , B ) is structurally controllableif and only if in the associated dynamic graph D N node sets U = { u i,t | i = 1 , , . . . , M ; t =0 , . . . , N − } (green nodes in Fig. 2b) and V N = { v i,t = N | i = 1 , . . . , N } (blue nodes)are connected by N disjoint paths (red links), i.e. there exists a set of disjoint paths C = { P , P , . . . , P N } such that U contains the set of starting points and V T is the set of endpoints.A path P of length l between node v i and v i l is a sequence of l consecutive links [( v i → v i ) , ( v i → v i ) , . . . , ( v i l − → v i l )] such that each node is traversed only once. Node v i is thestarting point and v i l is the endpoint of P . Two paths P and P are disjoint if no node istraversed by both P and P , a set of paths is disjoint if all paths in the set are pairwisedisjoint.A possible interpretation of this result is that if a P i path has starting point u j,t andendpoint v k,t , we say that the signal u j ( t ) is assigned to set x k ( t ), the state of node v k attime t , through path P i . Therefore we refer to path P i as a control path. The clear meaningof the dynamic graph and the control paths makes this condition useful to formulate proofsand to interpret results. However, it is rarely implemented to test controllability of largenetworks, because the size of the dynamical graph grows as N , rendering such algorithmstoo slow. In the following, we provide a condition that only requires the dynamic graph D as input; therefore it is more suitable for practical purposes.It was shown in Refs. [10, 30, 33] that a linear system ( A , B ) is structurally controllableif and only if (i) in D we can connect nodes U ∪ V = { u i,t =0 | i = 1 , , . . . , M } ∪ { v i,t =0 | i =1 , . . . , N } (green nodes in Fig. 2c) and nodes V = { v i,t =1 | i = 1 , . . . , N } (blue nodes) via N disjoint paths (red links) and (ii) all nodes are accessible from the inputs. This result canbe understood as a self-consistent version of the previous condition involving D N : Insteadof keeping track of the entire control paths as we previously did, we concentrate on a singletime step. Consider the dynamic graph D representing the time evolution of the system from8 = 0 to t = 1, and assume that the system is controllable. By definition we can set the stateof each node independently at t = 0; therefore we can treat them as control signals to controlthe system at a later time step. Now let us aim to control the system at t = 1, according toour previous condition, it is necessary that N disjoint paths exist between nodes U ∪ V = { u i,t =0 | i = 1 , , . . . , M } ∪ { v i,t =0 | i = 1 , . . . , N } and nodes V = { v i,t =1 | i = 1 , . . . , N } . Thisis exactly requirement (i), together with the accessibility requirement (ii) it is a sufficient andnecessary condition. Note that D is a bipartite network (each link is connected to exactlyone node in U ∪ V and one node in V ) and each disjoint path in D is a single link.The minimum input problem aims to identify the minimum number of inputs that guar-antee controllability for a given network, in other words, the goal is to design a B ∈ R N × M for a given A such that M is minimized. For this we consider the dynamic graph D withoutnodes representing control signals. We then find a maximum cardinality matching, where amatching is a set of links that do not share an endpoint. The matching is a set of disjointpaths connecting node sets V and V . Controllability requires N disjoint paths between U ∪ V and V ; therefore N i = N − N match , where N match is the size of the maximum matching(if N match = N , N i = 1). Allowing the inputs to be connected to multiple nodes we canguarantee that all nodes are accessible from the inputs. Thus we recovered the MinimumInput Theorem of Liu et al. [10].In summary, by relying on a self-consistent condition for structural controllability we re-derived the known result that identifying N i is equivalent to finding a maximum matchingin D . In the next section we show that this new self-consistent approach lends itself to beextended to the multiplex, multi-timescale model defined by Eq. (1), allowing us to deriveanalogous method to identify N i . B. Multiplex networks
To find the minimum number of inputs N i for multiplex, multi-timescale networks, wefirst extend the definition of the dynamic graph. We define the dynamic graph D τ II such thatit captures the time evolution of a multiplex system defined by ( A I , A II , D , B ) and Eq. (1)from t = 0 to t = τ II . For sake of brevity, we assume that τ I = 1 and τ II ≥
1, the case9f τ I > τ II = 1 is treated similarly (Fig. 1d). Each node v I i in Layer I is split into τ II + 1 copies { v I i, , v I i, , . . . , v I i,τ II } ; each node v II i in Layer II is split into two copies { v II i, , v II i,τ II } ,because Layer II does not update during the intermediate time steps. We draw a link from v I i,t to v I j,t +1 ( t = 0 , , . . . , τ II −
1) if they are connected in Layer I by a directed link ( v I i → v I j ),and similarly we connect v II i, to v II j,τ II if they are connected in Layer II. In addition we drawa link between each pair v I i, and v II i,τ II to account for the interconnectedness.As a natural extension of self-consistent approach introduced in Sec. III A, assume thatthe system is controllable. If the system is controllable, we can set the state of each node in-dependently at t = 0. To control the system at t = τ II , all nodes at t = τ II in D τ II (blue nodesin Fig. 1) have to be connected to a node at t = 0 or to a control signal (green nodes) via adisjoint path (red links). In other words, a linear two-layer system ( A I , A II , D , B ) is struc-turally controllable only if there exists 2 N disjoint paths in the dynamic graph connectingnode set U ∪ V = { u i,t | i = 1 , , . . . , M ; t = 0 , , . . . , τ II − } ∪ { v I i, | i = 1 , , . . . , N } ∪ { v II i, | i =1 , , . . . , N } and node set V τ II = { v I i,τ II | i = 1 , , . . . , N } ∪ { v II i,τ II | i = 1 , , . . . , N } . In otherwords, a linear two-layer system ( A I , A II , D , B ) is structurally controllable only if thereexists 2 N disjoint paths in the dynamic graph connecting node set U ∪ V = { u i,t | i =1 , , . . . , M ; t = 0 , , . . . , τ II − } ∪ { v I i, | i = 1 , , . . . , N } ∪ { v II i, | i = 1 , , . . . , N } and node set V τ II = { v I i,τ II | i = 1 , , . . . , N } ∪ { v II i,τ II | i = 1 , , . . . , N } .To test whether the system is controllable by M independent inputs, we need to find a B ∈ R N × M such that the system is controllable. We do not have to check all possibilities,because if such B exists, then the system is also controllable for B (cid:48) ∈ R N × M where B (cid:48) has nozero elements; therefore, we only check the case when each input is connected to each nodein Layer I. Given matrices ( A I , A II , D , B (cid:48) ), we now have to count the number of disjointpaths connecting U ∪ V and V τ II in the corresponding dynamic graph D II . We find thesepaths using maximum flow: We set the capacity of each link and each node to 1, we thenfind the maximum flow connecting source node set U ∪ V to target node set V τ II using anymaximum flow algorithm of choice. If the system is structurally controllable, the maximumflow equals to 2 N ; if it is less than 2 N , additional inputs are needed.We can now identify the minimum number of inputs N i by systematically scanning possiblevalues of M . A simple approach is to first set M = 1, and test if the system is controllable.10f not, increase M by one. Repeat this until the smallest M yielding full control is found.Significant increase in speed is possible if we find the minimum value of M using bisection.We initially know that N upperi = N ≥ N i ≥ N loweri = 1. We set M = ( N upperi + N loweri ) /
2, andtest if the system is controllable. If yes, we set N upperi = M ; if no, we set N loweri = M . Werepeat this until N upperi = N loweri , which provides N i . For implementation, we used GoogleOR-tools and igraph python packages [34, 35].The one-to-one coupling between Layer I and Layer II guarantees that full control ispossible with at most N independent inputs; therefore we often normalize N i by N , i.e. n i = N i /N .Note that in the above argument we rely on the test of structural controllability based onthe dynamic graph, which was originally introduced for single-timescale networks [30]. Thesufficiency of the condition relies on the fact that the zero is the only degenerate eigenvalueof a matrix A if the nonzero elements of A are uncorrelated. However, this might not remaintrue for the spectrum of A τ , where τ >
1, due to correlations arising in the nonzero elementsof A τ . If a λ (cid:54) = 0 eigenvalue has larger geometric multiplicity than the multiplicity of 0, N i would be larger than predicted by the dynamic graph; if a λ (cid:54) = 0 eigenvalue has largergeometric multiplicity than 1 but smaller than the multiplicity of zero, it does not affect N i , but may require connecting an input to multiple nodes [12]. In the τ I > τ II = 1case, a control signal is only injected into Layer II every τ I time step (Fig. 1d); therefore, thespectrum of A τ I II becomes relevant. However, we are interested in large and sparse complexnetworks whose spectra is dominated by the zero eigenvalue [12]. Therefore it is reasonableto expect that the spectrum of A τ will be dominated by zero eigenvalues as well. Meaningthat the minimum number of inputs is correctly given by this graph combinatorial condition.Furthermore the one-to-one coupling between the layers guarantees that control is possibleby only interacting with Layer I directly.So far, we developed a method to characterize controllability of a multiplex, multi-timescale system based on the underlying network structure and the timescale of each ofits layers. In the next section, we rely on these tools to systematic study how network char-acteristics and timescales affect N i . 11 V. RESULTS
In this section we investigate how different timescales and the degree distribution of eachlayer affect controllability. For timescales, we consider three scenarios: (i) the subsystemsoperate on the same timescale, i.e. τ I = τ II = 1; (ii) Layer I updates faster, i.e. τ I = 1 and τ II >
1; and (iii) Layer II updates faster τ I > τ II = 1. To uncover the effect of degreedistribution, we consider layers with Poisson (ER) or scale-free (SF) degree distribution, thelatter meaning that the distribution has a power-law tail.We generate scale-free layers using the static model [36]: We start with N unconnectednodes. Each node v i is assigned two hidden parameters w in ( i ) = i − ζ out and w out ( i ) = i − ζ out ,where i = 1 , , . . . , N . The weights are then shuffled to eliminate any correlations of the in-and out-degree of individual nodes and between layers. We then randomly place L directedlinks by choosing the start- and endpoint of the link with probability proportional to w in ( i )and w out ( i ), respectively. For large N this yields the degree distribution P SFin/out ( k ) = (cid:2) c (1 − ζ in/out ) /ζ in/out (cid:3) ζ in/out Γ( k − /ζ in/out , c [1 − ζ in/out ])Γ( k + 1) , (4)where c = L/N is equal to the average degree, and Γ( n, x ) is the upper incomplete gammafunction. For large k , P SFin/out ( k ) ∼ k − (1+1 /ζ in/out ) = k − γ in/out , where γ in/out = 1 + 1 /ζ in/out is theexponent characterizing the tail of the distribution.To reduce the number of parameters we only study layers with symmetric degree distri-bution, e.g. P ( k ) = P in ( k ) = P out ( k ); however, the in- and out-degree of a specific node canbe different. A. No timescale separation ( τ I = τ II = 1 ) In the special case when both layers operate on the same timescale, i.e. τ I = τ II = 1(Fig. 1b), there is no qualitative difference between the dynamics of the layers. The reasonwhy the system cannot be treated as a single large network is that we are only allowed todirectly interact with Layer I. Recently Iudice et al. developed methodology to identify N i ifthe control signals can only be connected to a subset of nodes [16]. However, the one-to-one12oupling between the layers enables us to find N i using a simpler approach.Finding N i for a single-layer network is equivalent to finding a maximum matching ofthe network [10]. A matching is a set of directed links that do not share starting or endpoints, and a node is unmatched if there is no link in the matching pointing at it. Liu et al.showed that full control of a network is possible if each unmatched node is controlled directlyby an independent input; therefore N i is provided by the minimum number of unmatchednodes. To determine N i for a two-layer network, we first find a maximum matching of thecombined network of Layer I and Layer II. If there are no unmatched nodes in Layer II, weonly have to interact with Layer I; therefore we are done. If a node v II i is unmatched in LayerII, v I i is necessarily matched by some node v I j , otherwise the size of the matching could beincreased by adding ( v I i → v II i ). By taking out the link ( v I j → v I i ) from the matching andincluding ( v I i → v II i ) the size of the maximum matching does not change, and we moved theunmatched node from Layer II to Layer I. We repeat this for all unmatched nodes in LayerII. (Note that it may be necessary to connect inputs to additional nodes so that all nodesare reached by the control signals. Due to the one-to-one coupling between the layers thistoo can be accomplished by interacting only with Layer I.) This simplified method allowsfaster identification of N i using the Hopcroft-Karp algorithm [37] and analytically solving n i = N i /N for random networks based on calculating the fraction of always matched nodesas described in Appendix A [38–41].First, we measure n i while fixing the average degree of Layer II ( c II ) and varying theaverage degree of Layer I ( c I ). For both ER-ER and SF-SF networks, we find that n i decreasesfor increasing values of c I and converges to n IIi = N IIi /N , the normalized number of inputsneeded to control Layer II in isolation (Fig. 3a). The latter observation is easily understood: n i is determined by the fraction of unmatched nodes in the combined network of the twolayers; if c I is high enough, Layer I is perfectly matched; therefore all unmatched nodes arein Layer II. Based on the same argument, n Ii also serves as a lower bound for n i .Varying both c I and c II for ER-ER and both γ I and γ II for SF-SF with constant averagedegrees c I = c II , we find that dense networks require less inputs than sparse networks (Fig. 3b)and degree heterogeneity makes control increasingly difficult (Fig. 3c) – in line with resultsfor single-layer networks [10]. We also observe that n i is invariant to exchanging Layer I and13ayer II. This is explained by the fact that the size of the maximum matching is invariantto flipping the direction of all links, and on the ensemble level this is the same as swappingthe two layers for networks with P ( k in ) = P ( k out ).In summary, for no timescale separation controllability is equally affected by the networkstructure of both layers, and n i is greater or equal to the number of inputs necessary tocontrol any of its layers in isolation. Similarly to single-layer networks, networks with lowaverage degree and high degree heterogeneity require more independent inputs than sparsehomogeneous networks. B. Layer I updates faster ( τ I = 1, τ II > ) In the previous section we found that the network structure of the two layers equallyaffect n i if τ I = τ II = 1. This is not the case if the timescales are different, for example ifLayer I updates faster than Layer II, we expect that we need fewer inputs than in the sametimescale case by the virtue of having more opportunity to interact with the faster system(Fig. 1b). In this section we systematically study this effect using the algorithm describedin Sec. III B and analytical arguments.By measuring n i for ER-ER and SF-SF networks as a function of τ II , we find that n i monotonically decreases with increasing τ II (Fig. 4a), confirming our expectations. For bothER-ER and SF-SF networks n i ( τ II ) converges to n Ii = N Ii /N which is the normalized numberof inputs needed to control Layer I in isolation. This can be understood by the followingargument: Suppose that τ II = N , the maximum number of time steps needed to imposecontrol on any network with N nodes [42]. We use the state of Layer I at t = 0 to set thestate of Layer II at t = N , and we have N time steps to impose control on Layer I as if itwas just by itself. For a given network we define the critical timescale parameter τ cII as theminimum value of τ II for which n i ( τ II ) = n Ii . Above the critical timescale separation, Layer Icompletely determines n i ( τ II ) independent of the structure of Layer II, in other words, themultiplex nature of the system no longer plays a role in determining n i .Measuring τ cII we find that for both ER-ER and SF-SF networks τ cII monotonically increaseswith increasing c I for fixed c II , and decreases with increasing c II for fixed c I (Fig 4b). That is14 cII is the highest if Layer I is dense and Layer II is sparse. SF-SF networks have significantlylower τ cII than ER-ER networks with the same average degree.To understand the observed pattern we provide an approximation to calculate τ cII . Wecall a node v I i externally controlled if in the dynamic graph v I i,τ II is connected to an externalsignal u j,t via a disjoint control path (e.g. nodes v I A and v I B in Fig. 1c), and the number ofsuch nodes is denoted by N e ( τ II ). We have previously shown that we require N Ii independentinputs at τ cII . For each independent input and each time step, we have one control signal u i,t ;therefore we need timescale parameter τ cII = (cid:100) N e ( τ cII ) /N Ii (cid:101) (5)to insert enough signals required by the N e ( τ cII ) externally controlled nodes, where (cid:100)·(cid:101) is theceiling function. Equation (5) is not yet useful as it contains τ cII on both side. Observing that N e ( τ II ) is a monotonically increasing function of τ II and N e ( τ II = 1) = N i ( τ II = 1), we canwrite N i ( τ II = 1) ≤ N e ( τ cII ) ≤ N. (6)In the special case when Layer II is fully connected, τ cII = 1 and N e ( τ cII = 1) = N i ( τ II =1). In the case when Layer II is entirely disconnected, i.e. is composed of isolated nodes, N e ( τ II ) = N = N i ( τ II = 1). These two opposite limiting cases suggest that it is reasonable toapproximate N e ( τ cII ) by its lower bound: τ cII ≈ (cid:100) N i ( τ II = 1) /N Ii (cid:101) , (7)which entirely depends on quantities that we can easily measure or analytically compute.We find that Eq. (7) preforms remarkably well: Figure 4b compares direct measurements of τ cII to approximations obtained by using measurements and analytically computed values of N i ( τ II = 1) and N Ii . The approximation based on measurements out performs the analyticalcalculations, because the analytical results provide the expectation value of the numeratorand denominator for ER and SF network ensembles; and therefore the ceiling function isapplied to the fraction of averages, instead of averaging after applying the ceiling function.To further test the Eq. (7), we fix n i ( τ II = 1) and n Ii and we analytically calculate c I and c II forSF-SF networks with varying degree exponent γ = γ I = γ II using the framework developed15n Appendix A. Then we generate SF-SF networks and measure τ cII as a function of γ . Theapproximation predicts that τ cII remains constant, in line with our observations (Fig. 4c).The good performance of Eq. (7) is partly due to the role of the ceiling function, as it isinsensitive to changes in the numerator that are small compared to N Ii . Indeed, errors aremore pronounced if N i ( τ II = 1) /N Ii is close to an integer (e.g. data point c I = 4 . c II = 1in Fig. 4b for ER-ER), or N i ( τ II = 1) (cid:29) N Ii (e.g. data points n Ii = 0 .
084 in Fig. 4c).What we learn from this approximation is that τ cII depends only indirectly on the degreedistribution of Layer I and Layer II through the control properties of the system withouttimescale separation – N i ( τ II = 1) and N Ii . In Sec. IV A, we showed that N i ( τ II = 1) ≥ N IIi ,therefore τ cII is expected to be large if Layer I is easy to control (e.g. it is dense and hashomogeneous degree distribution) and Layer II is hard to control (e.g. it is sparse and hasheterogeneous degree distribution).In summary, if Layer I updates faster, timescale separation enhances controllability upto a critical timescale parameter τ cII , above which n i ( τ II ) = n Ii and is completely determinedby Layer I. The critical timescale parameter τ cII largely depends on the controllability of thesystem without timescale separation, it is expected to be large if Layer I is easy and LayerII is hard to control. C. Layer II updates faster ( τ I > τ II = 1 ) Finally we investigate the case when Layer II operates faster than Layer I, i.e. τ I > τ II = 1 (Fig. 1d). Measurements show that n i monotonically increases in function of τ I for both ER-ER and SF-SF networks, and n i remains constant if τ I ≥ τ cI , where τ cI isdefined for a single network (Fig. 5a). To understand these results consider the followingargument: Some nodes of Layer II are internally controlled, meaning that the state of thesenodes at t = τ I is set by the state of nodes within Layer II at t = 0 connected to them viadisjoint control paths (node v II C in Fig. 1d); while the rest of the nodes of Layer II have to becontrolled by nodes of Layer I. The maximum number of internally controlled nodes is setby the number of disjoint paths of length τ I . A directed open path traversing l links in LayerII yields a path in the dynamic graph of at most length l ; therefore if τ I > l the path can16o longer be used for control. For example, in Fig. 1a path ( v II B → v II A ) consists of a singlelink; therefore, we can use it for control if τ I = 1 (Fig. 1b) and it is no longer useful if τ I > v II C → v II C ) in Fig. 1. This predicts that n i ( τ I = ∞ ) ≥ − n cycle , (8)where n cycle = N cycle /N is the maximum fraction of nodes that can be covered with cyclesin Layer II. Furthermore, it also means that τ cI ≤ l max + 1 , (9)where l max is the maximum length of a control path that does not involve cycles, a quantitythat only depends on the structure of Layer II. We provide the formal definition l max andalgorithms to measure n cycle and l max in Appendix B.Both l max and n cycle only depend on Layer II, furthermore both strongly depend onwhether Layer II contains a strongly connected component (SCC) or not. Uncorrelated ran-dom directed networks – both ER and SF – undergo a percolation transition at c = 1 [43].If c <
1, the network is composed of small tree components, meaning the n cycle = 0 and l max is equal to the diameter D of the network. If the system is in the critical point c = 1,the size of the largest component S diverges as N → ∞ , but the relative size S/N remainszero. The largest component contains a small number of cycles; therefore D is only approx-imately equal to l max . If c >
1, a unique giant SCC emerges which contains cycles; therefore n cycle > l max is no longer directly connected to the diameter. Rigorous mathemati-cal results show that the diameter of the ER model scales as D ∼ log( N ) for c (cid:54) = 1, and D ∼ N / for c = 1, the latter corresponding to percolation transition point [44], suggestingthat the critical timescale parameter τ cI also depends on N . Indeed, Figure 6 shows that τ cI monotonically increases with N for both ER-ER and SF-SF networks.We now scan possible values of c I while keeping c II and N fixed, we find that n i ( c I )and τ cI ( c I ) quickly converges to its respective lower and upper bound provided by Eqs. (8)and (9) (Fig. 5b-c). Varying c II and keeping c I fixed shows more intricate behavior: τ cI ( c II )increases, peaks and decreases again (Fig. 5d). This is explained by changes in the structure17f Layer II: For small c II the network is composed of small components with tree structure,increasing c II agglomerates these components, thus increasing l max . For large c II , a giant SCCexists supporting many cycles, as c II increases more and more nodes can be covered withcycles reducing l max . At the critical point c ∗ II = 1 the giant SCC emerges, and the largestcomponent consists of N α nodes (0 < α <
1) with only few cycles, providing the peak of τ cI ( c II ). Although c ∗ II = 1 for both ER and SF networks in the N → ∞ limit, finite sizeeffects delay the peak of τ cI for SF-SF networks. Below the transition point, τ cI is smaller forER-ER networks than for SF-SF networks with the same average degree. In contrast, abovethe transition point SF-SF networks have larger τ cI . A likely explanation is that the cyclecover of SF networks is smaller than the cycle cover of ER networks with the same averagedegree, thus more nodes can potentially participate in the longest control path that does notinvolve cycles.The number of inputs above the critical timescale parameter n i ( τ I = ∞ ) is also affectedby the cycle cover of Layer II (Fig. 5e): For c II <
1, Layer II does not contain cycles yielding n i ( τ I = ∞ ) = 1; for large c II , Layer II can be completely covered with cycles, and n i ( τ I = ∞ )is determined by n I i , the number of inputs needed to control Layer I in isolation.In summary, if Layer II updates faster, timescale separation reduces controllability upto a critical timescale parameter τ cI . For the model networks, the value of τ cI depends onwhether Layer II has a giant SCC; τ cI has the highest value at the percolation threshold ofLayer II. If Layer II does not contain a giant SCC, degree heterogeneity decreases τ cI ; abovethe percolation threshold homogeneous networks have lower τ cI . For all timescale parameters,it remains true that ER-ER networks require less independent inputs than SF-SF networkswith the same average degree. V. CONCLUSIONS
Here we explored controllability of interconnected complex systems with a model thatincorporates common properties of these systems: (i) it consists of two layers each describedby a complex network; (ii) the operation of each layer is characterized by a different, butoften comparable timescale and (iii) the external controller only interacts with one layer18irectly. We focused on two-layer multiplex networks, meaning that we assume one-to-onecoupling between the nodes of the two layers. Our motivation for this choice was to ensureanalytical tractability and to isolate the specific role of timescales from the effect of morecomplex multilayer network structure. Results obtained for more general multilayer networkswill ultimately be shaped by a variety of features such as complex interconnectivity structure,correlations in network structure and details of dynamics. However, even by studying mul-tiplex networks, we uncovered nontrivial phenomena, attesting that without understandingeach individual effect, it is impossible to fully understand a system as a whole.Using structural controllability we were able to solve the model, thereby directly link-ing controllability to a graph combinatorial problem. We investigated the effect of networkstructure and timescales by measuring the minimum number of independent inputs neededfor control, N i . Overall we found that dense networks with homogeneous degree distribu-tion require less inputs than sparse heterogeneous networks, in line with previous resultsfor single-layer networks [10]. We showed that if we control the faster layer directly, N i de-creases with increasing timescale difference, but only up to a critical value. Above the criticaltimescale difference, N i is completely determined by the faster layer and we do not have totake into account the multiplex structure of the system. This critical timescale separation isexpected to be large if the faster layer would be easy to control and the slower layer wouldbe hard to control in isolation. If we interact with the slower layer, control is increasinglydifficult for increasing timescale difference, again up to a critical value, above which N i stilldepends on the structure of both layers. In this case the critical timescale difference largelydepends on the longest control path that does not involve cycles in the faster layer.Although our model offers only a stylized description of real systems, it is a tractable firststep towards understanding the role of timescales in control of interconnected networks. Byidentifying the network characteristics that affect important measures of controllability, suchas minimum number of inputs needed for control and critical timescale difference, our resultsserve as a starting point for future work that aims to relax some of the model’s assumptions.Some of these extensions are relatively straightforward using the tool set developed here, forexample, the effect of higher order network structures can be studied by adding correlationsto the underlying networks. Other extensions are more challenging, e.g. if the interconnection19etween the layers is incomplete or the layers contain different number of nodes, the minimuminput problem is computationally more difficult; therefore investigating such systems wouldrequire development of efficient approximation schemes. Structural control theory does nottake the link weights into account; therefore answering questions that depend on the specificstrength of the connections require the development of different tools. For example, forcontinuous-time systems the timescales are encoded in the strength of the interactions; orthe minimum control energy also depends on value of the link weights. ACKNOWLEDGEMENTS
We thank Yang-Yu Liu, Philipp H¨ovel and Zs´ofia P´enzv´alt´o for useful discussions. Wegratefully acknowledge support from the US Army Research Office Cooperative AgreementNo. W911NF-09-2-0053 and MURI Award No. W911NF-13-1-0340, and the Defense ThreatReduction Agency Basic Research Awards HDTRA1-10-1-0088 and HDTRA1-10-1-00100.
Appendix A: Analytical solution for τ I = τ II = 1 In this section we derive an analytical solution of n i = N i /N in case of τ I = τ II = 1for two-layer random networks with predefined degree distribution as defined in Sec. IV.This network model is treelike in the N → ∞ limit; therefore it lends itself to the generatingfunction formalism. The approach described here is based on calculating the fraction of nodesthat are matched in all possible maximum matchings [39]. This solution is substantiallysimpler than the one described in Ref. [10]; however, it only applies to bipartite networks (orto bipartite representations of directed networks), and cannot be generalized to unipartitenetworks.We aim to calculate the expected size of the maximum matching of the following undi-rected bipartite network B . Layer I L I and Layer II L II are generated independently eitherusing the ER or the SF model; V I and E I are the node and link sets of L I and V II and E II are the node and link sets of L II . Each node in v I i ∈ V I is split into two copies v I i, ∈ V I0 and v I i, ∈ V I1 , we draw a link ( v I i, − v I j, ) if there exists a link ( v I i → v I j ) in L I . We treat L II v I i, − v II i, ) for all i . That is all links in B connect exactly onenode in V I0 ∪ V II0 to one node in V I1 ∪ V II1 . Nodes in V I0 ∪ V I1 belong to Layer I, and nodesin V II0 ∪ V II1 belong to Layer II. The network B is the undirected version of the dynamicalgraph D without control signals.In general, multiple possible maximum matchings may exist in a network. We first cal-culate the fraction of nodes that are matched in all possible maximum matchings. It wasshown in Ref. [39] that in any network G a node v is always matched if and only if at leastone of its neighbors is not always matched in G \ v , where G \ v is the network obtainedby removing node v from G . We translate this rule to a set of self-consistent equations tocalculate the expected fraction of always matched nodes in our random network model inthe N → ∞ limit. We provide comments on the issues of applying the rule proven for finitenetworks to infinite ones at the end of this section.To proceed we define a few probabilities. We randomly select a link e connecting twonodes v I i, ∈ V I0 and v I j, ∈ V I1 . Let θ I0 be the probability that v I i, is always matched in B \ e ,and θ I1 be the probability that v I j, is always matched in B \ e . Similarly we randomly selecta link e connecting a node v I i, ∈ V I0 with a node v II i, ∈ V II1 . Let θ I,II0 be the probability thatnode v I i, is always matched in B \ e , and θ I,II1 be the probability that node v II i, is alwaysmatched in B \ e . The probabilities θ II0 and θ II1 are defined similarly. According to the ruledescribed above these quantities can be determined by the following set of equations: θ I0 = 1 − H I ( θ I1 ) θ I,II1 ,θ I1 = 1 − H I ( θ I0 ) ,θ I,II0 = 1 − G I ( θ I1 ) ,θ I,II1 = 1 − G II ( θ II0 ) ,θ II0 = 1 − H II ( θ II1 ) ,θ II1 = 1 − H II ( θ II0 ) θ I,II0 , (A1)where G I/II ( x ) = (cid:80) ∞ k =0 P I/II ( k ) x k are the generating functions of the degree distributionsand H I/II ( x ) = (cid:80) ∞ k =1 k/ (cid:104) k (cid:105) P I/II ( k ) x k − are the generating functions of the excess degreedistributions.If we remove a node v which is not always matched, the size of the maximum matching21oes not decrease. However, if v is matched in all maximum matchings, the number ofmatched nodes will decrease by two. Therefore to count the size of the maximum matching,we first count the number of nodes that are always matched. By doing so, we have doublecounted the case when an always matched node is matched by another always matched one.This case occurs for each link e that connects two nodes that are not always matched in G \ e . Combining these two contributions, the expected number of links in the matching is N match = N [1 − G I ( θ I1 ) θ I,II1 ] + N [1 − G I ( θ II0 )] + N [1 − G II ( θ I1 )] + N [1 − G II ( θ II0 ) θ I,II0 ] −− c I N (1 − θ I0 )(1 − θ I1 ) − N (1 − θ I,II0 )(1 − θ I,II1 ) − c II N (1 − θ II0 )(1 − θ II1 ) , (A2)where the first four terms count the number of nodes that are always matched in V I0 , V I1 , V II0 and V II1 , respectively; and the last three terms correct the double counting. The expectednumber of independent inputs needed is determined by the number of unmatched nodes in V I1 and V II1 : N i = 2 N − N match . (A3)Due to the links between Layer I and Layer II, the size of the maximum matching is at least N , meaning that N i ≤ N . Therefore we normalize N i by N , yielding n i = G I ( θ I1 ) θ I,II1 + G I ( θ II0 ) + G II ( θ I1 ) + G II ( θ II0 ) θ I,II0 − c I (1 − θ I0 )(1 − θ I1 ) + (1 − θ I,II0 )(1 − θ I,II1 ) + c II (1 − θ II0 )(1 − θ II1 ) . (A4) Comments on matchings in the configuration model
The method we described to calculate the expected size of the maximum matching doesnot work for unipartite ER or SF networks generally. The reason for this is that above acritical average degree c ∗ a densely connected subgraph forms, which is referred to as thecore of the network (sometimes leaf removal core or computational core) [45–47]. To deriveEq. (A1), we assume that the neighbors of a randomly selected node v are independent ofeach other in B \ v and removing a single node does not influence macroscopic properties,e.g. θ . The effect of the core is that these assumptions no longer hold and removing justa few nodes may drastically change the number of always matched nodes. Possible way ofcircumventing this problem is to introduce a new category of nodes: in addition to keeping22rack of nodes that are sometimes matched and always matched, we separately account fornodes that are almost always matched [38].The reason why the calculation works for bipartite networks is that a core in the bipartitenetwork will have two sides: all nodes on one side will be always matched and all nodes onother will be some times matched [39–41]. If the expected size of the core on the two sidesis different, finite removal of nodes will not change macroscopic properties. If the expectedsize of the two sides of the core is the same, removal of finite nodes may change which sideis always matched and which side is sometimes matched [39]. However, this does not changeexpected fraction of matched nodes; therefore does not interfere with the calculations. Appendix B: Algorithms1. Cycle cover ( N cycle ) To find the maximum cycle cover of a directed network L , we assign weight 0 to eachlink in L ; and we add a self-loop with weight 1 to each node that does not already have aself-loop. Then we find the minimum weight maximum directed matching in L augmentedwith self-loops by converting the problem to a minimum cost maximum flow problem. Themaximum matching is guaranteed to be perfect, because each node has a self-loop. Theminimum weight perfect matching in the directed network corresponds to a perfect cyclecover where the number of self-loops with weight 1 is minimized. Therefore the maximumcycle cover in L without extra self-loops is N cycle = N − W, (B1)where W is the sum of the weights of the links in the minimum weight perfect matching.
2. Longest control path not involving cycles ( l max ) In this section we provide the algorithm to measure the longest control path not involvingcycles l max of Layer II of a two-layer network for the case τ I ≥ τ II = 1. The algorithmitself serves as the precise definition of l max . 23iven a two-layer directed network M , let N cycle be the maximum number of nodes thatcan be covered by node disjoint cycles in Layer II. To measure l max , first we construct thedynamical graph D II l representing the time evolution of the Layer II between time t = 0 and t = l as if it would be isolated as defined in Sec. III A. We search for disjoint control pathsconnecting nodes at time step t = 0 with nodes at time step t = l , e.g. each control pathconnects a node v II i, with v II j,l . The maximum number of such paths N path ( l ) provides themaximum number of internally controlled nodes if τ I = l . To determine N path ( l ) we convertthe problem to a maximum flow problem: We set the capacity of each link and each nodein D II l to 1. We then find the maximum flow connecting source node set V II0 = { v II i, | i =1 , , . . . , N } to target node set V II l = { v II i,l | i = 1 , , . . . , N } using a maximum flow algorithmof choice. The maximum flow provides N path ( l ). And l max is defined as one less than thesmallest value of l such that N path ( l ) = N cycle . (B2)Figures 7 and 8 provide two examples to illustrate the calculation of l max . [1] R´eka Albert and Albert-L´aszl´o Barab´asi, “Statistical mechanics of complex networks,” Reviewsof Modern Physics , 47 (2002).[2] Mark EJ Newman, “The structure and function of complex networks,” SIAM Review ,167–256 (2003).[3] Mikko Kivel¨a, Alex Arenas, Marc Barthelemy, James P Gleeson, Yamir Moreno, and Mason APorter, “Multilayer networks,” Journal of Complex Networks , 203–271 (2014).[4] Stefano Boccaletti, G Bianconi, R Criado, Charo I Del Genio, J G´omez-Garde˜nes, M Romance,I Sendina-Nadal, Z Wang, and M Zanin, “The structure and dynamics of multilayer networks,”Physics Reports , 1–122 (2014).[5] Sergey V Buldyrev, Roni Parshani, Gerald Paul, H Eugene Stanley, and Shlomo Havlin,“Catastrophic cascade of failures in interdependent networks,” Nature , 1025–1028 (2010).[6] Charles D Brummitt, Raissa M D’Souza, and EA Leicht, “Suppressing cascades of load in nterdependent networks,” Proceedings of the National Academy of Sciences , E680–E689(2012).[7] Sergio Gomez, Albert Diaz-Guilera, Jesus Gomez-Garde˜nes, Conrad J Perez-Vicente, YamirMoreno, and Alex Arenas, “Diffusion dynamics on multiplex networks,” Phys. Rev. Lett. ,028701 (2013).[8] Xiao Fan Wang and Guanrong Chen, “Pinning control of scale-free dynamical networks,”Physica A: Statistical Mechanics and its Applications , 521–531 (2002).[9] Francesco Sorrentino, Mario di Bernardo, Franco Garofalo, and Guanrong Chen, “Controlla-bility of complex networks via pinning,” Phys. Rev. E , 046103 (2007).[10] Yang-Yu Liu, Jean-Jacques Slotine, and Albert-L´aszl´o Barab´asi, “Controllability of complexnetworks,” Nature , 167–173 (2011).[11] Wen-Xu Wang, Xuan Ni, Ying-Cheng Lai, and Celso Grebogi, “Optimizing controllability ofcomplex networks by minimum structural perturbations,” Phys. Rev. E , 026115 (2012).[12] Zhengzhong Yuan, Chen Zhao, Zengru Di, Wen-Xu Wang, and Ying-Cheng Lai, “Exact con-trollability of complex networks,” Nature Communications (2013).[13] Sean P Cornelius, William L Kath, and Adilson E Motter, “Realistic control of networkdynamics,” Nature Communications (2013).[14] M´arton P´osfai, Yang-Yu Liu, Jean-Jacques Slotine, and Albert-L´aszl´o Barab´asi, “Effect ofcorrelations on network controllability,” Scientific Reports (2013).[15] Jianxi Gao, Yang-Yu Liu, Raissa M D’Souza, and Albert-L´aszl´o Barab´asi, “Target control ofcomplex networks,” Nature Communications (2014).[16] Francesco Lo Iudice, Franco Garofalo, and Francesco Sorrentino, “Structural permeability ofcomplex networks to control signals,” Nature Communications (2015).[17] Airlie Chapman, Marzieh Nabi-Abdolyousefi, and Mehran Mesbahi, “Controllability and ob-servability of network-of-networks via cartesian products,” IEEE Trans. on Automatic Control , 2668–2679 (2014).[18] Giulia Menichetti, Luca Dall’Asta, and Ginestra Bianconi, “Control of multilayer networks,”Scientific Reports (2016).[19] Zhengzhong Yuan, Chen Zhao, Wen-Xu Wang, Zengru Di, and Ying-Cheng Lai, “Exact con- rollability of multiplex networks,” New Journal of Physics , 103036 (2014).[20] Yan Zhang, Antonios Garas, and Frank Schweitzer, “Value of peripheral nodes in controllingmultilayer scale-free networks,” Phys. Rev. E , 012309 (2016).[21] Srilata Zaheer, Stuart Albert, and Akbar Zaheer, “Time scales and organizational theory,”Academy of Management Review , 725–741 (1999).[22] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and PeterWalter, Molecular Biology of the Cell , 4th ed. (Garland Science, New York, 2002).[23] Kristina Lerman and Rumi Ghosh, “Information contagion: An empirical study of the spreadof news on digg and twitter social networks,”
Proceedings of 4th International Conference onWeblogs and Social Media , , 90–97 (2010).[24] Haewoon Kwak, Changhyun Lee, Hosung Park, and Sue Moon, “What is twitter, a socialnetwork or a news media?” in Proceedings of the 19th International Conference on WorldWide Web (ACM, 2010) pp. 591–600.[25] Eytan Bakshy, Itamar Rosenn, Cameron Marlow, and Lada Adamic, “The role of socialnetworks in information diffusion,” in
Proceedings of the 21st International Conference onWorld Wide Web (ACM, 2012) pp. 519–528.[26] Bitly Science Team, “You just shared a link. how long willpeople pay attention?” http://blog.bitly.com/post/9887686919/you-just-shared-a-link-how-long-will-people-pay (2011), [Online; accessed 16-November-2015].[27] Ailing Qia, Xuewei Jua, Qing Zhanga, and Zengqiang Chenb, “Structural controllabilityof discrete-time linear control systems with time-delay: A delay node inserting approach,”Mathematical Problens in Engineering , 1429164 (2016).[28] Rudolf Emil Kalman, “Contributions to the theory of optimal control,” Boletin SociedadMatematica Mexicana , 102–119 (1960).[29] Ching Tai Lin, “Structural controllability,” IEEE Trans. on Automatic Control , 201–208(1974).[30] Kazuo Murota, “Systems analysis by graphs and matroids,” in Algorithms and Combinatorics ,Vol. 3 (Springer Verlag Berlin, 1987).
31] Ren´e Pfitzner, Ingo Scholtes, Antonios Garas, Claudio J Tessone, and Frank Schweitzer,“Betweenness preference: Quantifying correlations in the topological dynamics of temporalnetworks,” Phys. Rev. Lett. , 198701 (2013).[32] M´arton P´osfai and Philipp H¨ovel, “Structural controllability of temporal networks,” New Jour-nal of Physics , 123055 (2014).[33] Christian Commault and Jean-Michel Dion, “Input addition and leader selection for the con-trollability of graph-based systems,” Automatica , 3322–3328 (2013).[34] G´abor Cs´ardi and Tam´as Nepusz, “The igraph software package for complex network research,”InterJournal Complex Systems , 1695 (2006).[35] Google Optimization Tools, https://developers.google.com/optimization/ (2015).[36] K-I Goh, B Kahng, and D Kim, “Universal behavior of load distribution in scale-free net-works,” Phys. Rev. Lett. , 278701 (2001).[37] John E Hopcroft and Richard M Karp, “An nˆ5/2 algorithm for maximum matchings inbipartite graphs,” SIAM Journal on Computing , 225–231 (1973).[38] Lenka Zdeborov´a and Marc M´ezard, “The number of matchings in random graphs,” Journalof Statistical Mechanics: Theory and Experiment , P05003 (2006).[39] Tao Jia, Yang-Yu Liu, Endre Cs´oka, M´arton P´osfai, Jean-Jacques Slotine, and Albert-L´aszl´oBarab´asi, “Emergence of bimodality in controlling complex networks,” Nature Communica-tions (2013).[40] Tao Jia and M´arton P´osfai, “Connecting core percolation and controllability of complex net-works,” Scientific Reports (2014).[41] M´arton P´osfai, Structure and controllability of complex networks , Ph.D. thesis, E¨otv¨os Lor´andUniversity, Budapest, Hungary (2014).[42] Rudolf Emil Kalman, “Mathematical description of linear dynamical systems,” Journal of theSociety for Industrial & Applied Mathematics, Series A: Control , 152–192 (1963).[43] N Schwartz, R Cohen, D Ben-Avraham, A-L Barab´asi, and S Havlin, “Percolation in directedscale-free networks,” Phys. Rev. E , 015104 (2002).[44] Asaf Nachmias and Yuval Peres, “Critical random graphs: diameter and mixing time,” Annalsof Probability , 1267–1286 (2008).
45] M Bauer and O Golinelli, “Core percolation in random graphs: a critical phenomena analy-sis,” The European Physical Journal B-Condensed Matter and Complex Systems , 339–352(2001).[46] L Correale, M Leone, A Pagnani, M Weigt, and Riccardo Zecchina, “Core percolation andonset of complexity in boolean networks,” Phys. Rev. Lett. , 018101 (2006).[47] Yang-Yu Liu, Endre Cs´oka, Haijun Zhou, and M´arton P´osfai, “Core percolation on complexnetworks,” Phys. Rev. Lett. , 205703 (2012). B A C A B C t=0t=1 a b
B A CB A C
Layer ILayer II Layer I
A B C
Layer IIControlSignals
A B C t=0t=1 c Layer I
A B C
Layer IIControlSignals t=2
A B C t=0t=1 d Layer I
A B C
Layer IIControlSignals t=2τ I =2, τ II =1τ I =1, τ II =1τ I =1, τ II =2 FIG. 1.
Structural controllability of two-layer multiplex networks. (a)
A two-layer network. (b-c)
To determine N i , we construct the dynamic graph representing the time evolution of thesystem from t = 0 to t = max( τ I , τ II ). The system is controllable only if all nodes at t (blue)are connected to nodes at t or nodes representing control signals (green) via disjoint paths (red). (b) In case of no timescale separation ( τ I = τ II = 1), each disjoint control path consists of a singlelink, yielding N i = 2. (c) If Layer I updates twice as frequently as Layer II ( τ I = 1, τ II = 2), weare allowed to inject control signals at time steps t = 0 and 1, reducing the number of inputs to N i = 1. (d) On the other hand, if Layer II is faster ( τ I = 2, τ II = 1), Layer II needs to supportlonger control paths, yielding N i = 3. A a b B C c A B C t=0t=1NetworkControlSignals t=2t=3
A B C t=0t=1NetworkControlSignals
FIG. 2.
Structural controllability of single-layer networks. (a)
A single-layer network, weapply inputs to nodes v A and v B . (b) The dynamic graph D N representing the time evolution ofthe dynamics from t = 0 to t = N . The system is controllable, because we can connect the set ofnodes representing control signals (green) to the set of nodes at t = N (blue) via disjoint paths(red). (c) The dynamic graph D representing the time evolution of the dynamics from t = 0 to t = 1. The system is controllable, because we can connect the control signals and nodes at t = 0(green) to the set of nodes at t = 1 (blue) via disjoint paths (red), and all nodes are accessible fromcontrol signals. n i ER-ER c II =0.5c II =1.0c II =4.0analyticaln iII n i c ISF-SF 0 1 2 3 4 5 0 1 2 3 4 5 c II c I 0 0.2 0.4 0.6 0.8 1 n i c II c I ER-ER γ II γ I a cb SF-SF
FIG. 3.
No timescale separation. (a)
Number of inputs n i in function of c I for ER-ER andSF-SF ( γ I = γ II = 2 .
5) networks. The circles represent simulations, the continuous line is theanalytical solution, and the dashed line is the analytical solution of n IIi , the number of independentinputs necessary to control Layer II in isolation [10]. (b) n i for ER-ER networks with varyingaverage degrees c I and c II . In both layers P ( k ) = P ( k in ) = P ( k out ), therefore the heatmap issymmetric with respect to the diagonal. Increasing c in either layer enhances controllability. (c) n i for SF-SF networks with c I = c II = 4 . γ I and γ II . Increasing degreeheterogeneity in either layer increases n i . Each data point is the average over 10 randomly generatednetworks with N = 10 , n i τ cII c I =1.0,c II =0.5c I =2.0,c II =1.0c I =5.0,c II =1.5n iI ER-ER 0 0.4 0.8 1 2 3 4 5 6 n i τ II τ cII SF-SF 1 25 50 τ c II c II =0.5c II =1.0c II =2.0c II =4.0 ER-ER 1 3 5 0 1 2 3 4 5 τ c II c ISF-SF 0 2 4 6 8 10 12 2 3 4 5 τ c II γ = γ I = γ II a b c n i ( τ II =1)=0.8n Ii =0.53 n Ii =0.18 n Ii =0.11 n Ii =0.084 FIG. 4.
Layer I updates faster. (a)
Number of inputs n i for single ER-ER and SF-SF ( γ I = γ II = 2 .
5) networks with N = 10 ,
000 and varying timescale parameter τ II . The number of inputs n i monotonically decreases with increasing τ II , and for τ II ≥ τ cII , n i = n Ii . (b) The critical timescaleparameter τ cII for ER-ER and SF-SF ( γ I = γ II = 2 .
5) networks with varying average degree c I and c II .The crosses represent direct measurements of τ cII ; the squares represent the approximation obtainedby applying Eq. (7) to measurements of n i ( τ II = 1) and n Ii ; and the dashed line is an approximationobtained using analytically calculated expectation values of n i ( τ II = 1) and n Ii . (c) We measure τ cII for SF-SF networks with the same n i ( τ II = 1) and n Ii as a function of γ = γ I = γ II . Equation (7)predicts that τ cII remains constant (dashed line), in line with our observations. For (b-c), each datapoint is the average over 10 randomly generated networks with N = 10 ,
000 and error bars representthe standard deviation. n i τ cI ER-ER c I =1.0,c II =0.5c I =2.0,c II =1.0c I =5.0,c II =1.5n i ( τ I = ∞ ) 0.6 0.8 1 0 10 20 30 n i τ I τ cI SF-SF 0 30 60 τ I c ER-ER c II =0.5c II =1.0c II =4.0l max +1 0 15 30 0 1 2 3 4 5 τ I c c ISF-SF 0 0.5 1 n i ( τ I = ∞ ) ER-ER 0.6 0.8 1 0 1 2 3 4 5 n i ( τ I = ∞ ) c ISF-SF c II =0.5c II =2.0c II =4.01-n cycle τ I c c II ER-ER γ =3.0 γ =2.5 γ =2.1 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 n i ( τ I = ∞ ) c II a b cd e ER-ER γ =3.0 γ =2.5 γ =2.1n Ii FIG. 5.
Layer II updates faster. (a)
Number of inputs n i for single ER-ER and SF-SF ( γ I = γ II = 2 .
5) networks with N = 10 ,
000 and varying timescale parameter τ I . The number of inputs n i monotonically increases with increasing τ I ; for τ I ≥ τ cI , n i = n i ( τ I = ∞ ). (b) τ cI as a function of c I .For c I ≤ τ cI quickly reaches its upper bound; for c I >
1, the convergence is somewhat delayed. (c) n i ( τ I = ∞ ) as a function of c I . Increasing c I facilitates control, until n i reaches its lower bound. (d) τ cI as a function of c II with fix c I = 4 .
0. The peak of τ cI corresponds to the critical point wherethe giant strongly connected component in Layer II emerges. (e) n i ( τ I = ∞ ) as a function of c II with fix c I = 4 .
0. For c II <
1, Layer II does not contain cycles, therefore n i ( τ I = ∞ ) = 1; for large c II , Layer II can be completely covered with cycles, and n i ( τ I = ∞ ) is determined by n Ii . For (b-e),each data point is the average over 10 randomly generated networks with N = 10 ,
000 and errorbars represent the standard deviation. τ c I N ER-ER γ =3.0 γ =2.5 γ =2.1~log N 10 τ c I N ER-ER γ =3.0 γ =2.5 γ =2.1~ N τ c I Na b c
ER-ER γ =3.0 γ =2.5 γ =2.1~log N FIG. 6.
Layer II updates faster – Network size effects.
Critical timescale parameter for ER-ER networks and SF-SF networks with varying network size N . (a) Layer II has no giant stronglyconnected component ( c II = 0 . < l max equals the diameter D of Layer II which scales as D ∼ log N for ER networks, and the diameter of SF networks is smaller than the diameter of ERnetworks with the same average degree. The fact that l max + 1 ≥ τ cI suggest that τ cI ∼ log( N ). (b) At the critical point c II = 1 . D ∼ N / , suggestingthat τ cI scales as a powerlaw of N . (c) Above the critical point ( c II = 4 . >
1) there is no directconnection between D and τ cI , nonetheless observations suggest τ cI ∼ log N . In contrast with the c II ≤ τ cI increases more rapidly for SF-SF networks than for ER-ER networks. Each datapoint is the average over 100 randomly generated networks with c I = 4 . a b BAC A B C l =0 l =1 l =2 DD l =3 N path ( l =1) = 2 N path ( l =2) = 1 N path ( l =3) = 0 N cycle = 0 l max = 2 FIG. 7. l max – Example 1. (a) A directed network with tree structure; therefore not containingcycles. The diameter D = 2 is the length of the longest path. (b) We count the maximum numberof disjoint control paths N path ( l ) which connect nodes at time step 0 with nodes at time step l . Wefind that l (cid:48) = 3 is the smallest value of l such that N path ( l ) = N cycle = 0; therefore l max = 2. Thereare no cycles; therefore l max = D . a b BA C A B C l =0 l =1 l =2 l =3 N path ( l =1) = 2 N path ( l =2) = 1 N path ( l =3) = 1 N cycle = 1 l max = 1 FIG. 8. l max – Example 2. (a) A directed network containing a cycle. The size of the maximumcycle cover is N cycle = 1. (b) We count the maximum number of disjoint control paths N path ( l )which connect nodes at time step 0 with nodes at time step l . We find that l (cid:48) = 2 is the smallestvalue of l such that N path ( l ) = N cycle = 1; therefore l max = 1. N path ( l ) remains non-zero for l > l max ,showing that cycles can support control paths of any length.,showing that cycles can support control paths of any length.