Dynamical Response of Networks under External Perturbations: Exact Results

Dynamical Resp onse of Net w orks und er External P erturbations: Exact Results David D. Chinellato 1 , Marcus A.M. de Aguia r 1 , 2 , Irving R. Epstein 2 , 3 , Dan Braha 2 , 4 and Y aneer Bar- Y am 2 1 Instituto de F ´ ısic a ‘Gleb Wataghin ’, Uni versidade Estadual de Campi nas, Unic amp 13083-970, Campinas, SP, Br asil 2 New England Complex Systems I nstitut e, Cambridge, Massachusett s 02138 3 Dep artment of Chemistry, MS015, Br andeis University, Wal tham, Massachusett s 02454, USA 4 University of Massachus etts, Dartmouth, Massachusetts 02747 W e introduce and solv e a general model of dynamic response under external p erturbations. This mod el captures a wide range of systems out of equ ilibrium including Ising mo dels of p hysi cal sys- tems, social opinions, and p opulation genetics. The distribut ion of states under p erturbation and relaxation process refl ects t w o reg imes — one driven by the external p erturbation, and one driven by internal ordering. These regimes parallel th e disordered and ordered regimes of equilibrium ph ysical systems driven by thermal p erturbations but here are shown to be relev an t for non-th ermal and non-equilibrium ex ternal influences on complex biological and so cial sy stems. W e ex t end our results to a wide range of netw ork top ologies by in tro ducing an effective strength of ex ternal p erturbation by analytic mean-fi eld appro ximation. Sim ulations sho w this generaliza tion is remark ably accurate for many top ologies of current interes t in describing real systems. P ACS n umbers : 89.75.-k,05.50.+q,05.45.Xt Net w orks have b ecome a standard mo del for a wealth of complex systems, from physics to so cial sciences to bi- ology [1, 2]. A lar g e b ody of w ork has in v estigated to p o- logical prop erties [1, 3, 4, 5]. The r aison d’ ˆ etr e , though, of complex net work studies is to understand the rela - tionship b et w een structure a nd dynamics - from disease spreading and so c ial influence [6, 7, 8, 9, 1 0] to sea rch[11 ]. Y et, dynamic resp onse o f netw o rks under external p ertur- bations has b een less thoroughly inv es tigated [3, 12]. In this pap er w e consider a simple dynamical pro cess as a general framew ork for the dynamic r espo ns e of a netw or k to an external environmen t. The environmen t is initially treated as part of the netw ork and then generalized as an external system. W e obta in complete and exact results for the simplest case of fully connected netw orks a nd find a no n trivial dy- namic be havior that can be divided into t w o regimes. F or large p erturbations the en vironmental influence extends int o the system with a distr ibutio n which, in the ther- mo dynamic limit, b ecomes a Gaussian ar ound a v alue that reflects a balance b etw een the external p erturba- tions driving the system in different directions. F or small per turbations the distribution of states has p eaks at the t wo ordered states. Order arises from in teractions within the system, a nd power law tails res ult fr o m the external per turbation aw ay from these or dered states. The b ound- ary b et ween these reg imes is characterized by a uniform distribution where all states ar e equa lly likely . The time scale of eq uilibration is small for large p erturbations and diverges in versely as the strength of the p erturbation for small perturbations. This characterizes the switching time b ehavior o f the tw o o rdered states. W e ge neralize the exact results to net works of different top ologies using a mean field treatment. Sim ulations show that this ge n- eralization, which inv o lv es renormalizing the constants in the distr ibutions, is v ery ac c urate. Our results reveal and g eneralize key features of r e laxation and dynamic resp onse of mo dels of a wide ra nge of physical systems in the Ising universality class, elec to ral and con tagion mo dels of s o cia l systems, and the W rig h t-Fisher mo de l of evolution in po pulation biology . Spec ific a lly , we consider netw o rks with N + N 0 + N 1 no des. Each no de has an internal state which can take only the v a lues 0 or 1. W e let the N 0 no des b e frozen in state 0, and N 1 in state 1, and the remaining no des change by a dopting the state of a connected node. A t each time step a r andom fr ee node is selected; with prob- ability 1 − p the no de copies the state of o ne of its con- nected neighbo rs, and with probability p the state re- mains unchanged. The frozen no des can be interpreted as exter nal p erturbations to the subnetw o r k of free no des. Analytically extending N 0 and N 1 to b e smaller than 1 enables mo deling the case of weak coupling. This mo del generalizes our previous effor ts to derive exa c t re sults of net work dynamics [1 3] (see also[14]). This system is similar to the Ising mo del, where N 1 + N 0 by ex plicitly representing the impact of thermal per turbations play the roles of the temp eratur e T , and N 1 − N 0 acts as an external magnetic field h . O ur dynam- ics a re equiv a len t to Gla uber dynamics [15] for weak fields and high temp eratur es, wher e the Is ing model para me- ters are J / k T → 1 / ( z + N 0 + N 1 ) and h/ J → ( N 1 − N 0 ), where z is the num ber of nearest neighbors a nd J the nearest-neig h bor interaction s trength. F or low temp era- tures our model is an a lter native dynamics that also cap- tures the key kinetic prop erties of this system. Relev ant net work structures include crystalline 3-D la ttices and random net works for a morphous spin-g lasses; fully con- nected net works corresp ond to long ra nge in ter actions or the mean field a pproximation. Despite the relev a nce to the extensively studied Ising mo del, we are not aw are o f any o ther exa ct solution of the resp onse dynamics of a fully connected s ystem or explicit representation of ther- mal or o ther p erturbatio n for dynamic respons e. Specific results are av ailable only for z e r o temp eratur e dynamics in one-dimensio na l or mea n field systems. [1 6, 17] 2 Our sys tem can als o mo del an election w ith tw o can- didates [18, 19] whe r e s ome o f the voters ha ve a fixed opinion while the rest change their inten tion a ccording to the opinion of others. Another application is to e pi- demics that spread up on contact betw een infected no des (e.g., individuals or computers). Finally , the mo de l can represent an ev o lving populatio n of sexually repro ducing (haploid) o rganisms where the internal state repre s en ts one o f tw o a lleles of a gene [20]. T aking p = 1 / 2, the upda te of a no de mimics the mating o f tw o individua ls, with one parent be ing replaced by the offspring, which can receive the allele of either the mother or the father with 50% probability . Since a free node can als o copy the state of a frozen no de, the ratios N 0 / ( N + N 0 + N 1 − 1) and N 1 / ( N + N 0 + N 1 − 1) g iv e the mut ation ra tes. F or a fully connected net w ork the nodes are indistin- guishable and the state of the net w ork is fully sp ecified by the n um ber of nodes with in ternal state 1 [13]. There- fore, there are o nly N + 1 g lobal states, whic h we denote σ k , k = 0 , 1 , ..., N . The state σ k has k free no des in sta te 1 and N − k free no des in state 0. If P t ( m ) is the proba- bilit y of finding the netw or k in the sta te σ m at the time t , then P t +1 ( m ) can dep end only on P t ( m ), P t ( m + 1 ) and P t ( m − 1). The probabilities P t ( m ) define a vector of N + 1 comp onents P t . In terms of P t the dynamics is describ ed b y the equa tion P t +1 = UP t ≡  1 − (1 − p ) N ( N + N 0 + N 1 − 1) A  P t where the ev olution matrix U , and also the auxilia ry ma- trix A , is tri-diagona l. The non-zero e le men ts of A a re independent of p and are given by A m,m = 2 m ( N − m ) + N 1 ( N − m ) + N 0 m A m,m +1 = − ( m + 1)( N + N 0 − m − 1) A m,m − 1 = − ( N − m + 1)( N 1 + m − 1) . The tr a nsition proba bility from s tate σ M to σ L after a time t can b e written as P ( L, t ; M , 0) = N X r =0 1 Γ r b r M a r L λ t r . (1) where a r L and b r M are the comp onents of the right and left r-th eigenv ecto rs of the evolution matrix, a r and b r , with Γ r = b r · a r . Thus, the dynamical pro blem has been reduced to finding the righ t a nd le ft eigen vectors and the eigenv a lues o f A . It is ea sy to chec k by inspectio n of small matric es that the eigenv alues µ r of A are given b y µ r = r ( r − 1 + N 0 + N 1 ) . This implies 0 ≤ p ≤ λ r ≤ 1, where λ r are the eigenv alues of U . Because o f Eq.(1), the unit eig en v alues co mpletely determine the a s ymptotic behavior of the system. The eigensystem Aa r = µ r a r leads to the following recursion relatio n for the co efficient s a r m m +1 X j = m − 1 A mj a r j = µ r a r m (2) with a r,N + 1 = a r, − 1 ≡ 0. T o solve this eq uation w e mul- tiply the who le expressio n by x m , sum ov er m and de- fine the generating function p r ( x ) = P N m =0 a r m x m . The recursion relation then yields the following differential equation for p r x (1 − x ) p ′′ r + [(1 − N − N 0 ) − (1 + N 1 − N ) x ] p ′ r + [ N N 1 − µ r / (1 − x )] p r = 0 . (3) T o understand the a s ymptotic b eha vior o f the sys tem ( µ r = 0) w e ha ve to consider t w o cases: (a) If N 0 = N 1 = 0 then µ r = 0 leads to r = 0 o r r = 1 [13]. In this ca s e the differen tial equa tio n s implifies to xp ′′ r + (1 − N ) p ′ r = 0, whos e tw o indepe nden t solutions are p 0 ( x ) = 1 and p 1 ( x ) = x N , c o rresp onding to the all–no des–0 or all–no des–1 states res p ectively . (b) If N 0 , N 1 6 = 0 then µ r = 0 implies r = 0 . In this case equation (3) is that o f a hyperg eometric function F and we find p 0 ( x ) = F ( − N , N 1 , 1 − N − N 0 , x ), which is a finite polynomial with known co efficients a 0 m . Normaliz- ing this eigenv e c to r, we o btain the pr obability of finding the netw o rk in state σ m at larg e times: ρ ( m ) = A ( N 1 + m − 1)! ( N + N 0 − m − 1)! ( N − m )! m ! (4) where A = A ( N , N 0 , N 1 ) is a normalization. Because of the frozen no des, the dynamics will never stabilize in any sta te, but will alwa y s move from one sta te to another, with mean occ upation num ber ¯ m = N N 1 / ( N 0 + N 1 ). The surprising feature o f this solution is that for N 0 = N 1 = 1 we obtain ρ ( m ) = 1 / ( N + 1), for all v alues of N . Thus all ma croscopic states are eq ua lly likely and the system executes a random walk through the state spa ce. The dynamics at long times is dominated by the sec- ond larg est eigenv e c tor, with eigenv alue λ 1 . F or lar g e net works λ t 1 ≈ e − t/τ where τ = N ( N + N 0 + N 1 − 1) (1 − p )( N 0 + N 1 ) . (5) W e obtain a complete descr iption of the dynamics by deriving all eigenv ectors with µ r 6 = 0. The differential equation for p r ( x ) yields p r ( x ) = F (1 − r − N 0 , 1 − r − N − N 0 − N 1 , 1 − N − N 0 ,x ) (1 − x ) r − 1+ N 0 + N 1 . (6) Expanding the numerator a nd denominator in T aylor series g iv es the co efficients a r m . Although they ca n easily b e written down e xplicitly , we do not do so here. Simila r ly , defining the generating function q r ( x ) = P N + N 0 − 1 m =1 − N 1 b r m x m we obtain a differential equation for q r whose solution is q r ( x ) = x 1 − N 1 F (1 − r − N 1 , 1 − r − N − N 0 − N 1 , 1 − N − N 1 ,x ) (1 − x ) r +1 . (7) If N 0 = N 1 = 0 this solution is not v alid for r = 0 o r r = 1 , since the matrix A T bec omes singula r. In this 3 0 20 40 60 80 100 0.00 0.01 0.02 0.03 0.04 0.05 0.06 N 0 =5 N 1 =1 N 0 =5 N 1 =2 ρ (m) m N 0 =N 1 =1 N 0 =N 1 =0.5 N 0 =N 1 =10 N 0 =N 1 =2 FIG. 1: (color online) Asymptotic probability distribution for a netw ork with N = 100 and several val ues of N 0 and N 1 . case the tw o left eigenv ec tors are g iv en by b 0 ,m = 1 and b 1 ,m = N − 2 m . F or other ca ses the solution is o btained from the p ow er series expansion of q r ( x ). Equations (6) and (7) complete the solution of the proble m. In the thermodyna mic limit N → ∞ we can define contin uous v ariables x = m/ N , n 0 = N 0 / N and n 1 = N 1 / N a nd a pproximate the asymptotic distr ibution by a Gaussian ρ ( x ) = ρ 0 exp [ − ( x − x 0 ) 2 / 2 δ 2 ] with x 0 = n 1 / ( n 0 + n 1 ), ρ 0 = 1 / √ 2 π δ 2 and δ =  n 0 n 1 (1 + n 0 + n 1 ) N ( n 0 + n 1 ) 3  1 / 2 . (8) In the limit wher e n 0 , n 1 >> 1 the width dep ends only on the ratio α = n 0 /n 1 and is giv en by p α/ N / (1 + α ). The problem we just s o lv ed ca n b e generalized to tre at an external reservoir weakly coupled to the net w ork of N no des. W e note that the differential equatio ns for the generating functions p r ( x ) and q r ( x ) remain well defined for r eal N 0 and N 1 . The so lutions for the g enerating functions remain the same, except that factorials must b e replaced by g amma functions. Since N 0 / ( N + N 0 + N 1 − 1) and N 1 / ( N + N 0 + N 1 − 1) repres en t the probabilities that a free no de copies one of the frozen no de s , small v alues of N 0 and N 1 can b e interpreted as repres e nting a weak connection b et w een the free nodes and an external system con taining the frozen nodes. The external system can b e tho ug h t of as a r eservoir that affects the net work but is not affected by it. Alternatively , we can suppo se that there is a single no de fixed at 0 that is on for o nly a fraction N 0 of the time and off for the fraction 1 − N 0 , and similarly for a sing le node fixed at 1. Figure 1 shows exa mples of the distributio n ρ ( m ) for a netw ork with N = 100 and v ar ious v alues of N 0 and N 1 . Numerical simulations displaying simila r results are describ ed in [21]. Figure 2 s ho ws a n example of the time evolution of the probability density for a fully connected net work com- pared to n umerical simulations. The evolution fro m the 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 0.12 t = 10 ρ( m ) m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 0.12 t = 100 ρ( m ) m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 0.12 ρ( m ) t = 500 m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 0.10 0.12 ρ( m ) t = 10000 m FIG. 2: (color online) Time evo lution of the probability distri- bution P t for a netw ork with N = 100 and N 0 = N 1 = 5. The initial state is P 0 i = 0 . 5( δ i, 20 + δ i, 80 ). The h istograms show the av erage o ver 50,000 actual realizations of th e dynamics and the solid (red) line shows the analytical result. initial to the asymptotic time-indep endent distribution is the analog of an equilibr ation pro cess pro moted by the external system. F or sma ll v alues of N 0 and N 1 ( << 1 / ln N ), we can obtain a simplified expre ssion fo r ρ ( m ): ρ ( m ) ≈ N 1 N 0 N 0 + N 1  1 − N 1 ln N m 1 − N 1 + 1 − N 0 ln N ( N − m ) 1 − N 0  . (9) Thu s ρ ( m ) displa ys a p ow er law be ha vior on b oth ends of the curve: 1 /m for m close to 0 and 1 / ( N − m ) for m close to N (see, for instance, the curve with N 0 = N 1 = 0 . 5 in Fig. 1). Since the relaxation time τ is prop ortional to 1 / ( N 0 + N 1 ), the equilibr ation pro ces s bec o mes very slow in this limit. F or net w orks with differe nt top olog ies the effect o f the frozen no des is amplified. T o see this we note that the probability that a free no de copies a frozen no de is P i = ( N 0 + N 1 ) / ( N 0 + N 1 + k i ) where k i is the degree of the no de. F or fully connected netw o rks k i = N − 1 and P i ≡ P F C . F or ge ner al netw o rks an a verage v alue P av can be calc ulated by replacing k i by the av er age degree k av . W e can then define effective num b ers of frozen nodes, N 0 ef and N 1 ef , as b eing the v a lues of N 0 and N 1 in P F C for which P av ≡ P F C . This leads to N 0 ef = f N 0 , N 1 ef = f N 1 (10) where f = ( N − 1) /k av . Correctio ns inv olving higher mo- men ts c a n b e obtained by in tegrating P i with the degree distribution and expa nding a round k av . Figure 3 shows examples o f the equilibrium distri- bution for four different net w orks with N = 100 and N 0 = N 1 = 5. Panel (a) shows a random net work with connection probability 0 . 3 ( N av = 30 , f = 3 . 3). The theor etical res ult w as o btained with Eq. (4) with N 0 ef = N 1 ef = 17. F or a scale-free netw o rk (pa nel (b)) 4 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 Random ρ (m) m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 Scale Free ρ (m) m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 Regular 2-D Lattice ρ (m) m 0 20 40 60 80 100 0.00 0.02 0.04 0.06 0.08 Small World Network ρ (m) m FIG. 3: (color online) Asymptotic probability distribution for netw orks with differen t t opologies. I n all cases N = 100 , N 0 = N 1 = 5, t = 10 , 000, and th e num b er of reali zations is 50 , 000. The theoretical (red) cu rv e is drawn with effective num b ers of frozen no des N 0 ef = f N 0 and N 1 ef = f N 1 : (a) random netw ork N 0 ef = N 1 ef = 17; ( b ) scale-free N 0 ef = N 1 ef = 82; (c) regular 2-D lattice N 0 ef = N 1 ef = 140; ( d) small w orld netw ork N 0 ef = N 1 ef = 140. grown fro m an initial cluster o f 6 no des adding no des with 3 connections each following the prefere ntial a ttac hmen t rule [1], f = 99 / 6 a nd the effective v alues of N 0 and N 1 are appr o ximately 82. Panel (c) shows the probability distribution for a 2- D regular lattice with 10 × 10 no des for which f = 99 / 3 . 6 ≈ 28 . Finally , pa nel (d) shows a small world version of the reg ular lattice [1], where 30 connections were r andomly re- c o nnected, cr e ating short- cuts betw een o therwise dista n t no des. These results show that the mean field generaliza tion is accura te for many net work top ologies . Still, extreme cases such as a sta r net work should b e different and this is confirmed b y sim- ulations and preliminary analytic results. The relaxation time (5), in units of net w ork siz e and for p = 0, b ecomes τ / N = ( k av + N 0 + N 1 ) / ( N 0 + N 1 ). It increases linearly with N for fully connected or r andom netw orks, but is independent of N for regular and scale-free top ologies. Our results hav e impo rtant implications for real sys- tems. In the so cial sciences they show the impor tance of opinion makers in stabilizing the outcome of elec tio ns: weak external influences r esult in an arbitra ry but seem- ingly strong o pinio n that can switc h at random (see als o [22]), due to the arbitrar y c hoice of the o rdered state in the weak p erturbation regime. Thus, a n elected candi- date winning a landslide election may have no solid sup- po rt. The s lo w dynamics can play a crucial role, since the time to switching might o ccur only after the election day , especia lly if the nu mber of voters is large. Stronge r external influences , counter intu itively , reduce the rela - tion time g iving r is e to improved internal eq uilibration. In theoretical bio lo gy our r esults are equiv alent to the ex- act dynamical so lutio n of the W right-Fisher mo del [20] for ar bitrary p opulation sizes and mutation rates. Our equations give not only the a s ymptotic e q uilibrium dis - tribution o f alleles (see [20] for approximate expr essions), but also its time evolution in t wo regimes, one wher e mu- tations ha ve difficult y ov er coming a n ex is ting do minan t allele (the low p erturbation regime) and one where ra n- dom mutations do minate (the hig h perturbatio n regime). Again this is crucial information, since the equilibratio n time ca n b e e x tremely long for the typically small m uta- tion rates obser v ed in nature. Finally , we emphasize that exact dynamica l solutions, describing sys tems out of equilibrium, are ra re and impo rtant for the study of many statistica l pr oper ties that can b e describ ed in future publica tions. M.A.M.A. and D.D.C. a cknowledge financial supp ort from CNPq and F APESP . I.R.E. was supp orted by Na- tional Science F oundation grant CHE-0 526866 . [1] R . Alb ert and A.-L. Barab´ asi, Rev . Mo d. Phys. 74 , 47 (2002). [2] S . Boccaletti, V. Latora, Y. Moreno, M. Cha vez and D.- U. Hwang, Phys. Rep. 424 , 175 (2006). [3] Y . Bar-Y am and I . Epstein, PN AS 101 , 4341 (2004). [4] R . A lbert, H . Jeong, and A.-L. Barabsi, N atu re Lon d on 406 , 378 (2000). [5] R . Cohen, K. Erez, D. b en-Avraham, and S. Havlin, Phys. Rev. Lett. 85 , 4626 (2000). [6] R . P astor-Satorras and A. V espignani, Phys. Rev. Lett. 86 , 3200 (2001). [7] M. Barahona and L.M. Peco ra, Phys. Rev. Lett . 89 , 054101 (2002). [8] T. Nishik aw a, A.E. Motter, Y.-C. Lai, and F.C. Hoppen - steadt, Phys. Rev. Lett. 91 , 014101 (2003). [9] Y . Moreno, M. Neko vee and A.F. Pac h eco, Phys. R ev. E 69 , 066130 (2004). [10] M.F. Laguna, A. Guillermo, and H. Zanette Damin, Physica A 329 , 459 (2003). [11] R. Guimera et al., Phys. Rev. Lett. 89 , 248701 (2002). [12] X. W ang, Y.-C. Lai, and C. H. Lai, Ph ys. Rev . E 74 , 066104 (2006). [13] M.A.M. de Aguiar, I.R. Epstein and Y. Bar-Y am, Phys. Rev. E 72 , 067102 (2005). [14] H. Zhou and R. Lip o wsky , J. Stat. Mech., P01009 (2007). [15] R.J. Glauber, J. Math. Phys. 4 , 294 (1963) [16] A. Prados and J. J. Brey , J. Phys. A: Math. Gen. 34 L453 (2001). [17] Dj. Spaso jevic, S . Janicevic and M. Knezevi, Europhys. Lett., 76 912 (2006). [18] D. Vilone and C. Castellano, Phys. Rev. E 69 , 016109 (2004). [19] V. So o d and S. Red ner, Phys. Rev. Lett. 94 ,178701 (2005). [20] W.J. Ewens, Mathematic al Population Genetics , Springer - NY ( 2004). 5 [21] N. Boccara, eprint, arXiv:0704.179 0 [nlin.AO] 2007. [22] M. Mobilia, A. P etersen and S R edner, J. S tat. Mech. P0802(200 7).