Taming out-of-equilibrium dynamics on interconnected networks

A wide variety of social, biological or technological systems can be described as processes taking place on networked structures in continuous interaction with other networks. We propose here a new methodology to describe, anticipate and manage, in real time, the out-of-equilibrium dynamics of processes that evolve on interconnected networks. This goal is achieved through the full analytical treatment of the phenomenology and its reduction to a two-dimensional flux diagram, allowing us to predict at every time step the dynamical consequences of modifying the links between the different ensembles. Our results are consistent with real data and the methodology can be translated to clustered networks and/or interconnected networks of any size, topology or origin, from the struggle for knowledge on innovation structures to international economic relations or disease spreading on social groups.


Supplementary Figures
Supplementary Figure 1. Evolution of a dynamical process on two collaboration networks of physicians working in Bloomington (A) and Peoria (B) (Illinois, US). The networks' properties are N A = 44 nodes, N B = 113 nodes, L A = 99, L B = 246 links, λ A,1 = 6.4144 and λ B,1 = 6.1566, from [1,2]. We simulate the spreading of knowledge following a model of innovation networks [3][4][5] where all knowledge is initially placed at the weak network B, the networks are initially disconnected and then connect through a single connector link following different connection strategies. (a) Evolution with time of the knowledge accumulated at strong network A. The connection strategies are CC (the most central node in A is connected to the most central node in B for all times), PP (the most peripheral node in A is connected to the most peripheral node in B for all times), CC-PP (the networks swap from CC to PP strategy when the CC-asymptotic equilibrium has been reached, indicated by the red arrow at t = 60), and an exhaustive method choosing the most favourable connection for A from all possible scenarios after 5 and 10 steps.  [6]. We simulate the spreading of knowledge following a model of innovation networks [3][4][5] where all knowledge is initially placed at the weak network B, the networks are initially disconnected and then connect through a single connector link following different connection strategies. (a) Evolution with time of the knowledge accumulated at strong network A. The connection strategies are CC (the most central node in A is connected to the most central node in B for all times), PP (the most peripheral node in A is connected to the most peripheral node in B for all times), CC-PP (the networks swap from CC to PP strategy when the CC-asymptotic equilibrium has been reached, indicated by the red arrow at t = 30), and an exhaustive method choosing the most favourable connection for A from all possible scenarios after 5 and 10 steps.

Supplementary Note 1. Dynamical equations of evolutionary processes on networks
The phenomenology here analysed can be applied to evolutionary processes described as where M is the transition matrix, n(t) is the vector that stands for the state of the system at time t and n(0) is the initial condition. Let us call λ i the set of eigenvalues of M, with λ i ≥ λ i+1 , and u i the corresponding eigenvectors. Its eigenvectors verify u i · u j = 0, ∀i = j and |u i | = 1, ∀i. Since M is a primitive matrix, the Perron-Frobenius theorem assures that the largest eigenvalue of M is positive, λ 1 > |λ i |, ∀ i > 1, and its associated eigenvector is positive (i.e., (u 1 ) i > 0, ∀ i). The dynamics of the system, Eq. (1), can be thus written as where α i = n(0) · u i . Furthermore, as λ 1 > |λ i |, ∀i > 1, the asymptotic state of the population is proportional to the eigenvector that corresponds to the largest eigenvalue, u 1 , and the largest eigenvalue λ 1 yields the growth rate of the population at equilibrium. If we normalise the population n(t) such that |n(t)| = 1 after each generation, n(t) → u 1 when t → ∞. A note on the mathematical notation used in this paper: when the system represents two networks A and B interconnected through a set of connector links, the first N A elements of n(t) correspond to the nodes of network A and the N B remaining elements correspond to the nodes of network B. Therefore, the length of u T,i , u A,i and u B,i is N A + N B . The first N A elements of u A,i are equal to the i-eigenvector of network A isolated and their last N B elements are equal to zero, while the first N A elements of u B,i are equal to zero and their last N B elements are equal to the i-eigenvector of network B isolated.

Supplementary Note 2. Evolution of the population distribution
The evolution with time of a generic process on two interconnected non-directed networks is given by where n(t) represents the state of the system at time t, M is the symmetric transition matrix describing the dynamical process, N T = N A + N B is the number of nodes in the network of (two) networks T that contains A, B and the set of connector links, and α i = n(0) · u T,i is the projection of the initial condition on the i-eigenvector of matrix M of network T (see Supplementary Notes 5 and 6 for the generalisation to directed networks and more than two networks respectively). We consider networks A and B such that λ A,1 > λ B,1 , with u A,i and λ A,i the i-eigenvector and the corresponding eigenvalue of matrix M A of network A respectively, and u B,i and λ B,i the i-eigenvector and the corresponding eigenvalue of matrix M B of network B. Note that M consists of four blocks of elements: two diagonal blocks M A and M B that represent the internal interactions between nodes in networks A and B respectively, and two offdiagonal blocks that will be equal to zero if networks A and B are totally disconnected, and will show the connector links from A to B and vice versa if they eventually connect. Both internal and connector links are of the same nature, being the weight of the connector links.
Applying matrix perturbation theory (assuming small ) we obtain that the first and second eigenvalues and eigenvectors of T can be expressed as quantities that are only dependent on the isolated networks A and B [8]: where and P is a matrix whose elements are P lm = P ml = 1 if nodes l of A and m of B are connected through a connector link, and P lm = P ml = 0 elsewhere. With no loss of generality, we can neglect the k > 1 elements of the 2 terms in Eqs. (6-7) because they are far less relevant than the first element of the summation as λ A,1 − λ B,k > λ A,1 − λ B,1 for k > 1. Furthermore, we approximate the initial condition n(0) to where C 2 A + C 2 B = 1, and u A,1 and u B,1 are as mentioned the first eigenvectors associated to networks A and B, meaning that initial populations on A and B before interconnecting both networks are already at equilibrium.
We introduce the distribution of population at time t between both networks x(t) as Note that (i) The evolution with time of the distribution of population x(t) in the network of networks T , combining Eqs. (3) and (9), is given by Dividing numerator and denominator in Eq. (10) by λ t T,1 and due to the exponentially fast suppression of the contributions due to higher-order terms (since λ T,i ≥ λ T,i+1 , ∀i), we obtain Normalising Eqs. (4)(5) such that |u T,1 | = |u T,2 | = 1, and introducing Eqs. (4), (5), (6) and (7) in Eq. (11), we finally obtain the temporal evolution of the system, where being cl the set of connector links. Therefore, F represents the sum of the products of the eigenvector centralities of all connector nodes measured prior to interconnecting A to B, and for this reason we will name connection strength to F and normalised connection strength to F/∆λ. Note that peripheral-peripheral (PP) connections between networks lead to low values of F , while central-central (CC) connections (or a very large number of PP connections) lead to large values of F . All k > 1 elements in the -terms of Eqs. (4)(5), as well as all elements in their 2 -terms, are not present in Eq. (12) because they are cancelled when being multiplied by u A,1 and u B,1 . This fact leads to an expression for x(t) that only depends on the first eigenvalues of the isolated networks λ A,1 and λ B,1 , the centrality of their connector nodes (through F , that depends on the first eigenvectors u A,1 and u B,1 ), the value of the weight of the connector links and the initial condition (through x 0 ).
In summary, Eq. (12) describes the dynamics of the system for all times and can be used to obtain how much population will be at each network at any time, and how the system would evolve if we changed the connector links following any connection strategy.
In case we used a generic initial condition n(0) different to Eq. (8), the analytic expression becomes more troublesome during the first steps of the process, but in practice, after some steps, the population would spread over both networks and Eq. (12) would be fully applicable. See Supplementary Note 4 for details.

Supplementary Note 3. The population flux and the asymptotic equilibria
The population fluxẋ(t) measures the speed at which the population of network A moves towards network B at time t. Making use of the expression for x(t) at Eq. (11) we obtaiṅ and introducing Eqs. (4), (5), (6) and (7) in Eq. (14) yields where K, L, ∆λ, and F are as in Eq. (12).
x(t) = 0 means that the system is in the equilibrium state. We can therefore calculate the asymptotic equilibria with an initial condition given by Eq. (8) equalising to zero the population flux shown in Eq. (14) (see Supplementary Note 7 for details):ẋ(t) = 0 for (i) t → ∞ and (ii) α 2 = n(0) · u T,2 = 0. Equation (12) leads for the former that while the latter implies that the only initial conditions that coincide with equilibrium states are also given by x 0,eq = F/(λ A,1 − λ B,1 ).

Supplementary Note 4. Population distribution and flux for arbitrary initial conditions
The dependence with time of the population distribution x(t) and its fluxẋ(t) on two interconnected networks is formally given by Eqs. (11) and (14) respectively, and they can be approximated by Eqs. (12) and (15) if the initial condition n(0) is , meaning that initial populations on A and B before connection are already in equilibrium.
But what happens if we start the process with an arbitrary initial condition n(0)? The only term of Eqs. (11)(12)(13)(14)(15)) that depends on the initial condition is α 2 /α 1 , where u T,1 and u T,2 follow Eqs. (4) and (5). Substituting L by Eq. (17) in Eq. (12) for x(t) and in Eq. (15) forẋ(t) allows to describe very precisely the evolution with time of the population distribution and its flux for all times and any initial condition n(0). The price to pay for this high precision is that these expressions depend on all the eigenvalues λ A,i and λ B,i and eigenvectors u A,i and u B,i of matrices associated to networks A and B respectively, while Eqs. (11)(12)(13)(14)(15) only depend on the first eigenvalues and eigenvectors of such matrices. On the other hand, extensive numerical work shows that in most real cases the 2 -terms in Eqs. (4)(5) are negligible, and the same applies for the k > 1 elements of the -terms because they are far less relevant than the first element of the summation as λ A,1 − λ B,k > λ A,1 − λ B,1 for k > 1 in a k , and λ B,1 −λ A,k > λ A,1 −λ B,1 for k > 1 in b k for the relevant cases of λ A,1 ≈ λ B,1 . After these approximations, we obtain Note that this expression is equal to the one obtained for networks A and B initially in equilibrium -L in Eqs. (12) and (15)-, although now the initial condition is arbitrary and x 0 = n(0)·u B,1 n(0)·u A,1 for any n(0). This result indicates that, while Eq. (12) for x(t) and (15) forẋ(t) were obtained supposing populations on networks A and B initially at equilibrium, they are also applicable for most real and general cases after few time steps.

Supplementary Note 5. Generalisation to directed networks
When the transition matrix M is primitive and diagonalisable but not symmetric, that is, when the network under study is weighted and directed, the transition matrices M A , M B and M are not symmetric. As a consequence, the leftū L k and rightū l eigenvectors of these matrices do not coincide but are biorthonormal, that is, u L k · u l = 1, ∀ k = l, and u L k · u l = 0, ∀ k = l. The formal expression of the dynamics depends on both the right and left eigenvectors of M such that where α L i is the projection of the initial condition n(0) on the i−th left eigenvector u L i of M, α L i = n(0) u L i . However, the asymptotic state of the population is again given by u 1 , which corresponds to the right eigenvector associated to the largest eigenvalue λ 1 .
The calculations developed in the Results of the main manuscript and Supplementary Notes 2 and 3 should be conveniently modified for non-symmetric matrices taking into account that Eq. (19) depends on the left eigenvectors of the transition matrix. However, the results are formally very similar, making use of the expressions obtained in [8] for the eigenvectors and eigenvalues of two interconnected directed networks (in substitution of Eqs. (4-7)). As an example, Eqs. (4) and (6) become where a k = (u L B,K Pu A,1 )/(λ A,1 − λ B,k ), is the weight of the connector links and P is the matrix that defines the connector links, and is not symmetric.
Note that to first order of only the connector links that go from the strong network to the weak network are relevant for the calculation of the eigenvector u T,1 . Furthermore, in a general case where each connector link has a different weight (as it happens in the OECD Inter-Country Input-Output (ICIO) database studied in the main manuscript), the above equations show = 1 and the weights of the links become the elements of matrix P.
Finally, let us remark that as the Perron-Frobenius Theorem holds for any irreducible non-negative square matrix, our methodology can be applied to any network (directed or undirected) as far as its associated transition matrix is also irreducible. In practice, this happens when the graph is ergodic, that is, when any pair of nodes are connected through a sequence of links. If the matrix is not irreducible, it could happen that the system presented more than one asymptotic behaviour, and the phenomenology would be more complex than the one presented in our work.

Supplementary Note 6. Generalisation to more than 2 networks
The study of the out-of-equilibrium dynamics among m > 2 networks is a very rich phenomenon that deserves a detailed study beyond the scope of this work. However, we can present here a simplified example to show that the main highlights of the methodology presented in this paper are also applicable to the more complex system of m networks. Let us suppose that a network A is connected with a network BC that consists of two networks B and C joined by several connector links (that in this case can change in time), and λ A,1 > λ B,1 > λ C,1 . The expressions for the eigenvector centrality and the first eigenvalue associated to network BC follow Eqs. (4) and (6): We introduce the distribution of population at time t between network A and BC, x A−BC (t), as If we study the connection of network A to network BC, in the way already shown in former calculations for 2 networks, we obtain the temporal evolution of the system, where being cl 1 the set of connector links that connect network A with networks B and C, and cl 2 the set of connector links that connect networks B and C.
The population flux diagram, in its simplest version, consists now of a 3-dimensional representation of the evolution with time of the distribution of population between network A and the rest of networks, x(t), where the X-axis represents 1 F 1 /(λ A,1 − λ BC,1 ), where 1 and F 1 represent the weight and strength of the connector links that connect network A with the union of networks B and C, Y-axis represents 2 F 2 (λ B,1 − λ C,1 ), where 2 and F 2 represent the weight and strength of the connector links that connect networks B and C, and Z-axis is x(t), the evolution with time of the distribution of population between network A and the rest of networks.
A similar calculation would be derived if we focused on analysing the flux of centrality between networks A and B in the same 3-networked system. In this case, the temporal evolution of the system would be given by the following expressions, . (30)

Supplementary Note 7. Equilibrium configurations
We equalise the population flux to zero to calculate the asymptotic equilibria of the system. Regarding Eq. (14), we obtain thatẋ = 0 for (i) λ T,2 = 0, (ii) λ T,1 = λ T,2 , (iii) α 2 = 0 and (iv) (u B,1 · u T,1 )(u A,1 · u T,2 ) = (u A,1 · u T,1 )(u B,1 · u T,2 ). Note that the first case, λ T,2 = 0, is not possible in real networks, Perron-Frobenius theorem implies λ T,1 > λ T,2 , and (u B,1 · u T,1 )(u A,1 · u T,2 ) < (u A,1 · u T,1 )(u B,1 · u T,2 ) for all . Therefore, the only real equilibria are given by α 2 = 0, which implies that an initial population distribution, given by x 0 = (u B,1 · n(0))/(u A,1 · n(0)), coincides with an equilibrium point if We observe here that x 0,eq grows indefinitely with the connection strength F , but this behaviour has no physical meaning. Even when Eq. (31) allows for F of any size, there is a maximum possible value of F in real systems (and in general F << 1), as it will represent the case in which both networks are connected through their most populated nodes (CC strategy). Also, the weight of the connector links should be small to make the Eqs. (4-7) applicable.
It is easy to see computationally that increasing indefinitely the weight of the links that connect two arbitrary real networks leads to an equivalent distribution of the population between both networks (P A = P B = 0.5, x(t → ∞) = 1, see how the equilibrium lines in Fig. 2(a-b) and Fig. 3 of the main manuscript tend to a horizontal line for large F/(λ A,1 − λ B,1 )). Therefore, the system can never show an equilibrium in x 0,eq >> 1, that is, when substantially more population remains on B than on A, which in our calculation leads to F/(λ A,1 − λ B,1 ) ≤ 1. That is the limit of applicability of Eqs. (12) for x(t) and (15) forẋ(t), but it is safely beyond most practical cases.
In summary, the existence of an upper bound for the connection strength F (and therefore for the normalised connection strength F/∆λ) such that x 0,eq is maximum reveals that there are no equilibrium distributions very close to u B,1 , while they can be as close to u A,1 as desired. This result coincides with what was already known: even in optimal connections for the weak network (central nodes of A connected to central nodes of B), the strong network obtains at equilibrium a large fraction of the total final population.

Supplementary Note 8. Alternating strategies
Our methodology can be applied in the context of game theory. From this viewpoint, the different networks compete for population following a certain amount of rules that should be fixed in each case. A simple example is that in which the networks alternate the choice of the connector links, trying to optimise their own population. We can solve analytically this situation, as follows.
With no loss of generality, we develop the population fluxẋ(t) -see Eq. (15)-in series of and neglect the second order terms to simplify the calculations, obtaining In the equilibrium, the population flux provoked by one competitor must counteract the competing strategy, that is,ẋ(t) = −ẋ(t) , which yieldṡ From here it is easy to see that the equilibrium distribution of alternating two different strategies coincides with that of applying both strategies simultaneously but with the link weights equal to the half of the original ones. Applying Eq. (31) we finally obtain that In summary, the distribution of population between both networks at the equilibrium due to alternating strategies is the average between the equilibria of each of the strategies applied independently. For example, if the strong network connects to the weak network through PP connections in the even turns, while the weak one connects to the strong one through CC connections in the odd turns, both trying to maximise their own final population, the equilibrium would in fact be an intermediate situation between what they would obtain if they were able to decide in all steps the best strategy for each of them.

Supplementary Note 9. Applicability of the framework
Our aim is to present a general methodology and a guide for its use, and therefore a complete analysis of a particular sample beyond the ones already studied in the main manuscript is far from the scope of this work. Anyway, we focus on several widely used mathematical models of complex systems and a very renowned network-based biotechnological dataset to remark the versatility of the methodology and help its applicability in different contexts.

Mathematical models
Our framework can be applied to the description of mutation-selection evolutionary processes. They model the evolution of heterogeneous populations (e.g. viruses or bacteria) on the space of genomes [9,10], where each node represents a genome and two nodes are connected if they differ in only one nucleotide. In fact, the space of genomes has been recently proved to behave as a network of networks [9,11]. In this context the transition matrix becomes where R is the replication rate, µ < 1 is the probability of mutating from one genome to another, I is the identity matrix, and S stands for the maximum number of neighbours of a node. A different field of application could be the control and optimisation of the growth and spreading of species on fragmented habitats, where Eq. (35) is also valid being µ the probability per unit time of an organism moving from one region to another.
We can also focus on disease spreading in social environments. For instance, the Susceptible-Infected (SI) model can be described by a transition matrix of the form M = βG, and the Susceptible-Infected-Recovered (SIR) model M = βG − γI, where β and γ are the probabilities per unit time that an healthy individual will be infected or an infected individual will recover, respectively (see [12], for details of the two different models). In these two cases, n i (t) represents the probability of individual i of being infective at time t. Note that in both models there are only two potential final states, the infection or recovery of the whole population. This fact yields that our methodology will be applicable if (i) the parameters place the system close to the epidemic threshold (where the infected population remains bounded because the largest eigenvalue of matrix M verifies λ 1 = γ/β for the SIR model), or (ii) for short times, before the whole population has either been infected or totally recovered.
Finally, the diffusion of rumours is strongly related to epidemic processes. The Maki-Thompson (MK) model [13] considers three different kinds of actors in a population: ignorants (who never heard about the rumour), spreaders (who spread the rumour to their contacts) and stiflers (who have heard the rumour but decide to stop its diffusion). One ignorant becomes a spreader of the rumour with a probability α and when two spreaders meet, one of them realises that the rumour does not have novelty and becomes a stifler with a probability β. The dynamics of the fraction of spreaders in the MK model can be described by a transition matrix of the form M = (α + β)G − βI, similar expression to that of the SIR model.

Real datasets
A fruitful application of our methodology to real data associated to networked systems could be the analysis of an extensive 12-year (1988-1999) monitoring of the evolution of the structure and dynamics of interorganisational collaboration among almost three thousand biotechnological institutions [14]. This massive study optimally fulfils the conditions to be analysed with the methodology here presented, as (i) the networks of collaboration in biotechnology are composed of many inter-connected subnetworks [15], (ii) they are permanently out of equilibrium because they are characterised by a high rate of creation and destruction of connections, as entities that do not either expand or renew their networks finally lose their central positions, and (iii) the nodes with the most central positions in the field are not the firms with larger market power or more capability of creating novelty, but those with a more diverse group of well-connected collaborators [14], supporting eigenvector centrality as the ideal quantity to measure the importance of an organism in this kind of networks.